--- trunk/Sources/phylmd/hbtm.f 2015/04/29 15:47:56 134 +++ trunk/Sources/phylmd/hbtm.f 2016/03/21 15:36:26 186 @@ -4,61 +4,73 @@ contains - SUBROUTINE HBTM(knon, paprs, pplay, t2m, t10m, q2m, q10m, ustar, flux_t, & - flux_q, u, v, t, q, pblh, cape, EauLiq, ctei, pblT, therm, trmb1, & - trmb2, trmb3, plcl) + SUBROUTINE HBTM(knon, paprs, pplay, t2m, q2m, ustar, flux_t, flux_q, u, v, & + t, q, pblh, cape, EauLiq, ctei, pblT, therm, trmb1, trmb2, trmb3, plcl) - ! D'après Holstag et Boville et Troen et Mahrt + ! D'apr\'es Holstag et Boville et Troen et Mahrt ! JAS 47 BLM - ! Algorithme thèse Anne Mathieu - ! Critère d'entraînement Peter Duynkerke (JAS 50) - ! written by: Anne MATHIEU and Alain LAHELLEC, 22nd November 1999 - ! features : implem. exces Mathieu - - ! modifications : decembre 99 passage th a niveau plus bas. voir fixer - ! la prise du th a z/Lambda = -.2 (max Ray) - ! Autre algo : entrainement ~ Theta+v =cste mais comment=>The? - ! on peut fixer q a .7 qsat (cf. non adiabatique) => T2 et The2 - ! voir aussi //KE pblh = niveau The_e ou l = env. - - ! fin therm a la HBTM passage a forme Mathieu 12/09/2001 - - ! Adaptation a LMDZ version couplee - ! Pour le moment on fait passer en argument les grandeurs de surface : - ! flux, t, q2m, t, q10m, on va utiliser systematiquement les grandeurs a 2m - ! mais on garde la possibilité de changer si besoin est (jusqu'à présent - ! la forme de HB avec le 1er niveau modele etait conservee) + + ! Algorithme th\'ese Anne Mathieu. Crit\'ere d'entra\^inement + ! Peter Duynkerke (JAS 50). Written by: Anne MATHIEU and Alain + ! LAHELLEC, 22nd November 1999. + + ! Modifications : d\'ecembre 99 passage th \`a niveau plus bas. Voir fixer + ! la prise du th \`a z/Lambda = -.2 (max Ray) + ! Autre algorithme : entra\^inement ~ Theta + v =constante + ! mais comment ? The ? + ! On peut fixer q \`a 0.7 qsat (cf. non adiabatique) d'où T2 et The2. + ! Voir aussi KE pblh = niveau The_e ou l = env. + + ! Adaptation \`a LMDZ version coupl\'ee. Pour le moment on fait + ! passer en argument les grandeurs de surface : flux, t, q2m. On + ! va utiliser syst\'ematiquement les grandeurs \`a 2 m mais on + ! garde la possibilit\'e de changer si besoin (jusqu'\`a pr\'esent + ! la forme de HB avec le premier niveau mod\`ele \'etait + ! conserv\'ee). USE dimphy, ONLY: klev, klon - USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rlvtt, rtt, rv + USE suphec_m, ONLY: rcpd, rd, retv, rg, rkappa, rtt USE yoethf_m, ONLY: r2es, rvtmp2 USE fcttre, ONLY: foeew - REAL RLvCp, REPS ! Arguments: ! nombre de points a calculer INTEGER, intent(in):: knon - REAL, intent(in):: t2m(klon) ! temperature a 2 m - real t10m(klon) ! temperature a 10 m - ! q a 2 et 10m - REAL q2m(klon), q10m(klon) - REAL ustar(klon) ! pression a inter-couche (Pa) - REAL paprs(klon, klev+1) + REAL, intent(in):: paprs(klon, klev+1) ! pression au milieu de couche (Pa) - REAL pplay(klon, klev) + REAL, intent(in):: pplay(klon, klev) + REAL, intent(in):: t2m(klon) ! temperature a 2 m + ! q a 2 et 10m + REAL, intent(in):: q2m(klon) + REAL, intent(in):: ustar(klon) ! Flux - REAL flux_t(klon, klev), flux_q(klon, klev) + REAL, intent(in):: flux_t(klon, klev), flux_q(klon, klev) ! vitesse U (m/s) - REAL u(klon, klev) + REAL, intent(in):: u(klon, klev) ! vitesse V (m/s) - REAL v(klon, klev) + REAL, intent(in):: v(klon, klev) ! temperature (K) - REAL t(klon, klev) + REAL, intent(in):: t(klon, klev) ! vapeur d'eau (kg/kg) - REAL q(klon, klev) + REAL, intent(in):: q(klon, klev) + + REAL, intent(out):: pblh(:) ! (knon) + ! Cape du thermique + REAL Cape(klon) + ! Eau liqu integr du thermique + REAL EauLiq(klon) + ! Critere d'instab d'entrainmt des nuages de + REAL ctei(klon) + REAL pblT(klon) + ! thermal virtual temperature excess + REAL therm(klon) + REAL trmb1(klon), trmb2(klon), trmb3(klon) + REAL plcl(klon) + + ! Local: INTEGER isommet ! limite max sommet pbl @@ -123,84 +135,34 @@ REAL rhino(klon, klev) ! pts w/unstbl pbl (positive virtual ht flx) LOGICAL unstbl(klon) - ! stable pbl with levels within pbl - LOGICAL stblev(klon) - ! unstbl pbl with levels within pbl - LOGICAL unslev(klon) - ! unstb pbl w/lvls within srf pbl lyr - LOGICAL unssrf(klon) - ! unstb pbl w/lvls in outer pbl lyr - LOGICAL unsout(klon) LOGICAL check(klon) ! Richardson number > critical ! flag de prolongerment cape pour pt Omega LOGICAL omegafl(klon) - REAL pblh(klon) - REAL pblT(klon) - REAL plcl(klon) ! Monin-Obukhov lengh REAL obklen(klon) REAL zdu2 - ! thermal virtual temperature excess - REAL therm(klon) - REAL trmb1(klon), trmb2(klon), trmb3(klon) ! Algorithme thermique REAL s(klon, klev) ! [P/Po]^Kappa milieux couches - ! equivalent potential temperature of therma - REAL The_th(klon) ! total water of thermal REAL qT_th(klon) ! T thermique niveau precedent - REAL Tbef(klon) REAL qsatbef(klon) ! le thermique est sature LOGICAL Zsat(klon) - ! Cape du thermique - REAL Cape(klon) - ! Cape locale - REAL Kape(klon) - ! Eau liqu integr du thermique - REAL EauLiq(klon) - ! Critere d'instab d'entrainmt des nuages de - REAL ctei(klon) - REAL the1, the2, aa, zthvd, zthvu, xintpos, qqsat - REAL a1, a2, a3 - REAL xhis, rnum, th1, thv1, thv2, ql2 - REAL qsat2, qT1, q2, t1, t2, xnull - REAL quadsat, spblh, reduc + REAL zthvd, zthvu, qqsat + REAL t2 ! inverse phi function for momentum REAL phiminv(klon) - ! inverse phi function for heat - REAL phihinv(klon) ! turbulent velocity scale for momentum REAL wm(klon) - ! k*ustar*pblh - REAL fak1(klon) - ! k*wm*pblh - REAL fak2(klon) - ! fakn*wstr/wm - REAL fak3(klon) - ! level eddy diffusivity for momentum - REAL pblk(klon) - ! Prandtl number for eddy diffusivities - REAL pr(klon) - ! zmzp / Obukhov length - REAL zl(klon) - ! zmzp / pblh - REAL zh(klon) - ! (1-(zmzp/pblh))**2 - REAL zzh(klon) - ! w*, convective velocity scale - REAL wstr(klon) - ! current level height - REAL zm(klon) ! current level height + one level up REAL zp(klon) - REAL zcor, zcvm5 + REAL zcor - REAL fac, pblmin, zmzp, term + REAL pblmin !----------------------------------------------------------------- @@ -211,10 +173,6 @@ b212=sqrt(b1*b2) b2sr=sqrt(b2) - ! Initialisation - RLvCp = RLVTT/RCPD - REPS = RD/RV - ! Calculer les hauteurs de chaque couche ! (geopotentielle Int_dp/ro = Int_[Rd.T.dp/p] z = geop/g) ! pourquoi ne pas utiliser Phi/RG ? @@ -246,7 +204,7 @@ ! convention >0 vers le bas ds lmdz khfs(i) = - flux_t(i, 1)*zxt*Rd / (RCPD*paprs(i, 1)) - kqfs(i) = - flux_q(i, 1)*zxt*Rd / (paprs(i, 1)) + kqfs(i) = - flux_q(i, 1)*zxt*Rd / paprs(i, 1) ! verifier que khfs et kqfs sont bien de la forme w'l' heatv(i) = khfs(i) + 0.608*zxt*kqfs(i) ! a comparer aussi aux sorties de clqh : flux_T/RoCp et flux_q/RoLv @@ -274,7 +232,6 @@ ! until the Richardson number between the first level and the ! current level exceeds the "critical" value. (bonne idee Nu de ! separer le Ric et l'exces de temp du thermique) - fac = 100. DO k = 2, isommet DO i = 1, knon IF (check(i)) THEN @@ -339,11 +296,6 @@ q_star = kqfs(i)/wm(i) t_star = khfs(i)/wm(i) - a1=b1*(1.+2.*RETV*qT_th(i))*t_star**2 - a2=(RETV*T2m(i))**2*b2*q_star**2 - a3=2.*RETV*T2m(i)*b212*q_star*t_star - aa=a1+a2+a3 - therm(i) = sqrt( b1*(1.+2.*RETV*qT_th(i))*t_star**2 & + (RETV*T2m(i))**2*b2*q_star**2 & + max(0., 2.*RETV*T2m(i)*b212*q_star*t_star)) @@ -419,11 +371,8 @@ ! omegafl utilise pour prolongement CAPE omegafl(i) = .FALSE. Cape(i) = 0. - Kape(i) = 0. EauLiq(i) = 0. CTEI(i) = 0. - pblk(i) = 0.0 - fak1(i) = ustar(i)*pblh(i)*vk ! Do additional preparation for unstable cases only, set temperature ! and moisture perturbations depending on stability. @@ -433,15 +382,8 @@ zxt=(T2m(i)-zref*0.5*RG/RCPD/(1.+RVTMP2*qT_th(i))) & *(1.+RETV*qT_th(i)) phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet - phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i)) wm(i) = ustar(i)*phiminv(i) - fak2(i) = wm(i)*pblh(i)*vk - wstr(i) = (heatv(i)*RG*pblh(i)/zxt)**onet - fak3(i) = fakn*wstr(i)/wm(i) ENDIF - ! Computes Theta_e for thermal (all cases : to be modified) - ! attention ajout therm(i) = virtuelle - The_th(i) = T2m(i) + therm(i) + RLvCp*qT_th(i) ENDDO ! Main level loop to compute the diffusivities and @@ -449,69 +391,8 @@ loop_level: DO k = 2, isommet ! Find levels within boundary layer: DO i = 1, knon - unslev(i) = .FALSE. - stblev(i) = .FALSE. - zm(i) = z(i, k-1) zp(i) = z(i, k) IF (zkmin == 0. .AND. zp(i) > pblh(i)) zp(i) = pblh(i) - IF (zm(i) < pblh(i)) THEN - zmzp = 0.5*(zm(i) + zp(i)) - zh(i) = zmzp/pblh(i) - zl(i) = zmzp/obklen(i) - zzh(i) = 0. - IF (zh(i) <= 1.) zzh(i) = (1. - zh(i))**2 - - ! stblev for points zm < plbh and stable and neutral - ! unslev for points zm < plbh and unstable - IF (unstbl(i)) THEN - unslev(i) = .TRUE. - ELSE - stblev(i) = .TRUE. - ENDIF - ENDIF - ENDDO - - ! Stable and neutral points; set diffusivities; counter-gradient - ! terms zero for stable case: - DO i = 1, knon - IF (stblev(i)) THEN - IF (zl(i) <= 1.) THEN - pblk(i) = fak1(i)*zh(i)*zzh(i)/(1. + betas*zl(i)) - ELSE - pblk(i) = fak1(i)*zh(i)*zzh(i)/(betas + zl(i)) - ENDIF - ENDIF - ENDDO - - ! unssrf, unstable within surface layer of pbl - ! unsout, unstable within outer layer of pbl - DO i = 1, knon - unssrf(i) = .FALSE. - unsout(i) = .FALSE. - IF (unslev(i)) THEN - IF (zh(i) < sffrac) THEN - unssrf(i) = .TRUE. - ELSE - unsout(i) = .TRUE. - ENDIF - ENDIF - ENDDO - - ! Unstable for surface layer; counter-gradient terms zero - DO i = 1, knon - IF (unssrf(i)) THEN - term = (1. - betam*zl(i))**onet - pblk(i) = fak1(i)*zh(i)*zzh(i)*term - pr(i) = term/sqrt(1. - betah*zl(i)) - ENDIF - ENDDO - - ! Unstable for outer layer; counter-gradient terms non-zero: - DO i = 1, knon - IF (unsout(i)) THEN - pblk(i) = fak2(i)*zh(i)*zzh(i) - pr(i) = phiminv(i)/phihinv(i) + ccon*fak3(i)/fak - ENDIF ENDDO ! For all layers, compute integral info and CTEI @@ -534,7 +415,6 @@ * (qT_th(i)-qsatbef(i)) / (qsatbef(i)-qqsat) endif Zsat(i) = .true. - Tbef(i) = T2 endif endif qsatbef(i) = qqsat