source: trunk/src/timemanager.f90 @ 20

Last change on this file since 20 was 20, checked in by igpis, 10 years ago

move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others

File size: 22.9 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  !*****************************************************************************
48  !                                                                            *
49  ! Variables:                                                                 *
50  ! DEP                .true. if either wet or dry deposition is switched on   *
51  ! decay(maxspec) [1/s] decay constant for radioactive decay                  *
52  ! DRYDEP             .true. if dry deposition is switched on                 *
53  ! ideltas [s]        modelling period                                        *
54  ! itime [s]          actual temporal position of calculation                 *
55  ! ldeltat [s]        time since computation of radioact. decay of depositions*
56  ! loutaver [s]       averaging period for concentration calculations         *
57  ! loutend [s]        end of averaging for concentration calculations         *
58  ! loutnext [s]       next time at which output fields shall be centered      *
59  ! loutsample [s]     sampling interval for averaging of concentrations       *
60  ! loutstart [s]      start of averaging for concentration calculations       *
61  ! loutstep [s]       time interval for which concentrations shall be         *
62  !                    calculated                                              *
63  ! npoint(maxpart)    index, which starting point the trajectory has          *
64  !                    starting positions of trajectories                      *
65  ! nstop              serves as indicator for fate of particles               *
66  !                    in the particle loop                                    *
67  ! nstop1             serves as indicator for wind fields (see getfields)     *
68  ! outnum             number of samples for each concentration calculation    *
69  ! outnum             number of samples for each concentration calculation    *
70  ! prob               probability of absorption at ground due to dry          *
71  !                    deposition                                              *
72  ! WETDEP             .true. if wet deposition is switched on                 *
73  ! weight             weight for each concentration sample (1/2 or 1)         *
74  ! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to          *
75  !                    turbulence                                              *
76  ! us(maxpart),vs(maxpart),ws(maxpart) = random velocities due to inter-      *
77  !                    polation                                                *
78  ! xtra1(maxpart), ytra1(maxpart), ztra1(maxpart) =                           *
79  !                    spatial positions of trajectories                       *
80  !                                                                            *
81  ! Constants:                                                                 *
82  ! maxpart            maximum number of trajectories                          *
83  !                                                                            *
84  !*****************************************************************************
85
86  use unc_mod
87  use point_mod
88  use xmass_mod
89  use flux_mod
90  use outg_mod
91  use oh_mod
92  use par_mod
93  use com_mod
94
95  implicit none
96
97  integer :: j,ks,kp,l,n,itime,nstop,nstop1
98! integer :: ksp
99  integer :: loutnext,loutstart,loutend
100  integer :: ix,jy,ldeltat,itage,nage
101  real :: outnum,weight,prob(maxspec)
102  real :: uap(maxpart),ucp(maxpart),uzp(maxpart),decfact
103  real :: us(maxpart),vs(maxpart),ws(maxpart)
104  integer(kind=2) :: cbt(maxpart)
105  real :: drydeposit(maxspec),gridtotalunc,wetgridtotalunc
106  real :: drygridtotalunc,xold,yold,zold,xmassfract
107  !double precision xm(maxspec,maxpointspec_act),
108  !    +                 xm_depw(maxspec,maxpointspec_act),
109  !    +                 xm_depd(maxspec,maxpointspec_act)
110
111
112  !open(88,file='TEST.dat')
113
114  ! First output for time 0
115  !************************
116
117  loutnext=loutstep/2
118  outnum=0.
119  loutstart=loutnext-loutaver/2
120  loutend=loutnext+loutaver/2
121
122  !  open(127,file=path(2)(1:length(2))//'depostat.dat'
123  !    +  ,form='unformatted')
124  !write (*,*) 'writing deposition statistics depostat.dat!'
125
126  !**********************************************************************
127  ! Loop over the whole modelling period in time steps of mintime seconds
128  !**********************************************************************
129
130
131  if (verbosity.gt.0) then
132    write (*,*) 'timemanager> starting simulation'
133    if (verbosity.gt.1) then
134      CALL SYSTEM_CLOCK(count_clock)
135      WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
136    endif     
137  endif
138
139  do itime=0,ideltas,lsynctime
140
141
142  ! Computation of wet deposition, OH reaction and mass transfer
143  ! between two species every lsynctime seconds
144  ! maybe wet depo frequency can be relaxed later but better be on safe side
145  ! wetdepo must be called BEFORE new fields are read in but should not
146  ! be called in the very beginning before any fields are loaded, or
147  ! before particles are in the system
148  ! Code may not be correct for decay of deposition
149  ! changed by Petra Seibert 9/02
150  !********************************************************************
151
152    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
153        if (verbosity.gt.0) then
154           write (*,*) 'timemanager> call wetdepo'
155        endif     
156         call wetdepo(itime,lsynctime,loutnext)
157    endif
158
159    if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
160         call ohreaction(itime,lsynctime,loutnext)
161
162    if (ASSSPEC .and. itime .ne. 0 .and. numpart .gt. 0) then
163       stop 'associated species not yet implemented!'
164  !     call transferspec(itime,lsynctime,loutnext)
165    endif
166
167  ! compute convection for backward runs
168  !*************************************
169
170   if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) then
171        if (verbosity.gt.0) then
172           write (*,*) 'timemanager> call convmix -- backward'
173        endif         
174      call convmix(itime)
175        if (verbosity.gt.1) then
176          !CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
177          CALL SYSTEM_CLOCK(count_clock)
178          WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
179        endif
180   endif
181
182  ! Get necessary wind fields if not available
183  !*******************************************
184    if (verbosity.gt.0) then
185           write (*,*) 'timemanager> call getfields'
186    endif
187    call getfields(itime,nstop1)
188        if (verbosity.gt.1) then
189          CALL SYSTEM_CLOCK(count_clock)
190          WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
191        endif
192    if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
193
194  ! Release particles
195  !******************
196
197    if (verbosity.gt.0) then
198           write (*,*) 'timemanager>  Release particles'
199    endif
200
201    if (mdomainfill.ge.1) then
202      if (itime.eq.0) then
203        if (verbosity.gt.0) then
204          write (*,*) 'timemanager>  call init_domainfill'
205        endif       
206        call init_domainfill
207      else
208        if (verbosity.gt.0) then
209          write (*,*) 'timemanager>  call boundcond_domainfill'
210        endif   
211        call boundcond_domainfill(itime,loutend)
212      endif
213    else
214      if (verbosity.gt.0) then
215        print*,'call releaseparticles' 
216      endif
217      call releaseparticles(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
225  ! Compute convective mixing for forward runs
226  ! for backward runs it is done before next windfield is read in
227  !**************************************************************
228
229   if ((ldirect.eq.1).and.(lconvection.eq.1)) then
230     if (verbosity.gt.0) then
231       write (*,*) 'timemanager> call convmix -- forward'
232     endif   
233     call convmix(itime)
234   endif
235
236  ! If middle of averaging period of output fields is reached, accumulated
237  ! deposited mass radioactively decays
238  !***********************************************************************
239
240    if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
241      do ks=1,nspec
242      do kp=1,maxpointspec_act
243        if (decay(ks).gt.0.) then
244          do nage=1,nageclass
245            do l=1,nclassunc
246  ! Mother output grid
247              do jy=0,numygrid-1
248                do ix=0,numxgrid-1
249                  wetgridunc(ix,jy,ks,kp,l,nage)= &
250                       wetgridunc(ix,jy,ks,kp,l,nage)* &
251                       exp(-1.*outstep*decay(ks))
252                  drygridunc(ix,jy,ks,kp,l,nage)= &
253                       drygridunc(ix,jy,ks,kp,l,nage)* &
254                       exp(-1.*outstep*decay(ks))
255                end do
256              end do
257  ! Nested output grid
258              if (nested_output.eq.1) then
259                do jy=0,numygridn-1
260                  do ix=0,numxgridn-1
261                    wetgriduncn(ix,jy,ks,kp,l,nage)= &
262                         wetgriduncn(ix,jy,ks,kp,l,nage)* &
263                         exp(-1.*outstep*decay(ks))
264                    drygriduncn(ix,jy,ks,kp,l,nage)= &
265                         drygriduncn(ix,jy,ks,kp,l,nage)* &
266                         exp(-1.*outstep*decay(ks))
267                  end do
268                end do
269              endif
270            end do
271          end do
272        endif
273      end do
274      end do
275    endif
276
277  !!! CHANGE: These lines may be switched on to check the conservation
278  !!! of mass within FLEXPART
279  !   if (itime.eq.loutnext) then
280  !   do 247 ksp=1, nspec
281  !   do 247 kp=1, maxpointspec_act
282  !47         xm(ksp,kp)=0.
283
284  !   do 249 ksp=1, nspec
285  !     do 249 j=1,numpart
286  !          if (ioutputforeachrelease.eq.1) then
287  !            kp=npoint(j)
288  !          else
289  !            kp=1
290  !          endif
291  !       if (itra1(j).eq.itime) then
292  !          xm(ksp,kp)=xm(ksp,kp)+xmass1(j,ksp)
293  !         write(*,*) 'xmass: ',xmass1(j,ksp),j,ksp,nspec
294  !       endif
295  !49     continue
296  !  do 248 ksp=1,nspec
297  !  do 248 kp=1,maxpointspec_act
298  !  xm_depw(ksp,kp)=0.
299  !  xm_depd(ksp,kp)=0.
300  !     do 248 nage=1,nageclass
301  !       do 248 ix=0,numxgrid-1
302  !         do 248 jy=0,numygrid-1
303  !           do 248 l=1,nclassunc
304  !              xm_depw(ksp,kp)=xm_depw(ksp,kp)
305  !    +                  +wetgridunc(ix,jy,ksp,kp,l,nage)
306  !48                 xm_depd(ksp,kp)=xm_depd(ksp,kp)
307  !    +                  +drygridunc(ix,jy,ksp,kp,l,nage)
308  !             do 246 ksp=1,nspec
309  !46                    write(88,'(2i10,3e12.3)')
310  !    +              itime,ksp,(xm(ksp,kp),kp=1,maxpointspec_act),
311  !    +                (xm_depw(ksp,kp),kp=1,maxpointspec_act),
312  !    +                (xm_depd(ksp,kp),kp=1,maxpointspec_act)
313  !  endif
314  !!! CHANGE
315
316
317
318  ! Check whether concentrations are to be calculated
319  !**************************************************
320
321    if ((ldirect*itime.ge.ldirect*loutstart).and. &
322         (ldirect*itime.le.ldirect*loutend)) then ! add to grid
323      if (mod(itime-loutstart,loutsample).eq.0) then
324
325  ! If we are exactly at the start or end of the concentration averaging interval,
326  ! give only half the weight to this sample
327  !*****************************************************************************
328
329        if ((itime.eq.loutstart).or.(itime.eq.loutend)) then
330          weight=0.5
331        else
332          weight=1.0
333        endif
334        outnum=outnum+weight
335        call conccalc(itime,weight)
336      endif
337
338
339      if ((mquasilag.eq.1).and.(itime.eq.(loutstart+loutend)/2)) &
340           call partoutput_short(itime)    ! dump particle positions in extremely compressed format
341
342
343  ! Output and reinitialization of grid
344  ! If necessary, first sample of new grid is also taken
345  !*****************************************************
346
347      if ((itime.eq.loutend).and.(outnum.gt.0.)) then
348        if ((iout.le.3.).or.(iout.eq.5)) then
349          if (surf_only.ne.1) then
350          call concoutput(itime,outnum,gridtotalunc, &
351               wetgridtotalunc,drygridtotalunc)
352          else 
353  if (verbosity.eq.1) then
354     print*,'call concoutput_surf '
355     CALL SYSTEM_CLOCK(count_clock)
356     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
357  endif
358          call concoutput_surf(itime,outnum,gridtotalunc, &
359               wetgridtotalunc,drygridtotalunc)
360  if (verbosity.eq.1) then
361     print*,'called concoutput_surf '
362     CALL SYSTEM_CLOCK(count_clock)
363     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
364  endif
365          endif
366
367          if ((nested_output.eq.1).and.(surf_only.ne.1)) call concoutput_nest(itime,outnum)
368          if ((nested_output.eq.1).and.(surf_only.eq.1)) call concoutput_surf_nest(itime,outnum)
369          outnum=0.
370        endif
371        if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
372        if (iflux.eq.1) call fluxoutput(itime)
373        write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc, &
374             drygridtotalunc
37545      format(i9,' SECONDS SIMULATED: ',i8, &
376             ' PARTICLES:    Uncertainty: ',3f7.3)
377        if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
378        loutnext=loutnext+loutstep
379        loutstart=loutnext-loutaver/2
380        loutend=loutnext+loutaver/2
381        if (itime.eq.loutstart) then
382          weight=0.5
383          outnum=outnum+weight
384          call conccalc(itime,weight)
385        endif
386
387
388  ! Check, whether particles are to be split:
389  ! If so, create new particles and attribute all information from the old
390  ! particles also to the new ones; old and new particles both get half the
391  ! mass of the old ones
392  !************************************************************************
393
394        if (ldirect*itime.ge.ldirect*itsplit) then
395          n=numpart
396          do j=1,numpart
397            if (ldirect*itime.ge.ldirect*itrasplit(j)) then
398              if (n.lt.maxpart) then
399                n=n+1
400                itrasplit(j)=2*(itrasplit(j)-itramem(j))+itramem(j)
401                itrasplit(n)=itrasplit(j)
402                itramem(n)=itramem(j)
403                itra1(n)=itra1(j)
404                idt(n)=idt(j)
405                npoint(n)=npoint(j)
406                nclass(n)=nclass(j)
407                xtra1(n)=xtra1(j)
408                ytra1(n)=ytra1(j)
409                ztra1(n)=ztra1(j)
410                uap(n)=uap(j)
411                ucp(n)=ucp(j)
412                uzp(n)=uzp(j)
413                us(n)=us(j)
414                vs(n)=vs(j)
415                ws(n)=ws(j)
416                cbt(n)=cbt(j)
417                do ks=1,nspec
418                  xmass1(j,ks)=xmass1(j,ks)/2.
419                  xmass1(n,ks)=xmass1(j,ks)
420                end do
421              endif
422            endif
423          end do
424          numpart=n
425        endif
426      endif
427    endif
428
429
430    if (itime.eq.ideltas) exit         ! almost finished
431
432  ! Compute interval since radioactive decay of deposited mass was computed
433  !************************************************************************
434
435    if (itime.lt.loutnext) then
436      ldeltat=itime-(loutnext-loutstep)
437    else                                  ! first half of next interval
438      ldeltat=itime-loutnext
439    endif
440
441
442  ! Loop over all particles
443  !************************
444
445    do j=1,numpart
446
447
448  ! If integration step is due, do it
449  !**********************************
450
451      if (itra1(j).eq.itime) then
452
453        if (ioutputforeachrelease.eq.1) then
454            kp=npoint(j)
455        else
456            kp=1
457        endif
458  ! Determine age class of the particle
459        itage=abs(itra1(j)-itramem(j))
460        do nage=1,nageclass
461          if (itage.lt.lage(nage)) exit
462        end do
463
464  ! Initialize newly released particle
465  !***********************************
466
467        if ((itramem(j).eq.itime).or.(itime.eq.0)) &
468             call initialize(itime,idt(j),uap(j),ucp(j),uzp(j), &
469             us(j),vs(j),ws(j),xtra1(j),ytra1(j),ztra1(j),cbt(j))
470
471  ! Memorize particle positions
472  !****************************
473
474        xold=xtra1(j)
475        yold=ytra1(j)
476        zold=ztra1(j)
477
478  ! Integrate Lagevin equation for lsynctime seconds
479  !*************************************************
480
481        call advance(itime,npoint(j),idt(j),uap(j),ucp(j),uzp(j), &
482             us(j),vs(j),ws(j),nstop,xtra1(j),ytra1(j),ztra1(j),prob, &
483             cbt(j))
484
485  ! Calculate the gross fluxes across layer interfaces
486  !***************************************************
487
488        if (iflux.eq.1) call calcfluxes(nage,j,xold,yold,zold)
489
490
491  ! Determine, when next time step is due
492  ! If trajectory is terminated, mark it
493  !**************************************
494
495        if (nstop.gt.1) then
496          if (linit_cond.ge.1) call initial_cond_calc(itime,j)
497          itra1(j)=-999999999
498        else
499          itra1(j)=itime+lsynctime
500
501
502  ! Dry deposition and radioactive decay for each species
503  ! Also check maximum (of all species) of initial mass remaining on the particle;
504  ! if it is below a threshold value, terminate particle
505  !*****************************************************************************
506
507          xmassfract=0.
508          do ks=1,nspec
509            if (decay(ks).gt.0.) then             ! radioactive decay
510              decfact=exp(-real(abs(lsynctime))*decay(ks))
511            else
512              decfact=1.
513            endif
514
515            if (DRYDEPSPEC(ks)) then        ! dry deposition
516              drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact
517              xmass1(j,ks)=xmass1(j,ks)*(1.-prob(ks))*decfact
518              if (decay(ks).gt.0.) then   ! correct for decay (see wetdepo)
519                drydeposit(ks)=drydeposit(ks)* &
520                     exp(real(abs(ldeltat))*decay(ks))
521              endif
522            else                           ! no dry deposition
523              xmass1(j,ks)=xmass1(j,ks)*decfact
524            endif
525
526
527            if (mdomainfill.eq.0) then
528              if (xmass(npoint(j),ks).gt.0.) &
529                   xmassfract=max(xmassfract,real(npart(npoint(j)))* &
530                   xmass1(j,ks)/xmass(npoint(j),ks))
531            else
532              xmassfract=1.
533            endif
534          end do
535
536          if (xmassfract.lt.0.0001) then   ! terminate all particles carrying less mass
537            itra1(j)=-999999999
538          endif
539
540  !        Sabine Eckhardt, June 2008
541  !        don't create depofield for backward runs
542          if (DRYDEP.AND.(ldirect.eq.1)) then
543            call drydepokernel(nclass(j),drydeposit,real(xtra1(j)), &
544                 real(ytra1(j)),nage,kp)
545            if (nested_output.eq.1) call drydepokernel_nest( &
546                 nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)), &
547                 nage,kp)
548          endif
549
550  ! Terminate trajectories that are older than maximum allowed age
551  !***************************************************************
552
553          if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
554            if (linit_cond.ge.1) &
555                 call initial_cond_calc(itime+lsynctime,j)
556            itra1(j)=-999999999
557          endif
558        endif
559
560      endif
561
562    end do
563
564  end do
565
566
567  ! Complete the calculation of initial conditions for particles not yet terminated
568  !*****************************************************************************
569
570  do j=1,numpart
571    if (linit_cond.ge.1) call initial_cond_calc(itime,j)
572  end do
573
574  if (ipout.eq.2) call partoutput(itime)     ! dump particle positions
575
576  if (linit_cond.ge.1) call initial_cond_output(itime)   ! dump initial cond. field
577
578  close(104)
579
580  ! De-allocate memory and end
581  !***************************
582
583  if (iflux.eq.1) then
584      deallocate(flux)
585  endif
586  if (OHREA.eqv..TRUE.) then
587      deallocate(OH_field,OH_field_height)
588  endif
589  if (ldirect.gt.0) then
590  deallocate(drygridunc,wetgridunc)
591  endif
592  deallocate(gridunc)
593  deallocate(xpoint1,xpoint2,ypoint1,ypoint2,zpoint1,zpoint2,xmass)
594  deallocate(ireleasestart,ireleaseend,npart,kindz)
595  deallocate(xmasssave)
596  if (nested_output.eq.1) then
597     deallocate(orooutn, arean, volumen)
598     if (ldirect.gt.0) then
599     deallocate(griduncn,drygriduncn,wetgriduncn)
600     endif
601  endif
602  deallocate(outheight,outheighthalf)
603  deallocate(oroout, area, volume)
604
605end subroutine timemanager
606
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG