source: flexpart.git/src/timemanager_mpi.f90 @ 38b7917

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 38b7917 was 38b7917, checked in by Espen Sollum ATMOS <eso@…>, 8 years ago

Parallelization of domain fill option (save/restart not implemented yet)

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