MODULE lmdz_cloudth IMPLICIT NONE CONTAINS SUBROUTINE cloudth(ngrid, klev, ind2, & ztv, po, zqta, fraca, & qcloud, ctot, zpspsk, paprs, pplay, ztla, zthl, & ratqs, zqs, t, & cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv) USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, iflag_ratqs USE lmdz_yoethf USE lmdz_yomcst IMPLICIT NONE INCLUDE "FCTTRE.h" !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== INTEGER itap, ind1, ind2 INTEGER ngrid, klev, klon, l, ig REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv REAL ztv(ngrid, klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid, klev) REAL fraca(ngrid, klev + 1) REAL zpspsk(ngrid, klev) REAL paprs(ngrid, klev + 1) REAL pplay(ngrid, klev) REAL ztla(ngrid, klev) REAL zthl(ngrid, klev) REAL zqsatth(ngrid, klev) REAL zqsatenv(ngrid, klev) REAL sigma1(ngrid, klev) REAL sigma2(ngrid, klev) REAL qlth(ngrid, klev) REAL qlenv(ngrid, klev) REAL qltot(ngrid, klev) REAL cth(ngrid, klev) REAL cenv(ngrid, klev) REAL ctot(ngrid, klev) REAL rneb(ngrid, klev) REAL t(ngrid, klev) REAL qsatmmussig1, qsatmmussig2, sqrt2pi, pi REAL rdd, cppd, Lv REAL alth, alenv, ath, aenv REAL sth, senv, sigma1s, sigma2s, xth, xenv REAL Tbef, zdelta, qsatbef, zcor REAL qlbef REAL ratqs(ngrid, klev) ! determine la largeur de distribution de vapeur REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid) REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) REAL erf !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Gestion de deux versions de cloudth !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (iflag_cloudth_vert>=1) THEN CALL cloudth_vert(ngrid, klev, ind2, & ztv, po, zqta, fraca, & qcloud, ctot, zpspsk, paprs, pplay, ztla, zthl, & ratqs, zqs, t) RETURN ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !------------------------------------------------------------------------------- ! Initialisation des variables r?elles !------------------------------------------------------------------------------- sigma1(:, :) = 0. sigma2(:, :) = 0. qlth(:, :) = 0. qlenv(:, :) = 0. qltot(:, :) = 0. rneb(:, :) = 0. qcloud(:) = 0. cth(:, :) = 0. cenv(:, :) = 0. ctot(:, :) = 0. qsatmmussig1 = 0. qsatmmussig2 = 0. rdd = 287.04 cppd = 1005.7 pi = 3.14159 Lv = 2.5e6 sqrt2pi = sqrt(2. * pi) !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ?cart-types des distributions !------------------------------------------------------------------------------- DO ind1 = 1, ngrid IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2)) ! zqenv(ind1)=po(ind1) Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatenv(ind1, ind2) = qsatbef alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) aenv = 1. / (1. + (alenv * Lv / cppd)) senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatth(ind1, ind2) = qsatbef alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2) ath = 1. / (1. + (alth * Lv / cppd)) sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2)) !------------------------------------------------------------------------------ ! Calcul des ?cart-types pour s !------------------------------------------------------------------------------ ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2) ! if (paprs(ind1,ind2).gt.90000) THEN ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif sigma1s = (1.1**0.5) * (fraca(ind1, ind2)**0.6) / (1 - fraca(ind1, ind2)) * ((sth - senv)**2)**0.5 + 0.002 * po(ind1) sigma2s = 0.11 * ((sth - senv)**2)**0.5 / (fraca(ind1, ind2) + 0.01)**0.4 + 0.002 * zqta(ind1, ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 !------------------------------------------------------------------------------ ! Calcul de l'eau condens?e et de la couverture nuageuse !------------------------------------------------------------------------------ sqrt2pi = sqrt(2. * pi) xth = sth / (sqrt(2.) * sigma2s) xenv = senv / (sqrt(2.) * sigma1s) cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth)) cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2) qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth**2) / sqrt2pi) + xth * sqrt(2.) * cth(ind1, ind2)) qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv**2) / sqrt2pi) + xenv * sqrt(2.) * cenv(ind1, ind2)) qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (ctot(ind1, ind2)<1.e-10) THEN ctot(ind1, ind2) = 0. qcloud(ind1) = zqsatenv(ind1, ind2) else ctot(ind1, ind2) = ctot(ind1, ind2) qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1) endif ! PRINT*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif' else ! gaussienne environnement seule zqenv(ind1) = po(ind1) Tbef = t(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatenv(ind1, ind2) = qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd) alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) aenv = 1. / (1. + (alenv * Lv / cppd)) senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) sigma1s = ratqs(ind1, ind2) * zqenv(ind1) sqrt2pi = sqrt(2. * pi) xenv = senv / (sqrt(2.) * sigma1s) ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv**2) / sqrt2pi) + xenv * sqrt(2.) * cenv(ind1, ind2)) IF (ctot(ind1, ind2)<1.e-3) THEN ctot(ind1, ind2) = 0. qcloud(ind1) = zqsatenv(ind1, ind2) else ctot(ind1, ind2) = ctot(ind1, ind2) qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2) endif endif enddo RETURN ! end END SUBROUTINE cloudth !=========================================================================== SUBROUTINE cloudth_vert(ngrid, klev, ind2, & ztv, po, zqta, fraca, & qcloud, ctot, zpspsk, paprs, pplay, ztla, zthl, & ratqs, zqs, t) !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, vert_alpha USE lmdz_yoethf USE lmdz_yomcst IMPLICIT NONE INCLUDE "FCTTRE.h" INTEGER itap, ind1, ind2 INTEGER ngrid, klev, klon, l, ig REAL ztv(ngrid, klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid, klev) REAL fraca(ngrid, klev + 1) REAL zpspsk(ngrid, klev) REAL paprs(ngrid, klev + 1) REAL pplay(ngrid, klev) REAL ztla(ngrid, klev) REAL zthl(ngrid, klev) REAL zqsatth(ngrid, klev) REAL zqsatenv(ngrid, klev) REAL sigma1(ngrid, klev) REAL sigma2(ngrid, klev) REAL qlth(ngrid, klev) REAL qlenv(ngrid, klev) REAL qltot(ngrid, klev) REAL cth(ngrid, klev) REAL cenv(ngrid, klev) REAL ctot(ngrid, klev) REAL rneb(ngrid, klev) REAL t(ngrid, klev) REAL qsatmmussig1, qsatmmussig2, sqrt2pi, pi REAL rdd, cppd, Lv, sqrt2, sqrtpi REAL alth, alenv, ath, aenv REAL sth, senv, sigma1s, sigma2s, xth, xenv REAL xth1, xth2, xenv1, xenv2, deltasth, deltasenv REAL IntJ, IntI1, IntI2, IntI3, coeffqlenv, coeffqlth REAL Tbef, zdelta, qsatbef, zcor REAL qlbef REAL ratqs(ngrid, klev) ! determine la largeur de distribution de vapeur ! Change the width of the PDF used for vertical subgrid scale heterogeneity ! (J Jouhaud, JL Dufresne, JB Madeleine) REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid) REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) REAL erf !------------------------------------------------------------------------------ ! Initialisation des variables r?elles !------------------------------------------------------------------------------ sigma1(:, :) = 0. sigma2(:, :) = 0. qlth(:, :) = 0. qlenv(:, :) = 0. qltot(:, :) = 0. rneb(:, :) = 0. qcloud(:) = 0. cth(:, :) = 0. cenv(:, :) = 0. ctot(:, :) = 0. qsatmmussig1 = 0. qsatmmussig2 = 0. rdd = 287.04 cppd = 1005.7 pi = 3.14159 Lv = 2.5e6 sqrt2pi = sqrt(2. * pi) sqrt2 = sqrt(2.) sqrtpi = sqrt(pi) !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ?cart-types des distributions !------------------------------------------------------------------------------- DO ind1 = 1, ngrid IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2)) ! zqenv(ind1)=po(ind1) Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatenv(ind1, ind2) = qsatbef alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) aenv = 1. / (1. + (alenv * Lv / cppd)) senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatth(ind1, ind2) = qsatbef alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2) ath = 1. / (1. + (alth * Lv / cppd)) sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2)) !------------------------------------------------------------------------------ ! Calcul des ?cart-types pour s !------------------------------------------------------------------------------ sigma1s = (0.92**0.5) * (fraca(ind1, ind2)**0.5) / (1 - fraca(ind1, ind2)) * ((sth - senv)**2)**0.5 + ratqs(ind1, ind2) * po(ind1) sigma2s = 0.09 * ((sth - senv)**2)**0.5 / (fraca(ind1, ind2) + 0.02)**0.5 + 0.002 * zqta(ind1, ind2) ! if (paprs(ind1,ind2).gt.90000) THEN ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 !------------------------------------------------------------------------------ ! Calcul de l'eau condens?e et de la couverture nuageuse !------------------------------------------------------------------------------ sqrt2pi = sqrt(2. * pi) xth = sth / (sqrt(2.) * sigma2s) xenv = senv / (sqrt(2.) * sigma1s) cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth)) cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2) qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth**2) / sqrt2pi) + xth * sqrt(2.) * cth(ind1, ind2)) qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv**2) / sqrt2pi) + xenv * sqrt(2.) * cenv(ind1, ind2)) qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2) IF (iflag_cloudth_vert == 1) THEN !------------------------------------------------------------------------------- ! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) deltasenv = aenv * ratqs(ind1, ind2) * zqsatenv(ind1, ind2) deltasth = ath * ratqs(ind1, ind2) * zqsatth(ind1, ind2) ! deltasenv=aenv*0.01*po(ind1) ! deltasth=ath*0.01*zqta(ind1,ind2) xenv1 = (senv - deltasenv) / (sqrt(2.) * sigma1s) xenv2 = (senv + deltasenv) / (sqrt(2.) * sigma1s) xth1 = (sth - deltasth) / (sqrt(2.) * sigma2s) xth2 = (sth + deltasth) / (sqrt(2.) * sigma2s) coeffqlenv = (sigma1s)**2 / (2 * sqrtpi * deltasenv) coeffqlth = (sigma2s)**2 / (2 * sqrtpi * deltasth) cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth2)) cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv2)) ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2) IntJ = sigma1s * (exp(-1. * xenv1**2) / sqrt2pi) + 0.5 * senv * (1 + erf(xenv1)) IntI1 = coeffqlenv * 0.5 * (0.5 * sqrtpi * (erf(xenv2) - erf(xenv1)) + xenv1 * exp(-1. * xenv1**2) - xenv2 * exp(-1. * xenv2**2)) IntI2 = coeffqlenv * xenv2 * (exp(-1. * xenv2**2) - exp(-1. * xenv1**2)) IntI3 = coeffqlenv * 0.5 * sqrtpi * xenv2**2 * (erf(xenv2) - erf(xenv1)) qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3 ! qlenv(ind1,ind2)=IntJ ! PRINT*, qlenv(ind1,ind2),'VERIF EAU' IntJ = sigma2s * (exp(-1. * xth1**2) / sqrt2pi) + 0.5 * sth * (1 + erf(xth1)) ! IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2)) ! IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1)) IntI1 = coeffqlth * 0.5 * (0.5 * sqrtpi * (erf(xth2) - erf(xth1)) + xth1 * exp(-1. * xth1**2) - xth2 * exp(-1. * xth2**2)) IntI2 = coeffqlth * xth2 * (exp(-1. * xth2**2) - exp(-1. * xth1**2)) IntI3 = coeffqlth * 0.5 * sqrtpi * xth2**2 * (erf(xth2) - erf(xth1)) qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3 ! qlth(ind1,ind2)=IntJ ! PRINT*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2' qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2) ELSE IF (iflag_cloudth_vert == 2) THEN !------------------------------------------------------------------------------- ! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) ! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) ! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) deltasenv = aenv * vert_alpha * sigma1s deltasth = ath * vert_alpha * sigma2s xenv1 = -(senv + deltasenv) / (sqrt(2.) * sigma1s) xenv2 = -(senv - deltasenv) / (sqrt(2.) * sigma1s) xth1 = -(sth + deltasth) / (sqrt(2.) * sigma2s) xth2 = -(sth - deltasth) / (sqrt(2.) * sigma2s) ! coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv) ! coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth) cth(ind1, ind2) = 0.5 * (1. - 1. * erf(xth1)) cenv(ind1, ind2) = 0.5 * (1. - 1. * erf(xenv1)) ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2) IntJ = 0.5 * senv * (1 - erf(xenv2)) + (sigma1s / sqrt2pi) * exp(-1. * xenv2**2) IntI1 = (((senv + deltasenv)**2 + (sigma1s)**2) / (8 * deltasenv)) * (erf(xenv2) - erf(xenv1)) IntI2 = (sigma1s**2 / (4 * deltasenv * sqrtpi)) * (xenv1 * exp(-1. * xenv1**2) - xenv2 * exp(-1. * xenv2**2)) IntI3 = ((sqrt2 * sigma1s * (senv + deltasenv)) / (4 * sqrtpi * deltasenv)) * (exp(-1. * xenv1**2) - exp(-1. * xenv2**2)) ! IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2)) ! IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2)) ! IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1)) qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3 ! qlenv(ind1,ind2)=IntJ ! PRINT*, qlenv(ind1,ind2),'VERIF EAU' IntJ = 0.5 * sth * (1 - erf(xth2)) + (sigma2s / sqrt2pi) * exp(-1. * xth2**2) IntI1 = (((sth + deltasth)**2 + (sigma2s)**2) / (8 * deltasth)) * (erf(xth2) - erf(xth1)) IntI2 = (sigma2s**2 / (4 * deltasth * sqrtpi)) * (xth1 * exp(-1. * xth1**2) - xth2 * exp(-1. * xth2**2)) IntI3 = ((sqrt2 * sigma2s * (sth + deltasth)) / (4 * sqrtpi * deltasth)) * (exp(-1. * xth1**2) - exp(-1. * xth2**2)) qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3 ! qlth(ind1,ind2)=IntJ ! PRINT*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2' qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2) ENDIF ! of if (iflag_cloudth_vert==1 or 2) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (cenv(ind1, ind2)<1.e-10.OR.cth(ind1, ind2)<1.e-10) THEN ctot(ind1, ind2) = 0. qcloud(ind1) = zqsatenv(ind1, ind2) else ctot(ind1, ind2) = ctot(ind1, ind2) qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1) ! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) & ! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) endif ! PRINT*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif' else ! gaussienne environnement seule zqenv(ind1) = po(ind1) Tbef = t(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatenv(ind1, ind2) = qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd) alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) aenv = 1. / (1. + (alenv * Lv / cppd)) senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) sigma1s = ratqs(ind1, ind2) * zqenv(ind1) sqrt2pi = sqrt(2. * pi) xenv = senv / (sqrt(2.) * sigma1s) ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv**2) / sqrt2pi) + xenv * sqrt(2.) * cenv(ind1, ind2)) IF (ctot(ind1, ind2)<1.e-3) THEN ctot(ind1, ind2) = 0. qcloud(ind1) = zqsatenv(ind1, ind2) else ctot(ind1, ind2) = ctot(ind1, ind2) qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2) endif endif enddo RETURN ! end END SUBROUTINE cloudth_vert SUBROUTINE cloudth_v3(ngrid, klev, ind2, & ztv, po, zqta, fraca, & qcloud, ctot, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, & ratqs, zqs, t, & cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv) USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert USE lmdz_yoethf USE lmdz_yomcst IMPLICIT NONE INCLUDE "FCTTRE.h" !=========================================================================== ! Author : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== INTEGER, INTENT(IN) :: ind2 INTEGER, INTENT(IN) :: ngrid, klev REAL, DIMENSION(ngrid, klev), INTENT(IN) :: ztv REAL, DIMENSION(ngrid), INTENT(IN) :: po REAL, DIMENSION(ngrid, klev), INTENT(IN) :: zqta REAL, DIMENSION(ngrid, klev + 1), INTENT(IN) :: fraca REAL, DIMENSION(ngrid), INTENT(OUT) :: qcloud REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: ctot REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: ctot_vol REAL, DIMENSION(ngrid, klev), INTENT(IN) :: zpspsk REAL, DIMENSION(ngrid, klev + 1), INTENT(IN) :: paprs REAL, DIMENSION(ngrid, klev), INTENT(IN) :: pplay REAL, DIMENSION(ngrid, klev), INTENT(IN) :: ztla REAL, DIMENSION(ngrid, klev), INTENT(INOUT) :: zthl REAL, DIMENSION(ngrid, klev), INTENT(IN) :: ratqs REAL, DIMENSION(ngrid), INTENT(IN) :: zqs REAL, DIMENSION(ngrid, klev), INTENT(IN) :: t REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv REAL zqenv(ngrid) REAL zqsatth(ngrid, klev) REAL zqsatenv(ngrid, klev) REAL sigma1(ngrid, klev) REAL sigma2(ngrid, klev) REAL qlth(ngrid, klev) REAL qlenv(ngrid, klev) REAL qltot(ngrid, klev) REAL cth(ngrid, klev) REAL cenv(ngrid, klev) REAL cth_vol(ngrid, klev) REAL cenv_vol(ngrid, klev) REAL rneb(ngrid, klev) REAL qsatmmussig1, qsatmmussig2, sqrt2pi, sqrt2, sqrtpi, pi REAL rdd, cppd, Lv REAL alth, alenv, ath, aenv REAL sth, senv, sigma1s, sigma2s, xth, xenv, exp_xenv1, exp_xenv2, exp_xth1, exp_xth2 REAL inverse_rho, beta, a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks REAL Tbef, zdelta, qsatbef, zcor REAL qlbef REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid) REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid) REAL erf INTEGER :: ind1, l, ig IF (iflag_cloudth_vert>=1) THEN CALL cloudth_vert_v3(ngrid, klev, ind2, & ztv, po, zqta, fraca, & qcloud, ctot, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, & ratqs, zqs, t, & cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv) RETURN ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !------------------------------------------------------------------------------- ! Initialisation des variables r?elles !------------------------------------------------------------------------------- sigma1(:, :) = 0. sigma2(:, :) = 0. qlth(:, :) = 0. qlenv(:, :) = 0. qltot(:, :) = 0. rneb(:, :) = 0. qcloud(:) = 0. cth(:, :) = 0. cenv(:, :) = 0. ctot(:, :) = 0. cth_vol(:, :) = 0. cenv_vol(:, :) = 0. ctot_vol(:, :) = 0. qsatmmussig1 = 0. qsatmmussig2 = 0. rdd = 287.04 cppd = 1005.7 pi = 3.14159 Lv = 2.5e6 sqrt2pi = sqrt(2. * pi) sqrt2 = sqrt(2.) sqrtpi = sqrt(pi) !------------------------------------------------------------------------------- ! Cloud fraction in the thermals and standard deviation of the PDFs !------------------------------------------------------------------------------- DO ind1 = 1, ngrid IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2)) Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatenv(ind1, ind2) = qsatbef alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84 aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84 senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84 !po = qt de l'environnement ET des thermique !zqenv = qt environnement !zqsatenv = qsat environnement !zthl = Tl environnement Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatth(ind1, ind2) = qsatbef alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2) !qsl, p84 ath = 1. / (1. + (alth * Lv / cppd)) !al, p84 sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2)) !s, p84 !zqta = qt thermals !zqsatth = qsat thermals !ztla = Tl thermals !------------------------------------------------------------------------------ ! s standard deviations !------------------------------------------------------------------------------ ! tests ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1) ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2) ! final option sigma1s = (1.1**0.5) * (fraca(ind1, ind2)**0.6) / (1 - fraca(ind1, ind2)) * ((sth - senv)**2)**0.5 + ratqs(ind1, ind2) * po(ind1) sigma2s = 0.11 * ((sth - senv)**2)**0.5 / (fraca(ind1, ind2) + 0.01)**0.4 + 0.002 * zqta(ind1, ind2) !------------------------------------------------------------------------------ ! Condensed water and cloud cover !------------------------------------------------------------------------------ xth = sth / (sqrt2 * sigma2s) xenv = senv / (sqrt2 * sigma1s) cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2) ctot_vol(ind1, ind2) = ctot(ind1, ind2) qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth**2) / sqrt2pi) + xth * sqrt2 * cth(ind1, ind2)) qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv**2) / sqrt2pi) + xenv * sqrt2 * cenv(ind1, ind2)) qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2) IF (ctot(ind1, ind2)<1.e-10) THEN ctot(ind1, ind2) = 0. qcloud(ind1) = zqsatenv(ind1, ind2) else qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1) endif else ! Environnement only, follow the if l.110 zqenv(ind1) = po(ind1) Tbef = t(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatenv(ind1, ind2) = qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd) alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) aenv = 1. / (1. + (alenv * Lv / cppd)) senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) sigma1s = ratqs(ind1, ind2) * zqenv(ind1) xenv = senv / (sqrt2 * sigma1s) ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) ctot_vol(ind1, ind2) = ctot(ind1, ind2) qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv**2) / sqrt2pi) + xenv * sqrt2 * cenv(ind1, ind2)) IF (ctot(ind1, ind2)<1.e-3) THEN ctot(ind1, ind2) = 0. qcloud(ind1) = zqsatenv(ind1, ind2) else qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2) endif endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183 enddo ! from the loop on ngrid l.108 RETURN ! end END SUBROUTINE cloudth_v3 !=========================================================================== SUBROUTINE cloudth_vert_v3(ngrid, klev, ind2, & ztv, po, zqta, fraca, & qcloud, ctot, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, & ratqs, zqs, t, & cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv) !=========================================================================== ! Auteur : Arnaud Octavio Jam (LMD/CNRS) ! Date : 25 Mai 2010 ! Objet : calcule les valeurs de qc et rneb dans les thermiques !=========================================================================== USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, iflag_ratqs USE lmdz_cloudth_ini, ONLY: vert_alpha, vert_alpha_th, sigma1s_factor, sigma1s_power, sigma2s_factor, sigma2s_power, cloudth_ratqsmin, iflag_cloudth_vert_noratqs USE lmdz_yoethf USE lmdz_yomcst IMPLICIT NONE INCLUDE "FCTTRE.h" INTEGER itap, ind1, ind2 INTEGER ngrid, klev, klon, l, ig REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv REAL ztv(ngrid, klev) REAL po(ngrid) REAL zqenv(ngrid) REAL zqta(ngrid, klev) REAL fraca(ngrid, klev + 1) REAL zpspsk(ngrid, klev) REAL paprs(ngrid, klev + 1) REAL pplay(ngrid, klev) REAL ztla(ngrid, klev) REAL zthl(ngrid, klev) REAL zqsatth(ngrid, klev) REAL zqsatenv(ngrid, klev) REAL sigma1(ngrid, klev) REAL sigma2(ngrid, klev) REAL qlth(ngrid, klev) REAL qlenv(ngrid, klev) REAL qltot(ngrid, klev) REAL cth(ngrid, klev) REAL cenv(ngrid, klev) REAL ctot(ngrid, klev) REAL cth_vol(ngrid, klev) REAL cenv_vol(ngrid, klev) REAL ctot_vol(ngrid, klev) REAL rneb(ngrid, klev) REAL t(ngrid, klev) REAL qsatmmussig1, qsatmmussig2, sqrtpi, sqrt2, sqrt2pi, pi REAL rdd, cppd, Lv REAL alth, alenv, ath, aenv REAL sth, senv, sigma1s, sigma2s, sigma1s_fraca, sigma1s_ratqs REAL inverse_rho, beta, a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks REAL xth, xenv, exp_xenv1, exp_xenv2, exp_xth1, exp_xth2 REAL xth1, xth2, xenv1, xenv2, deltasth, deltasenv REAL IntJ, IntI1, IntI2, IntI3, IntJ_CF, IntI1_CF, IntI3_CF, coeffqlenv, coeffqlth REAL Tbef, zdelta, qsatbef, zcor REAL qlbef REAL ratqs(ngrid, klev) ! determine la largeur de distribution de vapeur ! Change the width of the PDF used for vertical subgrid scale heterogeneity ! (J Jouhaud, JL Dufresne, JB Madeleine) REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid) REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid) REAL zqs(ngrid), qcloud(ngrid) REAL erf REAL rhodz(ngrid, klev) REAL zrho(ngrid, klev) REAL dz(ngrid, klev) DO ind1 = 1, ngrid !Layer calculation rhodz(ind1, ind2) = (paprs(ind1, ind2) - paprs(ind1, ind2 + 1)) / rg !kg/m2 zrho(ind1, ind2) = pplay(ind1, ind2) / t(ind1, ind2) / rd !kg/m3 dz(ind1, ind2) = rhodz(ind1, ind2) / zrho(ind1, ind2) !m : epaisseur de la couche en metre END DO !------------------------------------------------------------------------------ ! Initialize !------------------------------------------------------------------------------ sigma1(:, :) = 0. sigma2(:, :) = 0. qlth(:, :) = 0. qlenv(:, :) = 0. qltot(:, :) = 0. rneb(:, :) = 0. qcloud(:) = 0. cth(:, :) = 0. cenv(:, :) = 0. ctot(:, :) = 0. cth_vol(:, :) = 0. cenv_vol(:, :) = 0. ctot_vol(:, :) = 0. qsatmmussig1 = 0. qsatmmussig2 = 0. rdd = 287.04 cppd = 1005.7 pi = 3.14159 Lv = 2.5e6 sqrt2pi = sqrt(2. * pi) sqrt2 = sqrt(2.) sqrtpi = sqrt(pi) !------------------------------------------------------------------------------- ! Calcul de la fraction du thermique et des ecart-types des distributions !------------------------------------------------------------------------------- DO ind1 = 1, ngrid IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) then !Thermal and environnement zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2)) !qt = a*qtth + (1-a)*qtenv Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatenv(ind1, ind2) = qsatbef alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84 aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84 senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84 !zqenv = qt environnement !zqsatenv = qsat environnement !zthl = Tl environnement Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatth(ind1, ind2) = qsatbef alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2) !qsl, p84 ath = 1. / (1. + (alth * Lv / cppd)) !al, p84 sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2)) !s, p84 !zqta = qt thermals !zqsatth = qsat thermals !ztla = Tl thermals !------------------------------------------------------------------------------ ! s standard deviation !------------------------------------------------------------------------------ sigma1s_fraca = (sigma1s_factor**0.5) * (fraca(ind1, ind2)**sigma1s_power) / & (1 - fraca(ind1, ind2)) * ((sth - senv)**2)**0.5 ! sigma1s_fraca = (1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5 IF (cloudth_ratqsmin>0.) THEN sigma1s_ratqs = cloudth_ratqsmin * po(ind1) ELSE sigma1s_ratqs = ratqs(ind1, ind2) * po(ind1) ENDIF sigma1s = sigma1s_fraca + sigma1s_ratqs IF (iflag_ratqs==11) THEN sigma1s = ratqs(ind1, ind2) * po(ind1) * aenv ENDIF sigma2s = (sigma2s_factor * (((sth - senv)**2)**0.5) / ((fraca(ind1, ind2) + 0.02)**sigma2s_power)) + 0.002 * zqta(ind1, ind2) ! tests ! sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1) ! sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1) ! sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2) ! sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2) ! if (paprs(ind1,ind2).gt.90000) THEN ! ratqs(ind1,ind2)=0.002 ! else ! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000 ! endif ! sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) ! sigma1s=ratqs(ind1,ind2)*po(ind1) ! sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 IF (iflag_cloudth_vert == 1) THEN !------------------------------------------------------------------------------- ! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs !------------------------------------------------------------------------------- deltasenv = aenv * ratqs(ind1, ind2) * zqsatenv(ind1, ind2) deltasth = ath * ratqs(ind1, ind2) * zqsatth(ind1, ind2) xenv1 = (senv - deltasenv) / (sqrt(2.) * sigma1s) xenv2 = (senv + deltasenv) / (sqrt(2.) * sigma1s) xth1 = (sth - deltasth) / (sqrt(2.) * sigma2s) xth2 = (sth + deltasth) / (sqrt(2.) * sigma2s) coeffqlenv = (sigma1s)**2 / (2 * sqrtpi * deltasenv) coeffqlth = (sigma2s)**2 / (2 * sqrtpi * deltasth) cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth2)) cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv2)) ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2) ! Environment IntJ = sigma1s * (exp(-1. * xenv1**2) / sqrt2pi) + 0.5 * senv * (1 + erf(xenv1)) IntI1 = coeffqlenv * 0.5 * (0.5 * sqrtpi * (erf(xenv2) - erf(xenv1)) + xenv1 * exp(-1. * xenv1**2) - xenv2 * exp(-1. * xenv2**2)) IntI2 = coeffqlenv * xenv2 * (exp(-1. * xenv2**2) - exp(-1. * xenv1**2)) IntI3 = coeffqlenv * 0.5 * sqrtpi * xenv2**2 * (erf(xenv2) - erf(xenv1)) qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3 ! Thermal IntJ = sigma2s * (exp(-1. * xth1**2) / sqrt2pi) + 0.5 * sth * (1 + erf(xth1)) IntI1 = coeffqlth * 0.5 * (0.5 * sqrtpi * (erf(xth2) - erf(xth1)) + xth1 * exp(-1. * xth1**2) - xth2 * exp(-1. * xth2**2)) IntI2 = coeffqlth * xth2 * (exp(-1. * xth2**2) - exp(-1. * xth1**2)) IntI3 = coeffqlth * 0.5 * sqrtpi * xth2**2 * (erf(xth2) - erf(xth1)) qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3 qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2) ELSE IF (iflag_cloudth_vert >= 3) THEN IF (iflag_cloudth_vert < 5) THEN !------------------------------------------------------------------------------- ! Version 3: Changes by J. Jouhaud; condensation for q > -delta s !------------------------------------------------------------------------------- ! deltasenv=aenv*ratqs(ind1,ind2)*po(ind1) ! deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2) ! deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2) ! deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2) IF (iflag_cloudth_vert == 3) THEN deltasenv = aenv * vert_alpha * sigma1s deltasth = ath * vert_alpha_th * sigma2s ELSE IF (iflag_cloudth_vert == 4) THEN IF (iflag_cloudth_vert_noratqs == 1) THEN deltasenv = vert_alpha * max(sigma1s_fraca, 1e-10) deltasth = vert_alpha_th * sigma2s ELSE deltasenv = vert_alpha * sigma1s deltasth = vert_alpha_th * sigma2s ENDIF ENDIF xenv1 = -(senv + deltasenv) / (sqrt(2.) * sigma1s) xenv2 = -(senv - deltasenv) / (sqrt(2.) * sigma1s) exp_xenv1 = exp(-1. * xenv1**2) exp_xenv2 = exp(-1. * xenv2**2) xth1 = -(sth + deltasth) / (sqrt(2.) * sigma2s) xth2 = -(sth - deltasth) / (sqrt(2.) * sigma2s) exp_xth1 = exp(-1. * xth1**2) exp_xth2 = exp(-1. * xth2**2) !CF_surfacique cth(ind1, ind2) = 0.5 * (1. - 1. * erf(xth1)) cenv(ind1, ind2) = 0.5 * (1. - 1. * erf(xenv1)) ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2) !CF_volumique & eau condense !environnement IntJ = 0.5 * senv * (1 - erf(xenv2)) + (sigma1s / sqrt2pi) * exp_xenv2 IntJ_CF = 0.5 * (1. - 1. * erf(xenv2)) IF (deltasenv < 1.e-10) THEN qlenv(ind1, ind2) = IntJ cenv_vol(ind1, ind2) = IntJ_CF else IntI1 = (((senv + deltasenv)**2 + (sigma1s)**2) / (8 * deltasenv)) * (erf(xenv2) - erf(xenv1)) IntI2 = (sigma1s**2 / (4 * deltasenv * sqrtpi)) * (xenv1 * exp_xenv1 - xenv2 * exp_xenv2) IntI3 = ((sqrt2 * sigma1s * (senv + deltasenv)) / (4 * sqrtpi * deltasenv)) * (exp_xenv1 - exp_xenv2) IntI1_CF = ((senv + deltasenv) * (erf(xenv2) - erf(xenv1))) / (4 * deltasenv) IntI3_CF = (sqrt2 * sigma1s * (exp_xenv1 - exp_xenv2)) / (4 * sqrtpi * deltasenv) qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3 cenv_vol(ind1, ind2) = IntJ_CF + IntI1_CF + IntI3_CF endif !thermique IntJ = 0.5 * sth * (1 - erf(xth2)) + (sigma2s / sqrt2pi) * exp_xth2 IntJ_CF = 0.5 * (1. - 1. * erf(xth2)) IF (deltasth < 1.e-10) THEN qlth(ind1, ind2) = IntJ cth_vol(ind1, ind2) = IntJ_CF else IntI1 = (((sth + deltasth)**2 + (sigma2s)**2) / (8 * deltasth)) * (erf(xth2) - erf(xth1)) IntI2 = (sigma2s**2 / (4 * deltasth * sqrtpi)) * (xth1 * exp_xth1 - xth2 * exp_xth2) IntI3 = ((sqrt2 * sigma2s * (sth + deltasth)) / (4 * sqrtpi * deltasth)) * (exp_xth1 - exp_xth2) IntI1_CF = ((sth + deltasth) * (erf(xth2) - erf(xth1))) / (4 * deltasth) IntI3_CF = (sqrt2 * sigma2s * (exp_xth1 - exp_xth2)) / (4 * sqrtpi * deltasth) qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3 cth_vol(ind1, ind2) = IntJ_CF + IntI1_CF + IntI3_CF endif qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2) ctot_vol(ind1, ind2) = fraca(ind1, ind2) * cth_vol(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv_vol(ind1, ind2) ELSE IF (iflag_cloudth_vert == 5) THEN sigma1s = (0.71794 + 0.000498239 * dz(ind1, ind2)) * (fraca(ind1, ind2)**0.5) & / (1 - fraca(ind1, ind2)) * (((sth - senv)**2)**0.5) & + ratqs(ind1, ind2) * po(ind1) !Environment sigma2s = (0.03218 + 0.000092655 * dz(ind1, ind2)) / ((fraca(ind1, ind2) + 0.02)**0.5) * (((sth - senv)**2)**0.5) + 0.002 * zqta(ind1, ind2) !Thermals !sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1) !sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2) xth = sth / (sqrt(2.) * sigma2s) xenv = senv / (sqrt(2.) * sigma1s) !Volumique cth_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xth)) cenv_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) ctot_vol(ind1, ind2) = fraca(ind1, ind2) * cth_vol(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv_vol(ind1, ind2) !PRINT *,'jeanjean_CV=',ctot_vol(ind1,ind2) qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth**2) / sqrt2pi) + xth * sqrt(2.) * cth_vol(ind1, ind2)) qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv**2) / sqrt2pi) + xenv * sqrt(2.) * cenv_vol(ind1, ind2)) qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2) !Surfacique !Neggers !beta=0.0044 !inverse_rho=1.+beta*dz(ind1,ind2) !PRINT *,'jeanjean : beta=',beta !cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho !cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho !ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) !Brooks a_Brooks = 0.6694 b_Brooks = 0.1882 A_Maj_Brooks = 0.1635 !-- sans shear !A_Maj_Brooks=0.17 !-- ARM LES !A_Maj_Brooks=0.18 !-- RICO LES !A_Maj_Brooks=0.19 !-- BOMEX LES Dx_Brooks = 200000. f_Brooks = A_Maj_Brooks * (dz(ind1, ind2)**(a_Brooks)) * (Dx_Brooks**(-b_Brooks)) !PRINT *,'jeanjean_f=',f_Brooks cth(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(cth_vol(ind1, ind2), 1.))) - 1.)) cenv(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(cenv_vol(ind1, ind2), 1.))) - 1.)) ctot(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(ctot_vol(ind1, ind2), 1.))) - 1.)) !PRINT *,'JJ_ctot_1',ctot(ind1,ind2) ENDIF ! of if (iflag_cloudth_vert<5) ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4) ! if (ctot(ind1,ind2).lt.1.e-10) THEN IF (cenv(ind1, ind2)<1.e-10.OR.cth(ind1, ind2)<1.e-10) THEN ctot(ind1, ind2) = 0. ctot_vol(ind1, ind2) = 0. qcloud(ind1) = zqsatenv(ind1, ind2) else qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1) ! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) & ! & +(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1) endif else ! gaussienne environnement seule zqenv(ind1) = po(ind1) Tbef = t(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatenv(ind1, ind2) = qsatbef ! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.) zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd) alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) aenv = 1. / (1. + (alenv * Lv / cppd)) senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) sth = 0. sigma1s = ratqs(ind1, ind2) * zqenv(ind1) sigma2s = 0. sqrt2pi = sqrt(2. * pi) xenv = senv / (sqrt(2.) * sigma1s) ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) ctot_vol(ind1, ind2) = ctot(ind1, ind2) qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv**2) / sqrt2pi) + xenv * sqrt(2.) * cenv(ind1, ind2)) IF (ctot(ind1, ind2)<1.e-3) THEN ctot(ind1, ind2) = 0. qcloud(ind1) = zqsatenv(ind1, ind2) else ! ctot(ind1,ind2)=ctot(ind1,ind2) qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2) endif endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492 ! Outputs used to check the PDFs cloudth_senv(ind1, ind2) = senv cloudth_sth(ind1, ind2) = sth cloudth_sigmaenv(ind1, ind2) = sigma1s cloudth_sigmath(ind1, ind2) = sigma2s enddo ! from the loop on ngrid l.333 RETURN ! end END SUBROUTINE cloudth_vert_v3 SUBROUTINE cloudth_v6(ngrid, klev, ind2, & ztv, po, zqta, fraca, & qcloud, ctot_surf, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, & ratqs, zqs, T, & cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv) USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert USE lmdz_yoethf USE lmdz_yomcst IMPLICIT NONE INCLUDE "FCTTRE.h" !Domain variables INTEGER ngrid !indice Max lat-lon INTEGER klev !indice Max alt REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv INTEGER ind1 !indice in [1:ngrid] INTEGER ind2 !indice in [1:klev] !thermal plume fraction REAL fraca(ngrid, klev + 1) !thermal plumes fraction in the gridbox !temperatures REAL T(ngrid, klev) !temperature REAL zpspsk(ngrid, klev) !factor (p/p0)**kappa (used for potential variables) REAL ztv(ngrid, klev) !potential temperature (voir thermcell_env.F90) REAL ztla(ngrid, klev) !liquid temperature in the thermals (Tl_th) REAL zthl(ngrid, klev) !liquid temperature in the environment (Tl_env) !pressure REAL paprs(ngrid, klev + 1) !pressure at the interface of levels REAL pplay(ngrid, klev) !pressure at the middle of the level !humidity REAL ratqs(ngrid, klev) !width of the total water subgrid-scale distribution REAL po(ngrid) !total water (qt) REAL zqenv(ngrid) !total water in the environment (qt_env) REAL zqta(ngrid, klev) !total water in the thermals (qt_th) REAL zqsatth(ngrid, klev) !water saturation level in the thermals (q_sat_th) REAL zqsatenv(ngrid, klev) !water saturation level in the environment (q_sat_env) REAL qlth(ngrid, klev) !condensed water in the thermals REAL qlenv(ngrid, klev) !condensed water in the environment REAL qltot(ngrid, klev) !condensed water in the gridbox !cloud fractions REAL cth_vol(ngrid, klev) !cloud fraction by volume in the thermals REAL cenv_vol(ngrid, klev) !cloud fraction by volume in the environment REAL ctot_vol(ngrid, klev) !cloud fraction by volume in the gridbox REAL cth_surf(ngrid, klev) !cloud fraction by surface in the thermals REAL cenv_surf(ngrid, klev) !cloud fraction by surface in the environment REAL ctot_surf(ngrid, klev) !cloud fraction by surface in the gridbox !PDF of saturation deficit variables REAL rdd, cppd, Lv REAL Tbef, zdelta, qsatbef, zcor REAL alth, alenv, ath, aenv REAL sth, senv !saturation deficits in the thermals and environment REAL sigma_env, sigma_th !standard deviations of the biGaussian PDF !cloud fraction variables REAL xth, xenv REAL inverse_rho, beta !Neggers et al. (2011) method REAL a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks !Brooks et al. (2005) method !Incloud total water variables REAL zqs(ngrid) !q_sat REAL qcloud(ngrid) !eau totale dans le nuage !Some arithmetic variables REAL erf, pi, sqrt2, sqrt2pi !Depth of the layer REAL dz(ngrid, klev) !epaisseur de la couche en metre REAL rhodz(ngrid, klev) REAL zrho(ngrid, klev) DO ind1 = 1, ngrid rhodz(ind1, ind2) = (paprs(ind1, ind2) - paprs(ind1, ind2 + 1)) / rg ![kg/m2] zrho(ind1, ind2) = pplay(ind1, ind2) / T(ind1, ind2) / rd ![kg/m3] dz(ind1, ind2) = rhodz(ind1, ind2) / zrho(ind1, ind2) ![m] END DO !------------------------------------------------------------------------------ ! Initialization !------------------------------------------------------------------------------ qlth(:, :) = 0. qlenv(:, :) = 0. qltot(:, :) = 0. cth_vol(:, :) = 0. cenv_vol(:, :) = 0. ctot_vol(:, :) = 0. cth_surf(:, :) = 0. cenv_surf(:, :) = 0. ctot_surf(:, :) = 0. qcloud(:) = 0. rdd = 287.04 cppd = 1005.7 pi = 3.14159 Lv = 2.5e6 sqrt2 = sqrt(2.) sqrt2pi = sqrt(2. * pi) DO ind1 = 1, ngrid !------------------------------------------------------------------------------- !Both thermal and environment in the gridbox !------------------------------------------------------------------------------- IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN !-------------------------------------------- !calcul de qsat_env !-------------------------------------------- Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatenv(ind1, ind2) = qsatbef !-------------------------------------------- !calcul de s_env !-------------------------------------------- alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84 these Arnaud Jam aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84 these Arnaud Jam senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul de qsat_th !-------------------------------------------- Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatth(ind1, ind2) = qsatbef !-------------------------------------------- !calcul de s_th !-------------------------------------------- alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2) !qsl, p84 these Arnaud Jam ath = 1. / (1. + (alth * Lv / cppd)) !al, p84 these Arnaud Jam sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul standard deviations bi-Gaussian PDF !-------------------------------------------- sigma_th = (0.03218 + 0.000092655 * dz(ind1, ind2)) / ((fraca(ind1, ind2) + 0.01)**0.5) * (((sth - senv)**2)**0.5) + 0.002 * zqta(ind1, ind2) sigma_env = (0.71794 + 0.000498239 * dz(ind1, ind2)) * (fraca(ind1, ind2)**0.5) & / (1 - fraca(ind1, ind2)) * (((sth - senv)**2)**0.5) & + ratqs(ind1, ind2) * po(ind1) xth = sth / (sqrt2 * sigma_th) xenv = senv / (sqrt2 * sigma_env) !-------------------------------------------- !Cloud fraction by volume CF_vol !-------------------------------------------- cth_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xth)) cenv_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) ctot_vol(ind1, ind2) = fraca(ind1, ind2) * cth_vol(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv_vol(ind1, ind2) !-------------------------------------------- !Condensed water qc !-------------------------------------------- qlth(ind1, ind2) = sigma_th * ((exp(-1. * xth**2) / sqrt2pi) + xth * sqrt2 * cth_vol(ind1, ind2)) qlenv(ind1, ind2) = sigma_env * ((exp(-1. * xenv**2) / sqrt2pi) + xenv * sqrt2 * cenv_vol(ind1, ind2)) qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2) !-------------------------------------------- !Cloud fraction by surface CF_surf !-------------------------------------------- !Method Neggers et al. (2011) : ok for cumulus clouds only !beta=0.0044 (Jouhaud et al.2018) !inverse_rho=1.+beta*dz(ind1,ind2) !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho !Method Brooks et al. (2005) : ok for all types of clouds a_Brooks = 0.6694 b_Brooks = 0.1882 A_Maj_Brooks = 0.1635 !-- sans dependence au cisaillement de vent Dx_Brooks = 200000. !-- si l'on considere des mailles de 200km de cote f_Brooks = A_Maj_Brooks * (dz(ind1, ind2)**(a_Brooks)) * (Dx_Brooks**(-b_Brooks)) ctot_surf(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(ctot_vol(ind1, ind2), 1.))) - 1.)) !-------------------------------------------- !Incloud Condensed water qcloud !-------------------------------------------- IF (ctot_surf(ind1, ind2) < 1.e-10) THEN ctot_vol(ind1, ind2) = 0. ctot_surf(ind1, ind2) = 0. qcloud(ind1) = zqsatenv(ind1, ind2) else qcloud(ind1) = qltot(ind1, ind2) / ctot_vol(ind1, ind2) + zqs(ind1) endif !------------------------------------------------------------------------------- !Environment only in the gridbox !------------------------------------------------------------------------------- ELSE !-------------------------------------------- !calcul de qsat_env !-------------------------------------------- Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2) zdelta = MAX(0., SIGN(1., RTT - Tbef)) qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2) qsatbef = MIN(0.5, qsatbef) zcor = 1. / (1. - retv * qsatbef) qsatbef = qsatbef * zcor zqsatenv(ind1, ind2) = qsatbef !-------------------------------------------- !calcul de s_env !-------------------------------------------- alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84 these Arnaud Jam aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84 these Arnaud Jam senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84 these Arnaud Jam !-------------------------------------------- !calcul standard deviations Gaussian PDF !-------------------------------------------- zqenv(ind1) = po(ind1) sigma_env = ratqs(ind1, ind2) * zqenv(ind1) xenv = senv / (sqrt2 * sigma_env) !-------------------------------------------- !Cloud fraction by volume CF_vol !-------------------------------------------- ctot_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) !-------------------------------------------- !Condensed water qc !-------------------------------------------- qltot(ind1, ind2) = sigma_env * ((exp(-1. * xenv**2) / sqrt2pi) + xenv * sqrt2 * ctot_vol(ind1, ind2)) !-------------------------------------------- !Cloud fraction by surface CF_surf !-------------------------------------------- !Method Neggers et al. (2011) : ok for cumulus clouds only !beta=0.0044 (Jouhaud et al.2018) !inverse_rho=1.+beta*dz(ind1,ind2) !ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho !Method Brooks et al. (2005) : ok for all types of clouds a_Brooks = 0.6694 b_Brooks = 0.1882 A_Maj_Brooks = 0.1635 !-- sans dependence au shear Dx_Brooks = 200000. f_Brooks = A_Maj_Brooks * (dz(ind1, ind2)**(a_Brooks)) * (Dx_Brooks**(-b_Brooks)) ctot_surf(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(ctot_vol(ind1, ind2), 1.))) - 1.)) !-------------------------------------------- !Incloud Condensed water qcloud !-------------------------------------------- IF (ctot_surf(ind1, ind2) < 1.e-8) THEN ctot_vol(ind1, ind2) = 0. ctot_surf(ind1, ind2) = 0. qcloud(ind1) = zqsatenv(ind1, ind2) else qcloud(ind1) = qltot(ind1, ind2) / ctot_vol(ind1, ind2) + zqsatenv(ind1, ind2) endif END IF ! From the separation (thermal/envrionnement) et (environnement only) ! Outputs used to check the PDFs cloudth_senv(ind1, ind2) = senv cloudth_sth(ind1, ind2) = sth cloudth_sigmaenv(ind1, ind2) = sigma_env cloudth_sigmath(ind1, ind2) = sigma_th END DO ! From the loop on ngrid END SUBROUTINE cloudth_v6 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE cloudth_mpc(klon, klev, ind2, mpc_bl_points, & & temp, qt, qt_th, frac_th, zpspsk, paprsup, paprsdn, play, thetal_th, & & ratqs, qcloud, qincloud, icefrac, ctot, ctot_vol, & & cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Author : Etienne Vignon (LMDZ/CNRS) ! Date: April 2024 ! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam ! IMPORTANT NOTE: we assume iflag_cloudth_vert=7 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, iflag_ratqs USE lmdz_cloudth_ini, ONLY: C_mpc, Ni, C_cap, Ei, d_top, vert_alpha, vert_alpha_th, sigma1s_factor, sigma1s_power, sigma2s_factor, sigma2s_power, cloudth_ratqsmin, iflag_cloudth_vert_noratqs USE lmdz_lscp_tools, only: CALC_QSAT_ECMWF USE lmdz_lscp_ini, only: RTT, RG, RPI, RD, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th IMPLICIT NONE !------------------------------------------------------------------------------ ! Declaration !------------------------------------------------------------------------------ ! INPUT/OUTPUT INTEGER, INTENT(IN) :: klon, klev, ind2 REAL, DIMENSION(klon), INTENT(IN) :: temp ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip REAL, DIMENSION(klon), INTENT(IN) :: qt ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip REAL, DIMENSION(klon), INTENT(IN) :: qt_th ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip REAL, DIMENSION(klon), INTENT(IN) :: thetal_th ! Liquid potential temperature in thermals [K]: has not seen the evap of precip REAL, DIMENSION(klon), INTENT(IN) :: frac_th ! Fraction of the mesh covered by thermals [0-1] REAL, DIMENSION(klon), INTENT(IN) :: zpspsk ! Exner potential REAL, DIMENSION(klon), INTENT(IN) :: paprsup ! Pressure at top layer interface [Pa] REAL, DIMENSION(klon), INTENT(IN) :: paprsdn ! Pressure at bottom layer interface [Pa] REAL, DIMENSION(klon), INTENT(IN) :: play ! Pressure of layers [Pa] REAL, DIMENSION(klon), INTENT(IN) :: ratqs ! Parameter that determines the width of the total water distrib. INTEGER, DIMENSION(klon, klev), INTENT(INOUT) :: mpc_bl_points ! grid points with convective (thermals) mixed phase clouds REAL, DIMENSION(klon), INTENT(OUT) :: ctot ! Cloud fraction [0-1] REAL, DIMENSION(klon), INTENT(OUT) :: ctot_vol ! Volume cloud fraction [0-1] REAL, DIMENSION(klon), INTENT(OUT) :: qcloud ! In cloud total water content [kg/kg] REAL, DIMENSION(klon), INTENT(OUT) :: qincloud ! In cloud condensed water content [kg/kg] REAL, DIMENSION(klon), INTENT(OUT) :: icefrac ! Fraction of ice in clouds [0-1] REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment ! LOCAL VARIABLES INTEGER itap, ind1, l, ig, iter, k INTEGER iflag_topthermals, niter REAL qcth(klon) REAL qcenv(klon) REAL qctot(klon) REAL cth(klon) REAL cenv(klon) REAL cth_vol(klon) REAL cenv_vol(klon) REAL qt_env(klon), thetal_env(klon) REAL icefracenv(klon), icefracth(klon) REAL sqrtpi, sqrt2, sqrt2pi REAL alth, alenv, ath, aenv REAL sth, senv, sigma1s, sigma2s, sigma1s_fraca, sigma1s_ratqs REAL inverse_rho, beta, a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks REAL xth, xenv, exp_xenv1, exp_xenv2, exp_xth1, exp_xth2 REAL xth1, xth2, xenv1, xenv2, deltasth, deltasenv REAL IntJ, IntI1, IntI2, IntI3, IntJ_CF, IntI1_CF, IntI3_CF, coeffqlenv, coeffqlth REAL zdelta, qsatbef, zcor REAL Tbefth(klon), Tbefenv(klon), zrho(klon), rhoth(klon) REAL qlbef REAL dqsatenv(klon), dqsatth(klon) REAL erf REAL zpdf_sig(klon), zpdf_k(klon), zpdf_delta(klon) REAL zpdf_a(klon), zpdf_b(klon), zpdf_e1(klon), zpdf_e2(klon) REAL rhodz(klon) REAL rho(klon) REAL dz(klon) REAL qslenv(klon), qslth(klon), qsienv(klon), qsith(klon) REAL alenvl, aenvl REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc REAL senvi, senvl, qbase, sbase, qliqth, qiceth REAL qmax, ttarget, stmp, cout, coutref, fraci REAL maxi, mini, pas ! Modifty the saturation deficit PDF in thermals ! in the presence of ice crystals CHARACTER (len = 80) :: abort_message CHARACTER (len = 20) :: routname = 'cloudth_mpc' !------------------------------------------------------------------------------ ! Initialisation !------------------------------------------------------------------------------ ! Few initial checksS DO k = 1, klev DO ind1 = 1, klon rhodz(ind1) = (paprsdn(ind1) - paprsup(ind1)) / rg !kg/m2 zrho(ind1) = play(ind1) / temp(ind1) / rd !kg/m3 dz(ind1) = rhodz(ind1) / zrho(ind1) !m : epaisseur de la couche en metre END DO END DO icefracth(:) = 0. icefracenv(:) = 0. sqrt2pi = sqrt(2. * rpi) sqrt2 = sqrt(2.) sqrtpi = sqrt(rpi) icefrac(:) = 0. !------------------------------------------------------------------------------- ! Identify grid points with potential mixed-phase conditions !------------------------------------------------------------------------------- DO ind1 = 1, klon IF ((temp(ind1) < RTT) .AND. (temp(ind1) > temp_nowater) & .AND. (ind2<=klev - 2) & .AND. (frac_th(ind1)>min_frac_th_cld)) THEN mpc_bl_points(ind1, ind2) = 1 ELSE mpc_bl_points(ind1, ind2) = 0 ENDIF ENDDO !------------------------------------------------------------------------------- ! Thermal fraction calculation and standard deviation of the distribution !------------------------------------------------------------------------------- ! initialisations and calculation of temperature, humidity and saturation specific humidity DO ind1 = 1, klon rhodz(ind1) = (paprsdn(ind1) - paprsup(ind1)) / rg ! kg/m2 rho(ind1) = play(ind1) / temp(ind1) / rd ! kg/m3 dz(ind1) = rhodz(ind1) / rho(ind1) ! m : epaisseur de la couche en metre Tbefenv(ind1) = temp(ind1) thetal_env(ind1) = Tbefenv(ind1) / zpspsk(ind1) Tbefth(ind1) = thetal_th(ind1) * zpspsk(ind1) rhoth(ind1) = play(ind1) / Tbefth(ind1) / RD qt_env(ind1) = (qt(ind1) - frac_th(ind1) * qt_th(ind1)) / (1. - frac_th(ind1)) !qt = a*qtth + (1-a)*qtenv ENDDO ! Calculation of saturation specific humidity CALL CALC_QSAT_ECMWF(klon, Tbefenv, qt_env, play, RTT, 1, .FALSE., qslenv, dqsatenv) CALL CALC_QSAT_ECMWF(klon, Tbefenv, qt_env, play, RTT, 2, .FALSE., qsienv, dqsatenv) CALL CALC_QSAT_ECMWF(klon, Tbefth, qt_th, play, RTT, 1, .FALSE., qslth, dqsatth) DO ind1 = 1, klon IF (frac_th(ind1)>min_frac_th_cld) THEN !Thermal and environnement ! unlike in the other cloudth routine, ! We consider distributions of the saturation deficit WRT liquid ! at positive AND negative celcius temperature ! subsequently, cloud fraction corresponds to the part of the pdf corresponding ! to superstauration with respect to liquid WHATEVER the temperature ! Environment: alenv = (0.622 * RLVTT * qslenv(ind1)) / (rd * thetal_env(ind1)**2) aenv = 1. / (1. + (alenv * RLVTT / rcpd)) senv = aenv * (qt_env(ind1) - qslenv(ind1)) ! Thermals: alth = (0.622 * RLVTT * qslth(ind1)) / (rd * thetal_th(ind1)**2) ath = 1. / (1. + (alth * RLVTT / rcpd)) sth = ath * (qt_th(ind1) - qslth(ind1)) ! IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC ! Standard deviation of the distributions sigma1s_fraca = (sigma1s_factor**0.5) * (frac_th(ind1)**sigma1s_power) / & (1 - frac_th(ind1)) * ((sth - senv)**2)**0.5 IF (cloudth_ratqsmin>0.) THEN sigma1s_ratqs = cloudth_ratqsmin * qt(ind1) ELSE sigma1s_ratqs = ratqs(ind1) * qt(ind1) ENDIF sigma1s = sigma1s_fraca + sigma1s_ratqs IF (iflag_ratqs==11) THEN sigma1s = ratqs(ind1) * qt(ind1) * aenv ENDIF sigma2s = (sigma2s_factor * (((sth - senv)**2)**0.5) / ((frac_th(ind1) + 0.02)**sigma2s_power)) + 0.002 * qt_th(ind1) deltasenv = aenv * vert_alpha * sigma1s deltasth = ath * vert_alpha_th * sigma2s xenv1 = -(senv + deltasenv) / (sqrt(2.) * sigma1s) xenv2 = -(senv - deltasenv) / (sqrt(2.) * sigma1s) exp_xenv1 = exp(-1. * xenv1**2) exp_xenv2 = exp(-1. * xenv2**2) xth1 = -(sth + deltasth) / (sqrt(2.) * sigma2s) xth2 = -(sth - deltasth) / (sqrt(2.) * sigma2s) exp_xth1 = exp(-1. * xth1**2) exp_xth2 = exp(-1. * xth2**2) !surface CF cth(ind1) = 0.5 * (1. - 1. * erf(xth1)) cenv(ind1) = 0.5 * (1. - 1. * erf(xenv1)) ctot(ind1) = frac_th(ind1) * cth(ind1) + (1. - 1. * frac_th(ind1)) * cenv(ind1) !volume CF, condensed water and ice fraction !environnement IntJ = 0.5 * senv * (1 - erf(xenv2)) + (sigma1s / sqrt2pi) * exp_xenv2 IntJ_CF = 0.5 * (1. - 1. * erf(xenv2)) IF (deltasenv < 1.e-10) THEN qcenv(ind1) = IntJ cenv_vol(ind1) = IntJ_CF ELSE IntI1 = (((senv + deltasenv)**2 + (sigma1s)**2) / (8 * deltasenv)) * (erf(xenv2) - erf(xenv1)) IntI2 = (sigma1s**2 / (4 * deltasenv * sqrtpi)) * (xenv1 * exp_xenv1 - xenv2 * exp_xenv2) IntI3 = ((sqrt2 * sigma1s * (senv + deltasenv)) / (4 * sqrtpi * deltasenv)) * (exp_xenv1 - exp_xenv2) IntI1_CF = ((senv + deltasenv) * (erf(xenv2) - erf(xenv1))) / (4 * deltasenv) IntI3_CF = (sqrt2 * sigma1s * (exp_xenv1 - exp_xenv2)) / (4 * sqrtpi * deltasenv) qcenv(ind1) = IntJ + IntI1 + IntI2 + IntI3 cenv_vol(ind1) = IntJ_CF + IntI1_CF + IntI3_CF IF (Tbefenv(ind1) < temp_nowater) THEN ! freeze all droplets in cirrus temperature regime icefracenv(ind1) = 1. ENDIF ENDIF !thermals IntJ = 0.5 * sth * (1 - erf(xth2)) + (sigma2s / sqrt2pi) * exp_xth2 IntJ_CF = 0.5 * (1. - 1. * erf(xth2)) IF (deltasth < 1.e-10) THEN qcth(ind1) = IntJ cth_vol(ind1) = IntJ_CF ELSE IntI1 = (((sth + deltasth)**2 + (sigma2s)**2) / (8 * deltasth)) * (erf(xth2) - erf(xth1)) IntI2 = (sigma2s**2 / (4 * deltasth * sqrtpi)) * (xth1 * exp_xth1 - xth2 * exp_xth2) IntI3 = ((sqrt2 * sigma2s * (sth + deltasth)) / (4 * sqrtpi * deltasth)) * (exp_xth1 - exp_xth2) IntI1_CF = ((sth + deltasth) * (erf(xth2) - erf(xth1))) / (4 * deltasth) IntI3_CF = (sqrt2 * sigma2s * (exp_xth1 - exp_xth2)) / (4 * sqrtpi * deltasth) qcth(ind1) = IntJ + IntI1 + IntI2 + IntI3 cth_vol(ind1) = IntJ_CF + IntI1_CF + IntI3_CF IF (Tbefth(ind1) < temp_nowater) THEN ! freeze all droplets in cirrus temperature regime icefracth(ind1) = 1. ENDIF ENDIF qctot(ind1) = frac_th(ind1) * qcth(ind1) + (1. - 1. * frac_th(ind1)) * qcenv(ind1) ctot_vol(ind1) = frac_th(ind1) * cth_vol(ind1) + (1. - 1. * frac_th(ind1)) * cenv_vol(ind1) IF (cenv(ind1) 0) THEN icefrac(ind1) = (frac_th(ind1) * qcth(ind1) * icefracth(ind1) + (1. - 1. * frac_th(ind1)) * qcenv(ind1) * icefracenv(ind1)) / qctot(ind1) icefrac(ind1) = max(min(1., icefrac(ind1)), 0.) ENDIF ENDIF ! ELSE ! mpc_bl_points>0 ELSE ! gaussian for environment only alenv = (0.622 * RLVTT * qslenv(ind1)) / (rd * thetal_env(ind1)**2) aenv = 1. / (1. + (alenv * RLVTT / rcpd)) senv = aenv * (qt_env(ind1) - qslenv(ind1)) sth = 0. sigma1s = ratqs(ind1) * qt_env(ind1) sigma2s = 0. xenv = senv / (sqrt(2.) * sigma1s) cenv(ind1) = 0.5 * (1. - 1. * erf(xenv)) ctot(ind1) = 0.5 * (1. + 1. * erf(xenv)) ctot_vol(ind1) = ctot(ind1) qctot(ind1) = sigma1s * ((exp(-1. * xenv**2) / sqrt2pi) + xenv * sqrt(2.) * cenv(ind1)) IF (ctot(ind1)0 ! ! Treat boundary layer mixed phase clouds ! ! thermals ! !========= ! ! ice phase ! !........... ! qiceth=0. ! deltazlev_mpc=dz(ind1,:) ! temp_mpc=ztla(ind1,:)*zpspsk(ind1,:) ! pres_mpc=pplay(ind1,:) ! fraca_mpc=fraca(ind1,:) ! snowf_mpc=snowflux(ind1,:) ! iflag_topthermals=0 ! IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN ! iflag_topthermals = 1 ! ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) & ! .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN ! iflag_topthermals = 2 ! ELSE ! iflag_topthermals = 0 ! ENDIF ! CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), & ! qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:)) ! ! qmax calculation ! sigma2s=(sigma2s_factor*((MAX((sthl-senvl),0.)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2) ! deltasth=athl*vert_alpha_th*sigma2s ! xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s) ! xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s) ! exp_xth1 = exp(-1.*xth1**2) ! exp_xth2 = exp(-1.*xth2**2) ! IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 ! IntJ_CF=0.5*(1.-1.*erf(xth2)) ! IntI1=(((sthl+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) ! IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) ! IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) ! IntI1_CF=((sthl+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) ! IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) ! qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.) ! ! Liquid phase ! !................ ! ! We account for the effect of ice crystals in thermals on sthl ! ! and on the width of the distribution ! sthlc=sthl*1./(1.+C_mpc*qith(ind1,ind2)) & ! + (1.-1./(1.+C_mpc*qith(ind1,ind2))) * athl*(qsith(ind1,ind2)-qslth(ind1)) ! sigma2sc=(sigma2s_factor*((MAX((sthlc-senvl),0.)**2)**0.5) & ! /((fraca(ind1,ind2)+0.02)**sigma2s_power)) & ! +0.002*zqta(ind1,ind2) ! deltasthc=athl*vert_alpha_th*sigma2sc ! xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc) ! xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc) ! exp_xth1 = exp(-1.*xth1**2) ! exp_xth2 = exp(-1.*xth2**2) ! IntJ=0.5*sthlc*(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2 ! IntJ_CF=0.5*(1.-1.*erf(xth2)) ! IntI1=(((sthlc+deltasthc)**2+(sigma2sc)**2)/(8*deltasthc))*(erf(xth2)-erf(xth1)) ! IntI2=(sigma2sc**2/(4*deltasthc*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) ! IntI3=((sqrt2*sigma2sc*(sthlc+deltasthc))/(4*sqrtpi*deltasthc))*(exp_xth1-exp_xth2) ! IntI1_CF=((sthlc+deltasthc)*(erf(xth2)-erf(xth1)))/(4*deltasthc) ! IntI3_CF=(sqrt2*sigma2sc*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasthc) ! qliqth=IntJ+IntI1+IntI2+IntI3 ! qlth(ind1,ind2)=MAX(0.,qliqth) ! ! Condensed water ! qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2) ! ! consistency with subgrid distribution ! IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN ! fraci=qith(ind1,ind2)/qcth(ind1,ind2) ! qcth(ind1,ind2)=qmax ! qith(ind1,ind2)=fraci*qmax ! qlth(ind1,ind2)=(1.-fraci)*qmax ! ENDIF ! ! Cloud Fraction ! !............... ! ! calculation of qbase which is the value of the water vapor within mixed phase clouds ! ! such that the total water in cloud = qbase+qliqth+qiceth ! ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction ! ! sbase and qbase calculation (note that sbase is wrt liq so negative) ! ! look for an approximate solution with iteration ! ttarget=qcth(ind1,ind2) ! mini= athl*(qsith(ind1,ind2)-qslth(ind1)) ! maxi= 0. !athl*(3.*sqrt(sigma2s)) ! niter=20 ! pas=(maxi-mini)/niter ! stmp=mini ! sbase=stmp ! coutref=1.E6 ! DO iter=1,niter ! cout=ABS(sigma2s/SQRT(2.*RPI)*EXP(-((sthl-stmp)/sigma2s)**2)+(sthl-stmp)/SQRT(2.)*(1.-erf(-(sthl-stmp)/sigma2s)) & ! + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget) ! IF (cout .LT. coutref) THEN ! sbase=stmp ! coutref=cout ! ELSE ! stmp=stmp+pas ! ENDIF ! ENDDO ! qbase=MAX(0., sbase/athl+qslth(ind1)) ! ! surface cloud fraction in thermals ! cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s)) ! cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.) ! !volume cloud fraction in thermals ! !to be checked ! xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s) ! xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s) ! exp_xth1 = exp(-1.*xth1**2) ! exp_xth2 = exp(-1.*xth2**2) ! IntJ=0.5*sthl*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 ! IntJ_CF=0.5*(1.-1.*erf(xth2)) ! IF (deltasth .LT. 1.e-10) THEN ! cth_vol(ind1,ind2)=IntJ_CF ! ELSE ! IntI1=(((sthl+deltasth-sbase)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) ! IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) ! IntI3=((sqrt2*sigma2s*(sthl+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) ! IntI1_CF=((sthl-sbase+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) ! IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) ! cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF ! ENDIF ! cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.) ! ! Environment ! !============= ! ! In the environment/downdrafts, ONLY liquid clouds ! ! See Shupe et al. 2008, JAS ! ! standard deviation of the distribution in the environment ! sth=sthl ! senv=senvl ! sigma1s_fraca = (sigma1s_factor**0.5)*(fraca(ind1,ind2)**sigma1s_power) / & ! & (1-fraca(ind1,ind2))*(MAX((sth-senv),0.)**2)**0.5 ! ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution ! ! in the environement ! sigma1s_ratqs=1E-10 ! IF (cloudth_ratqsmin>0.) THEN ! sigma1s_ratqs = cloudth_ratqsmin*po(ind1) ! ENDIF ! sigma1s = sigma1s_fraca + sigma1s_ratqs ! IF (iflag_ratqs.EQ.11) THEN ! sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv ! ENDIF ! IF (iflag_ratqs.EQ.11) THEN ! sigma1s = ratqs(ind1,ind2)*po(ind1)*aenvl ! ENDIF ! deltasenv=aenvl*vert_alpha*sigma1s ! xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s) ! xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s) ! exp_xenv1 = exp(-1.*xenv1**2) ! exp_xenv2 = exp(-1.*xenv2**2) ! !surface CF ! cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) ! !volume CF and condensed water ! IntJ=0.5*senvl*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 ! IntJ_CF=0.5*(1.-1.*erf(xenv2)) ! IF (deltasenv .LT. 1.e-10) THEN ! qcenv(ind1,ind2)=IntJ ! cenv_vol(ind1,ind2)=IntJ_CF ! ELSE ! IntI1=(((senvl+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) ! IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) ! IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) ! IntI1_CF=((senvl+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) ! IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) ! qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment ! cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF ! ENDIF ! qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.) ! cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.) ! ! Thermals + environment ! ! ======================= ! ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) ! qctot(ind1,ind2)=fraca(ind1,ind2)*qcth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2) ! ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) ! IF (qcth(ind1,ind2) .GT. 0) THEN ! icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) & ! /(fraca(ind1,ind2)*qcth(ind1,ind2) & ! +(1.-1.*fraca(ind1,ind2))*qcenv(ind1,ind2)) ! icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.) ! ELSE ! icefrac(ind1,ind2)=0. ! ENDIF ! IF (cenv(ind1,ind2).LT.1.e-10.OR.cth(ind1,ind2).LT.1.e-10) THEN ! ctot(ind1,ind2)=0. ! ctot_vol(ind1,ind2)=0. ! qincloud(ind1)=0. ! qcloud(ind1)=zqsatenv(ind1) ! ELSE ! qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) & ! +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1)) ! qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) & ! +(1.-1.*fraca(ind1,ind2))*(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.) ! ENDIF ! ENDIF ! mpc_bl_points !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE ICE_MPC_BL_CLOUDS(ind1, ind2, klev, Ni, Ei, C_cap, d_top, iflag_topthermals, temp, pres, qth, qsith, qlth, deltazlev, vith, fraca, qith) ! parameterization of ice for boundary ! layer mixed-phase clouds assuming a stationary system ! Note that vapor deposition on ice crystals and riming of liquid droplets ! depend on the ice number concentration Ni ! One could assume that Ni depends on qi, e.g., Ni=beta*(qi*rho)**xi ! and use values from Hong et al. 2004, MWR for instance ! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962 ! One could also think of a more complex expression of Ni; ! function of qi, T, the concentration in aerosols or INP .. ! Here we prefer fixing Ni to a tuning parameter ! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard ! in Mioche et al. 2017 ! References: !------------ ! This parameterization is thoroughly described in Vignon et al. ! More specifically ! for the Water vapor deposition process: ! Rotstayn, L. D., 1997: A physically based scheme for the treat- ! ment of stratiform cloudfs and precipitation in large-scale ! models. I: Description and evaluation of the microphysical ! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282. ! Morrison, H., and A. Gettelman, 2008: A new two-moment bulk ! stratiform cloud microphysics scheme in the NCAR Com- ! munity Atmosphere Model (CAM3). Part I: Description and ! numerical tests. J. Climate, 21, 3642–3659 ! for the Riming process: ! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro- ! scale structure and organization of clouds and precipitation in ! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’ ! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206 ! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit ! forecasts of winter precipitation using an improved bulk ! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit ! forecasts of winter precipitation using an improved bulk ! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542 ! For the formation of clouds by thermals: ! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of ! the Atmospheric Sciences, 65, 407–425. ! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a ! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3 ! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr !============================================================================= USE phys_state_var_mod, ONLY: fm_therm, detr_therm, entr_therm USE lmdz_yomcst IMPLICIT NONE INTEGER, INTENT(IN) :: ind1, ind2, klev ! horizontal and vertical indices and dimensions INTEGER, INTENT(IN) :: iflag_topthermals ! uppermost layer of thermals ? REAL, INTENT(IN) :: Ni ! ice number concentration [m-3] REAL, INTENT(IN) :: Ei ! ice-droplet collision efficiency REAL, INTENT(IN) :: C_cap ! ice crystal capacitance REAL, INTENT(IN) :: d_top ! cloud-top ice crystal mixing parameter REAL, DIMENSION(klev), INTENT(IN) :: temp ! temperature [K] within thermals REAL, DIMENSION(klev), INTENT(IN) :: pres ! pressure [Pa] REAL, DIMENSION(klev), INTENT(IN) :: qth ! mean specific water content in thermals [kg/kg] REAL, DIMENSION(klev), INTENT(IN) :: qsith ! saturation specific humidity wrt ice in thermals [kg/kg] REAL, DIMENSION(klev), INTENT(IN) :: qlth ! condensed liquid water in thermals, approximated value [kg/kg] REAL, DIMENSION(klev), INTENT(IN) :: deltazlev ! layer thickness [m] REAL, DIMENSION(klev), INTENT(IN) :: vith ! ice crystal fall velocity [m/s] REAL, DIMENSION(klev + 1), INTENT(IN) :: fraca ! fraction of the mesh covered by thermals REAL, DIMENSION(klev), INTENT(INOUT) :: qith ! condensed ice water , thermals [kg/kg] INTEGER ind2p1, ind2p2 REAL rho(klev) REAL unsurtaudet, unsurtaustardep, unsurtaurim REAL qi, AA, BB, Ka, Dv, rhoi REAL p0, t0, fp1, fp2 REAL alpha, flux_term REAL det_term, precip_term, rim_term, dep_term ind2p1 = ind2 + 1 ind2p2 = ind2 + 2 rho = pres / temp / RD ! air density kg/m3 Ka = 2.4e-2 ! thermal conductivity of the air, SI p0 = 101325.0 ! ref pressure T0 = 273.15 ! ref temp rhoi = 500.0 ! cloud ice density following Reisner et al. 1998 alpha = 700. ! fallvelocity param IF (iflag_topthermals > 0) THEN ! uppermost thermals levels Dv = 0.0001 * 0.211 * (p0 / pres(ind2)) * ((temp(ind2) / T0)**1.94) ! water vapor diffusivity in air, SI ! Detrainment term: unsurtaudet = detr_therm(ind1, ind2) / rho(ind2) / deltazlev(ind2) ! Deposition term AA = RLSTT / Ka / temp(ind2) * (RLSTT / RV / temp(ind2) - 1.) BB = 1. / (rho(ind2) * Dv * qsith(ind2)) unsurtaustardep = C_cap * (Ni**0.66) * (qth(ind2) - qsith(ind2)) / qsith(ind2) & * 4. * RPI / (AA + BB) * (6. * rho(ind2) / rhoi / RPI / Gamma(4.))**(0.33) ! Riming term neglected at this level !unsurtaurim=rho(ind2)*alpha*3./rhoi/2.*Ei*qlth(ind2)*((p0/pres(ind2))**0.4) qi = fraca(ind2) * unsurtaustardep / MAX((d_top * unsurtaudet), 1E-12) qi = MAX(qi, 0.)**(3. / 2.) ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2 Dv = 0.0001 * 0.211 * (p0 / pres(ind2p1)) * ((temp(ind2p1) / T0)**1.94) ! water vapor diffusivity in air, SI ! Detrainment term: unsurtaudet = detr_therm(ind1, ind2p1) / rho(ind2p1) / deltazlev(ind2p1) det_term = -unsurtaudet * qith(ind2p1) * rho(ind2p1) ! Deposition term AA = RLSTT / Ka / temp(ind2p1) * (RLSTT / RV / temp(ind2p1) - 1.) BB = 1. / (rho(ind2p1) * Dv * qsith(ind2p1)) unsurtaustardep = C_cap * (Ni**0.66) * (qth(ind2p1) - qsith(ind2p1)) & / qsith(ind2p1) * 4. * RPI / (AA + BB) & * (6. * rho(ind2p1) / rhoi / RPI / Gamma(4.))**(0.33) dep_term = rho(ind2p1) * fraca(ind2p1) * (qith(ind2p1)**0.33) * unsurtaustardep ! Riming term unsurtaurim = rho(ind2p1) * alpha * 3. / rhoi / 2. * Ei * qlth(ind2p1) * ((p0 / pres(ind2p1))**0.4) rim_term = rho(ind2p1) * fraca(ind2p1) * qith(ind2p1) * unsurtaurim ! Precip term ! We assume that there is no solid precipitation outside thermals ! so the precipitation flux within thermals is equal to the precipitation flux ! at mesh-scale divided by thermals fraction fp2 = 0. fp1 = 0. IF (fraca(ind2p1) > 0.) THEN fp2 = -qith(ind2p2) * rho(ind2p2) * vith(ind2p2) * fraca(ind2p2)! flux defined positive upward fp1 = -qith(ind2p1) * rho(ind2p1) * vith(ind2p1) * fraca(ind2p1) ENDIF precip_term = -1. / deltazlev(ind2p1) * (fp2 - fp1) ! Calculation in a top-to-bottom loop IF (fm_therm(ind1, ind2p1) > 0.) THEN qi = 1. / fm_therm(ind1, ind2p1) * & (deltazlev(ind2p1) * (-rim_term - dep_term - det_term - precip_term) + & fm_therm(ind1, ind2p2) * (qith(ind2p1))) END IF ENDIF ! top thermals qith(ind2) = MAX(0., qi) END SUBROUTINE ICE_MPC_BL_CLOUDS END MODULE lmdz_cloudth