source: branches/petra/src/advance.f90 @ 36

Last change on this file since 36 was 36, checked in by pesei, 7 years ago

Implement switch for incremental deposition, see ticket:113 and many small changes, see changelog.txt

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