source: trunk/src/advance.f90 @ 29

Last change on this file since 29 was 29, checked in by hasod, 7 years ago

BUGFIX: implemented bugfixes related to global domain-filling runs with
escaping particles at the poles and above the highest model level
(Tickets # 54 and 55)

File size: 29.8 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 advance(itime,nrelpoint,ldt,up,vp,wp, &
23       usigold,vsigold,wsigold,nstop,xt,yt,zt,prob,icbt)
24  !                     i    i  i/oi/oi/o
25  !  i/o     i/o     i/o     o  i/oi/oi/o i/o  i/o
26  !*****************************************************************************
27  !                                                                            *
28  !  Calculation of turbulent particle trajectories utilizing a                *
29  !  zero-acceleration scheme, which is corrected by a numerically more        *
30  !  accurate Petterssen scheme whenever possible.                             *
31  !                                                                            *
32  !  Particle positions are read in, incremented, and returned to the calling  *
33  !  program.                                                                  *
34  !                                                                            *
35  !  In different regions of the atmosphere (PBL vs. free troposphere),        *
36  !  different parameters are needed for advection, parameterizing turbulent   *
37  !  velocities, etc. For efficiency, different interpolation routines have    *
38  !  been written for these different cases, with the disadvantage that there  *
39  !  exist several routines doing almost the same. They all share the          *
40  !  included file 'interpol_mod'. The following                               *
41  !  interpolation routines are used:                                          *
42  !                                                                            *
43  !  interpol_all(_nests)     interpolates everything (called inside the PBL)  *
44  !  interpol_misslev(_nests) if a particle moves vertically in the PBL,       *
45  !                           additional parameters are interpolated if it     *
46  !                           crosses a model level                            *
47  !  interpol_wind(_nests)    interpolates the wind and determines the         *
48  !                           standard deviation of the wind (called outside   *
49  !                           PBL) also interpolates potential vorticity       *
50  !  interpol_wind_short(_nests) only interpolates the wind (needed for the    *
51  !                           Petterssen scheme)                               *
52  !  interpol_vdep(_nests)    interpolates deposition velocities               *
53  !                                                                            *
54  !                                                                            *
55  !     Author: A. Stohl                                                       *
56  !                                                                            *
57  !     16 December 1997                                                       *
58  !                                                                            *
59  !  Changes:                                                                  *
60  !                                                                            *
61  !  8 April 2000: Deep convection parameterization                            *
62  !                                                                            *
63  !  May 2002: Petterssen scheme introduced                                    *
64  !                                                                            *
65  !*****************************************************************************
66  !                                                                            *
67  ! Variables:                                                                 *
68  ! icbt               1 if particle not transferred to forbidden state,       *
69  !                    else -1                                                 *
70  ! dawsave            accumulated displacement in along-wind direction        *
71  ! dcwsave            accumulated displacement in cross-wind direction        *
72  ! dxsave             accumulated displacement in longitude                   *
73  ! dysave             accumulated displacement in latitude                    *
74  ! h [m]              Mixing height                                           *
75  ! lwindinterv [s]    time interval between two wind fields                   *
76  ! itime [s]          time at which this subroutine is entered                *
77  ! itimec [s]         actual time, which is incremented in this subroutine    *
78  ! href [m]           height for which dry deposition velocity is calculated  *
79  ! ladvance [s]       Total integration time period                           *
80  ! ldirect            1 forward, -1 backward                                  *
81  ! ldt [s]            Time step for the next integration                      *
82  ! lsynctime [s]      Synchronisation interval of FLEXPART                    *
83  ! ngrid              index which grid is to be used                          *
84  ! nrand              index for a variable to be picked from rannumb          *
85  ! nstop              if > 1 particle has left domain and must be stopped     *
86  ! prob               probability of absorption due to dry deposition         *
87  ! rannumb(maxrand)   normally distributed random variables                   *
88  ! rhoa               air density                                             *
89  ! rhograd            vertical gradient of the air density                    *
90  ! up,vp,wp           random velocities due to turbulence (along wind, cross  *
91  !                    wind, vertical wind                                     *
92  ! usig,vsig,wsig     mesoscale wind fluctuations                             *
93  ! usigold,vsigold,wsigold  like usig, etc., but for the last time step       *
94  ! vdepo              Deposition velocities for all species                   *
95  ! xt,yt,zt           Particle position                                       *
96  !                                                                            *
97  !*****************************************************************************
98
99  use point_mod
100  use par_mod
101  use com_mod
102  use interpol_mod
103  use hanna_mod
104  use cmapf_mod
105
106  implicit none
107
108  real(kind=dp) :: xt,yt
109  real :: zt,xts,yts,weight
110  integer :: itime,itimec,nstop,ldt,i,j,k,nrand,loop,memindnext
111  integer :: ngr,nix,njy,ks,nsp,nrelpoint
112  real :: dz,dz1,dz2,xlon,ylat,xpol,ypol,gridsize
113  real :: ru,rv,rw,dt,ux,vy,cosfact,xtn,ytn,tropop
114  real :: prob(maxspec),up,vp,wp,dxsave,dysave,dawsave
115  real :: dcwsave
116  real :: usigold,vsigold,wsigold,r,rs
117  real :: uold,vold,wold,vdepo(maxspec)
118  !real uprof(nzmax),vprof(nzmax),wprof(nzmax)
119  !real usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax)
120  !real rhoprof(nzmax),rhogradprof(nzmax)
121  real :: rhoa,rhograd,ran3,delz,dtf,rhoaux,dtftlw,uxscale,wpscale
122  integer(kind=2) :: icbt
123  real,parameter :: eps=nxmax/3.e5,eps2=1.e-9
124
125
126  !!! CHANGE: TEST OF THE WELL-MIXED CRITERION
127  !  integer,parameter :: iclass=10
128  !  real(kind=dp) :: zacc,tacc,t(iclass),th(0:iclass),hsave
129  !  logical dump
130  !  save zacc,tacc,t,th,hsave,dump
131  !!! CHANGE
132
133  integer :: idummy = -7
134  real    :: settling = 0.
135
136
137  !!! CHANGE: TEST OF THE WELL-MIXED CRITERION
138  !if (idummy.eq.-7) then
139  !open(550,file='WELLMIXEDTEST')
140  !do 17 i=0,iclass
141  !7      th(i)=real(i)/real(iclass)
142  !endif
143  !!! CHANGE
144
145
146  nstop=0
147  do i=1,nmixz
148    indzindicator(i)=.true.
149  end do
150
151
152  if (DRYDEP) then    ! reset probability for deposition
153    do ks=1,nspec
154      depoindicator(ks)=.true.
155      prob(ks)=0.
156    end do
157  endif
158
159  dxsave=0.           ! reset position displacements
160  dysave=0.           ! due to mean wind
161  dawsave=0.          ! and turbulent wind
162  dcwsave=0.
163
164  itimec=itime
165
166  nrand=int(ran3(idummy)*real(maxrand-1))+1
167
168
169  ! Determine whether lat/long grid or polarstereographic projection
170  ! is to be used
171  ! Furthermore, determine which nesting level to be used
172  !*****************************************************************
173
174  if (nglobal.and.(yt.gt.switchnorthg)) then
175    ngrid=-1
176  else if (sglobal.and.(yt.lt.switchsouthg)) then
177    ngrid=-2
178  else
179    ngrid=0
180    do j=numbnests,1,-1
181      if ((xt.gt.xln(j)+eps).and.(xt.lt.xrn(j)-eps).and. &
182           (yt.gt.yln(j)+eps).and.(yt.lt.yrn(j)-eps)) then
183        ngrid=j
184        goto 23
185      endif
186    end do
18723   continue
188  endif
189
190
191  !***************************
192  ! Interpolate necessary data
193  !***************************
194
195  if (abs(itime-memtime(1)).lt.abs(itime-memtime(2))) then
196    memindnext=1
197  else
198    memindnext=2
199  endif
200
201  ! Determine nested grid coordinates
202  !**********************************
203
204  if (ngrid.gt.0) then
205    xtn=(xt-xln(ngrid))*xresoln(ngrid)
206    ytn=(yt-yln(ngrid))*yresoln(ngrid)
207    ix=int(xtn)
208    jy=int(ytn)
209    nix=nint(xtn)
210    njy=nint(ytn)
211  else
212    ix=int(xt)
213    jy=int(yt)
214    nix=nint(xt)
215    njy=nint(yt)
216  endif
217  ixp=ix+1
218  jyp=jy+1
219
220
221  ! Compute maximum mixing height around particle position
222  !*******************************************************
223
224  h=0.
225  if (ngrid.le.0) then
226    do k=1,2
227      do j=jy,jyp
228        do i=ix,ixp
229          if (hmix(i,j,1,k).gt.h) h=hmix(i,j,1,k)
230        end do
231      end do
232    end do
233    tropop=tropopause(nix,njy,1,1)
234  else
235    do k=1,2
236      do j=jy,jyp
237        do i=ix,ixp
238          if (hmixn(i,j,1,k,ngrid).gt.h) h=hmixn(i,j,1,k,ngrid)
239        end do
240      end do
241    end do
242    tropop=tropopausen(nix,njy,1,1,ngrid)
243  endif
244
245  zeta=zt/h
246
247
248
249  !*************************************************************
250  ! If particle is in the PBL, interpolate once and then make a
251  ! time loop until end of interval is reached
252  !*************************************************************
253
254  if (zeta.le.1.) then
255
256  ! BEGIN TIME LOOP
257  !================
258
259    loop=0
260100   loop=loop+1
261      if (method.eq.1) then
262        ldt=min(ldt,abs(lsynctime-itimec+itime))
263        itimec=itimec+ldt*ldirect
264      else
265        ldt=abs(lsynctime)
266        itimec=itime+lsynctime
267      endif
268      dt=real(ldt)
269
270      zeta=zt/h
271
272
273      if (loop.eq.1) then
274        if (ngrid.le.0) then
275          xts=real(xt)
276          yts=real(yt)
277          call interpol_all(itime,xts,yts,zt)
278        else
279          call interpol_all_nests(itime,xtn,ytn,zt)
280        endif
281
282      else
283
284
285  ! Determine the level below the current position for u,v,rho
286  !***********************************************************
287
288        do i=2,nz
289          if (height(i).gt.zt) then
290            indz=i-1
291            indzp=i
292            goto 6
293          endif
294        end do
2956       continue
296
297  ! If one of the levels necessary is not yet available,
298  ! calculate it
299  !*****************************************************
300
301        do i=indz,indzp
302          if (indzindicator(i)) then
303            if (ngrid.le.0) then
304              call interpol_misslev(i)
305            else
306              call interpol_misslev_nests(i)
307            endif
308          endif
309        end do
310      endif
311
312
313  ! Vertical interpolation of u,v,w,rho and drhodz
314  !***********************************************
315
316  ! Vertical distance to the level below and above current position
317  ! both in terms of (u,v) and (w) fields
318  !****************************************************************
319
320      dz=1./(height(indzp)-height(indz))
321      dz1=(zt-height(indz))*dz
322      dz2=(height(indzp)-zt)*dz
323
324      u=dz1*uprof(indzp)+dz2*uprof(indz)
325      v=dz1*vprof(indzp)+dz2*vprof(indz)
326      w=dz1*wprof(indzp)+dz2*wprof(indz)
327      rhoa=dz1*rhoprof(indzp)+dz2*rhoprof(indz)
328      rhograd=dz1*rhogradprof(indzp)+dz2*rhogradprof(indz)
329
330
331  ! Compute the turbulent disturbances
332  ! Determine the sigmas and the timescales
333  !****************************************
334
335      if (turbswitch) then
336        call hanna(zt)
337      else
338        call hanna1(zt)
339      endif
340
341
342  !*****************************************
343  ! Determine the new diffusivity velocities
344  !*****************************************
345
346  ! Horizontal components
347  !**********************
348
349      if (nrand+1.gt.maxrand) nrand=1
350      if (dt/tlu.lt..5) then
351        up=(1.-dt/tlu)*up+rannumb(nrand)*sigu*sqrt(2.*dt/tlu)
352      else
353        ru=exp(-dt/tlu)
354        up=ru*up+rannumb(nrand)*sigu*sqrt(1.-ru**2)
355      endif
356      if (dt/tlv.lt..5) then
357        vp=(1.-dt/tlv)*vp+rannumb(nrand+1)*sigv*sqrt(2.*dt/tlv)
358      else
359        rv=exp(-dt/tlv)
360        vp=rv*vp+rannumb(nrand+1)*sigv*sqrt(1.-rv**2)
361      endif
362      nrand=nrand+2
363
364
365      if (nrand+ifine.gt.maxrand) nrand=1
366      rhoaux=rhograd/rhoa
367      dtf=dt*fine
368
369      dtftlw=dtf/tlw
370
371  ! Loop over ifine short time steps for vertical component
372  !********************************************************
373
374      do i=1,ifine
375
376  ! Determine the drift velocity and density correction velocity
377  !*************************************************************
378
379        if (turbswitch) then
380          if (dtftlw.lt..5) then
381            wp=((1.-dtftlw)*wp+rannumb(nrand+i)*sqrt(2.*dtftlw) &
382                 +dtf*(dsigwdz+rhoaux*sigw))*real(icbt)
383          else
384            rw=exp(-dtftlw)
385            wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2) &
386                 +tlw*(1.-rw)*(dsigwdz+rhoaux*sigw))*real(icbt)
387          endif
388          delz=wp*sigw*dtf
389        else
390          rw=exp(-dtftlw)
391          wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2)*sigw &
392               +tlw*(1.-rw)*(dsigw2dz+rhoaux*sigw**2))*real(icbt)
393          delz=wp*dtf
394        endif
395
396  !****************************************************
397  ! Compute turbulent vertical displacement of particle
398  !****************************************************
399
400        if (abs(delz).gt.h) delz=mod(delz,h)
401
402  ! Determine if particle transfers to a "forbidden state" below the ground
403  ! or above the mixing height
404  !************************************************************************
405
406        if (delz.lt.-zt) then         ! reflection at ground
407          icbt=-1
408          zt=-zt-delz
409        else if (delz.gt.(h-zt)) then ! reflection at h
410          icbt=-1
411          zt=-zt-delz+2.*h
412        else                         ! no reflection
413          icbt=1
414          zt=zt+delz
415        endif
416
417        if (i.ne.ifine) then
418          zeta=zt/h
419          call hanna_short(zt)
420        endif
421
422      end do
423      nrand=nrand+i
424
425  ! Determine time step for next integration
426  !*****************************************
427
428      if (turbswitch) then
429        ldt=int(min(tlw,h/max(2.*abs(wp*sigw),1.e-5), &
430             0.5/abs(dsigwdz))*ctl)
431      else
432        ldt=int(min(tlw,h/max(2.*abs(wp),1.e-5))*ctl)
433      endif
434      ldt=max(ldt,mintime)
435
436
437  ! If particle represents only a single species, add gravitational settling
438  ! velocity. The settling velocity is zero for gases, or if particle
439  ! represents more than one species
440  !*************************************************************************
441
442      if (mdomainfill.eq.0) then
443        do nsp=1,nspec
444          if (xmass(nrelpoint,nsp).gt.eps2) goto 887
445        end do
446887     nsp=min(nsp,nspec)
447!!$        if (density(nsp).gt.0.) &
448!!$             call get_settling(itime,xts,yts,zt,nsp,settling)    !old
449        if (density(nsp).gt.0.) &
450             call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
451        w=w+settling
452      endif
453
454  ! Horizontal displacements during time step dt are small real values compared
455  ! to the position; adding the two, would result in large numerical errors.
456  ! Thus, displacements are accumulated during lsynctime and are added to the
457  ! position at the end
458  !****************************************************************************
459
460      dxsave=dxsave+u*dt
461      dysave=dysave+v*dt
462      dawsave=dawsave+up*dt
463      dcwsave=dcwsave+vp*dt
464      zt=zt+w*dt*real(ldirect)
465
466      ! HSO/AL: Particle managed to go over highest level -> interpolation error in goto 700
467      !          alias interpol_wind (division by zero)
468      if (zt.ge.height(nz)) zt=height(nz)-100.*eps
469
470      if (zt.gt.h) then
471        if (itimec.eq.itime+lsynctime) goto 99
472        goto 700    ! complete the current interval above PBL
473      endif
474
475
476  !!! CHANGE: TEST OF THE WELL-MIXED CRITERION
477  !!! These lines may be switched on to test the well-mixed criterion
478  !if (zt.le.h) then
479  !  zacc=zacc+zt/h*dt
480  !  hsave=hsave+h*dt
481  !  tacc=tacc+dt
482  !  do 67 i=1,iclass
483  !    if ((zt/h.gt.th(i-1)).and.(zt/h.le.th(i)))
484  !    +    t(i)=t(i)+dt
485  !7        continue
486  !endif
487  !if ((mod(itime,10800).eq.0).and.dump) then
488  ! dump=.false.
489  ! write(550,'(i6,12f10.3)') itime,hsave/tacc,zacc/tacc,
490  !    + (t(i)/tacc*real(iclass),i=1,iclass)
491  !  zacc=0.
492  !  tacc=0.
493  !  do 68 i=1,iclass
494  !8        t(i)=0.
495  !  hsave=0.
496  !endif
497  !if (mod(itime,10800).ne.0) dump=.true.
498  !!! CHANGE
499
500
501  ! Determine probability of deposition
502  !************************************
503
504      if ((DRYDEP).and.(zt.lt.2.*href)) then
505        do ks=1,nspec
506          if (DRYDEPSPEC(ks)) then
507            if (depoindicator(ks)) then
508              if (ngrid.le.0) then
509                call interpol_vdep(ks,vdepo(ks))
510              else
511                call interpol_vdep_nests(ks,vdepo(ks))
512              endif
513            endif
514  ! correction by Petra Seibert, 10 April 2001
515  !   this formulation means that prob(n) = 1 - f(0)*...*f(n)
516  !   where f(n) is the exponential term
517               prob(ks)=1.+(prob(ks)-1.)* &
518                    exp(-vdepo(ks)*abs(dt)/(2.*href))
519          endif
520        end do
521      endif
522
523      if (zt.lt.0.) zt=min(h-eps2,-1.*zt)    ! if particle below ground -> reflection
524
525      if (itimec.eq.(itime+lsynctime)) then
526        usig=0.5*(usigprof(indzp)+usigprof(indz))
527        vsig=0.5*(vsigprof(indzp)+vsigprof(indz))
528        wsig=0.5*(wsigprof(indzp)+wsigprof(indz))
529        goto 99  ! finished
530      endif
531      goto 100
532
533  ! END TIME LOOP
534  !==============
535
536
537  endif
538
539
540
541  !**********************************************************
542  ! For all particles that are outside the PBL, make a single
543  ! time step. Only horizontal turbulent disturbances are
544  ! calculated. Vertical disturbances are reset.
545  !**********************************************************
546
547
548  ! Interpolate the wind
549  !*********************
550
551700   continue
552  if (ngrid.le.0) then
553    xts=real(xt)
554    yts=real(yt)
555    call interpol_wind(itime,xts,yts,zt)
556  else
557    call interpol_wind_nests(itime,xtn,ytn,zt)
558  endif
559
560
561  ! Compute everything for above the PBL
562
563  ! Assume constant, uncorrelated, turbulent perturbations
564  ! In the stratosphere, use a small vertical diffusivity d_strat,
565  ! in the troposphere, use a larger horizontal diffusivity d_trop.
566  ! Turbulent velocity scales are determined based on sqrt(d_trop/dt)
567  !******************************************************************
568
569  ldt=abs(lsynctime-itimec+itime)
570  dt=real(ldt)
571
572  if (zt.lt.tropop) then  ! in the troposphere
573    uxscale=sqrt(2.*d_trop/dt)
574    if (nrand+1.gt.maxrand) nrand=1
575    ux=rannumb(nrand)*uxscale
576    vy=rannumb(nrand+1)*uxscale
577    nrand=nrand+2
578    wp=0.
579  else if (zt.lt.tropop+1000.) then     ! just above the tropopause: make transition
580    weight=(zt-tropop)/1000.
581    uxscale=sqrt(2.*d_trop/dt*(1.-weight))
582    if (nrand+2.gt.maxrand) nrand=1
583    ux=rannumb(nrand)*uxscale
584    vy=rannumb(nrand+1)*uxscale
585    wpscale=sqrt(2.*d_strat/dt*weight)
586    wp=rannumb(nrand+2)*wpscale+d_strat/1000.
587    nrand=nrand+3
588  else                 ! in the stratosphere
589    if (nrand.gt.maxrand) nrand=1
590    ux=0.
591    vy=0.
592    wpscale=sqrt(2.*d_strat/dt)
593    wp=rannumb(nrand)*wpscale
594    nrand=nrand+1
595  endif
596
597
598  ! If particle represents only a single species, add gravitational settling
599  ! velocity. The settling velocity is zero for gases
600  !*************************************************************************
601
602
603
604    if (mdomainfill.eq.0) then
605      do nsp=1,nspec
606        if (xmass(nrelpoint,nsp).gt.eps2) goto 888
607      end do
608888   nsp=min(nsp,nspec)
609!!$      if (density(nsp).gt.0.) &
610!!$           call get_settling(itime,xts,yts,zt,nsp,settling)    !old
611      if (density(nsp).gt.0.) &
612           call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
613      w=w+settling
614    endif
615
616  ! Calculate position at time step itime+lsynctime
617  !************************************************
618
619  dxsave=dxsave+(u+ux)*dt
620  dysave=dysave+(v+vy)*dt
621  zt=zt+(w+wp)*dt*real(ldirect)
622  if (zt.lt.0.) zt=min(h-eps2,-1.*zt)    ! if particle below ground -> reflection
623
62499   continue
625
626
627
628  !****************************************************************
629  ! Add mesoscale random disturbances
630  ! This is done only once for the whole lsynctime interval to save
631  ! computation time
632  !****************************************************************
633
634
635  ! Mesoscale wind velocity fluctuations are obtained by scaling
636  ! with the standard deviation of the grid-scale winds surrounding
637  ! the particle location, multiplied by a factor turbmesoscale.
638  ! The autocorrelation time constant is taken as half the
639  ! time interval between wind fields
640  !****************************************************************
641
642  r=exp(-2.*real(abs(lsynctime))/real(lwindinterv))
643  rs=sqrt(1.-r**2)
644  if (nrand+2.gt.maxrand) nrand=1
645  usigold=r*usigold+rs*rannumb(nrand)*usig*turbmesoscale
646  vsigold=r*vsigold+rs*rannumb(nrand+1)*vsig*turbmesoscale
647  wsigold=r*wsigold+rs*rannumb(nrand+2)*wsig*turbmesoscale
648
649  dxsave=dxsave+usigold*real(lsynctime)
650  dysave=dysave+vsigold*real(lsynctime)
651
652  zt=zt+wsigold*real(lsynctime)
653  if (zt.lt.0.) zt=-1.*zt    ! if particle below ground -> refletion
654
655  !*************************************************************
656  ! Transform along and cross wind components to xy coordinates,
657  ! add them to u and v, transform u,v to grid units/second
658  ! and calculate new position
659  !*************************************************************
660
661  call windalign(dxsave,dysave,dawsave,dcwsave,ux,vy)
662  dxsave=dxsave+ux
663  dysave=dysave+vy
664  if (ngrid.ge.0) then
665    cosfact=dxconst/cos((yt*dy+ylat0)*pi180)
666    xt=xt+real(dxsave*cosfact*real(ldirect),kind=dp)
667    yt=yt+real(dysave*dyconst*real(ldirect),kind=dp)
668  else if (ngrid.eq.-1) then      ! around north pole
669    xlon=xlon0+xt*dx
670    ylat=ylat0+yt*dy
671    call cll2xy(northpolemap,ylat,xlon,xpol,ypol)
672    gridsize=1000.*cgszll(northpolemap,ylat,xlon)
673    dxsave=dxsave/gridsize
674    dysave=dysave/gridsize
675    xpol=xpol+dxsave*real(ldirect)
676    ypol=ypol+dysave*real(ldirect)
677    call cxy2ll(northpolemap,xpol,ypol,ylat,xlon)
678    xt=(xlon-xlon0)/dx
679    yt=(ylat-ylat0)/dy
680  else if (ngrid.eq.-2) then    ! around south pole
681    xlon=xlon0+xt*dx
682    ylat=ylat0+yt*dy
683    call cll2xy(southpolemap,ylat,xlon,xpol,ypol)
684    gridsize=1000.*cgszll(southpolemap,ylat,xlon)
685    dxsave=dxsave/gridsize
686    dysave=dysave/gridsize
687    xpol=xpol+dxsave*real(ldirect)
688    ypol=ypol+dysave*real(ldirect)
689    call cxy2ll(southpolemap,xpol,ypol,ylat,xlon)
690    xt=(xlon-xlon0)/dx
691    yt=(ylat-ylat0)/dy
692  endif
693
694
695  ! If global data are available, use cyclic boundary condition
696  !************************************************************
697
698  if (xglobal) then
699    if (xt.ge.real(nxmin1)) xt=xt-real(nxmin1)
700    if (xt.lt.0.) xt=xt+real(nxmin1)
701    if (xt.le.eps) xt=eps
702    if (abs(xt-real(nxmin1)).le.eps) xt=real(nxmin1)-eps
703  endif
704
705  ! HSO/AL: Prevent particles from disappearing at the pole
706  !******************************************************************
707
708  if ( yt.lt.0. ) then
709    xt=mod(xt+180.,360.)
710    yt=-yt
711  else if ( yt.gt.real(nymin1) ) then
712    xt=mod(xt+180.,360.)
713    yt=2*real(nymin1)-yt
714  endif
715
716  ! Check position: If trajectory outside model domain, terminate it
717  !*****************************************************************
718
719  if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. &
720       (yt.gt.real(nymin1))) then
721    nstop=3
722    return
723  endif
724
725  ! If particle above highest model level, set it back into the domain
726  !*******************************************************************
727
728  if (zt.ge.height(nz)) zt=height(nz)-100.*eps
729
730
731  !************************************************************************
732  ! Now we could finish, as this was done in FLEXPART versions up to 4.0.
733  ! However, truncation errors of the advection can be significantly
734  ! reduced by doing one iteration of the Petterssen scheme, if this is
735  ! possible.
736  ! Note that this is applied only to the grid-scale winds, not to
737  ! the turbulent winds.
738  !************************************************************************
739
740  ! The Petterssen scheme can only applied with long time steps (only then u
741  ! is the "old" wind as required by the scheme); otherwise do nothing
742  !*************************************************************************
743
744  if (ldt.ne.abs(lsynctime)) return
745
746  ! The Petterssen scheme can only be applied if the ending time of the time step
747  ! (itime+ldt*ldirect) is still between the two wind fields held in memory;
748  ! otherwise do nothing
749  !******************************************************************************
750
751  if (abs(itime+ldt*ldirect).gt.abs(memtime(2))) return
752
753  ! Apply it also only if starting and ending point of current time step are on
754  ! the same grid; otherwise do nothing
755  !*****************************************************************************
756  if (nglobal.and.(yt.gt.switchnorthg)) then
757    ngr=-1
758  else if (sglobal.and.(yt.lt.switchsouthg)) then
759    ngr=-2
760  else
761    ngr=0
762    do j=numbnests,1,-1
763      if ((xt.gt.xln(j)+eps).and.(xt.lt.xrn(j)-eps).and. &
764           (yt.gt.yln(j)+eps).and.(yt.lt.yrn(j)-eps)) then
765        ngr=j
766        goto 43
767      endif
768    end do
76943   continue
770  endif
771
772  if (ngr.ne.ngrid) return
773
774  ! Determine nested grid coordinates
775  !**********************************
776
777  if (ngrid.gt.0) then
778    xtn=(xt-xln(ngrid))*xresoln(ngrid)
779    ytn=(yt-yln(ngrid))*yresoln(ngrid)
780    ix=int(xtn)
781    jy=int(ytn)
782  else
783    ix=int(xt)
784    jy=int(yt)
785  endif
786  ixp=ix+1
787  jyp=jy+1
788
789
790  ! Memorize the old wind
791  !**********************
792
793  uold=u
794  vold=v
795  wold=w
796
797  ! Interpolate wind at new position and time
798  !******************************************
799
800  if (ngrid.le.0) then
801    xts=real(xt)
802    yts=real(yt)
803    call interpol_wind_short(itime+ldt*ldirect,xts,yts,zt)
804  else
805    call interpol_wind_short_nests(itime+ldt*ldirect,xtn,ytn,zt)
806  endif
807
808  if (mdomainfill.eq.0) then
809    do nsp=1,nspec
810      if (xmass(nrelpoint,nsp).gt.eps2) goto 889
811    end do
812889   nsp=min(nsp,nspec)
813!!$    if (density(nsp).gt.0.) &
814!!$         call get_settling(itime+ldt,xts,yts,zt,nsp,settling)    !old
815    if (density(nsp).gt.0.) &
816         call get_settling(itime+ldt,real(xt),real(yt),zt,nsp,settling)  !bugfix
817    w=w+settling
818  endif
819
820
821  ! Determine the difference vector between new and old wind
822  ! (use half of it to correct position according to Petterssen)
823  !*************************************************************
824
825  u=(u-uold)/2.
826  v=(v-vold)/2.
827  w=(w-wold)/2.
828
829
830  ! Finally, correct the old position
831  !**********************************
832
833  zt=zt+w*real(ldt*ldirect)
834  if (zt.lt.0.) zt=min(h-eps2,-1.*zt)    ! if particle below ground -> reflection
835  if (ngrid.ge.0) then
836    cosfact=dxconst/cos((yt*dy+ylat0)*pi180)
837    xt=xt+real(u*cosfact*real(ldt*ldirect),kind=dp)
838    yt=yt+real(v*dyconst*real(ldt*ldirect),kind=dp)
839  else if (ngrid.eq.-1) then      ! around north pole
840    xlon=xlon0+xt*dx
841    ylat=ylat0+yt*dy
842    call cll2xy(northpolemap,ylat,xlon,xpol,ypol)
843    gridsize=1000.*cgszll(northpolemap,ylat,xlon)
844    u=u/gridsize
845    v=v/gridsize
846    xpol=xpol+u*real(ldt*ldirect)
847    ypol=ypol+v*real(ldt*ldirect)
848    call cxy2ll(northpolemap,xpol,ypol,ylat,xlon)
849    xt=(xlon-xlon0)/dx
850    yt=(ylat-ylat0)/dy
851  else if (ngrid.eq.-2) then    ! around south pole
852    xlon=xlon0+xt*dx
853    ylat=ylat0+yt*dy
854    call cll2xy(southpolemap,ylat,xlon,xpol,ypol)
855    gridsize=1000.*cgszll(southpolemap,ylat,xlon)
856    u=u/gridsize
857    v=v/gridsize
858    xpol=xpol+u*real(ldt*ldirect)
859    ypol=ypol+v*real(ldt*ldirect)
860    call cxy2ll(southpolemap,xpol,ypol,ylat,xlon)
861    xt=(xlon-xlon0)/dx
862    yt=(ylat-ylat0)/dy
863  endif
864
865  ! If global data are available, use cyclic boundary condition
866  !************************************************************
867
868  if (xglobal) then
869    if (xt.ge.real(nxmin1)) xt=xt-real(nxmin1)
870    if (xt.lt.0.) xt=xt+real(nxmin1)
871    if (xt.le.eps) xt=eps
872    if (abs(xt-real(nxmin1)).le.eps) xt=real(nxmin1)-eps
873  endif
874
875  ! HSO/AL: Prevent particles from disappearing at the pole
876  !******************************************************************
877
878  if ( yt.lt.0. ) then
879    xt=mod(xt+180.,360.)
880    yt=-yt
881  else if ( yt.gt.real(nymin1) ) then
882    xt=mod(xt+180.,360.)
883    yt=2*real(nymin1)-yt
884  endif
885
886  ! Check position: If trajectory outside model domain, terminate it
887  !*****************************************************************
888
889  if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. &
890       (yt.gt.real(nymin1))) then
891    nstop=3
892    return
893  endif
894
895  ! If particle above highest model level, set it back into the domain
896  !*******************************************************************
897
898  if (zt.ge.height(nz)) zt=height(nz)-100.*eps
899
900
901end subroutine advance
902
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG