source: flexpart.git/src/timemanager_mpi.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was d8eed02, checked in by Espen Sollum ATMOS <eso@…>, 2 years ago

Minor edits

  • Property mode set to 100644
File size: 33.4 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine timemanager(metdata_format)
23
24!*****************************************************************************
25!                                                                            *
26! Handles the computation of trajectories, i.e. determines which             *
27! trajectories have to be computed at what time.                             *
28! Manages dry+wet deposition routines, radioactive decay and the computation *
29! of concentrations.                                                         *
30!                                                                            *
31!     Author: A. Stohl                                                       *
32!                                                                            *
33!     20 May 1996                                                            *
34!                                                                            *
35!*****************************************************************************
36!  Changes, Bernd C. Krueger, Feb. 2001:                                     *
37!        Call of convmix when new windfield is read                          *
38!------------------------------------                                        *
39!  Changes Petra Seibert, Sept 2002                                          *
40!     fix wet scavenging problem                                             *
41!     Code may not be correct for decay of deposition!                       *
42!  Changes Petra Seibert, Nov 2002                                           *
43!     call convection BEFORE new fields are read in BWD mode                 *
44!  Changes Caroline Forster, Feb 2005                                        *
45!   new interface between flexpart and convection scheme                     *
46!   Emanuel's latest subroutine convect43c.f is used                         *
47!  Changes Stefan Henne, Harald Sodemann, 2013-2014                          *
48!   added netcdf output code                                                 *
49!  Changes Espen Sollum 2014                                                 *
50!   MPI version                                                              *
51!   Variables uap,ucp,uzp,us,vs,ws,cbt now in module com_mod                 *
52!  Unified ECMWF and GFS builds                                              *
53!   Marian Harustak, 12.5.2017                                               *
54!   - Added passing of metdata_format as it was needed by called routines    *
55!*****************************************************************************
56!                                                                            *
57! Variables:                                                                 *
58! dep                .true. if either wet or dry deposition is switched on   *
59! decay(maxspec) [1/s] decay constant for radioactive decay                  *
60! drydep             .true. if dry deposition is switched on                 *
61! ideltas [s]        modelling period                                        *
62! itime [s]          actual temporal position of calculation                 *
63! ldeltat [s]        time since computation of radioact. decay of depositions*
64! loutaver [s]       averaging period for concentration calculations         *
65! loutend [s]        end of averaging for concentration calculations         *
66! loutnext [s]       next time at which output fields shall be centered      *
67! loutsample [s]     sampling interval for averaging of concentrations       *
68! loutstart [s]      start of averaging for concentration calculations       *
69! loutstep [s]       time interval for which concentrations shall be         *
70!                    calculated                                              *
71! npoint(maxpart)    index, which starting point the trajectory has          *
72!                    starting positions of trajectories                      *
73! nstop              serves as indicator for fate of particles               *
74!                    in the particle loop                                    *
75! nstop1             serves as indicator for wind fields (see getfields)     *
76! memstat            additional indicator for wind fields (see getfields)    *
77! outnum             number of samples for each concentration calculation    *
78! prob               probability of absorption at ground due to dry          *
79!                    deposition                                              *
80! wetdep             .true. if wet deposition is switched on                 *
81! weight             weight for each concentration sample (1/2 or 1)         *
82! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to          *
83!                    turbulence                                              *
84! us(maxpart),vs(maxpart),ws(maxpart) = random velocities due to inter-      *
85!                    polation                                                *
86! xtra1(maxpart), ytra1(maxpart), ztra1(maxpart) =                           *
87!                    spatial positions of trajectories                       *
88! metdata_format     format of metdata (ecmwf/gfs)                           *
89!                                                                            *
90! Constants:                                                                 *
91! maxpart            maximum number of trajectories                          *
92!                                                                            *
93!*****************************************************************************
94
95  use unc_mod
96  use point_mod
97  use xmass_mod
98  use flux_mod
99  use outg_mod
100  use oh_mod
101  use par_mod
102  use com_mod
103  use mpi_mod
104  use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
105       &concoutput_surf_netcdf,concoutput_surf_nest_netcdf
106
107  implicit none
108
109  integer :: metdata_format
110  logical :: reqv_state=.false. ! .true. if waiting for a MPI_Irecv to complete
111  integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0
112! integer :: ksp
113  integer :: ip
114  integer :: loutnext,loutstart,loutend
115  integer :: ix,jy,ldeltat,itage,nage
116  integer :: i_nan=0,ii_nan,total_nan_intl=0  !added by mc to check instability in CBL scheme
117  integer :: numpart_tot_mpi ! for summing particles on all processes
118  real :: outnum,weight,prob(maxspec)
119  real :: decfact
120
121  real(sp) :: gridtotalunc
122  real(dep_prec) :: drygridtotalunc=0_dep_prec,wetgridtotalunc=0_dep_prec,&
123       & drydeposit(maxspec)=0_dep_prec
124  real :: xold,yold,zold,xmassfract
125  real, parameter :: e_inv = 1.0/exp(1.0)
126
127! Measure time spent in timemanager
128  if (mp_measure_time) call mpif_mtime('timemanager',0)
129
130! First output for time 0
131!************************
132
133  loutnext=loutstep/2
134  outnum=0.
135  loutstart=loutnext-loutaver/2
136  loutend=loutnext+loutaver/2
137
138
139!**********************************************************************
140! Loop over the whole modelling period in time steps of mintime seconds
141!**********************************************************************
142
143
144  if (lroot.or.mp_dev_mode) then
145  !  write(*,45) itime,numpart*mp_partgroup_np,gridtotalunc,wetgridtotalunc,drygridtotalunc
146    write(*,46) float(itime)/3600,itime,numpart*mp_partgroup_np
147   
148    if (verbosity.gt.0) then
149      write (*,*) 'timemanager> starting simulation'
150    end if
151  end if ! (lroot)
152
153!CGZ-lifetime: set lifetime to 0
154  ! checklifetime(:,:)=0
155  ! species_lifetime(:,:)=0
156  ! print*, 'Initialized lifetime'
157!CGZ-lifetime: set lifetime to 0
158
159
160
161  do itime=0,ideltas,lsynctime
162   
163
164! Computation of wet deposition, OH reaction and mass transfer
165! between two species every lsynctime seconds
166! maybe wet depo frequency can be relaxed later but better be on safe side
167! wetdepo must be called BEFORE new fields are read in but should not
168! be called in the very beginning before any fields are loaded, or
169! before particles are in the system
170! Code may not be correct for decay of deposition
171! changed by Petra Seibert 9/02
172!********************************************************************
173
174    if (mp_dbg_mode) write(*,*) 'myid, itime: ',mp_pid,itime
175   
176    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
177      if (verbosity.gt.0) then
178        write (*,*) 'timemanager> call wetdepo'
179      endif
180
181      if (mp_measure_time) call mpif_mtime('wetdepo',0)
182       
183! readwind process skips this step
184      if (.not.(lmpreader.and.lmp_use_reader)) then
185        call wetdepo(itime,lsynctime,loutnext)
186      end if
187
188      if (mp_measure_time) call mpif_mtime('wetdepo',1)
189    end if
190
191    if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) then
192! readwind process skips this step
193      if (.not.(lmpreader.and.lmp_use_reader)) then
194        call ohreaction(itime,lsynctime,loutnext)
195      endif
196    end if
197
198    if (assspec .and. itime .ne. 0 .and. numpart .gt. 0) then
199      stop 'associated species not yet implemented!'
200!     call transferspec(itime,lsynctime,loutnext)
201    endif
202
203
204! compute convection for backward runs
205!*************************************
206
207    if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) then
208      if (verbosity.gt.0) then
209        write (*,*) 'timemanager> call convmix -- backward'
210      endif
211! readwind process skips this step
212      if (.not.(lmpreader.and.lmp_use_reader)) call convmix(itime,metdata_format)
213    endif
214
215! Get necessary wind fields if not available
216!*******************************************
217    if (verbosity.gt.0 .and. lmpreader) then
218      write (*,*) 'timemanager> call getfields'
219    endif
220
221! This time measure includes reading/MPI communication (for the reader process),
222! or MPI communication time only (for other processes)
223    if (mp_measure_time) call mpif_mtime('getfields',0)
224
225    call getfields(itime,nstop1,memstat,metdata_format)
226
227    if (mp_measure_time) call mpif_mtime('getfields',1)
228
229
230! Broadcast fields to all MPI processes
231! Skip if all processes have called getfields or if no new fields
232!*****************************************************************
233
234    if (mp_measure_time.and..not.(lmpreader.and.lmp_use_reader)) call mpif_mtime('getfields',0)
235
236! Two approaches to MPI getfields is implemented:
237! Version 1 (lmp_sync=.true.) uses a read-ahead process where send/recv is done
238! in sync at start of each new field time interval
239!
240! Version 2 (lmp_sync=.false.) is for holding three fields in memory. Uses a
241! read-ahead process where sending/receiving of the 3rd fields is done in
242! the background in parallel with performing computations with fields 1&2
243!********************************************************************************
244
245    if (lmp_sync.and.lmp_use_reader.and.memstat.gt.0) then
246      call mpif_gf_send_vars(memstat)
247      if (numbnests>0) call mpif_gf_send_vars_nest(memstat)
248! Version 2  (lmp_sync=.false.) is also used whenever 2 new fields are
249! read (as at first time step), in which case async send/recv is impossible.
250    else if (.not.lmp_sync.and.lmp_use_reader.and.memstat.ge.32) then
251      call mpif_gf_send_vars(memstat)
252      if (numbnests>0) call mpif_gf_send_vars_nest(memstat)
253    end if
254
255    if (.not.lmp_sync) then
256   
257! Reader process:
258      if (memstat.gt.0..and.memstat.lt.32.and.lmp_use_reader.and.lmpreader) then
259        if (mp_dev_mode) write(*,*) 'Reader process: calling mpif_gf_send_vars_async'
260        call mpif_gf_send_vars_async(memstat)
261        if (numbnests>0) call mpif_gf_send_vars_nest_async(memstat)
262      end if
263
264! Completion check:
265! Issued at start of each new field period.
266      if (memstat.ne.0.and.memstat.lt.32.and.lmp_use_reader) then
267        call mpif_gf_request
268      end if
269
270! Recveiving process(es):
271! eso TODO: at this point we do not know if clwc/ciwc will be available
272! at next time step. Issue receive request anyway, cancel at mpif_gf_request
273      if (memstat.gt.0.and.lmp_use_reader.and..not.lmpreader) then
274        if (mp_dev_mode) write(*,*) 'Receiving process: calling mpif_gf_send_vars_async. PID: ', mp_pid
275        call mpif_gf_recv_vars_async(memstat)
276        if (numbnests>0) call mpif_gf_recv_vars_nest_async(memstat)
277      end if
278
279    end if
280
281    if (mp_measure_time.and..not.(lmpreader.and.lmp_use_reader)) call mpif_mtime('getfields',1)
282
283! For validation and tests: call the function below to set all fields to simple
284! homogeneous values
285!    if (mp_dev_mode) call set_fields_synthetic
286
287!*******************************************************************************
288
289    if (lmpreader.and.nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
290
291! Reader process goes back to top of time loop (unless simulation end)
292!*********************************************************************
293
294    if (lmpreader.and.lmp_use_reader) then
295      if (itime.lt.ideltas*ldirect) then
296        cycle
297      else
298        goto 999
299      end if
300    end if
301
302
303! Get hourly OH fields if not available
304!****************************************************
305    if (OHREA) then
306      if (verbosity.gt.0) then
307        write (*,*) 'timemanager> call gethourlyOH'
308      endif
309      call gethourlyOH(itime)
310    endif
311
312
313! Release particles
314!******************
315
316    if (verbosity.gt.0.and.lroot) then
317      write (*,*) 'timemanager>  Release particles'
318    endif
319
320    if (mdomainfill.ge.1) then
321      if (itime.eq.0) then
322        if (verbosity.gt.0) then
323          write (*,*) 'timemanager>  call init_domainfill'
324        endif
325        call init_domainfill
326      else
327        if (verbosity.gt.0.and.lroot) then
328          write (*,*) 'timemanager>  call boundcond_domainfill'
329        endif
330        call boundcond_domainfill(itime,loutend)
331      endif
332    else
333      if (verbosity.gt.0.and.lroot) then
334        print*,'call releaseparticles' 
335      endif
336      call releaseparticles(itime)
337    endif
338
339
340! Check if particles should be redistributed among processes
341!***********************************************************
342    call mpif_calculate_part_redist(itime)
343
344
345! Compute convective mixing for forward runs
346! for backward runs it is done before next windfield is read in
347!**************************************************************
348
349    if ((ldirect.eq.1).and.(lconvection.eq.1)) then
350      if (verbosity.gt.0) then
351        write (*,*) 'timemanager> call convmix -- forward'
352      endif
353      call convmix(itime,metdata_format)
354    endif
355
356! If middle of averaging period of output fields is reached, accumulated
357! deposited mass radioactively decays
358!***********************************************************************
359
360    if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
361      do ks=1,nspec
362        do kp=1,maxpointspec_act
363          if (decay(ks).gt.0.) then ! TODO move this statement up 2 levels
364            do nage=1,nageclass
365              do l=1,nclassunc
366! Mother output grid
367                do jy=0,numygrid-1
368                  do ix=0,numxgrid-1
369                    wetgridunc(ix,jy,ks,kp,l,nage)= &
370                         wetgridunc(ix,jy,ks,kp,l,nage)* &
371                         exp(-1.*outstep*decay(ks))
372                    drygridunc(ix,jy,ks,kp,l,nage)= &
373                         drygridunc(ix,jy,ks,kp,l,nage)* &
374                         exp(-1.*outstep*decay(ks))
375                  end do
376                end do
377! Nested output grid
378                if (nested_output.eq.1) then
379                  do jy=0,numygridn-1
380                    do ix=0,numxgridn-1
381                      wetgriduncn(ix,jy,ks,kp,l,nage)= &
382                           wetgriduncn(ix,jy,ks,kp,l,nage)* &
383                           exp(-1.*outstep*decay(ks))
384                      drygriduncn(ix,jy,ks,kp,l,nage)= &
385                           drygriduncn(ix,jy,ks,kp,l,nage)* &
386                           exp(-1.*outstep*decay(ks))
387                    end do
388                  end do
389                endif
390              end do
391            end do
392          endif
393        end do
394      end do
395    endif
396
397
398!!! CHANGE: These lines may be switched on to check the conservation
399!!! of mass within FLEXPART
400!   if (itime.eq.loutnext) then
401!   do 247 ksp=1, nspec
402!   do 247 kp=1, maxpointspec_act
403!47         xm(ksp,kp)=0.
404
405!   do 249 ksp=1, nspec
406!     do 249 j=1,numpart
407!          if (ioutputforeachrelease.eq.1) then
408!            kp=npoint(j)
409!          else
410!            kp=1
411!          endif
412!       if (itra1(j).eq.itime) then
413!          xm(ksp,kp)=xm(ksp,kp)+xmass1(j,ksp)
414!         write(*,*) 'xmass: ',xmass1(j,ksp),j,ksp,nspec
415!       endif
416!49     continue
417!  do 248 ksp=1,nspec
418!  do 248 kp=1,maxpointspec_act
419!  xm_depw(ksp,kp)=0.
420!  xm_depd(ksp,kp)=0.
421!     do 248 nage=1,nageclass
422!       do 248 ix=0,numxgrid-1
423!         do 248 jy=0,numygrid-1
424!           do 248 l=1,nclassunc
425!              xm_depw(ksp,kp)=xm_depw(ksp,kp)
426!    +                  +wetgridunc(ix,jy,ksp,kp,l,nage)
427!48                 xm_depd(ksp,kp)=xm_depd(ksp,kp)
428!    +                  +drygridunc(ix,jy,ksp,kp,l,nage)
429!             do 246 ksp=1,nspec
430!46                    write(88,'(2i10,3e12.3)')
431!    +              itime,ksp,(xm(ksp,kp),kp=1,maxpointspec_act),
432!    +                (xm_depw(ksp,kp),kp=1,maxpointspec_act),
433!    +                (xm_depd(ksp,kp),kp=1,maxpointspec_act)
434!  endif
435!!! CHANGE
436
437
438
439
440! Check whether concentrations are to be calculated
441!**************************************************
442
443    if ((ldirect*itime.ge.ldirect*loutstart).and. &
444         (ldirect*itime.le.ldirect*loutend)) then ! add to grid
445      if (mod(itime-loutstart,loutsample).eq.0) then
446
447! If we are exactly at the start or end of the concentration averaging interval,
448! give only half the weight to this sample
449!*****************************************************************************
450
451        if ((itime.eq.loutstart).or.(itime.eq.loutend)) then
452          weight=0.5
453        else
454          weight=1.0
455        endif
456        outnum=outnum+weight
457
458        call conccalc(itime,weight)
459
460      endif
461
462! :TODO: MPI output of particle positions;  each process sequentially
463!   access the same file
464      if ((mquasilag.eq.1).and.(itime.eq.(loutstart+loutend)/2)) &
465           call partoutput_short(itime)    ! dump particle positions in extremely compressed format
466
467
468! Output and reinitialization of grid
469! If necessary, first sample of new grid is also taken
470!*****************************************************
471
472      if ((itime.eq.loutend).and.(outnum.gt.0.)) then
473        if ((iout.le.3.).or.(iout.eq.5)) then
474
475! MPI: Root process collects/sums grids
476!**************************************
477          call mpif_tm_reduce_grid
478
479          if (mp_measure_time) call mpif_mtime('iotime',0)
480          if (surf_only.ne.1) then
481            if (lroot) then
482              if (lnetcdfout.eq.1) then
483                call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,&
484                     &drygridtotalunc)
485              else
486                call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
487              endif
488            else
489! zero arrays on non-root processes
490              gridunc(:,:,:,:,:,:,:)=0.
491              creceptor(:,:)=0.
492            end if
493          else
494            if (lroot) then
495              if (lnetcdfout.eq.1) then
496                call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,&
497                     &drygridtotalunc)
498              else
499                call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
500              end if
501            else
502! zero arrays on non-root processes
503              gridunc(:,:,:,:,:,:,:)=0.
504              creceptor(:,:)=0.
505            endif
506          endif
507          if (mp_measure_time) call mpif_mtime('iotime',1)
508
509          if (nested_output.eq.1) then
510
511! MPI: Root process collects/sums nested grids
512!*********************************************
513            call mpif_tm_reduce_grid_nest
514 
515           if (mp_measure_time) call mpif_mtime('iotime',0)
516
517            if (lnetcdfout.eq.0) then
518              if (surf_only.ne.1) then
519
520                if (lroot) then
521                  call concoutput_nest(itime,outnum)
522                else
523                  griduncn(:,:,:,:,:,:,:)=0.
524                end if
525
526              else  ! :TODO: check for zeroing in the netcdf module
527                call concoutput_surf_nest(itime,outnum)
528
529              end if
530
531            else
532
533              if (surf_only.ne.1) then
534                if (lroot) then             
535                  call concoutput_nest_netcdf(itime,outnum)
536                else
537                  griduncn(:,:,:,:,:,:,:)=0.
538                end if
539              else
540                if (lroot) then
541                  call concoutput_surf_nest_netcdf(itime,outnum)
542                else
543                  griduncn(:,:,:,:,:,:,:)=0.
544                end if
545              endif
546
547
548            end if
549          end if
550         
551
552          outnum=0.
553        endif
554        if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
555        if (iflux.eq.1) call fluxoutput(itime)
556        if (mp_measure_time) call mpif_mtime('iotime',1)
557
558
559! Decide whether to write an estimate of the number of particles released,
560! or exact number (require MPI reduce operation)
561        if (mp_dbg_mode) then
562          numpart_tot_mpi = numpart
563        else
564          numpart_tot_mpi = numpart*mp_partgroup_np
565        end if
566
567        if (mp_exact_numpart.and..not.(lmpreader.and.lmp_use_reader).and.&
568             &.not.mp_dbg_mode) then
569          call MPI_Reduce(numpart, numpart_tot_mpi, 1, MPI_INTEGER, MPI_SUM, id_root, &
570               & mp_comm_used, mp_ierr)
571        endif
572       
573        !CGZ-lifetime: output species lifetime
574        if (lroot.or.mp_dbg_mode) then
575        !   write(*,*) 'Overview species lifetime in days', &
576        !        real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
577        !   write(*,*) 'all info:',species_lifetime
578          write(*,45) itime,numpart_tot_mpi,gridtotalunc,&
579               &wetgridtotalunc,drygridtotalunc
580        !   if (verbosity.gt.0) then
581        !     write (*,*) 'timemanager> starting simulation'
582        !   end if
583        end if
584
585        ! Write number of particles for all processes
586        if (mp_dev_mode) write(*,*) "PID, itime, numpart", mp_pid,itime,numpart
587
588
58945      format(i13,' SECONDS SIMULATED: ',i13, ' PARTICLES:    Uncertainty: ',3f7.3)
59046      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
591        if (ipout.ge.1) then
592          do ip=0, mp_partgroup_np-1
593            if (ip.eq.mp_partid) call partoutput(itime) ! dump particle positions
594            call mpif_mpi_barrier
595          end do
596        end if
597
598        loutnext=loutnext+loutstep
599        loutstart=loutnext-loutaver/2
600        loutend=loutnext+loutaver/2
601        if (itime.eq.loutstart) then
602          weight=0.5
603          outnum=outnum+weight
604          call conccalc(itime,weight)
605        endif
606
607
608! Check, whether particles are to be split:
609! If so, create new particles and attribute all information from the old
610! particles also to the new ones; old and new particles both get half the
611! mass of the old ones
612!************************************************************************
613       
614        if (ldirect*itime.ge.ldirect*itsplit) then
615          n=numpart
616          do j=1,numpart
617            if (ldirect*itime.ge.ldirect*itrasplit(j)) then
618!                if (n.lt.maxpart) then
619              if (n.lt.maxpart_mpi) then
620                n=n+1
621                itrasplit(j)=2*(itrasplit(j)-itramem(j))+itramem(j)
622                itrasplit(n)=itrasplit(j)
623                itramem(n)=itramem(j)
624                itra1(n)=itra1(j)
625                idt(n)=idt(j)
626                npoint(n)=npoint(j)
627                nclass(n)=nclass(j)
628                xtra1(n)=xtra1(j)
629                ytra1(n)=ytra1(j)
630                ztra1(n)=ztra1(j)
631                uap(n)=uap(j)
632                ucp(n)=ucp(j)
633                uzp(n)=uzp(j)
634                us(n)=us(j)
635                vs(n)=vs(j)
636                ws(n)=ws(j)
637                cbt(n)=cbt(j)
638                do ks=1,nspec
639                  xmass1(j,ks)=xmass1(j,ks)/2.
640                  xmass1(n,ks)=xmass1(j,ks)
641                end do
642              endif
643            endif
644          end do
645          numpart=n
646        endif
647      endif
648    endif
649   
650
651    if (itime.eq.ideltas) exit         ! almost finished
652
653! Compute interval since radioactive decay of deposited mass was computed
654!************************************************************************
655
656    if (itime.lt.loutnext) then
657      ldeltat=itime-(loutnext-loutstep)
658    else                                  ! first half of next interval
659      ldeltat=itime-loutnext
660    endif
661
662
663! Loop over all particles
664!************************
665    if (mp_measure_time) call mpif_mtime('partloop1',0)
666
667!--------------------------------------------------------------------------------
668! various variables for testing reason of CBL scheme, by mc
669    well_mixed_vector=0. !erase vector to test well mixed condition: modified by mc
670    well_mixed_norm=0.   !erase normalization to test well mixed condition: modified by mc
671    avg_ol=0.
672    avg_wst=0.
673    avg_h=0.
674    avg_air_dens=0.  !erase vector to obtain air density at particle positions: modified by mc
675!--------------------------------------------------------------------------------
676
677    do j=1,numpart
678
679! If integration step is due, do it
680!**********************************
681
682      if (itra1(j).eq.itime) then
683
684        if (ioutputforeachrelease.eq.1) then
685          kp=npoint(j)
686        else
687          kp=1
688        endif
689! Determine age class of the particle
690        itage=abs(itra1(j)-itramem(j))
691        do nage=1,nageclass
692          if (itage.lt.lage(nage)) exit
693        end do
694
695! Initialize newly released particle
696!***********************************
697
698        if ((itramem(j).eq.itime).or.(itime.eq.0)) &
699             call initialize(itime,idt(j),uap(j),ucp(j),uzp(j), &
700             us(j),vs(j),ws(j),xtra1(j),ytra1(j),ztra1(j),cbt(j))
701
702! Memorize particle positions
703!****************************
704
705        xold=xtra1(j)
706        yold=ytra1(j)
707        zold=ztra1(j)
708
709! Integrate Lagevin equation for lsynctime seconds
710!*************************************************
711
712        if (mp_measure_time) call mpif_mtime('advance',0)
713
714        call advance(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), &
715             us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
716             cbt(j))
717
718        if (mp_measure_time) call mpif_mtime('advance',1)
719
720
721! Calculate the gross fluxes across layer interfaces
722!***************************************************
723
724        if (iflux.eq.1) call calcfluxes(nage,j,xold,yold,zold)
725
726
727! Determine, when next time step is due
728! If trajectory is terminated, mark it
729!**************************************
730
731        if (nstop.gt.1) then
732          if (linit_cond.ge.1) call initial_cond_calc(itime,j)
733          itra1(j)=-999999999
734        else
735          itra1(j)=itime+lsynctime
736
737
738! Dry deposition and radioactive decay for each species
739! Also check maximum (of all species) of initial mass remaining on the particle;
740! if it is below a threshold value, terminate particle
741!*****************************************************************************
742
743          xmassfract=0.
744          do ks=1,nspec
745            if (decay(ks).gt.0.) then             ! radioactive decay
746              decfact=exp(-real(abs(lsynctime))*decay(ks))
747            else
748              decfact=1.
749            endif
750
751            if (DRYDEPSPEC(ks)) then        ! dry deposition
752              drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact
753              xmass1(j,ks)=xmass1(j,ks)*(1.-prob(ks))*decfact
754              if (decay(ks).gt.0.) then   ! correct for decay (see wetdepo)
755                drydeposit(ks)=drydeposit(ks)* &
756                     exp(real(abs(ldeltat))*decay(ks))
757              endif
758            else                           ! no dry deposition
759              xmass1(j,ks)=xmass1(j,ks)*decfact
760            endif
761
762! Skip check on mass fraction when npoint represents particle number
763            if (mdomainfill.eq.0.and.mquasilag.eq.0) then
764              if (xmass(npoint(j),ks).gt.0.)then
765                   xmassfract=max(xmassfract,real(npart(npoint(j)))* &
766                   xmass1(j,ks)/xmass(npoint(j),ks))
767                   
768                   !CGZ-lifetime: Check mass fraction left/save lifetime
769                   ! if(lroot.and.real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.e_inv.and.checklifetime(j,ks).eq.0.)then
770                       !Mass below 1% of initial >register lifetime
771                   !     checklifetime(j,ks)=abs(itra1(j)-itramem(j))
772
773                   !     species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
774                   !     species_lifetime(ks,2)= species_lifetime(ks,2)+1
775                   ! endif
776                   !CGZ-lifetime: Check mass fraction left/save lifetime
777                   
778              endif
779            else
780              xmassfract=1.0
781            endif
782          end do
783
784          if (xmassfract.lt.minmass) then ! .and. sum(real(npart(npoint(j)))*xmass1(j,:)).lt.1.0) then   ! terminate all particles carrying less mass
785          !            print*,'terminated particle ',j,' for small mass (', sum(real(npart(npoint(j)))* &
786          !         xmass1(j,:)), ' of ', sum(xmass(npoint(j),:)),')'
787            itra1(j)=-999999999
788            if (verbosity.gt.0) then
789              print*,'terminated particle ',j,' for small mass'
790            endif
791          endif
792
793!        Sabine Eckhardt, June 2008
794!        don't create depofield for backward runs
795          if (DRYDEP.AND.(ldirect.eq.1)) then
796            call drydepokernel(nclass(j),drydeposit,real(xtra1(j)), &
797                 real(ytra1(j)),nage,kp)
798            if (nested_output.eq.1) call drydepokernel_nest( &
799                 nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)), &
800                 nage,kp)
801          endif
802
803! Terminate trajectories that are older than maximum allowed age
804!***************************************************************
805
806          if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
807            if (linit_cond.ge.1) &
808                 call initial_cond_calc(itime+lsynctime,j)
809            itra1(j)=-999999999
810            !print*, 'terminated particle ',j,'for age'
811          endif
812        endif
813
814      endif
815
816    end do ! j=1, numpart
817
818    if(mp_measure_time) call mpif_mtime('partloop1',1)
819
820
821! Added by mc: counter of "unstable" particle velocity during a time scale
822!   of maximumtl=20 minutes (defined in com_mod)
823
824    total_nan_intl=0
825    i_nan=i_nan+1 ! added by mc to count nan during a time of maxtl (i.e. maximum tl fixed here to 20 minutes, see com_mod)
826    sum_nan_count(i_nan)=nan_count
827    if (i_nan > maxtl/lsynctime) i_nan=1 !lsynctime must be <= maxtl
828    do ii_nan=1, (maxtl/lsynctime)
829      total_nan_intl=total_nan_intl+sum_nan_count(ii_nan)
830    end do
831
832! Output to keep track of the numerical instabilities in CBL simulation
833! and if they are compromising the final result (or not):
834    if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl 
835
836
837  end do ! itime=0,ideltas,lsynctime
838
839
840! Complete the calculation of initial conditions for particles not yet terminated
841!*****************************************************************************
842
843! eso :TODO: this not implemented yet (transfer particles to PID 0 or rewrite)
844! the tools to do this are already in mpi_mod.f90
845  if (lroot) then
846    do j=1,numpart
847      if (linit_cond.ge.1) call initial_cond_calc(itime,j)
848    end do
849  end if
850
851
852  if (ipout.eq.2) then
853! MPI process 0 creates the file, the other processes append to it
854    do ip=0, mp_partgroup_np-1
855      if (ip.eq.mp_partid) then
856        !if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
857        call partoutput(itime)    ! dump particle positions
858      end if
859      call mpif_mpi_barrier
860    end do
861  end if
862
863! eso :TODO: MPI
864  if (linit_cond.ge.1.and.lroot) call initial_cond_output(itime)   ! dump initial cond. field
865
866
867! De-allocate memory and end
868!***************************
869
870999 if (iflux.eq.1) then
871    deallocate(flux)
872  endif
873  if (OHREA) then
874    deallocate(OH_field,OH_hourly,lonOH,latOH,altOH)
875  endif
876  if (ldirect.gt.0) then
877    deallocate(drygridunc,wetgridunc)
878  endif
879  deallocate(gridunc)
880  deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass)
881  if (allocated(checklifetime)) deallocate(checklifetime)
882  deallocate(ireleasestart,ireleaseend,npart,kindz)
883  deallocate(xmasssave)
884  if (nested_output.eq.1) then
885    deallocate(orooutn, arean, volumen)
886    if (ldirect.gt.0) then
887      deallocate(griduncn,drygriduncn,wetgriduncn)
888    endif
889  endif
890  deallocate(outheight,outheighthalf)
891  deallocate(oroout, area, volume)
892
893  if (mp_measure_time) call mpif_mtime('timemanager',1)
894
895end subroutine timemanager
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG