!********************************************************************** ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 * ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * ! * ! This file is part of FLEXPART. * ! * ! FLEXPART is free software: you can redistribute it and/or modify * ! it under the terms of the GNU General Public License as published by* ! the Free Software Foundation, either version 3 of the License, or * ! (at your option) any later version. * ! * ! FLEXPART is distributed in the hope that it will be useful, * ! but WITHOUT ANY WARRANTY; without even the implied warranty of * ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * ! GNU General Public License for more details. * ! * ! You should have received a copy of the GNU General Public License * ! along with FLEXPART. If not, see . * !********************************************************************** subroutine advance(itime,nrelpoint,ldt,up,vp,wp, & usigold,vsigold,wsigold,nstop,xt,yt,zt,prob,icbt) ! i i i/oi/oi/o ! i/o i/o i/o o i/oi/oi/o i/o i/o !***************************************************************************** ! * ! Calculation of turbulent particle trajectories utilizing a * ! zero-acceleration scheme, which is corrected by a numerically more * ! accurate Petterssen scheme whenever possible. * ! * ! Particle positions are read in, incremented, and returned to the calling * ! program. * ! * ! In different regions of the atmosphere (PBL vs. free troposphere), * ! different parameters are needed for advection, parameterizing turbulent * ! velocities, etc. For efficiency, different interpolation routines have * ! been written for these different cases, with the disadvantage that there * ! exist several routines doing almost the same. They all share the * ! included file 'interpol_mod'. The following * ! interpolation routines are used: * ! * ! interpol_all(_nests) interpolates everything (called inside the PBL) * ! interpol_misslev(_nests) if a particle moves vertically in the PBL, * ! additional parameters are interpolated if it * ! crosses a model level * ! interpol_wind(_nests) interpolates the wind and determines the * ! standard deviation of the wind (called outside * ! PBL) also interpolates potential vorticity * ! interpol_wind_short(_nests) only interpolates the wind (needed for the * ! Petterssen scheme) * ! interpol_vdep(_nests) interpolates deposition velocities * ! * ! * ! Author: A. Stohl * ! * ! 16 December 1997 * ! * ! Changes: * ! * ! 8 April 2000: Deep convection parameterization * ! * ! May 2002: Petterssen scheme introduced * ! * !***************************************************************************** ! * ! Variables: * ! icbt 1 if particle not transferred to forbidden state, * ! else -1 * ! dawsave accumulated displacement in along-wind direction * ! dcwsave accumulated displacement in cross-wind direction * ! dxsave accumulated displacement in longitude * ! dysave accumulated displacement in latitude * ! h [m] Mixing height * ! lwindinterv [s] time interval between two wind fields * ! itime [s] time at which this subroutine is entered * ! itimec [s] actual time, which is incremented in this subroutine * ! href [m] height for which dry deposition velocity is calculated * ! ladvance [s] Total integration time period * ! ldirect 1 forward, -1 backward * ! ldt [s] Time step for the next integration * ! lsynctime [s] Synchronisation interval of FLEXPART * ! ngrid index which grid is to be used * ! nrand index for a variable to be picked from rannumb * ! nstop if > 1 particle has left domain and must be stopped * ! prob probability of absorption due to dry deposition * ! rannumb(maxrand) normally distributed random variables * ! rhoa air density * ! rhograd vertical gradient of the air density * ! up,vp,wp random velocities due to turbulence (along wind, cross * ! wind, vertical wind * ! usig,vsig,wsig mesoscale wind fluctuations * ! usigold,vsigold,wsigold like usig, etc., but for the last time step * ! vdepo Deposition velocities for all species * ! xt,yt,zt Particle position * ! * !***************************************************************************** use point_mod use par_mod use com_mod use interpol_mod use hanna_mod use cmapf_mod implicit none real(kind=dp) :: xt,yt real :: zt,xts,yts,weight integer :: itime,itimec,nstop,ldt,i,j,k,nrand,loop,memindnext integer :: ngr,nix,njy,ks,nsp,nrelpoint real :: dz,dz1,dz2,xlon,ylat,xpol,ypol,gridsize real :: ru,rv,rw,dt,ux,vy,cosfact,xtn,ytn,tropop real :: prob(maxspec),up,vp,wp,dxsave,dysave,dawsave real :: dcwsave real :: usigold,vsigold,wsigold,r,rs real :: uold,vold,wold,vdepo(maxspec) !real uprof(nzmax),vprof(nzmax),wprof(nzmax) !real usigprof(nzmax),vsigprof(nzmax),wsigprof(nzmax) !real rhoprof(nzmax),rhogradprof(nzmax) real :: rhoa,rhograd,ran3,delz,dtf,rhoaux,dtftlw,uxscale,wpscale integer(kind=2) :: icbt real,parameter :: eps=nxmax/3.e5,eps2=1.e-9 !!! CHANGE: TEST OF THE WELL-MIXED CRITERION ! integer,parameter :: iclass=10 ! real(kind=dp) :: zacc,tacc,t(iclass),th(0:iclass),hsave ! logical dump ! save zacc,tacc,t,th,hsave,dump !!! CHANGE integer :: idummy = -7 real :: settling = 0. !!! CHANGE: TEST OF THE WELL-MIXED CRITERION !if (idummy.eq.-7) then !open(550,file='WELLMIXEDTEST') !do 17 i=0,iclass !7 th(i)=real(i)/real(iclass) !endif !!! CHANGE nstop=0 do i=1,nmixz indzindicator(i)=.true. end do if (DRYDEP) then ! reset probability for deposition do ks=1,nspec depoindicator(ks)=.true. prob(ks)=0. end do endif dxsave=0. ! reset position displacements dysave=0. ! due to mean wind dawsave=0. ! and turbulent wind dcwsave=0. itimec=itime nrand=int(ran3(idummy)*real(maxrand-1))+1 ! Determine whether lat/long grid or polarstereographic projection ! is to be used ! Furthermore, determine which nesting level to be used !***************************************************************** if (nglobal.and.(yt.gt.switchnorthg)) then ngrid=-1 else if (sglobal.and.(yt.lt.switchsouthg)) then ngrid=-2 else ngrid=0 do j=numbnests,1,-1 if ((xt.gt.xln(j)+eps).and.(xt.lt.xrn(j)-eps).and. & (yt.gt.yln(j)+eps).and.(yt.lt.yrn(j)-eps)) then ngrid=j goto 23 endif end do 23 continue endif !*************************** ! Interpolate necessary data !*************************** if (abs(itime-memtime(1)).lt.abs(itime-memtime(2))) then memindnext=1 else memindnext=2 endif ! Determine nested grid coordinates !********************************** if (ngrid.gt.0) then xtn=(xt-xln(ngrid))*xresoln(ngrid) ytn=(yt-yln(ngrid))*yresoln(ngrid) ix=int(xtn) jy=int(ytn) nix=nint(xtn) njy=nint(ytn) else ix=int(xt) jy=int(yt) nix=nint(xt) njy=nint(yt) endif ixp=ix+1 jyp=jy+1 ! Compute maximum mixing height around particle position !******************************************************* h=0. if (ngrid.le.0) then do k=1,2 do j=jy,jyp do i=ix,ixp if (hmix(i,j,1,k).gt.h) h=hmix(i,j,1,k) end do end do end do tropop=tropopause(nix,njy,1,1) else do k=1,2 do j=jy,jyp do i=ix,ixp if (hmixn(i,j,1,k,ngrid).gt.h) h=hmixn(i,j,1,k,ngrid) end do end do end do tropop=tropopausen(nix,njy,1,1,ngrid) endif zeta=zt/h !************************************************************* ! If particle is in the PBL, interpolate once and then make a ! time loop until end of interval is reached !************************************************************* if (zeta.le.1.) then ! BEGIN TIME LOOP !================ loop=0 100 loop=loop+1 if (method.eq.1) then ldt=min(ldt,abs(lsynctime-itimec+itime)) itimec=itimec+ldt*ldirect else ldt=abs(lsynctime) itimec=itime+lsynctime endif dt=real(ldt) zeta=zt/h if (loop.eq.1) then if (ngrid.le.0) then xts=real(xt) yts=real(yt) call interpol_all(itime,xts,yts,zt) else call interpol_all_nests(itime,xtn,ytn,zt) endif else ! Determine the level below the current position for u,v,rho !*********************************************************** do i=2,nz if (height(i).gt.zt) then indz=i-1 indzp=i goto 6 endif end do 6 continue ! If one of the levels necessary is not yet available, ! calculate it !***************************************************** do i=indz,indzp if (indzindicator(i)) then if (ngrid.le.0) then call interpol_misslev(i) else call interpol_misslev_nests(i) endif endif end do endif ! Vertical interpolation of u,v,w,rho and drhodz !*********************************************** ! Vertical distance to the level below and above current position ! both in terms of (u,v) and (w) fields !**************************************************************** dz=1./(height(indzp)-height(indz)) dz1=(zt-height(indz))*dz dz2=(height(indzp)-zt)*dz u=dz1*uprof(indzp)+dz2*uprof(indz) v=dz1*vprof(indzp)+dz2*vprof(indz) w=dz1*wprof(indzp)+dz2*wprof(indz) rhoa=dz1*rhoprof(indzp)+dz2*rhoprof(indz) rhograd=dz1*rhogradprof(indzp)+dz2*rhogradprof(indz) ! Compute the turbulent disturbances ! Determine the sigmas and the timescales !**************************************** if (turbswitch) then call hanna(zt) else call hanna1(zt) endif !***************************************** ! Determine the new diffusivity velocities !***************************************** ! Horizontal components !********************** if (nrand+1.gt.maxrand) nrand=1 if (dt/tlu.lt..5) then up=(1.-dt/tlu)*up+rannumb(nrand)*sigu*sqrt(2.*dt/tlu) else ru=exp(-dt/tlu) up=ru*up+rannumb(nrand)*sigu*sqrt(1.-ru**2) endif if (dt/tlv.lt..5) then vp=(1.-dt/tlv)*vp+rannumb(nrand+1)*sigv*sqrt(2.*dt/tlv) else rv=exp(-dt/tlv) vp=rv*vp+rannumb(nrand+1)*sigv*sqrt(1.-rv**2) endif nrand=nrand+2 if (nrand+ifine.gt.maxrand) nrand=1 rhoaux=rhograd/rhoa dtf=dt*fine dtftlw=dtf/tlw ! Loop over ifine short time steps for vertical component !******************************************************** do i=1,ifine ! Determine the drift velocity and density correction velocity !************************************************************* if (turbswitch) then if (dtftlw.lt..5) then wp=((1.-dtftlw)*wp+rannumb(nrand+i)*sqrt(2.*dtftlw) & +dtf*(dsigwdz+rhoaux*sigw))*real(icbt) else rw=exp(-dtftlw) wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2) & +tlw*(1.-rw)*(dsigwdz+rhoaux*sigw))*real(icbt) endif delz=wp*sigw*dtf else rw=exp(-dtftlw) wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2)*sigw & +tlw*(1.-rw)*(dsigw2dz+rhoaux*sigw**2))*real(icbt) delz=wp*dtf endif !**************************************************** ! Compute turbulent vertical displacement of particle !**************************************************** if (abs(delz).gt.h) delz=mod(delz,h) ! Determine if particle transfers to a "forbidden state" below the ground ! or above the mixing height !************************************************************************ if (delz.lt.-zt) then ! reflection at ground icbt=-1 zt=-zt-delz else if (delz.gt.(h-zt)) then ! reflection at h icbt=-1 zt=-zt-delz+2.*h else ! no reflection icbt=1 zt=zt+delz endif if (i.ne.ifine) then zeta=zt/h call hanna_short(zt) endif end do nrand=nrand+i ! Determine time step for next integration !***************************************** if (turbswitch) then ldt=int(min(tlw,h/max(2.*abs(wp*sigw),1.e-5), & 0.5/abs(dsigwdz))*ctl) else ldt=int(min(tlw,h/max(2.*abs(wp),1.e-5))*ctl) endif ldt=max(ldt,mintime) ! If particle represents only a single species, add gravitational settling ! velocity. The settling velocity is zero for gases, or if particle ! represents more than one species !************************************************************************* if (mdomainfill.eq.0) then do nsp=1,nspec if (xmass(nrelpoint,nsp).gt.eps2) goto 887 end do 887 nsp=min(nsp,nspec) !!$ if (density(nsp).gt.0.) & !!$ call get_settling(itime,xts,yts,zt,nsp,settling) !old if (density(nsp).gt.0.) & call get_settling(itime,real(xt),real(yt),zt,nsp,settling) !bugfix w=w+settling endif ! Horizontal displacements during time step dt are small real values compared ! to the position; adding the two, would result in large numerical errors. ! Thus, displacements are accumulated during lsynctime and are added to the ! position at the end !**************************************************************************** dxsave=dxsave+u*dt dysave=dysave+v*dt dawsave=dawsave+up*dt dcwsave=dcwsave+vp*dt zt=zt+w*dt*real(ldirect) if (zt.gt.h) then if (itimec.eq.itime+lsynctime) goto 99 goto 700 ! complete the current interval above PBL endif !!! CHANGE: TEST OF THE WELL-MIXED CRITERION !!! These lines may be switched on to test the well-mixed criterion !if (zt.le.h) then ! zacc=zacc+zt/h*dt ! hsave=hsave+h*dt ! tacc=tacc+dt ! do 67 i=1,iclass ! if ((zt/h.gt.th(i-1)).and.(zt/h.le.th(i))) ! + t(i)=t(i)+dt !7 continue !endif !if ((mod(itime,10800).eq.0).and.dump) then ! dump=.false. ! write(550,'(i6,12f10.3)') itime,hsave/tacc,zacc/tacc, ! + (t(i)/tacc*real(iclass),i=1,iclass) ! zacc=0. ! tacc=0. ! do 68 i=1,iclass !8 t(i)=0. ! hsave=0. !endif !if (mod(itime,10800).ne.0) dump=.true. !!! CHANGE ! Determine probability of deposition !************************************ if ((DRYDEP).and.(zt.lt.2.*href)) then do ks=1,nspec if (DRYDEPSPEC(ks)) then if (depoindicator(ks)) then if (ngrid.le.0) then call interpol_vdep(ks,vdepo(ks)) else call interpol_vdep_nests(ks,vdepo(ks)) endif endif ! correction by Petra Seibert, 10 April 2001 ! this formulation means that prob(n) = 1 - f(0)*...*f(n) ! where f(n) is the exponential term prob(ks)=1.+(prob(ks)-1.)* & exp(-vdepo(ks)*abs(dt)/(2.*href)) endif end do endif if (zt.lt.0.) zt=min(h-eps2,-1.*zt) ! if particle below ground -> reflection if (itimec.eq.(itime+lsynctime)) then usig=0.5*(usigprof(indzp)+usigprof(indz)) vsig=0.5*(vsigprof(indzp)+vsigprof(indz)) wsig=0.5*(wsigprof(indzp)+wsigprof(indz)) goto 99 ! finished endif goto 100 ! END TIME LOOP !============== endif !********************************************************** ! For all particles that are outside the PBL, make a single ! time step. Only horizontal turbulent disturbances are ! calculated. Vertical disturbances are reset. !********************************************************** ! Interpolate the wind !********************* 700 continue if (ngrid.le.0) then xts=real(xt) yts=real(yt) call interpol_wind(itime,xts,yts,zt) else call interpol_wind_nests(itime,xtn,ytn,zt) endif ! Compute everything for above the PBL ! Assume constant, uncorrelated, turbulent perturbations ! In the stratosphere, use a small vertical diffusivity d_strat, ! in the troposphere, use a larger horizontal diffusivity d_trop. ! Turbulent velocity scales are determined based on sqrt(d_trop/dt) !****************************************************************** ldt=abs(lsynctime-itimec+itime) dt=real(ldt) if (zt.lt.tropop) then ! in the troposphere uxscale=sqrt(2.*d_trop/dt) if (nrand+1.gt.maxrand) nrand=1 ux=rannumb(nrand)*uxscale vy=rannumb(nrand+1)*uxscale nrand=nrand+2 wp=0. else if (zt.lt.tropop+1000.) then ! just above the tropopause: make transition weight=(zt-tropop)/1000. uxscale=sqrt(2.*d_trop/dt*(1.-weight)) if (nrand+2.gt.maxrand) nrand=1 ux=rannumb(nrand)*uxscale vy=rannumb(nrand+1)*uxscale wpscale=sqrt(2.*d_strat/dt*weight) wp=rannumb(nrand+2)*wpscale+d_strat/1000. nrand=nrand+3 else ! in the stratosphere if (nrand.gt.maxrand) nrand=1 ux=0. vy=0. wpscale=sqrt(2.*d_strat/dt) wp=rannumb(nrand)*wpscale nrand=nrand+1 endif ! If particle represents only a single species, add gravitational settling ! velocity. The settling velocity is zero for gases !************************************************************************* if (mdomainfill.eq.0) then do nsp=1,nspec if (xmass(nrelpoint,nsp).gt.eps2) goto 888 end do 888 nsp=min(nsp,nspec) !!$ if (density(nsp).gt.0.) & !!$ call get_settling(itime,xts,yts,zt,nsp,settling) !old if (density(nsp).gt.0.) & call get_settling(itime,real(xt),real(yt),zt,nsp,settling) !bugfix w=w+settling endif ! Calculate position at time step itime+lsynctime !************************************************ dxsave=dxsave+(u+ux)*dt dysave=dysave+(v+vy)*dt zt=zt+(w+wp)*dt*real(ldirect) if (zt.lt.0.) zt=min(h-eps2,-1.*zt) ! if particle below ground -> reflection 99 continue !**************************************************************** ! Add mesoscale random disturbances ! This is done only once for the whole lsynctime interval to save ! computation time !**************************************************************** ! Mesoscale wind velocity fluctuations are obtained by scaling ! with the standard deviation of the grid-scale winds surrounding ! the particle location, multiplied by a factor turbmesoscale. ! The autocorrelation time constant is taken as half the ! time interval between wind fields !**************************************************************** r=exp(-2.*real(abs(lsynctime))/real(lwindinterv)) rs=sqrt(1.-r**2) if (nrand+2.gt.maxrand) nrand=1 usigold=r*usigold+rs*rannumb(nrand)*usig*turbmesoscale vsigold=r*vsigold+rs*rannumb(nrand+1)*vsig*turbmesoscale wsigold=r*wsigold+rs*rannumb(nrand+2)*wsig*turbmesoscale dxsave=dxsave+usigold*real(lsynctime) dysave=dysave+vsigold*real(lsynctime) zt=zt+wsigold*real(lsynctime) if (zt.lt.0.) zt=-1.*zt ! if particle below ground -> refletion !************************************************************* ! Transform along and cross wind components to xy coordinates, ! add them to u and v, transform u,v to grid units/second ! and calculate new position !************************************************************* call windalign(dxsave,dysave,dawsave,dcwsave,ux,vy) dxsave=dxsave+ux dysave=dysave+vy if (ngrid.ge.0) then cosfact=dxconst/cos((yt*dy+ylat0)*pi180) xt=xt+real(dxsave*cosfact*real(ldirect),kind=dp) yt=yt+real(dysave*dyconst*real(ldirect),kind=dp) else if (ngrid.eq.-1) then ! around north pole xlon=xlon0+xt*dx ylat=ylat0+yt*dy call cll2xy(northpolemap,ylat,xlon,xpol,ypol) gridsize=1000.*cgszll(northpolemap,ylat,xlon) dxsave=dxsave/gridsize dysave=dysave/gridsize xpol=xpol+dxsave*real(ldirect) ypol=ypol+dysave*real(ldirect) call cxy2ll(northpolemap,xpol,ypol,ylat,xlon) xt=(xlon-xlon0)/dx yt=(ylat-ylat0)/dy else if (ngrid.eq.-2) then ! around south pole xlon=xlon0+xt*dx ylat=ylat0+yt*dy call cll2xy(southpolemap,ylat,xlon,xpol,ypol) gridsize=1000.*cgszll(southpolemap,ylat,xlon) dxsave=dxsave/gridsize dysave=dysave/gridsize xpol=xpol+dxsave*real(ldirect) ypol=ypol+dysave*real(ldirect) call cxy2ll(southpolemap,xpol,ypol,ylat,xlon) xt=(xlon-xlon0)/dx yt=(ylat-ylat0)/dy endif ! If global data are available, use cyclic boundary condition !************************************************************ if (xglobal) then if (xt.ge.real(nxmin1)) xt=xt-real(nxmin1) if (xt.lt.0.) xt=xt+real(nxmin1) if (xt.le.eps) xt=eps if (abs(xt-real(nxmin1)).le.eps) xt=real(nxmin1)-eps endif ! Check position: If trajectory outside model domain, terminate it !***************************************************************** if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. & (yt.ge.real(nymin1))) then nstop=3 return endif ! If particle above highest model level, set it back into the domain !******************************************************************* if (zt.ge.height(nz)) zt=height(nz)-100.*eps !************************************************************************ ! Now we could finish, as this was done in FLEXPART versions up to 4.0. ! However, truncation errors of the advection can be significantly ! reduced by doing one iteration of the Petterssen scheme, if this is ! possible. ! Note that this is applied only to the grid-scale winds, not to ! the turbulent winds. !************************************************************************ ! The Petterssen scheme can only applied with long time steps (only then u ! is the "old" wind as required by the scheme); otherwise do nothing !************************************************************************* if (ldt.ne.abs(lsynctime)) return ! The Petterssen scheme can only be applied if the ending time of the time step ! (itime+ldt*ldirect) is still between the two wind fields held in memory; ! otherwise do nothing !****************************************************************************** if (abs(itime+ldt*ldirect).gt.abs(memtime(2))) return ! Apply it also only if starting and ending point of current time step are on ! the same grid; otherwise do nothing !***************************************************************************** if (nglobal.and.(yt.gt.switchnorthg)) then ngr=-1 else if (sglobal.and.(yt.lt.switchsouthg)) then ngr=-2 else ngr=0 do j=numbnests,1,-1 if ((xt.gt.xln(j)+eps).and.(xt.lt.xrn(j)-eps).and. & (yt.gt.yln(j)+eps).and.(yt.lt.yrn(j)-eps)) then ngr=j goto 43 endif end do 43 continue endif if (ngr.ne.ngrid) return ! Determine nested grid coordinates !********************************** if (ngrid.gt.0) then xtn=(xt-xln(ngrid))*xresoln(ngrid) ytn=(yt-yln(ngrid))*yresoln(ngrid) ix=int(xtn) jy=int(ytn) else ix=int(xt) jy=int(yt) endif ixp=ix+1 jyp=jy+1 ! Memorize the old wind !********************** uold=u vold=v wold=w ! Interpolate wind at new position and time !****************************************** if (ngrid.le.0) then xts=real(xt) yts=real(yt) call interpol_wind_short(itime+ldt*ldirect,xts,yts,zt) else call interpol_wind_short_nests(itime+ldt*ldirect,xtn,ytn,zt) endif if (mdomainfill.eq.0) then do nsp=1,nspec if (xmass(nrelpoint,nsp).gt.eps2) goto 889 end do 889 nsp=min(nsp,nspec) !!$ if (density(nsp).gt.0.) & !!$ call get_settling(itime+ldt,xts,yts,zt,nsp,settling) !old if (density(nsp).gt.0.) & call get_settling(itime+ldt,real(xt),real(yt),zt,nsp,settling) !bugfix w=w+settling endif ! Determine the difference vector between new and old wind ! (use half of it to correct position according to Petterssen) !************************************************************* u=(u-uold)/2. v=(v-vold)/2. w=(w-wold)/2. ! Finally, correct the old position !********************************** zt=zt+w*real(ldt*ldirect) if (zt.lt.0.) zt=min(h-eps2,-1.*zt) ! if particle below ground -> reflection if (ngrid.ge.0) then cosfact=dxconst/cos((yt*dy+ylat0)*pi180) xt=xt+real(u*cosfact*real(ldt*ldirect),kind=dp) yt=yt+real(v*dyconst*real(ldt*ldirect),kind=dp) else if (ngrid.eq.-1) then ! around north pole xlon=xlon0+xt*dx ylat=ylat0+yt*dy call cll2xy(northpolemap,ylat,xlon,xpol,ypol) gridsize=1000.*cgszll(northpolemap,ylat,xlon) u=u/gridsize v=v/gridsize xpol=xpol+u*real(ldt*ldirect) ypol=ypol+v*real(ldt*ldirect) call cxy2ll(northpolemap,xpol,ypol,ylat,xlon) xt=(xlon-xlon0)/dx yt=(ylat-ylat0)/dy else if (ngrid.eq.-2) then ! around south pole xlon=xlon0+xt*dx ylat=ylat0+yt*dy call cll2xy(southpolemap,ylat,xlon,xpol,ypol) gridsize=1000.*cgszll(southpolemap,ylat,xlon) u=u/gridsize v=v/gridsize xpol=xpol+u*real(ldt*ldirect) ypol=ypol+v*real(ldt*ldirect) call cxy2ll(southpolemap,xpol,ypol,ylat,xlon) xt=(xlon-xlon0)/dx yt=(ylat-ylat0)/dy endif ! If global data are available, use cyclic boundary condition !************************************************************ if (xglobal) then if (xt.ge.real(nxmin1)) xt=xt-real(nxmin1) if (xt.lt.0.) xt=xt+real(nxmin1) if (xt.le.eps) xt=eps if (abs(xt-real(nxmin1)).le.eps) xt=real(nxmin1)-eps endif ! Check position: If trajectory outside model domain, terminate it !***************************************************************** if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. & (yt.ge.real(nymin1))) then nstop=3 return endif ! If particle above highest model level, set it back into the domain !******************************************************************* if (zt.ge.height(nz)) zt=height(nz)-100.*eps end subroutine advance