source: flexpart.git/src/timemanager_mpi.f90 @ 7e52e2e

devrelease-10univie
Last change on this file since 7e52e2e was 7e52e2e, checked in by Espen Sollum ATMOS <eso@…>, 4 years ago

Fixed issue with mpi backward runs

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