SUBROUTINE drag_noro (nlon,nlev,dtime,paprs,pplay, e pmea,pstd, psig, pgam, pthe,ppic,pval, e kgwd,kgwdim,kdx,ktest, e t, u, v, s pulow, pvlow, pustr, pvstr, s d_t, d_u, d_v) c IMPLICIT none c====================================================================== c Auteur(s): F.Lott (LMD/CNRS) date: 19950201 c Objet: Frottement de la montagne Interface c====================================================================== c Arguments: c dtime---input-R- pas d'integration (s) c paprs---input-R-pression pour chaque inter-couche (en Pa) c pplay---input-R-pression pour le mileu de chaque couche (en Pa) c t-------input-R-temperature (K) c u-------input-R-vitesse horizontale (m/s) c v-------input-R-vitesse horizontale (m/s) c c d_t-----output-R-increment de la temperature c d_u-----output-R-increment de la vitesse u c d_v-----output-R-increment de la vitesse v c====================================================================== #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" c c ARGUMENTS c INTEGER nlon,nlev REAL dtime REAL paprs(klon,klev+1) REAL pplay(klon,klev) REAL pmea(nlon),pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon) REAL ppic(nlon),pval(nlon) REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon) REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev) REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev) c INTEGER i, k, kgwd, kgwdim, kdx(nlon), ktest(nlon) c c Variables locales: c REAL zgeom(klon,klev) REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev) REAL pt(klon,klev), pu(klon,klev), pv(klon,klev) REAL papmf(klon,klev),papmh(klon,klev+1) c c initialiser les variables de sortie (pour securite) c DO i = 1,klon pulow(i) = 0.0 pvlow(i) = 0.0 pustr(i) = 0.0 pvstr(i) = 0.0 ENDDO DO k = 1, klev DO i = 1, klon d_t(i,k) = 0.0 d_u(i,k) = 0.0 d_v(i,k) = 0.0 pdudt(i,k)=0.0 pdvdt(i,k)=0.0 pdtdt(i,k)=0.0 ENDDO ENDDO c c preparer les variables d'entree (attention: l'ordre des niveaux c verticaux augmente du haut vers le bas) c DO k = 1, klev DO i = 1, klon pt(i,k) = t(i,klev-k+1) pu(i,k) = u(i,klev-k+1) pv(i,k) = v(i,klev-k+1) papmf(i,k) = pplay(i,klev-k+1) ENDDO ENDDO DO k = 1, klev+1 DO i = 1, klon papmh(i,k) = paprs(i,klev-k+2) ENDDO ENDDO DO i = 1, klon zgeom(i,klev) = RD * pt(i,klev) . * LOG(papmh(i,klev+1)/papmf(i,klev)) ENDDO DO k = klev-1, 1, -1 DO i = 1, klon zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0 . * LOG(papmf(i,k+1)/papmf(i,k)) ENDDO ENDDO c c appeler la routine principale c CALL orodrag(klon,klev,kgwd,kgwdim,kdx,ktest, . dtime, . papmh, papmf, zgeom, . pt, pu, pv, . pmea, pstd, psig, pgam, pthe, ppic,pval, . pulow,pvlow, . pdudt,pdvdt,pdtdt) C DO k = 1, klev DO i = 1, klon d_u(i,klev+1-k) = dtime*pdudt(i,k) d_v(i,klev+1-k) = dtime*pdvdt(i,k) d_t(i,klev+1-k) = dtime*pdtdt(i,k) pustr(i) = pustr(i) . +rg*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k)) pvstr(i) = pvstr(i) . +rg*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k)) ENDDO ENDDO c RETURN END SUBROUTINE orodrag( nlon,nlev i , kgwd, kgwdim, kdx, ktest r , ptsphy r , paphm1,papm1,pgeom1,ptm1,pum1,pvm1 r , pmea, pstd, psig, pgamma, ptheta, ppic, pval c outputs r , pulow,pvlow r , pvom,pvol,pte ) c c c**** *gwdrag* - does the gravity wave parametrization. c c purpose. c -------- c c this routine computes the physical tendencies of the c prognostic variables u,v and t due to vertical transports by c subgridscale orographically excited gravity waves c c** interface. c ---------- c called from *callpar*. c c the routine takes its input from the long-term storage: c u,v,t and p at t-1. c c explicit arguments : c -------------------- c ==== inputs === c ==== outputs === c c implicit arguments : none c -------------------- c implicit logical (l) c c method. c ------- c c externals. c ---------- integer ismin, ismax external ismin, ismax c c reference. c ---------- c c author. c ------- c m.miller + b.ritter e.c.m.w.f. 15/06/86. c c f.lott + m. miller e.c.m.w.f. 22/11/94 c----------------------------------------------------------------------- c c #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" #include "YOEGWD.h" c----------------------------------------------------------------------- c c* 0.1 arguments c --------- c c real pte(nlon,nlev), * pvol(nlon,nlev), * pvom(nlon,nlev), * pulow(klon), * pvlow(klon) real pum1(nlon,nlev), * pvm1(nlon,nlev), * ptm1(nlon,nlev), * pmea(nlon),pstd(nlon),psig(nlon), * pgamma(nlon),ptheta(nlon),ppic(nlon),pval(nlon), * pgeom1(nlon,nlev), * papm1(nlon,nlev), * paphm1(nlon,nlev+1) c integer kdx(nlon),ktest(nlon) c----------------------------------------------------------------------- c c* 0.2 local arrays c ------------ integer isect(klon), * icrit(klon), * ikcrith(klon), * ikenvh(klon), * iknu(klon), * iknu2(klon), * ikcrit(klon), * ikhlim(klon) c real ztau(klon,klev+1), * zstab(klon,klev+1), * zvph(klon,klev+1), * zrho(klon,klev+1), * zri(klon,klev+1), * zpsi(klon,klev+1), * zzdep(klon,klev) real zdudt(klon), * zdvdt(klon), * zdtdt(klon), * zdedt(klon), * zvidis(klon), * ztfr(klon), * znu(klon), * zd1(klon), * zd2(klon), * zdmod(klon) c c------------------------------------------------------------------ c c* 1. initialization c -------------- c 100 continue c c ------------------------------------------------------------------ c c* 1.1 computational constants c ----------------------- c 110 continue c c ztmst=twodt c if(nstep.eq.nstart) ztmst=0.5*twodt klevm1=klev-1 ztmst=ptsphy zrtmst=1./ztmst c ------------------------------------------------------------------ c 120 continue c c ------------------------------------------------------------------ c c* 1.3 check whether row contains point for printing c --------------------------------------------- c 130 continue c c ------------------------------------------------------------------ c c* 2. precompute basic state variables. c* ---------- ----- ----- ---------- c* define low level wind, project winds in plane of c* low level wind, determine sector in which to take c* the variance and set indicator for critical levels. c 200 continue c c c call orosetup * ( nlon, nlev , ktest * , ikcrit, ikcrith, icrit, isect, ikhlim, ikenvh,iknu,iknu2 * , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd * , zrho , zri , zstab , ztau , zvph , zpsi, zzdep * , pulow, pvlow * , ptheta,pgamma,pmea,ppic,pval,znu ,zd1, zd2, zdmod ) c c c c*********************************************************** c c c* 3. compute low level stresses using subcritical and c* supercritical forms.computes anisotropy coefficient c* as measure of orographic twodimensionality. c 300 continue c call gwstress * ( nlon , nlev * , ikcrit, isect, ikhlim, ktest, ikcrith, icrit, ikenvh, iknu * , zrho , zstab, zvph , pstd, psig, pmea, ppic, pval * , ztfr , ztau * , pgeom1,pgamma,zd1,zd2,zdmod,znu) c c c* 4. compute stress profile. c* ------- ------ -------- c 400 continue c c call gwprofil * ( nlon , nlev * , kgwd , kdx , ktest * , ikcrit, ikcrith, icrit , ikenvh, iknu * ,iknu2 , paphm1, zrho , zstab , ztfr , zvph * , zri , ztau , ztauf c * , zdmod , znu , psig , pgamma , pstd , ppic , pval) c c c* 5. compute tendencies. c* ------------------- c 500 continue c c explicit solution at all levels for the gravity wave c implicit solution for the blocked levels do 510 jl=kidia,kfdia zvidis(jl)=0.0 zdudt(jl)=0.0 zdvdt(jl)=0.0 zdtdt(jl)=0.0 510 continue c ilevp1=klev+1 c c do 524 jk=1,klev c c do 523 jl=1,kgwd ji=kdx(jl) zdelp=paphm1(ji,jk+1)-paphm1(ji,jk) ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp) zdudt(ji)=(pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji) zdvdt(ji)=(pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji) c c controle des overshoots: c zforc=sqrt(zdudt(ji)**2+zdvdt(ji)**2)+1.E-12 ztend=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst+1.E-12 rover=0.25 if(zforc.ge.rover*ztend)then zdudt(ji)=rover*ztend/zforc*zdudt(ji) zdvdt(ji)=rover*ztend/zforc*zdvdt(ji) endif c c fin du controle des overshoots c if(jk.ge.ikenvh(ji)) then zb=1.0-0.18*pgamma(ji)-0.04*pgamma(ji)**2 zc=0.48*pgamma(ji)+0.3*pgamma(ji)**2 zconb=2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji)) zabsv=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2. zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2 ratio=(cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji,jk))**2)/ * (pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2) zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv c c simplement oppose au vent c zdudt(ji)=-pum1(ji,jk)/ztmst zdvdt(ji)=-pvm1(ji,jk)/ztmst c c projection dans la direction de l'axe principal de l'orographie cmod zdudt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.) cmod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.)) cmod * *cos(ptheta(ji)*rpi/180.)/ztmst cmod zdvdt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.) cmod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.)) cmod * *sin(ptheta(ji)*rpi/180.)/ztmst zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet)) zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet)) end if pvom(ji,jk)=zdudt(ji) pvol(ji,jk)=zdvdt(ji) zust=pum1(ji,jk)+ztmst*zdudt(ji) zvst=pvm1(ji,jk)+ztmst*zdvdt(ji) zdis=0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2) zdedt(ji)=zdis/ztmst zvidis(ji)=zvidis(ji)+zdis*zdelp zdtdt(ji)=zdedt(ji)/rcpd c pte(ji,jk)=zdtdt(ji) c c ENCORE UN TRUC POUR EVITER LES EXPLOSIONS c pte(ji,jk)=0.0 523 continue 524 continue c c return end SUBROUTINE orosetup * ( nlon , nlev , ktest * , kkcrit, kkcrith, kcrit, ksect , kkhlim * , kkenvh, kknu , kknu2 * , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd * , prho , pri , pstab , ptau , pvph ,ppsi, pzdep * , pulow , pvlow * , ptheta, pgamma, pmea, ppic, pval * , pnu , pd1 , pd2 ,pdmod ) c c**** *gwsetup* c c purpose. c -------- c c** interface. c ---------- c from *orodrag* c c explicit arguments : c -------------------- c ==== inputs === c ==== outputs === c c implicit arguments : none c -------------------- c c method. c ------- c c c externals. c ---------- c c c reference. c ---------- c c see ecmwf research department documentation of the "i.f.s." c c author. c ------- c c modifications. c -------------- c f.lott for the new-gwdrag scheme november 1993 c c----------------------------------------------------------------------- implicit logical (l) c #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" #include "YOEGWD.h" c----------------------------------------------------------------------- c c* 0.1 arguments c --------- c integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),ksect(nlon), * kkhlim(nlon),ktest(nlon),kkenvh(nlon) c real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev), * pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev), * prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1), * ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1), * pzdep(nlon,klev) real pulow(nlon),pvlow(nlon),ptheta(nlon),pgamma(nlon),pnu(nlon), * pd1(nlon),pd2(nlon),pdmod(nlon) real pstd(nlon),pmea(nlon),ppic(nlon),pval(nlon) c c----------------------------------------------------------------------- c c* 0.2 local arrays c ------------ c c logical ll1(klon,klev+1) integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon), * kentp(klon),ncount(klon) c real zhcrit(klon,klev),zvpf(klon,klev), * zdp(klon,klev) real znorm(klon),zpsi(klon),zb(klon),zc(klon), * zulow(klon),zvlow(klon),znup(klon),znum(klon) c c ------------------------------------------------------------------ c c* 1. initialization c -------------- c c print *,' entree gwsetup' 100 continue c c ------------------------------------------------------------------ c c* 1.1 computational constants c ----------------------- c 110 continue c ilevm1=klev-1 ilevm2=klev-2 ilevh =klev/3 c zcons1=1./rd cold zcons2=g**2/cpd zcons2=rg**2/rcpd cold zcons3=1.5*api zcons3=1.5*rpi c c c ------------------------------------------------------------------ c c* 2. c -------------- c 200 continue c c ------------------------------------------------------------------ c c* 2.1 define low level wind, project winds in plane of c* low level wind, determine sector in which to take c* the variance and set indicator for critical levels. c c c do 2001 jl=kidia,kfdia kknu(jl) =klev kknu2(jl) =klev kknub(jl) =klev kknul(jl) =klev pgamma(jl) =max(pgamma(jl),gtsec) ll1(jl,klev+1)=.false. 2001 continue c c Ajouter une initialisation (L. Li, le 23fev99): c do jk=klev,ilevh,-1 do jl=kidia,kfdia ll1(jl,jk)= .FALSE. ENDDO ENDDO c c* define top of low level flow c ---------------------------- do 2002 jk=klev,ilevh,-1 do 2003 jl=kidia,kfdia lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr if(lo) then kkcrit(jl)=jk endif zhcrit(jl,jk)=ppic(jl) zhgeo=pgeom1(jl,jk)/rg ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk)) if(ll1(jl,jk).xor.ll1(jl,jk+1)) then kknu(jl)=jk endif if(.not.ll1(jl,ilevh))kknu(jl)=ilevh 2003 continue 2002 continue do 2004 jk=klev,ilevh,-1 do 2005 jl=kidia,kfdia zhcrit(jl,jk)=ppic(jl)-pval(jl) zhgeo=pgeom1(jl,jk)/rg ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk)) if(ll1(jl,jk).xor.ll1(jl,jk+1)) then kknu2(jl)=jk endif if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh 2005 continue 2004 continue do 2006 jk=klev,ilevh,-1 do 2007 jl=kidia,kfdia zhcrit(jl,jk)=amax1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl)) zhgeo=pgeom1(jl,jk)/rg ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk)) if(ll1(jl,jk).xor.ll1(jl,jk+1)) then kknub(jl)=jk endif if(.not.ll1(jl,ilevh))kknub(jl)=ilevh 2007 continue 2006 continue c do 2010 jl=kidia,kfdia kknu(jl)=min(kknu(jl),nktopg) kknu2(jl)=min(kknu2(jl),nktopg) kknub(jl)=min(kknub(jl),nktopg) kknul(jl)=klev 2010 continue c 210 continue c c cc* initialize various arrays c do 2107 jl=kidia,kfdia prho(jl,klev+1) =0.0 pstab(jl,klev+1) =0.0 pstab(jl,1) =0.0 pri(jl,klev+1) =9999.0 ppsi(jl,klev+1) =0.0 pri(jl,1) =0.0 pvph(jl,1) =0.0 pulow(jl) =0.0 pvlow(jl) =0.0 zulow(jl) =0.0 zvlow(jl) =0.0 kkcrith(jl) =klev kkenvh(jl) =klev kentp(jl) =klev kcrit(jl) =1 ncount(jl) =0 ll1(jl,klev+1) =.false. 2107 continue c c* define low-level flow c --------------------- c do 223 jk=klev,2,-1 do 222 jl=kidia,kfdia if(ktest(jl).eq.1) then zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1) prho(jl,jk)=2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1)) pstab(jl,jk)=2.*zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1))* * (1.-rcpd*prho(jl,jk)*(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk)) pstab(jl,jk)=max(pstab(jl,jk),gssec) endif 222 continue 223 continue c c******************************************************************** c c* define blocked flow c ------------------- do 2115 jk=klev,ilevh,-1 do 2116 jl=kidia,kfdia if(jk.ge.kknub(jl).and.jk.le.kknul(jl)) then pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk)) pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk)) end if 2116 continue 2115 continue do 2110 jl=kidia,kfdia pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl))) pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl))) znorm(jl)=max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec) pvph(jl,klev+1)=znorm(jl) 2110 continue c c******* setup orography axes and define plane of profiles ******* c do 2112 jl=kidia,kfdia lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec) if(lo) then zu=pulow(jl)+2.*gvsec else zu=pulow(jl) endif zphi=atan(pvlow(jl)/zu) ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi zb(jl)=1.-0.18*pgamma(jl)-0.04*pgamma(jl)**2 zc(jl)=0.48*pgamma(jl)+0.3*pgamma(jl)**2 pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2) pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))*cos(ppsi(jl,klev+1)) pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2) 2112 continue c c ************ define flow in plane of lowlevel stress ************* c do 213 jk=1,klev do 212 jl=kidia,kfdia if(ktest(jl).eq.1) then zvt1 =pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk) zvt2 =-pvlow(jl)*pum1(jl,jk)+pulow(jl)*pvm1(jl,jk) zvpf(jl,jk)=(zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl)) endif ptau(jl,jk) =0.0 pzdep(jl,jk) =0.0 ppsi(jl,jk) =0.0 ll1(jl,jk) =.false. 212 continue 213 continue do 215 jk=2,klev do 214 jl=kidia,kfdia if(ktest(jl).eq.1) then zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1) pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+ * (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1)) * /zdp(jl,jk) if(pvph(jl,jk).lt.gvsec) then pvph(jl,jk)=gvsec kcrit(jl)=jk endif endif 214 continue 215 continue c c c* 2.2 brunt-vaisala frequency and density at half levels. c 220 continue c do 2211 jk=ilevh,klev do 221 jl=kidia,kfdia if(ktest(jl).eq.1) then if(jk.ge.(kknub(jl)+1).and.jk.le.kknul(jl)) then zst=zcons2/ptm1(jl,jk)*(1.-rcpd*prho(jl,jk)* * (ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk)) pstab(jl,klev+1)=pstab(jl,klev+1)+zst*zdp(jl,jk) pstab(jl,klev+1)=max(pstab(jl,klev+1),gssec) prho(jl,klev+1)=prho(jl,klev+1)+paphm1(jl,jk)*2.*zdp(jl,jk) * *zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1)) endif endif 221 continue 2211 continue c do 2212 jl=kidia,kfdia pstab(jl,klev+1)=pstab(jl,klev+1)/(papm1(jl,kknul(jl)) * -papm1(jl,kknub(jl))) prho(jl,klev+1)=prho(jl,klev+1)/(papm1(jl,kknul(jl)) * -papm1(jl,kknub(jl))) zvar=pstd(jl) 2212 continue c c* 2.3 mean flow richardson number. c* and critical height for froude layer c 230 continue c do 232 jk=2,klev do 231 jl=kidia,kfdia if(ktest(jl).eq.1) then zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec) pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk) * /(rg*prho(jl,jk)*zdwind))**2 pri(jl,jk)=max(pri(jl,jk),grcrit) endif 231 continue 232 continue c c c* define top of 'envelope' layer c ---------------------------- do 233 jl=kidia,kfdia pnu (jl)=0.0 znum(jl)=0.0 233 continue do 234 jk=2,klev-1 do 234 jl=kidia,kfdia if(ktest(jl).eq.1) then if (jk.ge.kknub(jl)) then znum(jl)=pnu(jl) zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/ * max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec) zwind=max(sqrt(zwind**2),gvsec) zdelp=paphm1(jl,jk+1)-paphm1(jl,jk) zstabm=sqrt(max(pstab(jl,jk ),gssec)) zstabp=sqrt(max(pstab(jl,jk+1),gssec)) zrhom=prho(jl,jk ) zrhop=prho(jl,jk+1) pnu(jl) = pnu(jl) + (zdelp/rg)* * ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit) * .and.(kkenvh(jl).eq.klev)) * kkenvh(jl)=jk endif endif 234 continue c calculation of a dynamical mixing height for the breaking c of gravity waves: do 235 jl=kidia,kfdia znup(jl)=0.0 znum(jl)=0.0 235 continue do 236 jk=klev-1,2,-1 do 236 jl=kidia,kfdia if(ktest(jl).eq.1) then znum(jl)=znup(jl) zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/ * max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec) zwind=max(sqrt(zwind**2),gvsec) zdelp=paphm1(jl,jk+1)-paphm1(jl,jk) zstabm=sqrt(max(pstab(jl,jk ),gssec)) zstabp=sqrt(max(pstab(jl,jk+1),gssec)) zrhom=prho(jl,jk ) zrhop=prho(jl,jk+1) znup(jl) = znup(jl) + (zdelp/rg)* * ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind if((znum(jl).le.rpi/2.).and.(znup(jl).gt.rpi/2.) * .and.(kkcrith(jl).eq.klev)) * kkcrith(jl)=jk endif 236 continue do 237 jl=kidia,kfdia kkcrith(jl)=min0(kkcrith(jl),kknu2(jl)) kkcrith(jl)=max0(kkcrith(jl),ilevh*2) 237 continue c c directional info for flow blocking ************************* c do 251 jk=ilevh,klev do 252 jl=kidia,kfdia if(jk.ge.kkenvh(jl)) then lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec) if(lo) then zu=pum1(jl,jk)+2.*gvsec else zu=pum1(jl,jk) endif zphi=atan(pvm1(jl,jk)/zu) ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi end if 252 continue 251 continue c forms the vertical 'leakiness' ************************** alpha=3. do 254 jk=ilevh,klev do 253 jl=kidia,kfdia if(jk.ge.kkenvh(jl)) then zggeenv=amax1(1., * (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)-1))/2.) zggeom1=amax1(pgeom1(jl,jk),1.) zgvar=amax1(pstd(jl)*rg,1.) cmod pzdep(jl,jk)=sqrt((zggeenv-zggeom1)/(zggeom1+zgvar)) pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl, jk))/ * (pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl,klev)) end if 253 continue 254 continue 260 continue return end SUBROUTINE gwstress * ( nlon , nlev * , kkcrit, ksect, kkhlim, ktest, kkcrith, kcrit, kkenvh * , kknu * , prho , pstab , pvph , pstd, psig * , pmea , ppic , pval , ptfr , ptau * , pgeom1 , pgamma , pd1 , pd2 , pdmod , pnu ) c c**** *gwstress* c c purpose. c -------- c c** interface. c ---------- c call *gwstress* from *gwdrag* c c explicit arguments : c -------------------- c ==== inputs === c ==== outputs === c c implicit arguments : none c -------------------- c c method. c ------- c c c externals. c ---------- c c c reference. c ---------- c c see ecmwf research department documentation of the "i.f.s." c c author. c ------- c c modifications. c -------------- c f. lott put the new gwd on ifs 22/11/93 c c----------------------------------------------------------------------- implicit logical (l) #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" #include "YOEGWD.h" c----------------------------------------------------------------------- c c* 0.1 arguments c --------- c integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),ksect(nlon), * kkhlim(nlon),ktest(nlon),kkenvh(nlon),kknu(nlon) c real prho(nlon,nlev+1),pstab(nlon,nlev+1),ptau(nlon,nlev+1), * pvph(nlon,nlev+1),ptfr(nlon), * pgeom1(nlon,nlev),pstd(nlon) c real pd1(nlon),pd2(nlon),pnu(nlon),psig(nlon),pgamma(nlon) real pmea(nlon),ppic(nlon),pval(nlon) real pdmod(nlon) c c----------------------------------------------------------------------- c c* 0.2 local arrays c ------------ c c----------------------------------------------------------------------- c c* 0.3 functions c --------- c ------------------------------------------------------------------ c c* 1. initialization c -------------- c 100 continue c c* 3.1 gravity wave stress. c 300 continue c c do 301 jl=kidia,kfdia if(ktest(jl).eq.1) then c effective mountain height above the blocked flow if(kkenvh(jl).eq.klev)then zblock=0.0 else zblock=(pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)+1))/2./rg endif zvar=ppic(jl)-pmea(jl) zeff=amax1(0.,zvar-zblock) ptau(jl,klev+1)=prho(jl,klev+1)*gkdrag*psig(jl)*zeff**2 * /4./pstd(jl)*pvph(jl,klev+1)*pdmod(jl)*sqrt(pstab(jl,klev+1)) c too small value of stress or low level flow include critical level c or low level flow: gravity wave stress nul. lo=(ptau(jl,klev+1).lt.gtsec).or.(kcrit(jl).ge.kknu(jl)) * .or.(pvph(jl,klev+1).lt.gvcrit) c if(lo) ptau(jl,klev+1)=0.0 else ptau(jl,klev+1)=0.0 endif 301 continue c return end SUBROUTINE GWPROFIL * ( NLON, NLEV * , kgwd ,kdx , ktest * , KKCRIT, KKCRITH, KCRIT , kkenvh, kknu,kknu2 * , PAPHM1, PRHO , PSTAB , PTFR , PVPH , PRI , PTAU * , ptauf ,pdmod , pnu , psig ,pgamma, pvar, ppic,pval) C**** *GWPROFIL* C C PURPOSE. C -------- C C** INTERFACE. C ---------- C FROM *GWDRAG* C C EXPLICIT ARGUMENTS : C -------------------- C ==== INPUTS === C ==== OUTPUTS === C C IMPLICIT ARGUMENTS : NONE C -------------------- C C METHOD: C ------- C THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS: C IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND C AND THE TOP OF THE BLOCKED LAYER (KKENVH). C IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE C BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR C NONLINEAR GRAVITY WAVE BREAKING. C ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL C LEVEL (KCRIT) OR WHEN IT BREAKS. C C C C EXTERNALS. C ---------- C C C REFERENCE. C ---------- C C SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S." C C AUTHOR. C ------- C C MODIFICATIONS. C -------------- C PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93) C----------------------------------------------------------------------- IMPLICIT LOGICAL (L) C C #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" #include "YOEGWD.h" C----------------------------------------------------------------------- C C* 0.1 ARGUMENTS C --------- C INTEGER KKCRIT(NLON),KKCRITH(NLON),KCRIT(NLON) * ,kdx(nlon),ktest(NLON) * ,kkenvh(NLON),kknu(NLON),kknu2(NLON) C REAL PAPHM1(NLON,NLEV+1), PSTAB(NLON,NLEV+1), * PRHO (NLON,NLEV+1), PVPH (NLON,NLEV+1), * PRI (NLON,NLEV+1), PTFR (NLON), PTAU(NLON,NLEV+1), * ptauf (NLON,nlev+1) REAL pdmod (NLON) , pnu (NLON) , psig(NLON), * pgamma(NLON) , pvar(NLON) , ppic(NLON), pval(NLON) C----------------------------------------------------------------------- C C* 0.2 LOCAL ARRAYS C ------------ C REAL ZDZ2 (KLON,KLEV) , ZNORM(KLON) , zoro(KLON) REAL ZTAU (KLON,KLEV+1) C C----------------------------------------------------------------------- C C* 1. INITIALIZATION C -------------- C c print *,' entree gwprofil' 100 CONTINUE C C C* COMPUTATIONAL CONSTANTS. C ------------- ---------- C ilevh=KLEV/3 C DO 400 ji=1,kgwd jl=kdx(ji) Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0) ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1) 400 CONTINUE C DO 430 JK=KLEV,2,-1 C C C* 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE C BLOCKING LAYER. 410 CONTINUE C DO 411 ji=1,kgwd jl=kdx(ji) IF(JK.GT.KKCRITH(JL)) THEN PTAU(JL,JK)=ZTAU(JL,KLEV+1) C ENDIF C IF(JK.EQ.KKCRITH(JL)) THEN ELSE PTAU(JL,JK)=GRAHILO*ZTAU(JL,KLEV+1) ENDIF 411 CONTINUE C C* 4.15 CONSTANT SHEAR STRESS UNTIL THE TOP OF THE C LOW LEVEL FLOW LAYER. 415 CONTINUE C C C* 4.2 WAVE DISPLACEMENT AT NEXT LEVEL. C 420 CONTINUE C DO 421 ji=1,kgwd jl=kdx(ji) IF(JK.LT.KKCRITH(JL)) THEN ZNORM(JL)=gkdrag*PRHO(JL,JK)*SQRT(PSTAB(JL,JK))*PVPH(JL,JK) * *zoro(jl) ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec) ENDIF 421 CONTINUE C C* 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT C* AND STRESS: BREAKING EVALUATION AND CRITICAL C LEVEL C DO 431 ji=1,kgwd jl=kdx(ji) IF(JK.LT.KKCRITH(JL)) THEN IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN PTAU(JL,JK)=0.0 ELSE ZSQR=SQRT(PRI(JL,JK)) ZALFA=SQRT(PSTAB(JL,JK)*ZDZ2(JL,JK))/PVPH(JL,JK) ZRIW=PRI(JL,JK)*(1.-ZALFA)/(1+ZALFA*ZSQR)**2 IF(ZRIW.LT.GRCRIT) THEN ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT ZB=1./GRCRIT+2./ZSQR ZALPHA=0.5*(-ZB+SQRT(ZDEL)) ZDZ2N=(PVPH(JL,JK)*ZALPHA)**2/PSTAB(JL,JK) PTAU(JL,JK)=ZNORM(JL)*ZDZ2N ELSE PTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK) ENDIF PTAU(JL,JK)=MIN(PTAU(JL,JK),PTAU(JL,JK+1)) ENDIF ENDIF 431 CONTINUE 430 CONTINUE 440 CONTINUE C REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL DO 530 ji=1,kgwd jl=kdx(ji) ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL)) ZTAU(JL,NSTRA)=PTAU(JL,NSTRA) 530 CONTINUE DO 531 JK=1,KLEV DO 532 ji=1,kgwd jl=kdx(ji) IF(JK.GT.KKCRITH(JL))THEN ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KLEV+1 ) ZDELPT=PAPHM1(JL,KKCRITH(JL))-PAPHM1(JL,KLEV+1 ) PTAU(JL,JK)=ZTAU(JL,KLEV+1 ) + . (ZTAU(JL,KKCRITH(JL))-ZTAU(JL,KLEV+1 ) )* . ZDELP/ZDELPT ENDIF 532 CONTINUE C REORGANISATION IN THE STRATOSPHERE DO 533 ji=1,kgwd jl=kdx(ji) IF(JK.LT.NSTRA)THEN ZDELP =PAPHM1(JL,NSTRA) ZDELPT=PAPHM1(JL,JK) PTAU(JL,JK)=ZTAU(JL,NSTRA)*ZDELPT/ZDELP ENDIF 533 CONTINUE C REORGANISATION IN THE TROPOSPHERE DO 534 ji=1,kgwd jl=kdx(ji) IF(JK.LT.KKCRITH(JL).AND.JK.GT.NSTRA)THEN ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KKCRITH(JL)) ZDELPT=PAPHM1(JL,NSTRA)-PAPHM1(JL,KKCRITH(JL)) PTAU(JL,JK)=ZTAU(JL,KKCRITH(JL)) + * (ZTAU(JL,NSTRA)-ZTAU(JL,KKCRITH(JL)))*ZDELP/ZDELPT ENDIF 534 CONTINUE 531 CONTINUE RETURN END SUBROUTINE lift_noro (nlon,nlev,dtime,paprs,pplay, e plat,pmea,pstd, psig, pgam, pthe, ppic,pval, e kgwd,kgwdim,kdx,ktest, e t, u, v, s pulow, pvlow, pustr, pvstr, s d_t, d_u, d_v) c IMPLICIT none c====================================================================== c Auteur(s): F.Lott (LMD/CNRS) date: 19950201 c Objet: Frottement de la montagne Interface c====================================================================== c Arguments: c dtime---input-R- pas d'integration (s) c paprs---input-R-pression pour chaque inter-couche (en Pa) c pplay---input-R-pression pour le mileu de chaque couche (en Pa) c t-------input-R-temperature (K) c u-------input-R-vitesse horizontale (m/s) c v-------input-R-vitesse horizontale (m/s) c c d_t-----output-R-increment de la temperature c d_u-----output-R-increment de la vitesse u c d_v-----output-R-increment de la vitesse v c====================================================================== #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" c c ARGUMENTS c INTEGER nlon,nlev REAL dtime REAL paprs(klon,klev+1) REAL pplay(klon,klev) REAL plat(nlon),pmea(nlon) REAL pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon) REAL ppic(nlon),pval(nlon) REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon) REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev) REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev) c INTEGER i, k, kgwd, kgwdim, kdx(nlon), ktest(nlon) c c Variables locales: c REAL zgeom(klon,klev) REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev) REAL pt(klon,klev), pu(klon,klev), pv(klon,klev) REAL papmf(klon,klev),papmh(klon,klev+1) c c initialiser les variables de sortie (pour securite) c DO i = 1,klon pulow(i) = 0.0 pvlow(i) = 0.0 pustr(i) = 0.0 pvstr(i) = 0.0 ENDDO DO k = 1, klev DO i = 1, klon d_t(i,k) = 0.0 d_u(i,k) = 0.0 d_v(i,k) = 0.0 pdudt(i,k)=0.0 pdvdt(i,k)=0.0 pdtdt(i,k)=0.0 ENDDO ENDDO c c preparer les variables d'entree (attention: l'ordre des niveaux c verticaux augmente du haut vers le bas) c DO k = 1, klev DO i = 1, klon pt(i,k) = t(i,klev-k+1) pu(i,k) = u(i,klev-k+1) pv(i,k) = v(i,klev-k+1) papmf(i,k) = pplay(i,klev-k+1) ENDDO ENDDO DO k = 1, klev+1 DO i = 1, klon papmh(i,k) = paprs(i,klev-k+2) ENDDO ENDDO DO i = 1, klon zgeom(i,klev) = RD * pt(i,klev) . * LOG(papmh(i,klev+1)/papmf(i,klev)) ENDDO DO k = klev-1, 1, -1 DO i = 1, klon zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0 . * LOG(papmf(i,k+1)/papmf(i,k)) ENDDO ENDDO c c appeler la routine principale c CALL OROLIFT(klon,klev,kgwd,kgwdim,kdx,ktest, . dtime, . papmh, papmf, zgeom, . pt, pu, pv, . plat,pmea, pstd, psig, pgam, pthe, ppic,pval, . pulow,pvlow, . pdudt,pdvdt,pdtdt) C DO k = 1, klev DO i = 1, klon d_u(i,klev+1-k) = dtime*pdudt(i,k) d_v(i,klev+1-k) = dtime*pdvdt(i,k) d_t(i,klev+1-k) = dtime*pdtdt(i,k) pustr(i) = pustr(i) . +RG*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k)) pvstr(i) = pvstr(i) . +RG*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k)) ENDDO ENDDO c RETURN END SUBROUTINE OROLIFT( NLON,NLEV I , KGWD, KGWDIM, KDX, KTEST R , PTSPHY R , PAPHM1,PAPM1,PGEOM1,PTM1,PUM1,PVM1 R , PLAT R , PMEA, PVAROR, PSIG, PGAMMA, PTHETA,ppic,pval C OUTPUTS R , PULOW,PVLOW R , PVOM,PVOL,PTE ) C C**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT. C C PURPOSE. C -------- C C** INTERFACE. C ---------- C CALLED FROM *lift_noro C ---------- C C AUTHOR. C ------- C F.LOTT LMD 22/11/95 C IMPLICIT LOGICAL (L) C C #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" #include "YOEGWD.h" C----------------------------------------------------------------------- C C* 0.1 ARGUMENTS C --------- C C REAL PTE(NLON,NLEV), * PVOL(NLON,NLEV), * PVOM(NLON,NLEV), * PULOW(NLON), * PVLOW(NLON) REAL PUM1(NLON,NLEV), * PVM1(NLON,NLEV), * PTM1(NLON,NLEV), * PLAT(NLON),PMEA(NLON), * PVAROR(NLON),PSIG(NLON),PGAMMA(NLON), * PTHETA(NLON),ppic(NLON),pval(NLON), * PGEOM1(NLON,NLEV), * PAPM1(NLON,NLEV), * PAPHM1(NLON,NLEV+1) C INTEGER KDX(NLON),KTEST(NLON) C----------------------------------------------------------------------- C C* 0.2 LOCAL ARRAYS C ------------ INTEGER ISECT(KLON), * ICRIT(KLON), * IKCRITH(KLON), * IKenvh(klon), * IKNUB(klon), * IKNUL(klon), * IKCRIT(KLON), * IKHLIM(KLON) LOGICAL LL1(KLON,KLEV+1) C REAL ZTAU(KLON,KLEV+1), * ZTAV(KLON,KLEV+1), * ZSTAB(KLON,KLEV+1), * ZVPH(KLON,KLEV+1), * ZRHO(KLON,KLEV+1), * ZRI(KLON,KLEV+1), * ZpsI(KLON,KLEV+1), * Zzdep(KLON,KLEV) REAL ZDUDT(KLON), * ZDVDT(KLON), * ZDTDT(KLON), * ZDEDT(KLON), * ZVIDIS(KLON), * ZTFR(KLON), * Znu(KLON), * Zd1(KLON), * Zd2(KLON), * Zdmod(KLON) REAL ZHCRIT(KLON,KLEV) C----------------------------------------------------------------------- C C* 1.1 INITIALIZATIONS C --------------- LIFTHIGH=.FALSE. IF(NLON.NE.KLON.OR.NLEV.NE.KLEV)STOP ZCONS1=1./RD KLEVM1=KLEV-1 ZTMST=PTSPHY ZRTMST=1./ZTMST ZPI=ACOS(-1.) C DO 1001 JL=kidia,kfdia ZRHO(JL,KLEV+1) =0.0 PULOW(JL) =0.0 PVLOW(JL) =0.0 iknub(JL) =klev iknul(JL) =klev ilevh=klev/3 ll1(jl,klev+1)=.false. DO 1000 JK=1,KLEV PVOM(JL,JK)=0.0 PVOL(JL,JK)=0.0 PTE (JL,JK)=0.0 1000 CONTINUE 1001 CONTINUE C C* 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF C* LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE C* THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS. C C C DO 2006 JK=KLEV,1,-1 DO 2007 JL=kidia,kfdia IF(KTEST(JL).EQ.1) THEN ZHCRIT(JL,JK)=amax1(Ppic(JL)-pmea(JL),100.) ZHGEO=PGEOM1(JL,JK)/RG ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK)) IF(ll1(JL,JK).XOR.ll1(JL,JK+1)) THEN iknub(JL)=JK ENDIF ENDIF 2007 CONTINUE 2006 CONTINUE C do 2010 jl=kidia,kfdia IF(KTEST(JL).EQ.1) THEN iknub(jl)=max(iknub(jl),klev/2) iknul(jl)=max(iknul(jl),2*klev/3) if(iknub(jl).gt.nktopg) iknub(jl)=nktopg if(iknub(jl).eq.nktopg) iknul(jl)=klev if(iknub(jl).eq.iknul(jl)) iknub(jl)=iknul(jl)-1 ENDIF 2010 continue C do 2011 jl=kidia,kfdia C IF(KTEST(JL).EQ.1) THEN C print *,' iknul= ',iknul(jl),' iknub=',iknub(jl) C ENDIF C2011 continue C PRINT *,' DANS OROLIFT: 2010' DO 223 JK=KLEV,2,-1 DO 222 JL=kidia,kfdia ZRHO(JL,JK)=2.*PAPHM1(JL,JK)*ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1)) 222 CONTINUE 223 CONTINUE C PRINT *,' DANS OROLIFT: 223' C******************************************************************** C C* DEFINE LOW LEVEL FLOW C ------------------- DO 2115 JK=klev,1,-1 DO 2116 JL=kidia,kfdia IF(KTEST(JL).EQ.1) THEN if(jk.ge.iknub(jl).and.jk.le.iknul(jl)) then pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK)) pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK)) zrho(JL,klev+1)=zrho(JL,klev+1) * +zrho(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK)) end if ENDIF 2116 CONTINUE 2115 CONTINUE DO 2110 JL=kidia,kfdia IF(KTEST(JL).EQ.1) THEN pulow(JL)=pulow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl))) pvlow(JL)=pvlow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl))) zrho(JL,klev+1)=zrho(JL,klev+1) * /(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl))) ENDIF 2110 CONTINUE 200 CONTINUE C*********************************************************** C C* 3. COMPUTE MOUNTAIN LIFT C 300 CONTINUE C DO 301 JL=kidia,kfdia IF(KTEST(JL).EQ.1) THEN ZTAU(JL,KLEV+1)= - GKLIFT*ZRHO(JL,KLEV+1)*2.*ROMEGA* C * (2*PVAROR(JL)+PMEA(JL))* * 2*PVAROR(JL)* * SIN(ZPI/180.*PLAT(JL))*PVLOW(JL) ZTAV(JL,KLEV+1)= GKLIFT*ZRHO(JL,KLEV+1)*2.*ROMEGA* C * (2*PVAROR(JL)+PMEA(JL))* * 2*PVAROR(JL)* * SIN(ZPI/180.*PLAT(JL))*PULOW(JL) ELSE ZTAU(JL,KLEV+1)=0.0 ZTAV(JL,KLEV+1)=0.0 ENDIF 301 CONTINUE C C* 4. COMPUTE LIFT PROFILE C* -------------------- C 400 CONTINUE DO 401 JK=1,KLEV DO 401 JL=kidia,kfdia IF(KTEST(JL).EQ.1) THEN ZTAU(JL,JK)=ZTAU(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1) ZTAV(JL,JK)=ZTAV(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1) ELSE ZTAU(JL,JK)=0.0 ZTAV(JL,JK)=0.0 ENDIF 401 CONTINUE C C C* 5. COMPUTE TENDENCIES. C* ------------------- IF(LIFTHIGH)THEN C 500 CONTINUE C PRINT *,' DANS OROLIFT: 500' C C EXPLICIT SOLUTION AT ALL LEVELS C DO 524 JK=1,klev DO 523 JL=KIDIA,KFDIA IF(KTEST(JL).EQ.1) THEN ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK) ZDUDT(JL)=-RG*(ZTAU(JL,JK+1)-ZTAU(JL,JK))/ZDELP ZDVDT(JL)=-RG*(ZTAV(JL,JK+1)-ZTAV(JL,JK))/ZDELP ENDIF 523 CONTINUE 524 CONTINUE C C PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY C DO 530 JK=1,klev DO 530 JL=KIDIA,KFDIA IF(KTEST(JL).EQ.1) THEN ZSLOW=SQRT(PULOW(JL)**2+PVLOW(JL)**2) ZSQUA=AMAX1(SQRT(PUM1(JL,JK)**2+PVM1(JL,JK)**2),GVSEC) ZSCAV=-ZDUDT(JL)*PVM1(JL,JK)+ZDVDT(JL)*PUM1(JL,JK) IF(ZSQUA.GT.GVSEC)THEN PVOM(JL,JK)=-ZSCAV*PVM1(JL,JK)/ZSQUA**2 PVOL(JL,JK)= ZSCAV*PUM1(JL,JK)/ZSQUA**2 ELSE PVOM(JL,JK)=0.0 PVOL(JL,JK)=0.0 ENDIF ZSQUA=SQRT(PUM1(JL,JK)**2+PUM1(JL,JK)**2) IF(ZSQUA.LT.ZSLOW)THEN PVOM(JL,JK)=ZSQUA/ZSLOW*PVOM(JL,JK) PVOL(JL,JK)=ZSQUA/ZSLOW*PVOL(JL,JK) ENDIF ENDIF 530 CONTINUE C C 6. LOW LEVEL LIFT, SEMI IMPLICIT: C ---------------------------------- ELSE DO 601 JL=KIDIA,KFDIA IF(KTEST(JL).EQ.1) THEN DO JK=KLEV,IKNUB(JL),-1 ZBET=GKLIFT*2.*ROMEGA*SIN(ZPI/180.*PLAT(JL))*ztmst* * (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL, JK))/ * (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL,KLEV)) ZDUDT(JL)=-PUM1(JL,JK)/ztmst/(1+ZBET**2) ZDVDT(JL)=-PVM1(JL,JK)/ztmst/(1+ZBET**2) PVOM(JL,JK)= ZBET**2*ZDUDT(JL) - ZBET *ZDVDT(JL) PVOL(JL,JK)= ZBET *ZDUDT(JL) + ZBET**2*ZDVDT(JL) ENDDO ENDIF 601 CONTINUE ENDIF RETURN END SUBROUTINE SUGWD(NLON,NLEV,paprs,pplay) C C**** *SUGWD* INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG C C PURPOSE. C -------- C INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE C GRAVITY WAVE DRAG PARAMETRIZATION. C C** INTERFACE. C ---------- C CALL *SUGWD* FROM *SUPHEC* C ----- ------ C C EXPLICIT ARGUMENTS : C -------------------- C PSIG : VERTICAL COORDINATE TABLE C NLEV : NUMBER OF MODEL LEVELS C C IMPLICIT ARGUMENTS : C -------------------- C COMMON YOEGWD C C METHOD. C ------- C SEE DOCUMENTATION C C EXTERNALS. C ---------- C NONE C C REFERENCE. C ---------- C ECMWF Research Department documentation of the IFS C C AUTHOR. C ------- C MARTIN MILLER *ECMWF* C C MODIFICATIONS. C -------------- C ORIGINAL : 90-01-01 C ------------------------------------------------------------------ IMPLICIT LOGICAL(L) C C ----------------------------------------------------------------- #include "YOEGWD.h" C ---------------------------------------------------------------- C integer nlon,nlev REAL paprs(nlon,nlev+1) REAL pplay(nlon,nlev) C C* 1. SET THE VALUES OF THE PARAMETERS C -------------------------------- C 100 CONTINUE C PRINT *,' DANS SUGWD NLEV=',NLEV GHMAX=10000. C ZPR=100000. ZSTRA=0.1 ZSIGT=0.94 cold ZPR=80000. cold ZSIGT=0.85 C DO 110 JK=1,NLEV ZPM1R=pplay(nlon/2,jk)/paprs(nlon/2,1) IF(ZPM1R.GE.ZSIGT)THEN nktopg=JK ENDIF ZPM1R=pplay(nlon/2,jk)/paprs(nlon/2,1) IF(ZPM1R.GE.ZSTRA)THEN NSTRA=JK ENDIF 110 CONTINUE c c inversion car dans orodrag on compte les niveaux a l'envers nktopg=nlev-nktopg+1 nstra=nlev-nstra print *,' DANS SUGWD nktopg=', nktopg print *,' DANS SUGWD nstra=', nstra C GSIGCR=0.80 C GKDRAG=0.2 GRAHILO=1 GRCRIT=0.01 GFRCRIT=1.0 GKWAKE=0.50 C GKLIFT=0.50 GVCRIT =0.0 C C C ---------------------------------------------------------------- C C* 2. SET VALUES OF SECURITY PARAMETERS C --------------------------------- C 200 CONTINUE C GVSEC=0.10 GSSEC=1.E-12 C GTSEC=1.E-07 C C ---------------------------------------------------------------- C RETURN END