source: trunk/src/timemanager.f90 @ 27

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