22 subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
23 usigold,vsigold,wsigold,nstop,xt,yt,zt,prob,icbt)
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
116 real :: usigold,vsigold,wsigold,r,rs
117 real :: uold,vold,wold,vdepo(maxspec)
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
133 integer :: idummy = -7
134 real :: settling = 0.
148 indzindicator(i)=.true.
154 depoindicator(ks)=.true.
166 nrand=int(
ran3(idummy)*
real(maxrand-1))+1
174 if (nglobal.and.(yt.gt.switchnorthg))
then
176 else if (sglobal.and.(yt.lt.switchsouthg))
then
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
195 if (abs(itime-memtime(1)).lt.abs(itime-memtime(2)))
then
205 xtn=(xt-xln(ngrid))*xresoln(ngrid)
206 ytn=(yt-yln(ngrid))*yresoln(ngrid)
229 if (hmix(i,j,1,k).gt.h) h=hmix(i,j,1,k)
233 tropop=tropopause(nix,njy,1,1)
238 if (hmixn(i,j,1,k,ngrid).gt.h) h=hmixn(i,j,1,k,ngrid)
242 tropop=tropopausen(nix,njy,1,1,ngrid)
261 if (method.eq.1)
then
262 ldt=min(ldt,abs(lsynctime-itimec+itime))
263 itimec=itimec+ldt*ldirect
266 itimec=itime+lsynctime
289 if (height(i).gt.zt)
then
302 if (indzindicator(i))
then
320 dz=1./(height(indzp)-height(indz))
321 dz1=(zt-height(indz))*dz
322 dz2=(height(indzp)-zt)*dz
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)
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)
354 up=ru*up+rannumb(nrand)*sigu*sqrt(1.-ru**2)
356 if (dt/tlv.lt..5)
then
357 vp=(1.-dt/tlv)*vp+rannumb(nrand+1)*sigv*sqrt(2.*dt/tlv)
360 vp=rv*vp+rannumb(nrand+1)*sigv*sqrt(1.-rv**2)
365 if (nrand+ifine.gt.maxrand) nrand=1
380 if (dtftlw.lt..5)
then
381 wp=((1.-dtftlw)*wp+rannumb(nrand+i)*sqrt(2.*dtftlw) &
382 +dtf*(dsigwdz+rhoaux*sigw))*
real(icbt)
385 wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2) &
386 +tlw*(1.-rw)*(dsigwdz+rhoaux*sigw))*
real(icbt)
391 wp=(rw*wp+rannumb(nrand+i)*sqrt(1.-rw**2)*sigw &
392 +tlw*(1.-rw)*(dsigw2dz+rhoaux*sigw**2))*
real(icbt)
400 if (abs(delz).gt.h) delz=mod(delz,h)
406 if (delz.lt.-zt)
then
409 else if (delz.gt.(h-zt))
then
429 ldt=int(min(tlw,h/max(2.*abs(wp*sigw),1.e-5), &
430 0.5/abs(dsigwdz))*ctl)
432 ldt=int(min(tlw,h/max(2.*abs(wp),1.e-5))*ctl)
442 if (mdomainfill.eq.0)
then
444 if (xmass(nrelpoint,nsp).gt.eps2) goto 887
446 887 nsp=min(nsp,nspec)
449 if (density(nsp).gt.0.) &
450 call
get_settling(itime,
real(xt),
real(yt),zt,nsp,settling)
462 dawsave=dawsave+up*dt
463 dcwsave=dcwsave+vp*dt
464 zt=zt+w*dt*
real(ldirect)
467 if (itimec.eq.itime+lsynctime) goto 99
500 if ((drydep).and.(zt.lt.2.*href))
then
502 if (drydepspec(ks))
then
503 if (depoindicator(ks))
then
513 prob(ks)=1.+(prob(ks)-1.)* &
514 exp(-vdepo(ks)*abs(dt)/(2.*href))
519 if (zt.lt.0.) zt=min(h-eps2,-1.*zt)
521 if (itimec.eq.(itime+lsynctime))
then
522 usig=0.5*(usigprof(indzp)+usigprof(indz))
523 vsig=0.5*(vsigprof(indzp)+vsigprof(indz))
524 wsig=0.5*(wsigprof(indzp)+wsigprof(indz))
565 ldt=abs(lsynctime-itimec+itime)
568 if (zt.lt.tropop)
then
569 uxscale=sqrt(2.*d_trop/dt)
570 if (nrand+1.gt.maxrand) nrand=1
571 ux=rannumb(nrand)*uxscale
572 vy=rannumb(nrand+1)*uxscale
575 else if (zt.lt.tropop+1000.)
then
576 weight=(zt-tropop)/1000.
577 uxscale=sqrt(2.*d_trop/dt*(1.-weight))
578 if (nrand+2.gt.maxrand) nrand=1
579 ux=rannumb(nrand)*uxscale
580 vy=rannumb(nrand+1)*uxscale
581 wpscale=sqrt(2.*d_strat/dt*weight)
582 wp=rannumb(nrand+2)*wpscale+d_strat/1000.
585 if (nrand.gt.maxrand) nrand=1
588 wpscale=sqrt(2.*d_strat/dt)
589 wp=rannumb(nrand)*wpscale
600 if (mdomainfill.eq.0)
then
602 if (xmass(nrelpoint,nsp).gt.eps2) goto 888
604 888 nsp=min(nsp,nspec)
607 if (density(nsp).gt.0.) &
608 call
get_settling(itime,
real(xt),
real(yt),zt,nsp,settling)
615 dxsave=dxsave+(u+ux)*dt
616 dysave=dysave+(v+vy)*dt
617 zt=zt+(w+wp)*dt*
real(ldirect)
618 if (zt.lt.0.) zt=min(h-eps2,-1.*zt)
638 r=exp(-2.*
real(abs(lsynctime))/
real(lwindinterv))
640 if (nrand+2.gt.maxrand) nrand=1
641 usigold=r*usigold+rs*rannumb(nrand)*usig*turbmesoscale
642 vsigold=r*vsigold+rs*rannumb(nrand+1)*vsig*turbmesoscale
643 wsigold=r*wsigold+rs*rannumb(nrand+2)*wsig*turbmesoscale
645 dxsave=dxsave+usigold*
real(lsynctime)
646 dysave=dysave+vsigold*
real(lsynctime)
648 zt=zt+wsigold*
real(lsynctime)
649 if (zt.lt.0.) zt=-1.*zt
657 call
windalign(dxsave,dysave,dawsave,dcwsave,ux,vy)
661 cosfact=dxconst/cos((yt*dy+ylat0)*pi180)
662 xt=xt+
real(dxsave*cosfact*real(ldirect),kind=dp)
663 yt=yt+
real(dysave*dyconst*real(ldirect),kind=dp)
664 else if (ngrid.eq.-1)
then
667 call
cll2xy(northpolemap,ylat,xlon,xpol,ypol)
668 gridsize=1000.*
cgszll(northpolemap,ylat,xlon)
669 dxsave=dxsave/gridsize
670 dysave=dysave/gridsize
671 xpol=xpol+dxsave*
real(ldirect)
672 ypol=ypol+dysave*
real(ldirect)
673 call
cxy2ll(northpolemap,xpol,ypol,ylat,xlon)
676 else if (ngrid.eq.-2)
then
679 call
cll2xy(southpolemap,ylat,xlon,xpol,ypol)
680 gridsize=1000.*
cgszll(southpolemap,ylat,xlon)
681 dxsave=dxsave/gridsize
682 dysave=dysave/gridsize
683 xpol=xpol+dxsave*
real(ldirect)
684 ypol=ypol+dysave*
real(ldirect)
685 call
cxy2ll(southpolemap,xpol,ypol,ylat,xlon)
695 if (xt.ge.
real(nxmin1)) xt=xt-
real(nxmin1)
696 if (xt.lt.0.) xt=xt+
real(nxmin1)
697 if (xt.le.eps) xt=eps
698 if (abs(xt-
real(nxmin1)).le.eps) xt=
real(nxmin1)-eps
705 if ((xt.lt.0.).or.(xt.ge.
real(nxmin1)).or.(yt.lt.0.).or. &
706 (yt.ge.
real(nymin1))) then
714 if (zt.ge.height(nz)) zt=height(nz)-100.*eps
730 if (ldt.ne.abs(lsynctime))
return
737 if (abs(itime+ldt*ldirect).gt.abs(memtime(2)))
return
742 if (nglobal.and.(yt.gt.switchnorthg))
then
744 else if (sglobal.and.(yt.lt.switchsouthg))
then
749 if ((xt.gt.xln(j)+eps).and.(xt.lt.xrn(j)-eps).and. &
750 (yt.gt.yln(j)+eps).and.(yt.lt.yrn(j)-eps))
then
758 if (ngr.ne.ngrid)
return
764 xtn=(xt-xln(ngrid))*xresoln(ngrid)
765 ytn=(yt-yln(ngrid))*yresoln(ngrid)
794 if (mdomainfill.eq.0)
then
796 if (xmass(nrelpoint,nsp).gt.eps2) goto 889
798 889 nsp=min(nsp,nspec)
801 if (density(nsp).gt.0.) &
802 call
get_settling(itime+ldt,
real(xt),
real(yt),zt,nsp,settling)
819 zt=zt+w*
real(ldt*ldirect)
820 if (zt.lt.0.) zt=min(h-eps2,-1.*zt)
822 cosfact=dxconst/cos((yt*dy+ylat0)*pi180)
823 xt=xt+
real(u*cosfact*real(ldt*ldirect),kind=dp)
824 yt=yt+
real(v*dyconst*real(ldt*ldirect),kind=dp)
825 else if (ngrid.eq.-1)
then
828 call
cll2xy(northpolemap,ylat,xlon,xpol,ypol)
829 gridsize=1000.*
cgszll(northpolemap,ylat,xlon)
832 xpol=xpol+u*
real(ldt*ldirect)
833 ypol=ypol+v*
real(ldt*ldirect)
834 call
cxy2ll(northpolemap,xpol,ypol,ylat,xlon)
837 else if (ngrid.eq.-2)
then
840 call
cll2xy(southpolemap,ylat,xlon,xpol,ypol)
841 gridsize=1000.*
cgszll(southpolemap,ylat,xlon)
844 xpol=xpol+u*
real(ldt*ldirect)
845 ypol=ypol+v*
real(ldt*ldirect)
846 call
cxy2ll(southpolemap,xpol,ypol,ylat,xlon)
855 if (xt.ge.
real(nxmin1)) xt=xt-
real(nxmin1)
856 if (xt.lt.0.) xt=xt+
real(nxmin1)
857 if (xt.le.eps) xt=eps
858 if (abs(xt-
real(nxmin1)).le.eps) xt=
real(nxmin1)-eps
864 if ((xt.lt.0.).or.(xt.ge.
real(nxmin1)).or.(yt.lt.0.).or. &
865 (yt.ge.
real(nymin1))) then
873 if (zt.ge.height(nz)) zt=height(nz)-100.*eps
real function, public cgszll(strcmp, xlat, xlong)
subroutine windalign(u, v, ffap, ffcp, ux, vy)
subroutine interpol_all_nests(itime, xt, yt, zt)
subroutine, public cxy2ll(strcmp, x, y, xlat, xlong)
subroutine hanna_short(z)
subroutine interpol_all(itime, xt, yt, zt)
subroutine interpol_misslev_nests(n)
subroutine interpol_vdep_nests(level, vdepo)
subroutine interpol_misslev(n)
subroutine advance(itime, nrelpoint, ldt, up, vp, wp, usigold, vsigold, wsigold, nstop, xt, yt, zt, prob, icbt)
subroutine interpol_vdep(level, vdepo)
subroutine get_settling(itime, xt, yt, zt, nsp, settling)
subroutine, public cll2xy(strcmp, xlat, xlong, x, y)
subroutine interpol_wind_nests(itime, xt, yt, zt)
subroutine interpol_wind(itime, xt, yt, zt)
subroutine interpol_wind_short_nests(itime, xt, yt, zt)
subroutine interpol_wind_short(itime, xt, yt, zt)