! $Id: orografi.F90 5160 2024-08-03 12:56:58Z evignon $ SUBROUTINE drag_noro(nlon, nlev, dtime, paprs, pplay, pmea, pstd, psig, pgam, & pthe, ppic, pval, kgwd, kdx, ktest, t, u, v, pulow, pvlow, pustr, pvstr, & d_t, d_u, d_v) USE dimphy USE lmdz_yomcst IMPLICIT NONE ! ====================================================================== ! Auteur(s): F.Lott (LMD/CNRS) date: 19950201 ! Objet: Frottement de la montagne Interface ! ====================================================================== ! Arguments: ! dtime---input-R- pas d'integration (s) ! paprs---input-R-pression pour chaque inter-couche (en Pa) ! pplay---input-R-pression pour le mileu de chaque couche (en Pa) ! t-------input-R-temperature (K) ! u-------input-R-vitesse horizontale (m/s) ! v-------input-R-vitesse horizontale (m/s) ! d_t-----output-R-increment de la temperature ! d_u-----output-R-increment de la vitesse u ! d_v-----output-R-increment de la vitesse v ! ====================================================================== ! ARGUMENTS 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) INTEGER i, k, kgwd, kdx(nlon), ktest(nlon) ! Variables locales: 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) ! initialiser les variables de sortie (pour securite) DO i = 1, klon pulow(i) = 0.0 pvlow(i) = 0.0 pustr(i) = 0.0 pvstr(i) = 0.0 END DO 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 END DO END DO ! preparer les variables d'entree (attention: l'ordre des niveaux ! verticaux augmente du haut vers le bas) 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) END DO END DO DO k = 1, klev + 1 DO i = 1, klon papmh(i, k) = paprs(i, klev - k + 2) END DO END DO DO i = 1, klon zgeom(i, klev) = rd * pt(i, klev) * log(papmh(i, klev + 1) / papmf(i, klev)) END DO 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)) END DO END DO ! appeler la routine principale CALL orodrag(klon, klev, kgwd, kdx, ktest, dtime, papmh, papmf, zgeom, pt, & pu, pv, pmea, pstd, psig, pgam, pthe, ppic, pval, pulow, pvlow, pdudt, & pdvdt, pdtdt) 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) & ! IM BUG . ! +rg*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k)) + pdudt(i, k) * (papmh(i, k + 1) - papmh(i, k)) / rg pvstr(i) = pvstr(i) & ! IM BUG . ! +rg*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k)) + pdvdt(i, k) * (papmh(i, k + 1) - papmh(i, k)) / rg END DO END DO END SUBROUTINE drag_noro SUBROUTINE orodrag(nlon, nlev, kgwd, kdx, ktest, ptsphy, paphm1, papm1, & pgeom1, ptm1, pum1, pvm1, pmea, pstd, psig, pgamma, ptheta, ppic, pval & ! outputs , pulow, pvlow, pvom, pvol, pte) USE dimphy USE lmdz_YOEGWD, ONLY: GFRCRIT, GKWAKE, GRCRIT, GVCRIT, GKDRAG, GKLIFT, GHMAX, GRAHILO, GSIGCR, NKTOPG, NSTRA, GSSEC, GTSEC, GVSEC, & GWD_RANDO_RUWMAX, gwd_rando_sat, GWD_FRONT_RUWMAX, gwd_front_sat USE lmdz_libmath, ONLY: ismax, ismin USE lmdz_yomcst IMPLICIT NONE ! **** *gwdrag* - does the gravity wave parametrization. ! purpose. ! -------- ! this routine computes the physical tendencies of the ! prognostic variables u,v and t due to vertical transports by ! subgridscale orographically excited gravity waves ! ** interface. ! ---------- ! called from *callpar*. ! the routine takes its input from the long-term storage: ! u,v,t and p at t-1. ! explicit arguments : ! -------------------- ! ==== inputs === ! ==== outputs === ! implicit arguments : none ! -------------------- ! implicit LOGICAL (l) ! author. ! ------- ! m.miller + b.ritter e.c.m.w.f. 15/06/86. ! f.lott + m. miller e.c.m.w.f. 22/11/94 ! ----------------------------------------------------------------------- ! * 0.1 arguments ! --------- ! ym integer nlon, nlev, klevm1 INTEGER nlon, nlev INTEGER kgwd, jl, ilevp1, jk, ji REAL zdelp, ztemp, zforc, ztend REAL rover, zb, zc, zconb, zabsv REAL zzd1, ratio, zbet, zust, zvst, zdis 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) INTEGER kdx(nlon), ktest(nlon) ! ----------------------------------------------------------------------- ! * 0.2 local arrays ! ------------ INTEGER isect(klon), icrit(klon), ikcrith(klon), ikenvh(klon), iknu(klon), & iknu2(klon), ikcrit(klon), ikhlim(klon) REAL ztau(klon, klev + 1), ztauf(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), & znu(klon), zd1(klon), zd2(klon), zdmod(klon) REAL ztmst, ptsphy, zrtmst ! ------------------------------------------------------------------ ! * 1. initialization ! -------------- ! ------------------------------------------------------------------ ! * 1.1 computational constants ! ----------------------- ! ztmst=twodt ! IF(nstep.EQ.nstart) ztmst=0.5*twodt ! ym klevm1=klev-1 ztmst = ptsphy zrtmst = 1. / ztmst ! ------------------------------------------------------------------ ! * 1.3 check whether row contains point for printing ! --------------------------------------------- ! ------------------------------------------------------------------ ! * 2. precompute basic state variables. ! * ---------- ----- ----- ---------- ! * define low level wind, project winds in plane of ! * low level wind, determine sector in which to take ! * the variance and set indicator for critical levels. CALL orosetup(nlon, ktest, ikcrit, ikcrith, icrit, 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) ! *********************************************************** ! * 3. compute low level stresses using subcritical and ! * supercritical forms.computes anisotropy coefficient ! * as measure of orographic twodimensionality. CALL gwstress(nlon, nlev, ktest, icrit, ikenvh, iknu, zrho, zstab, zvph, & pstd, psig, pmea, ppic, ztau, pgeom1, zdmod) ! * 4. compute stress profile. ! * ------- ------ -------- CALL gwprofil(nlon, nlev, kgwd, kdx, ktest, ikcrith, icrit, paphm1, zrho, & zstab, zvph, zri, ztau, zdmod, psig, pstd) ! * 5. compute tendencies. ! * ------------------- ! explicit solution at all levels for the gravity wave ! implicit solution for the blocked levels DO jl = kidia, kfdia zvidis(jl) = 0.0 zdudt(jl) = 0.0 zdvdt(jl) = 0.0 zdtdt(jl) = 0.0 END DO ilevp1 = klev + 1 DO jk = 1, klev ! do 523 jl=1,kgwd ! ji=kdx(jl) ! Modif vectorisation 02/04/2004 DO ji = kidia, kfdia IF (ktest(ji)==1) THEN 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) ! controle des overshoots: 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>=rover * ztend) THEN zdudt(ji) = rover * ztend / zforc * zdudt(ji) zdvdt(ji) = rover * ztend / zforc * zdvdt(ji) END IF ! fin du controle des overshoots IF (jk>=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 ! simplement oppose au vent zdudt(ji) = -pum1(ji, jk) / ztmst zdvdt(ji) = -pvm1(ji, jk) / ztmst ! projection dans la direction de l'axe principal de l'orographie ! mod zdudt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.) ! mod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.)) ! mod * *cos(ptheta(ji)*rpi/180.)/ztmst ! mod zdvdt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.) ! mod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.)) ! mod * *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 ! pte(ji,jk)=zdtdt(ji) ! ENCORE UN TRUC POUR EVITER LES EXPLOSIONS pte(ji, jk) = 0.0 END IF END DO END DO END SUBROUTINE orodrag SUBROUTINE orosetup(nlon, ktest, kkcrit, kkcrith, kcrit, 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) ! **** *gwsetup* ! purpose. ! -------- ! ** interface. ! ---------- ! from *orodrag* ! explicit arguments : ! -------------------- ! ==== inputs === ! ==== outputs === ! implicit arguments : none ! -------------------- ! method. ! ------- ! externals. ! ---------- ! reference. ! ---------- ! see ecmwf research department documentation of the "i.f.s." ! author. ! ------- ! modifications. ! -------------- ! f.lott for the new-gwdrag scheme november 1993 ! ----------------------------------------------------------------------- USE dimphy USE lmdz_YOEGWD, ONLY: GFRCRIT, GKWAKE, GRCRIT, GVCRIT, GKDRAG, GKLIFT, GHMAX, GRAHILO, GSIGCR, NKTOPG, NSTRA, GSSEC, GTSEC, GVSEC, & GWD_RANDO_RUWMAX, gwd_rando_sat, GWD_FRONT_RUWMAX, gwd_front_sat USE lmdz_yomcst IMPLICIT NONE ! ----------------------------------------------------------------------- ! * 0.1 arguments ! --------- INTEGER nlon INTEGER jl, jk REAL zdelp INTEGER kkcrit(nlon), kkcrith(nlon), kcrit(nlon), ktest(nlon), kkenvh(nlon) 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) ! ----------------------------------------------------------------------- ! * 0.2 local arrays ! ------------ INTEGER ilevm1, ilevm2, ilevh REAL zcons1, zcons2, zcons3, zhgeo REAL zu, zphi, zvt1, zvt2, zst, zvar, zdwind, zwind REAL zstabm, zstabp, zrhom, zrhop, alpha REAL zggeenv, zggeom1, zgvar LOGICAL lo LOGICAL ll1(klon, klev + 1) INTEGER kknu(klon), kknu2(klon), kknub(klon), kknul(klon), kentp(klon), & ncount(klon) REAL zhcrit(klon, klev), zvpf(klon, klev), zdp(klon, klev) REAL znorm(klon), zb(klon), zc(klon), zulow(klon), zvlow(klon), znup(klon), & znum(klon) ! ------------------------------------------------------------------ ! * 1. initialization ! -------------- ! PRINT *,' entree gwsetup' ! ------------------------------------------------------------------ ! * 1.1 computational constants ! ----------------------- ilevm1 = klev - 1 ilevm2 = klev - 2 ilevh = klev / 3 zcons1 = 1. / rd ! old zcons2=g**2/cpd zcons2 = rg**2 / rcpd ! old zcons3=1.5*api zcons3 = 1.5 * rpi ! ------------------------------------------------------------------ ! * 2. ! -------------- ! ------------------------------------------------------------------ ! * 2.1 define low level wind, project winds in plane of ! * low level wind, determine sector in which to take ! * the variance and set indicator for critical levels. DO 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. END DO ! Ajouter une initialisation (L. Li, le 23fev99): DO jk = klev, ilevh, -1 DO jl = kidia, kfdia ll1(jl, jk) = .FALSE. END DO END DO ! * define top of low level flow ! ---------------------------- DO jk = klev, ilevh, -1 DO jl = kidia, kfdia lo = (paphm1(jl, jk) / paphm1(jl, klev + 1)) >= gsigcr IF (lo) THEN kkcrit(jl) = jk END IF zhcrit(jl, jk) = ppic(jl) zhgeo = pgeom1(jl, jk) / rg ll1(jl, jk) = (zhgeo>zhcrit(jl, jk)) IF (ll1(jl, jk) .NEQV. ll1(jl, jk + 1)) THEN kknu(jl) = jk END IF IF (.NOT. ll1(jl, ilevh)) kknu(jl) = ilevh END DO END DO DO jk = klev, ilevh, -1 DO jl = kidia, kfdia zhcrit(jl, jk) = ppic(jl) - pval(jl) zhgeo = pgeom1(jl, jk) / rg ll1(jl, jk) = (zhgeo>zhcrit(jl, jk)) IF (ll1(jl, jk) .NEQV. ll1(jl, jk + 1)) THEN kknu2(jl) = jk END IF IF (.NOT. ll1(jl, ilevh)) kknu2(jl) = ilevh END DO END DO DO jk = klev, ilevh, -1 DO jl = kidia, kfdia zhcrit(jl, jk) = amax1(ppic(jl) - pmea(jl), pmea(jl) - pval(jl)) zhgeo = pgeom1(jl, jk) / rg ll1(jl, jk) = (zhgeo>zhcrit(jl, jk)) IF (ll1(jl, jk) .NEQV. ll1(jl, jk + 1)) THEN kknub(jl) = jk END IF IF (.NOT. ll1(jl, ilevh)) kknub(jl) = ilevh END DO END DO DO jl = kidia, kfdia kknu(jl) = min(kknu(jl), nktopg) kknu2(jl) = min(kknu2(jl), nktopg) kknub(jl) = min(kknub(jl), nktopg) kknul(jl) = klev END DO ! c* initialize various arrays DO 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. END DO ! * define low-level flow ! --------------------- DO jk = klev, 2, -1 DO jl = kidia, kfdia IF (ktest(jl)==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) END IF END DO END DO ! ******************************************************************** ! * define blocked flow ! ------------------- DO jk = klev, ilevh, -1 DO jl = kidia, kfdia IF (jk>=kknub(jl) .AND. jk<=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 END DO END DO DO 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) END DO ! ******* setup orography axes and define plane of profiles ******* DO jl = kidia, kfdia lo = (pulow(jl)=-gvsec) IF (lo) THEN zu = pulow(jl) + 2. * gvsec ELSE zu = pulow(jl) END IF 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) END DO ! ************ define flow in plane of lowlevel stress ************* DO jk = 1, klev DO jl = kidia, kfdia IF (ktest(jl)==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)) END IF ptau(jl, jk) = 0.0 pzdep(jl, jk) = 0.0 ppsi(jl, jk) = 0.0 ll1(jl, jk) = .FALSE. END DO END DO DO jk = 2, klev DO jl = kidia, kfdia IF (ktest(jl)==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)=(kknub(jl) + 1) .AND. jk<=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)) END IF END IF END DO END DO DO 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) END DO ! * 2.3 mean flow richardson number. ! * and critical height for froude layer DO jk = 2, klev DO jl = kidia, kfdia IF (ktest(jl)==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) END IF END DO END DO ! * define top of 'envelope' layer ! ---------------------------- DO jl = kidia, kfdia pnu(jl) = 0.0 znum(jl) = 0.0 END DO DO jk = 2, klev - 1 DO jl = kidia, kfdia IF (ktest(jl)==1) THEN IF (jk>=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)<=gfrcrit) .AND. (pnu(jl)>gfrcrit) .AND. (kkenvh(& jl)==klev)) kkenvh(jl) = jk END IF END IF END DO END DO ! calculation of a dynamical mixing height for the breaking ! of gravity waves: DO jl = kidia, kfdia znup(jl) = 0.0 znum(jl) = 0.0 END DO DO jk = klev - 1, 2, -1 DO jl = kidia, kfdia IF (ktest(jl)==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)<=rpi / 2.) .AND. (znup(jl)>rpi / 2.) .AND. (kkcrith(& jl)==klev)) kkcrith(jl) = jk END IF END DO END DO DO jl = kidia, kfdia kkcrith(jl) = min0(kkcrith(jl), kknu2(jl)) kkcrith(jl) = max0(kkcrith(jl), ilevh * 2) END DO ! directional info for flow blocking ************************* DO jk = ilevh, klev DO jl = kidia, kfdia IF (jk>=kkenvh(jl)) THEN lo = (pum1(jl, jk)=-gvsec) IF (lo) THEN zu = pum1(jl, jk) + 2. * gvsec ELSE zu = pum1(jl, jk) END IF zphi = atan(pvm1(jl, jk) / zu) ppsi(jl, jk) = ptheta(jl) * rpi / 180. - zphi END IF END DO END DO ! forms the vertical 'leakiness' ************************** alpha = 3. DO jk = ilevh, klev DO jl = kidia, kfdia IF (jk>=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.) ! mod 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 END DO END DO END SUBROUTINE orosetup SUBROUTINE gwstress(nlon, nlev, ktest, kcrit, kkenvh, kknu, prho, pstab, & pvph, pstd, psig, pmea, ppic, ptau, pgeom1, pdmod) ! **** *gwstress* ! purpose. ! -------- ! ** interface. ! ---------- ! CALL *gwstress* from *gwdrag* ! explicit arguments : ! -------------------- ! ==== inputs === ! ==== outputs === ! implicit arguments : none ! -------------------- ! method. ! ------- ! externals. ! ---------- ! reference. ! ---------- ! see ecmwf research department documentation of the "i.f.s." ! author. ! ------- ! modifications. ! -------------- ! f. lott put the new gwd on ifs 22/11/93 ! ----------------------------------------------------------------------- USE dimphy USE lmdz_YOEGWD, ONLY: GFRCRIT, GKWAKE, GRCRIT, GVCRIT, GKDRAG, GKLIFT, GHMAX, GRAHILO, GSIGCR, NKTOPG, NSTRA, GSSEC, GTSEC, GVSEC, & GWD_RANDO_RUWMAX, gwd_rando_sat, GWD_FRONT_RUWMAX, gwd_front_sat USE lmdz_yomcst IMPLICIT NONE ! ----------------------------------------------------------------------- ! * 0.1 arguments ! --------- INTEGER nlon, nlev INTEGER kcrit(nlon), ktest(nlon), kkenvh(nlon), kknu(nlon) REAL prho(nlon, nlev + 1), pstab(nlon, nlev + 1), ptau(nlon, nlev + 1), & pvph(nlon, nlev + 1), pgeom1(nlon, nlev), pstd(nlon) REAL psig(nlon) REAL pmea(nlon), ppic(nlon) REAL pdmod(nlon) ! ----------------------------------------------------------------------- ! * 0.2 local arrays ! ------------ INTEGER jl REAL zblock, zvar, zeff LOGICAL lo ! ----------------------------------------------------------------------- ! * 0.3 functions ! --------- ! ------------------------------------------------------------------ ! * 1. initialization ! -------------- ! * 3.1 gravity wave stress. DO jl = kidia, kfdia IF (ktest(jl)==1) THEN ! effective mountain height above the blocked flow IF (kkenvh(jl)==klev) THEN zblock = 0.0 ELSE zblock = (pgeom1(jl, kkenvh(jl)) + pgeom1(jl, kkenvh(jl) + 1)) / 2. / rg END IF 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)) ! too small value of stress or low level flow include critical level ! or low level flow: gravity wave stress nul. lo = (ptau(jl, klev + 1)=kknu(jl)) .OR. & (pvph(jl, klev + 1)kkcrith(jl)) THEN ptau(jl, jk) = ztau(jl, klev + 1) ! ENDIF ! IF(JK.EQ.KKCRITH(JL)) THEN ELSE ptau(jl, jk) = grahilo * ztau(jl, klev + 1) END IF END IF END DO ! * 4.15 CONSTANT SHEAR STRESS UNTIL THE TOP OF THE ! LOW LEVEL FLOW LAYER. ! * 4.2 WAVE DISPLACEMENT AT NEXT LEVEL. ! DO 421 ji=1,kgwd ! jl=kdx(ji) ! Modif vectorisation 02/04/2004 DO jl = kidia, kfdia IF (ktest(jl)==1) THEN IF (jkkkcrith(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 END IF END IF END DO ! REORGANISATION IN THE STRATOSPHERE ! DO 533 ji=1,kgwd ! jl=kdx(ji) ! Modif vectorisation 02/04/2004 DO jl = kidia, kfdia IF (ktest(jl)==1) THEN IF (jknstra) 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 END IF END IF END DO END DO END SUBROUTINE gwprofil SUBROUTINE lift_noro(nlon, nlev, dtime, paprs, pplay, plat, pmea, pstd, ppic, & ktest, t, u, v, pulow, pvlow, pustr, pvstr, d_t, d_u, d_v) USE dimphy USE lmdz_yomcst IMPLICIT NONE ! ====================================================================== ! Auteur(s): F.Lott (LMD/CNRS) date: 19950201 ! Objet: Frottement de la montagne Interface ! ====================================================================== ! Arguments: ! dtime---input-R- pas d'integration (s) ! paprs---input-R-pression pour chaque inter-couche (en Pa) ! pplay---input-R-pression pour le mileu de chaque couche (en Pa) ! t-------input-R-temperature (K) ! u-------input-R-vitesse horizontale (m/s) ! v-------input-R-vitesse horizontale (m/s) ! d_t-----output-R-increment de la temperature ! d_u-----output-R-increment de la vitesse u ! d_v-----output-R-increment de la vitesse v ! ====================================================================== ! ARGUMENTS INTEGER nlon, nlev REAL dtime REAL paprs(klon, klev + 1) REAL pplay(klon, klev) REAL plat(nlon), pmea(nlon) REAL pstd(nlon) REAL ppic(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) INTEGER i, k, ktest(nlon) ! Variables locales: 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) ! initialiser les variables de sortie (pour securite) DO i = 1, klon pulow(i) = 0.0 pvlow(i) = 0.0 pustr(i) = 0.0 pvstr(i) = 0.0 END DO 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 END DO END DO ! preparer les variables d'entree (attention: l'ordre des niveaux ! verticaux augmente du haut vers le bas) 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) END DO END DO DO k = 1, klev + 1 DO i = 1, klon papmh(i, k) = paprs(i, klev - k + 2) END DO END DO DO i = 1, klon zgeom(i, klev) = rd * pt(i, klev) * log(papmh(i, klev + 1) / papmf(i, klev)) END DO 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)) END DO END DO ! appeler la routine principale CALL orolift(klon, klev, ktest, dtime, papmh, zgeom, pt, pu, pv, plat, & pmea, pstd, ppic, pulow, pvlow, pdudt, pdvdt, pdtdt) 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) & ! IM BUG . ! +RG*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k)) + pdudt(i, k) * (papmh(i, k + 1) - papmh(i, k)) / rg pvstr(i) = pvstr(i) & ! IM BUG . ! +RG*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k)) + pdvdt(i, k) * (papmh(i, k + 1) - papmh(i, k)) / rg END DO END DO END SUBROUTINE lift_noro SUBROUTINE orolift(nlon, nlev, ktest, ptsphy, paphm1, pgeom1, ptm1, pum1, & pvm1, plat, pmea, pvaror, ppic & ! OUTPUTS , pulow, pvlow, pvom, pvol, pte) ! **** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT. ! PURPOSE. ! -------- ! ** INTERFACE. ! ---------- ! CALLED FROM *lift_noro ! ---------- ! AUTHOR. ! ------- ! F.LOTT LMD 22/11/95 USE dimphy USE lmdz_abort_physic, ONLY: abort_physic USE lmdz_YOEGWD, ONLY: GFRCRIT, GKWAKE, GRCRIT, GVCRIT, GKDRAG, GKLIFT, GHMAX, GRAHILO, GSIGCR, NKTOPG, NSTRA, GSSEC, GTSEC, GVSEC, & GWD_RANDO_RUWMAX, gwd_rando_sat, GWD_FRONT_RUWMAX, gwd_front_sat USE lmdz_yomcst IMPLICIT NONE ! * 0.1 ARGUMENTS ! --------- INTEGER nlon, nlev 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), ppic(nlon), pgeom1(nlon, nlev), & paphm1(nlon, nlev + 1) INTEGER ktest(nlon) REAL ptsphy ! ----------------------------------------------------------------------- ! * 0.2 LOCAL ARRAYS ! ------------ LOGICAL lifthigh ! ym integer klevm1, jl, ilevh, jk INTEGER jl, ilevh, jk REAL zcons1, ztmst, zrtmst, zpi, zhgeo REAL zdelp, zslow, zsqua, zscav, zbet INTEGER iknub(klon), iknul(klon) LOGICAL ll1(klon, klev + 1) REAL ztau(klon, klev + 1), ztav(klon, klev + 1), zrho(klon, klev + 1) REAL zdudt(klon), zdvdt(klon) REAL zhcrit(klon, klev) CHARACTER (LEN = 20) :: modname = 'orografi' CHARACTER (LEN = 80) :: abort_message ! ----------------------------------------------------------------------- ! * 1.1 INITIALIZATIONS ! --------------- lifthigh = .FALSE. IF (nlon/=klon .OR. nlev/=klev) THEN abort_message = 'pb dimension' CALL abort_physic(modname, abort_message, 1) END IF zcons1 = 1. / rd ! ym KLEVM1=KLEV-1 ztmst = ptsphy zrtmst = 1. / ztmst zpi = acos(-1.) DO 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 jk = 1, klev pvom(jl, jk) = 0.0 pvol(jl, jk) = 0.0 pte(jl, jk) = 0.0 END DO END DO ! * 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF ! * LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE ! * THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS. DO jk = klev, 1, -1 DO jl = kidia, kfdia IF (ktest(jl)==1) THEN zhcrit(jl, jk) = amax1(ppic(jl) - pmea(jl), 100.) zhgeo = pgeom1(jl, jk) / rg ll1(jl, jk) = (zhgeo>zhcrit(jl, jk)) IF (ll1(jl, jk) .NEQV. ll1(jl, jk + 1)) THEN iknub(jl) = jk END IF END IF END DO END DO DO jl = kidia, kfdia IF (ktest(jl)==1) THEN iknub(jl) = max(iknub(jl), klev / 2) iknul(jl) = max(iknul(jl), 2 * klev / 3) IF (iknub(jl)>nktopg) iknub(jl) = nktopg IF (iknub(jl)==nktopg) iknul(jl) = klev IF (iknub(jl)==iknul(jl)) iknub(jl) = iknul(jl) - 1 END IF END DO ! do 2011 jl=kidia,kfdia ! IF(KTEST(JL).EQ.1) THEN ! PRINT *,' iknul= ',iknul(jl),' iknub=',iknub(jl) ! ENDIF ! 2011 continue ! PRINT *,' DANS OROLIFT: 2010' DO jk = klev, 2, -1 DO jl = kidia, kfdia zrho(jl, jk) = 2. * paphm1(jl, jk) * zcons1 / (ptm1(jl, jk) + ptm1(jl, jk - 1)) END DO END DO ! PRINT *,' DANS OROLIFT: 223' ! ******************************************************************** ! * DEFINE LOW LEVEL FLOW ! ------------------- DO jk = klev, 1, -1 DO jl = kidia, kfdia IF (ktest(jl)==1) THEN IF (jk>=iknub(jl) .AND. jk<=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 END IF END DO END DO DO jl = kidia, kfdia IF (ktest(jl)==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))) END IF END DO ! *********************************************************** ! * 3. COMPUTE MOUNTAIN LIFT DO jl = kidia, kfdia IF (ktest(jl)==1) THEN ztau(jl, klev + 1) = -gklift * zrho(jl, klev + 1) * 2. * romega * & ! * ! (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 * & ! * ! (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 END IF END DO ! * 4. COMPUTE LIFT PROFILE ! * -------------------- DO jk = 1, klev DO jl = kidia, kfdia IF (ktest(jl)==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 END IF END DO END DO ! * 5. COMPUTE TENDENCIES. ! * ------------------- IF (lifthigh) THEN ! PRINT *,' DANS OROLIFT: 500' ! EXPLICIT SOLUTION AT ALL LEVELS DO jk = 1, klev DO jl = kidia, kfdia IF (ktest(jl)==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 END IF END DO END DO ! PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY DO jk = 1, klev DO jl = kidia, kfdia IF (ktest(jl)==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>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 END IF zsqua = sqrt(pum1(jl, jk)**2 + pum1(jl, jk)**2) IF (zsqua