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

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

Added option to use double precision for calculating the deposition fields

  • Property mode set to 100644
File size: 27.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
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  !   For compatibility with 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  ! outnum             number of samples for each concentration calculation    *
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 netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
100       &concoutput_surf_netcdf,concoutput_surf_nest_netcdf
101
102  implicit none
103
104  integer :: j,ks,kp,l,n,itime,nstop,nstop1
105! integer :: ksp
106  integer :: loutnext,loutstart,loutend
107  integer :: ix,jy,ldeltat,itage,nage
108  integer :: i_nan=0,ii_nan,total_nan_intl=0  !added by mc to check instability in CBL scheme
109  real :: outnum,weight,prob(maxspec),decfact
110  ! real :: uap(maxpart),ucp(maxpart),uzp(maxpart)
111  ! real :: us(maxpart),vs(maxpart),ws(maxpart)
112  ! integer(kind=2) :: cbt(maxpart)
113  real(sp) :: gridtotalunc
114  real(dep_prec) :: drydeposit(maxspec),wetgridtotalunc,drygridtotalunc
115  real :: xold,yold,zold,xmassfract
116  real, parameter :: e_inv = 1.0/exp(1.0)
117  !double precision xm(maxspec,maxpointspec_act),
118  !    +                 xm_depw(maxspec,maxpointspec_act),
119  !    +                 xm_depd(maxspec,maxpointspec_act)
120
121
122  !open(88,file='TEST.dat')
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  !  open(127,file=path(2)(1:length(2))//'depostat.dat'
133  !    +  ,form='unformatted')
134  !write (*,*) 'writing deposition statistics depostat.dat!'
135
136  !**********************************************************************
137  ! Loop over the whole modelling period in time steps of mintime seconds
138  !**********************************************************************
139
140!ZHG 2015
141!CGZ-lifetime: set lifetime to 0
142  ! checklifetime(:,:)=0
143  ! species_lifetime(:,:)=0
144  ! print*, 'Initialized lifetime'
145!CGZ-lifetime: set lifetime to 0
146 
147
148
149  if (verbosity.gt.0) then
150    write (*,*) 'timemanager> starting simulation'
151    if (verbosity.gt.1) then
152      CALL SYSTEM_CLOCK(count_clock)
153      WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
154    endif     
155  endif
156
157  do itime=0,ideltas,lsynctime
158
159  ! Computation of wet deposition, OH reaction and mass transfer
160  ! between two species every lsynctime seconds
161  ! maybe wet depo frequency can be relaxed later but better be on safe side
162  ! wetdepo must be called BEFORE new fields are read in but should not
163  ! be called in the very beginning before any fields are loaded, or
164  ! before particles are in the system
165  ! Code may not be correct for decay of deposition
166  ! changed by Petra Seibert 9/02
167  !********************************************************************
168
169    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
170        if (verbosity.gt.0) then
171           write (*,*) 'timemanager> call wetdepo'
172        endif     
173         call wetdepo(itime,lsynctime,loutnext)
174    endif
175
176    if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
177         call ohreaction(itime,lsynctime,loutnext)
178
179    if (ASSSPEC .and. itime .ne. 0 .and. numpart .gt. 0) then
180       stop 'associated species not yet implemented!'
181  !     call transferspec(itime,lsynctime,loutnext)
182    endif
183
184  ! compute convection for backward runs
185  !*************************************
186
187   if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) then
188        if (verbosity.gt.0) then
189           write (*,*) 'timemanager> call convmix -- backward'
190        endif         
191      call convmix(itime)
192        if (verbosity.gt.1) then
193          !CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
194          CALL SYSTEM_CLOCK(count_clock)
195          WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
196        endif
197   endif
198
199  ! Get necessary wind fields if not available
200  !*******************************************
201    if (verbosity.gt.0) then
202           write (*,*) 'timemanager> call getfields'
203    endif
204    call getfields(itime,nstop1)
205        if (verbosity.gt.1) then
206          CALL SYSTEM_CLOCK(count_clock)
207          WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
208        endif
209    if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
210
211  ! Get hourly OH fields if not available
212  !****************************************************
213    if (OHREA) then
214      if (verbosity.gt.0) then
215             write (*,*) 'timemanager> call gethourlyOH'
216      endif
217      call gethourlyOH(itime)
218          if (verbosity.gt.1) then
219            CALL SYSTEM_CLOCK(count_clock)
220            WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
221          endif
222    endif
223       
224  ! Release particles
225  !******************
226
227    if (verbosity.gt.0) then
228           write (*,*) 'timemanager>  Release particles'
229    endif
230
231    if (mdomainfill.ge.1) then
232      if (itime.eq.0) then
233        if (verbosity.gt.0) then
234          write (*,*) 'timemanager>  call init_domainfill'
235        endif       
236        call init_domainfill
237      else
238        if (verbosity.gt.0) then
239          write (*,*) 'timemanager>  call boundcond_domainfill'
240        endif   
241        call boundcond_domainfill(itime,loutend)
242      endif
243    else
244      if (verbosity.gt.0) then
245        print*,'call releaseparticles' 
246      endif
247      call releaseparticles(itime)
248      if (verbosity.gt.1) then
249        CALL SYSTEM_CLOCK(count_clock)
250        WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
251      endif
252    endif
253
254
255  ! Compute convective mixing for forward runs
256  ! for backward runs it is done before next windfield is read in
257  !**************************************************************
258
259   if ((ldirect.eq.1).and.(lconvection.eq.1)) then
260     if (verbosity.gt.0) then
261       write (*,*) 'timemanager> call convmix -- forward'
262     endif   
263     call convmix(itime)
264   endif
265
266  ! If middle of averaging period of output fields is reached, accumulated
267  ! deposited mass radioactively decays
268  !***********************************************************************
269
270    if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
271      do ks=1,nspec
272      do kp=1,maxpointspec_act
273        if (decay(ks).gt.0.) then
274          do nage=1,nageclass
275            do l=1,nclassunc
276  ! Mother output grid
277              do jy=0,numygrid-1
278                do ix=0,numxgrid-1
279                  wetgridunc(ix,jy,ks,kp,l,nage)= &
280                       wetgridunc(ix,jy,ks,kp,l,nage)* &
281                       exp(-1.*outstep*decay(ks))
282                  drygridunc(ix,jy,ks,kp,l,nage)= &
283                       drygridunc(ix,jy,ks,kp,l,nage)* &
284                       exp(-1.*outstep*decay(ks))
285                end do
286              end do
287  ! Nested output grid
288              if (nested_output.eq.1) then
289                do jy=0,numygridn-1
290                  do ix=0,numxgridn-1
291                    wetgriduncn(ix,jy,ks,kp,l,nage)= &
292                         wetgriduncn(ix,jy,ks,kp,l,nage)* &
293                         exp(-1.*outstep*decay(ks))
294                    drygriduncn(ix,jy,ks,kp,l,nage)= &
295                         drygriduncn(ix,jy,ks,kp,l,nage)* &
296                         exp(-1.*outstep*decay(ks))
297                  end do
298                end do
299              endif
300            end do
301          end do
302        endif
303      end do
304      end do
305    endif
306
307  !!! CHANGE: These lines may be switched on to check the conservation
308  !!! of mass within FLEXPART
309  !   if (itime.eq.loutnext) then
310  !   do 247 ksp=1, nspec
311  !   do 247 kp=1, maxpointspec_act
312  !47         xm(ksp,kp)=0.
313
314  !   do 249 ksp=1, nspec
315  !     do 249 j=1,numpart
316  !          if (ioutputforeachrelease.eq.1) then
317  !            kp=npoint(j)
318  !          else
319  !            kp=1
320  !          endif
321  !       if (itra1(j).eq.itime) then
322  !          xm(ksp,kp)=xm(ksp,kp)+xmass1(j,ksp)
323  !         write(*,*) 'xmass: ',xmass1(j,ksp),j,ksp,nspec
324  !       endif
325  !49     continue
326  !  do 248 ksp=1,nspec
327  !  do 248 kp=1,maxpointspec_act
328  !  xm_depw(ksp,kp)=0.
329  !  xm_depd(ksp,kp)=0.
330  !     do 248 nage=1,nageclass
331  !       do 248 ix=0,numxgrid-1
332  !         do 248 jy=0,numygrid-1
333  !           do 248 l=1,nclassunc
334  !              xm_depw(ksp,kp)=xm_depw(ksp,kp)
335  !    +                  +wetgridunc(ix,jy,ksp,kp,l,nage)
336  !48                 xm_depd(ksp,kp)=xm_depd(ksp,kp)
337  !    +                  +drygridunc(ix,jy,ksp,kp,l,nage)
338  !             do 246 ksp=1,nspec
339  !46                    write(88,'(2i10,3e12.3)')
340  !    +              itime,ksp,(xm(ksp,kp),kp=1,maxpointspec_act),
341  !    +                (xm_depw(ksp,kp),kp=1,maxpointspec_act),
342  !    +                (xm_depd(ksp,kp),kp=1,maxpointspec_act)
343  !  endif
344  !!! CHANGE
345
346
347
348  ! Check whether concentrations are to be calculated
349  !**************************************************
350
351    if ((ldirect*itime.ge.ldirect*loutstart).and. &
352         (ldirect*itime.le.ldirect*loutend)) then ! add to grid
353      if (mod(itime-loutstart,loutsample).eq.0) then
354
355  ! If we are exactly at the start or end of the concentration averaging interval,
356  ! give only half the weight to this sample
357  !*****************************************************************************
358
359        if ((itime.eq.loutstart).or.(itime.eq.loutend)) then
360          weight=0.5
361        else
362          weight=1.0
363        endif
364        outnum=outnum+weight
365        call conccalc(itime,weight)
366      endif
367
368
369      if ((mquasilag.eq.1).and.(itime.eq.(loutstart+loutend)/2)) &
370           call partoutput_short(itime)    ! dump particle positions in extremely compressed format
371
372
373  ! Output and reinitialization of grid
374  ! If necessary, first sample of new grid is also taken
375  !*****************************************************
376
377      if ((itime.eq.loutend).and.(outnum.gt.0.)) then
378        if ((iout.le.3.).or.(iout.eq.5)) then
379          if (surf_only.ne.1) then
380            if (lnetcdfout.eq.1) then
381              call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
382            else
383              call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
384            endif
385          else 
386            if (verbosity.eq.1) then
387             print*,'call concoutput_surf '
388             call system_clock(count_clock)
389             write(*,*) 'system clock',count_clock - count_clock0   
390            endif
391            if (lnetcdfout.eq.1) then
392              call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
393            else
394              call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
395              if (verbosity.eq.1) then
396                print*,'called concoutput_surf '
397                call system_clock(count_clock)
398                write(*,*) 'system clock',count_clock - count_clock0   
399              endif
400            endif
401          endif
402
403          if (nested_output .eq. 1) then
404            if (lnetcdfout.eq.0) then
405              if (surf_only.ne.1) then
406                call concoutput_nest(itime,outnum)
407              else
408                call concoutput_surf_nest(itime,outnum)
409              endif
410            else
411              if (surf_only.ne.1) then
412                call concoutput_nest_netcdf(itime,outnum)
413              else
414                call concoutput_surf_nest_netcdf(itime,outnum)
415              endif
416            endif
417          endif
418          outnum=0.
419        endif
420        if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
421        if (iflux.eq.1) call fluxoutput(itime)
422        write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
423 
424        !CGZ-lifetime: output species lifetime
425!ZHG
426        ! write(*,*) 'Overview species lifetime in days', &
427        !      real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
428        ! write(*,*) 'all info:',species_lifetime
429!ZHG
430        !CGZ-lifetime: output species lifetime
431
432        !write(*,46) float(itime)/3600,itime,numpart
43345      format(i9,' SECONDS SIMULATED: ',i8, ' PARTICLES:    Uncertainty: ',3f7.3)
43446      format(' Simulated ',f7.1,' hours (',i9,' s), ',i8, ' particles')
435        if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
436        loutnext=loutnext+loutstep
437        loutstart=loutnext-loutaver/2
438        loutend=loutnext+loutaver/2
439        if (itime.eq.loutstart) then
440          weight=0.5
441          outnum=outnum+weight
442          call conccalc(itime,weight)
443        endif
444
445
446  ! Check, whether particles are to be split:
447  ! If so, create new particles and attribute all information from the old
448  ! particles also to the new ones; old and new particles both get half the
449  ! mass of the old ones
450  !************************************************************************
451
452        if (ldirect*itime.ge.ldirect*itsplit) then
453          n=numpart
454          do j=1,numpart
455            if (ldirect*itime.ge.ldirect*itrasplit(j)) then
456              if (n.lt.maxpart) then
457                n=n+1
458                itrasplit(j)=2*(itrasplit(j)-itramem(j))+itramem(j)
459                itrasplit(n)=itrasplit(j)
460                itramem(n)=itramem(j)
461                itra1(n)=itra1(j)
462                idt(n)=idt(j)
463                npoint(n)=npoint(j)
464                nclass(n)=nclass(j)
465                xtra1(n)=xtra1(j)
466                ytra1(n)=ytra1(j)
467                ztra1(n)=ztra1(j)
468                uap(n)=uap(j)
469                ucp(n)=ucp(j)
470                uzp(n)=uzp(j)
471                us(n)=us(j)
472                vs(n)=vs(j)
473                ws(n)=ws(j)
474                cbt(n)=cbt(j)
475                do ks=1,nspec
476                  xmass1(j,ks)=xmass1(j,ks)/2.
477                  xmass1(n,ks)=xmass1(j,ks)
478                end do
479              endif
480            endif
481          end do
482          numpart=n
483        endif
484      endif
485    endif
486
487
488    if (itime.eq.ideltas) exit         ! almost finished
489
490  ! Compute interval since radioactive decay of deposited mass was computed
491  !************************************************************************
492
493    if (itime.lt.loutnext) then
494      ldeltat=itime-(loutnext-loutstep)
495    else                                  ! first half of next interval
496      ldeltat=itime-loutnext
497    endif
498
499
500  ! Loop over all particles
501  !************************
502  ! Various variables for testing reason of CBL scheme, by mc
503    well_mixed_vector=0. !erase vector to test well mixed condition: modified by mc
504    well_mixed_norm=0.   !erase normalization to test well mixed condition: modified by mc
505    avg_ol=0.
506    avg_wst=0.
507    avg_h=0.
508    avg_air_dens=0.  !erase vector to obtain air density at particle positions: modified by mc
509  !-----------------------------------------------------------------------------
510    do j=1,numpart
511
512
513  ! If integration step is due, do it
514  !**********************************
515
516      if (itra1(j).eq.itime) then
517
518        if (ioutputforeachrelease.eq.1) then
519            kp=npoint(j)
520        else
521            kp=1
522        endif
523  ! Determine age class of the particle
524        itage=abs(itra1(j)-itramem(j))
525        do nage=1,nageclass
526          if (itage.lt.lage(nage)) exit
527        end do
528
529  ! Initialize newly released particle
530  !***********************************
531
532        if ((itramem(j).eq.itime).or.(itime.eq.0)) &
533             call initialize(itime,idt(j),uap(j),ucp(j),uzp(j), &
534             us(j),vs(j),ws(j),xtra1(j),ytra1(j),ztra1(j),cbt(j))
535
536  ! Memorize particle positions
537  !****************************
538
539        xold=xtra1(j)
540        yold=ytra1(j)
541        zold=ztra1(j)
542
543  ! Integrate Lagevin equation for lsynctime seconds
544  !*************************************************
545
546        call advance(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), &
547             us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
548             cbt(j))
549
550  ! Calculate the gross fluxes across layer interfaces
551  !***************************************************
552
553        if (iflux.eq.1) call calcfluxes(nage,j,xold,yold,zold)
554
555
556  ! Determine, when next time step is due
557  ! If trajectory is terminated, mark it
558  !**************************************
559
560        if (nstop.gt.1) then
561          if (linit_cond.ge.1) call initial_cond_calc(itime,j)
562          itra1(j)=-999999999
563        else
564          itra1(j)=itime+lsynctime
565
566
567  ! Dry deposition and radioactive decay for each species
568  ! Also check maximum (of all species) of initial mass remaining on the particle;
569  ! if it is below a threshold value, terminate particle
570  !*****************************************************************************
571
572          xmassfract=0.
573          do ks=1,nspec
574            if (decay(ks).gt.0.) then             ! radioactive decay
575              decfact=exp(-real(abs(lsynctime))*decay(ks))
576            else
577              decfact=1.
578            endif
579
580            if (DRYDEPSPEC(ks)) then        ! dry deposition
581              drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact
582              xmass1(j,ks)=xmass1(j,ks)*(1.-prob(ks))*decfact
583              if (decay(ks).gt.0.) then   ! correct for decay (see wetdepo)
584                drydeposit(ks)=drydeposit(ks)* &
585                     exp(real(abs(ldeltat))*decay(ks))
586              endif
587            else                           ! no dry deposition
588              xmass1(j,ks)=xmass1(j,ks)*decfact
589            endif
590
591
592            if (mdomainfill.eq.0) then
593              if (xmass(npoint(j),ks).gt.0.) &
594                   xmassfract=max(xmassfract,real(npart(npoint(j)))* &
595                   xmass1(j,ks)/xmass(npoint(j),ks))
596!ZHG 2015
597                  !CGZ-lifetime: Check mass fraction left/save lifetime
598                   ! if(real(npart(npoint(j)))*xmass1(j,ks)/xmass(npoint(j),ks).lt.e_inv.and.checklifetime(j,ks).eq.0.)then
599                       !Mass below 1% of initial >register lifetime
600                       ! checklifetime(j,ks)=abs(itra1(j)-itramem(j))
601                       ! species_lifetime(ks,1)=species_lifetime(ks,1)+abs(itra1(j)-itramem(j))
602                       ! species_lifetime(ks,2)= species_lifetime(ks,2)+1
603                   ! endif
604                   !CGZ-lifetime: Check mass fraction left/save lifetime
605!ZHG 2015
606            else
607              xmassfract=1.
608            endif
609          end do
610
611          if (xmassfract.lt.minmass) then   ! terminate all particles carrying less mass
612            itra1(j)=-999999999
613            if (verbosity.gt.0) then
614              print*,'terminated particle ',j,' for small mass'
615            endif
616          endif
617
618  !        Sabine Eckhardt, June 2008
619  !        don't create depofield for backward runs
620          if (DRYDEP.AND.(ldirect.eq.1)) then
621            call drydepokernel(nclass(j),drydeposit,real(xtra1(j)), &
622                 real(ytra1(j)),nage,kp)
623            if (nested_output.eq.1) call drydepokernel_nest( &
624                 nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)), &
625                 nage,kp)
626          endif
627
628  ! Terminate trajectories that are older than maximum allowed age
629  !***************************************************************
630
631          if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
632            if (linit_cond.ge.1) call initial_cond_calc(itime+lsynctime,j)
633            itra1(j)=-999999999
634            if (verbosity.gt.0) then
635              print*,'terminated particle ',j,' for age'
636            endif
637          endif
638        endif
639
640      endif
641
642    end do !loop over particles
643   
644  ! Counter of "unstable" particle velocity during a time scale of
645  ! maximumtl=20 minutes (defined in com_mod)
646  !***************************************************************
647   
648    total_nan_intl=0
649    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)
650    sum_nan_count(i_nan)=nan_count
651    if (i_nan > maxtl/lsynctime) i_nan=1 !lsynctime must be <= maxtl
652    do ii_nan=1, (maxtl/lsynctime)
653      total_nan_intl=total_nan_intl+sum_nan_count(ii_nan)
654    end do
655  ! Output to keep track of the numerical instabilities in CBL simulation and if
656  ! they are compromising the final result (or not)
657    if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl 
658         
659  end do
660
661
662  ! Complete the calculation of initial conditions for particles not yet terminated
663  !*****************************************************************************
664
665  do j=1,numpart
666    if (linit_cond.ge.1) call initial_cond_calc(itime,j)
667  end do
668
669  if (ipout.eq.2) call partoutput(itime)     ! dump particle positions
670
671  if (linit_cond.ge.1) call initial_cond_output(itime)   ! dump initial cond. field
672
673  !close(104)
674
675  ! De-allocate memory and end
676  !***************************
677
678  if (iflux.eq.1) then
679      deallocate(flux)
680  endif
681  if (OHREA) then
682      deallocate(OH_field,OH_hourly,lonOH,latOH,altOH)
683  endif
684  if (ldirect.gt.0) then
685  deallocate(drygridunc,wetgridunc)
686  endif
687  deallocate(gridunc)
688  deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass)
689  deallocate(ireleasestart,ireleaseend,npart,kindz)
690  deallocate(xmasssave)
691  if (nested_output.eq.1) then
692     deallocate(orooutn, arean, volumen)
693     if (ldirect.gt.0) then
694     deallocate(griduncn,drygriduncn,wetgriduncn)
695     endif
696  endif
697  deallocate(outheight,outheighthalf)
698  deallocate(oroout, area, volume)
699
700end subroutine timemanager
701
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG