! $Id: newmicro.F 1750 2013-04-25 15:27:27Z fairhead $ ! SUBROUTINE newmicro (ok_cdnc, bl95_b0, bl95_b1, & & paprs, pplay, & & t, pqlwp, pclc, pcltau, pclemi, & & pch, pcl, pcm, pct, pctlwp, & & xflwp, xfiwp, xflwc, xfiwc, & & mass_solu_aero, mass_solu_aero_pi, & & pcldtaupi, re, fl, reliq, reice) !c USE dimphy USE phys_local_var_mod, only: scdnc,cldncl,reffclwtop,lcc, & & reffclws,reffclwc,cldnvi,lcc3d, & & lcc3dcon,lcc3dstra USE phys_state_var_mod, only: rnebcon,clwcon IMPLICIT none !c====================================================================== !c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910 !c O. Boucher (LMD/CNRS) mise a jour en 201212 !c Objet: Calculer epaisseur optique et emmissivite des nuages !c====================================================================== !c Arguments: !c ok_cdnc-input-L-flag pour calculer les rayons a partir des aerosols !c !c t-------input-R-temperature !c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere dans la partie nuageuse (kg/kg) !c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1) !c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3] !c mass_solu_aero_pi--input-R-ditto, pre-industrial value !c !c bl95_b0-input-R-a PARAMETER, may be varied for tests (s-sea, l-land) !c bl95_b1-input-R-a PARAMETER, may be varied for tests ( -"- ) !c !c re------output-R-Cloud droplet effective radius multiplied by fl [um] !c fl------output-R-Denominator to re, introduced to avoid problems in !c the averaging of the output. fl is the fraction of liquid !c water clouds within a grid cell !c !c pcltau--output-R-epaisseur optique des nuages !c pclemi--output-R-emissivite des nuages (0 a 1) !c pcldtaupi-output-R-pre-industrial value of cloud optical thickness, !c !c pcl-output-R-2D low-level cloud cover !c pcm-output-R-2D mid-level cloud cover !c pch-output-R-2D high-level cloud cover !c pct-output-R-2D total cloud cover !c====================================================================== !C #include "YOMCST.h" #include "nuage.h" #include "radepsi.h" #include "radopt.h" !c choix de l'hypothese de recouvrememnt nuageuse LOGICAL RANDOM, MAXIMUM_RANDOM, MAXIMUM PARAMETER (RANDOM=.FALSE., MAXIMUM_RANDOM=.TRUE., MAXIMUM=.FALSE.) !c LOGICAL, SAVE :: FIRST=.TRUE. !$OMP THREADPRIVATE(FIRST) INTEGER flag_max !c !c threshold PARAMETERs REAL thres_tau,thres_neb PARAMETER (thres_tau=0.3, thres_neb=0.001) !c REAL phase3d(klon, klev) REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon) !c REAL paprs(klon,klev+1) REAL pplay(klon,klev) REAL t(klon,klev) REAL pclc(klon,klev) REAL pqlwp(klon,klev) REAL pcltau(klon,klev) REAL pclemi(klon,klev) REAL pcldtaupi(klon, klev) !c REAL pct(klon) REAL pcl(klon) REAL pcm(klon) REAL pch(klon) REAL pctlwp(klon) !c LOGICAL lo !c !!Abderr modif JL mail du 19.01.2011 18:31 ! REAL cetahb, cetamb ! PARAMETER (cetahb = 0.45, cetamb = 0.80) ! Remplacer ! cetahb*paprs(i,1) par prmhc ! cetamb*paprs(i,1) par prlmc REAL prmhc ! Pressure between medium and high level cloud in Pa REAL prlmc ! Pressure between low and medium level cloud in Pa PARAMETER (prmhc = 440.*100., prlmc = 680.*100.) !C INTEGER i, k REAL xflwp(klon), xfiwp(klon) REAL xflwc(klon,klev), xfiwc(klon,klev) !c REAL radius !c REAL coef_froi, coef_chau PARAMETER (coef_chau=0.13, coef_froi=0.09) !c REAL seuil_neb PARAMETER (seuil_neb=0.001) !c INTEGER nexpo ! exponentiel pour glace/eau PARAMETER (nexpo=6) !c PARAMETER (nexpo=1) REAL rel, tc, rei REAL k_ice0, k_ice, DF PARAMETER (k_ice0=0.005) ! units=m2/g PARAMETER (DF=1.66) ! diffusivity factor !c !cjq for the aerosol indirect effect !cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003 !cjq REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols [ug m-3] REAL mass_solu_aero_pi(klon, klev) ! - " - (pre-industrial value) REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3] REAL re(klon, klev) ! cloud droplet effective radius [um] REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value) REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value) REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell) LOGICAL ok_cdnc REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula !cjq-end !IM cf. CR:parametres supplementaires REAL zclear(klon) REAL zcloud(klon) REAL zcloudh(klon) REAL zcloudm(klon) REAL zcloudl(klon) REAL rhodz(klon, klev) !--rho*dz pour la couche REAL zrho(klon, klev) !--rho pour la couche REAL dh(klon, klev) !--dz pour la couche REAL zfice(klon, klev) REAL rad_chaud(klon, klev) !--rayon pour les nuages chauds REAL zflwp_var, zfiwp_var REAL d_rei_dt ! Abderrahmane oct 2009 Real reliq(klon, klev), reice(klon, klev) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! FH : 2011/05/24 ! ! rei = ( rei_max - rei_min ) * T(°C) / 81.4 + rei_max ! to be used for a temperature in celcius T(°C) < 0 ! rei=rei_min for T(°C) < -81.4 ! ! Calcul de la pente de la relation entre rayon effective des cristaux ! et la température. ! Pour retrouver les résultats numériques de la version d'origine, ! on impose 0.71 quand on est proche de 0.71 !c d_rei_dt=(rei_max-rei_min)/81.4 if (abs(d_rei_dt-0.71)<1.e-4) d_rei_dt=0.71 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !c !c Calculer l'epaisseur optique et l'emmissivite des nuages !c IM inversion des DO !c xflwp = 0.d0 xfiwp = 0.d0 xflwc = 0.d0 xfiwc = 0.d0 !c reliq=0. reice=0. !c DO k = 1, klev DO i = 1, klon !c-layer calculation rhodz(i,k) = (paprs(i,k)-paprs(i,k+1))/RG ! kg/m2 zrho(i,k)=pplay(i,k)/t(i,k)/RD ! kg/m3 dh(i,k)=rhodz(i,k)/zrho(i,k) ! m !c-Fraction of ice in cloud using a linear transition zfice(i,k) = 1.0 - (t(i,k)-t_glace_min) / & & (t_glace_max-t_glace_min) zfice(i,k) = MIN(MAX(zfice(i,k),0.0),1.0) !c-IM Total Liquid/Ice water content xflwc(i,k) = (1.-zfice(i,k))*pqlwp(i,k) xfiwc(i,k) = zfice(i,k)*pqlwp(i,k) ENDDO ENDDO IF (ok_cdnc) THEN !c !c--we compute cloud properties as a function of the aerosol load !c DO k = 1, klev DO i = 1, klon !c !c Formula "D" of Boucher and Lohmann, Tellus, 1995 !c Cloud droplet number concentration (CDNC) is restricted !c to be within [20, 1000 cm^3] !c !c--present-day case cdnc(i,k) = 10.**(bl95_b0+bl95_b1* & & log(MAX(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3 cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k))) !c !c--pre-industrial case cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1* & & log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.)) & & *1.e6 !-m-3 cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k))) !c !c--present-day case rad_chaud(i,k) = & & 1.1*((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) & & /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.) rad_chaud(i,k) = MAX(rad_chaud(i,k) * 1.e6, 5.) !c !c--pre-industrial case radius = & & 1.1*((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) & & /(4./3.*RPI*1000.*cdnc_pi(i,k)))**(1./3.) radius = MAX(radius*1.e6, 5.) !c !c--pre-industrial case !c--liquid/ice cloud water paths: IF (pclc(i,k) .LE. seuil_neb) THEN !c pcldtaupi(i,k) = 0.0 !c ELSE !c zflwp_var= 1000.*(1.-zfice(i,k))*pqlwp(i,k)/pclc(i,k) & & *rhodz(i,k) zfiwp_var= 1000.*zfice(i,k)*pqlwp(i,k)/pclc(i,k) & & *rhodz(i,k) tc = t(i,k)-273.15 rei = d_rei_dt*tc + rei_max if (tc.le.-81.4) rei = rei_min !c !c-- cloud optical thickness : !c [for liquid clouds, traditional formula, !c for ice clouds, Ebert & Curry (1992)] !c if (zflwp_var.eq.0.) radius = 1. if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1. pcldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius & & + zfiwp_var * (3.448e-03 + 2.431/rei) !c ENDIF !c ENDDO ENDDO !c ELSE !--not ok_cdnc !c !c-prescribed cloud droplet radius !c DO k = 1, MIN(3,klev) DO i = 1, klon rad_chaud(i,k) = rad_chau2 ENDDO ENDDO DO k = MIN(3,klev)+1, klev DO i = 1, klon rad_chaud(i,k) = rad_chau1 ENDDO ENDDO ENDIF !--ok_cdnc !c !c--computation of cloud optical depth and emissivity !c--in the general case !c DO k = 1, klev DO i = 1, klon !c IF (pclc(i,k) .LE. seuil_neb) THEN !c !c effective cloud droplet radius (microns) for liquid water clouds: !c For output diagnostics cloud droplet effective radius [um] !c we multiply here with f * xl (fraction of liquid water !c clouds in the grid cell) to avoid problems in the averaging of the output. !c In the output of IOIPSL, derive the REAL cloud droplet !c effective radius as re/fl !c fl(i,k) = seuil_neb*(1.-zfice(i,k)) re(i,k) = rad_chaud(i,k)*fl(i,k) rel = 0. rei = 0. pclc(i,k) = 0.0 pcltau(i,k) = 0.0 pclemi(i,k) = 0.0 !c ELSE !c !c -- liquid/ice cloud water paths: zflwp_var= 1000.*(1.-zfice(i,k))*pqlwp(i,k)/pclc(i,k) & & *rhodz(i,k) zfiwp_var= 1000.*zfice(i,k)*pqlwp(i,k)/pclc(i,k) & & *rhodz(i,k) !c !c effective cloud droplet radius (microns) for liquid water clouds: !c For output diagnostics cloud droplet effective radius [um] !c we multiply here with f * xl (fraction of liquid water !c clouds in the grid cell) to avoid problems in the averaging of the output. !c In the output of IOIPSL, derive the REAL cloud droplet !c effective radius as re/fl !c fl(i,k) = pclc(i,k)*(1.-zfice(i,k)) re(i,k) = rad_chaud(i,k)*fl(i,k) !c rel = rad_chaud(i,k) !c !c for ice clouds: as a function of the ambiant temperature !c [formula used by Iacobellis and Somerville (2000), with an !c asymptotical value of 3.5 microns at T<-81.4 C added to be !c consistent with observations of Heymsfield et al. 1986]: !c 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as rei_max=61.29 !c tc = t(i,k)-273.15 rei = d_rei_dt*tc + rei_max if (tc.le.-81.4) rei = rei_min !c !c-- cloud optical thickness : !c [for liquid clouds, traditional formula, !c for ice clouds, Ebert & Curry (1992)] !c if (zflwp_var.eq.0.) rel = 1. if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1. pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel ) & & + zfiwp_var * (3.448e-03 + 2.431/rei) !c !c -- cloud infrared emissivity: !c [the broadband infrared absorption coefficient is PARAMETERized !c as a function of the effective cld droplet radius] !c Ebert and Curry (1992) formula as used by Kiehl & Zender (1995): !c k_ice = k_ice0 + 1.0/rei !c pclemi(i,k) = 1.0 & & - EXP( -coef_chau*zflwp_var - DF*k_ice*zfiwp_var) !c ENDIF !c reliq(i,k)=rel reice(i,k)=rei !c xflwp(i) = xflwp(i)+ xflwc(i,k) * rhodz(i,k) xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * rhodz(i,k) !c ENDDO ENDDO !c !c--if cloud droplet radius is fixed, then pcldtaupi=pcltau !c IF (.NOT.ok_cdnc) THEN DO k = 1, klev DO i = 1, klon pcldtaupi(i,k)=pcltau(i,k) ENDDO ENDDO ENDIF !C !C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS !c IM cf. CR:test: calcul prenant ou non en compte le recouvrement !c initialisations !c DO i=1,klon zclear(i)=1. zcloud(i)=0. zcloudh(i)=0. zcloudm(i)=0. zcloudl(i)=0. pch(i)=1.0 pcm(i) = 1.0 pcl(i) = 1.0 pctlwp(i) = 0.0 ENDDO !C !c--calculation of liquid water path !c DO k = klev, 1, -1 DO i = 1, klon pctlwp(i) = pctlwp(i)+ pqlwp(i,k)*rhodz(i,k) ENDDO ENDDO !c !c--calculation of cloud properties with cloud overlap !c IF (NOVLP.EQ.1) THEN DO k = klev, 1, -1 DO i = 1, klon zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i))) & & /(1.-MIN(REAL(zcloud(i), kind=8),1.-ZEPSEC)) pct(i)=1.-zclear(i) IF (paprs(i,k).LT.prmhc) THEN pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloudh(i))) & & /(1.-MIN(REAL(zcloudh(i), kind=8),1.-ZEPSEC)) zcloudh(i)=pclc(i,k) ELSE IF (paprs(i,k).GE.prmhc .AND. & & paprs(i,k).LT.prlmc) THEN pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloudm(i))) & & /(1.-MIN(REAL(zcloudm(i), kind=8),1.-ZEPSEC)) zcloudm(i)=pclc(i,k) ELSE IF (paprs(i,k).GE.prlmc) THEN pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloudl(i))) & & /(1.-MIN(REAL(zcloudl(i), kind=8),1.-ZEPSEC)) zcloudl(i)=pclc(i,k) endif zcloud(i)=pclc(i,k) ENDDO ENDDO ELSE IF (NOVLP.EQ.2) THEN DO k = klev, 1, -1 DO i = 1, klon zcloud(i)=MAX(pclc(i,k),zcloud(i)) pct(i)=zcloud(i) IF (paprs(i,k).LT.prmhc) THEN pch(i) = MIN(pclc(i,k),pch(i)) ELSE IF (paprs(i,k).GE.prmhc .AND. & & paprs(i,k).LT.prlmc) THEN pcm(i) = MIN(pclc(i,k),pcm(i)) ELSE IF (paprs(i,k).GE.prlmc) THEN pcl(i) = MIN(pclc(i,k),pcl(i)) endif ENDDO ENDDO ELSE IF (NOVLP.EQ.3) THEN DO k = klev, 1, -1 DO i = 1, klon zclear(i)=zclear(i)*(1.-pclc(i,k)) pct(i)=1-zclear(i) IF (paprs(i,k).LT.prmhc) THEN pch(i) = pch(i)*(1.0-pclc(i,k)) ELSE IF (paprs(i,k).GE.prmhc .AND. & & paprs(i,k).LT.prlmc) THEN pcm(i) = pcm(i)*(1.0-pclc(i,k)) ELSE IF (paprs(i,k).GE.prlmc) THEN pcl(i) = pcl(i)*(1.0-pclc(i,k)) endif ENDDO ENDDO ENDIF !C DO i = 1, klon pch(i)=1.-pch(i) pcm(i)=1.-pcm(i) pcl(i)=1.-pcl(i) ENDDO !c !c ======================================================== !c DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL !c ======================================================== !c change by Nicolas Yan (LSCE) !c Cloud Droplet Number Concentration (CDNC) : 3D variable !c Fractionnal cover by liquid water cloud (LCC3D) : 3D variable !c Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable !c Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable !c Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable !c IF (ok_cdnc) THEN !c DO k = 1, klev DO i = 1, klon phase3d(i,k)=1-zfice(i,k) IF (pclc(i,k) .LE. seuil_neb) THEN lcc3d(i,k)=seuil_neb*phase3d(i,k) ELSE lcc3d(i,k)=pclc(i,k)*phase3d(i,k) ENDIF scdnc(i,k)=lcc3d(i,k)*cdnc(i,k) ! m-3 ENDDO ENDDO !c DO i=1,klon lcc(i)=0. reffclwtop(i)=0. cldncl(i)=0. IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i) = 1. IF(MAXIMUM) tcc(i) = 0. ENDDO !c DO i=1,klon DO k=klev-1,1,-1 !From TOA down !c ! Test, if the cloud optical depth exceeds the necessary ! threshold: IF (pcltau(i,k).GT.thres_tau & & .AND. pclc(i,k).GT.thres_neb) THEN IF (MAXIMUM) THEN IF (FIRST) THEN write(*,*)'Hypothese de recouvrement: MAXIMUM' FIRST=.FALSE. ENDIF flag_max= -1. ftmp(i) = MAX(tcc(i),pclc(i,k)) ENDIF IF (RANDOM) THEN IF (FIRST) THEN write(*,*)'Hypothese de recouvrement: RANDOM' FIRST=.FALSE. ENDIF flag_max= 1. ftmp(i) = tcc(i) * (1-pclc(i,k)) ENDIF IF (MAXIMUM_RANDOM) THEN IF (FIRST) THEN write(*,*)'Hypothese de recouvrement: MAXIMUM_ & & RANDOM' FIRST=.FALSE. ENDIF flag_max= 1. ftmp(i) = tcc(i) * & & (1. - MAX(pclc(i,k),pclc(i,k+1))) / & & (1. - MIN(pclc(i,k+1),1.-thres_neb)) ENDIF !c Effective radius of cloud droplet at top of cloud (m) reffclwtop(i) = reffclwtop(i) + rad_chaud(i,k) * & & 1.0E-06 * phase3d(i,k) * ( tcc(i) - ftmp(i))*flag_max !c CDNC at top of cloud (m-3) cldncl(i) = cldncl(i) + cdnc(i,k) * phase3d(i,k) * & & (tcc(i) - ftmp(i))*flag_max !c Liquid Cloud Content at top of cloud lcc(i) = lcc(i) + phase3d(i,k) * (tcc(i)-ftmp(i))* & & flag_max !c Total Cloud Content at top of cloud tcc(i)=ftmp(i) ENDIF ! is there a visible, not-too-small cloud? ENDDO ! loop over k !c IF (RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i) !c ENDDO ! loop over i !! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC REFFCLWS) DO i = 1, klon DO k = 1, klev ! Weight to be used for outputs: eau_liquide*couverture nuageuse lcc3dcon(i,k) =rnebcon(i,k)*phase3d(i,k)*clwcon(i,k) ! eau liquide convective lcc3dstra(i,k)=pclc(i,k)*pqlwp(i,k)*phase3d(i,k) lcc3dstra(i,k)=lcc3dstra(i,k)-lcc3dcon(i,k) ! eau liquide stratiforme lcc3dstra(i,k)=MAX(lcc3dstra(i,k),0.0) ! Compute cloud droplet radius as above in meter radius=1.1*((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) & & /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.) radius=MAX(radius, 5.e-6) ! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D reffclwc(i,k)=radius reffclwc(i,k)=reffclwc(i,k)*lcc3dcon(i,k) ! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D reffclws(i,k)=radius reffclws(i,k)=reffclws(i,k)*lcc3dstra(i,k) ENDDO !klev ENDDO !klon !c !c Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D !c DO i = 1, klon cldnvi(i)=0. lcc_integrat(i)=0. height(i)=0. DO k = 1, klev cldnvi(i)=cldnvi(i)+cdnc(i,k)*lcc3d(i,k)*dh(i,k) lcc_integrat(i)=lcc_integrat(i)+lcc3d(i,k)*dh(i,k) height(i)=height(i)+dh(i,k) ENDDO ! klev lcc_integrat(i)=lcc_integrat(i)/height(i) IF (lcc_integrat(i) .LE. 1.0E-03) THEN cldnvi(i)=cldnvi(i)*lcc(i)/seuil_neb ELSE cldnvi(i)=cldnvi(i)*lcc(i)/lcc_integrat(i) ENDIF ENDDO ! klon DO i = 1, klon DO k = 1, klev IF (scdnc(i,k) .LE. 0.0) scdnc(i,k)=0.0 IF (reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0 IF (reffclwc(i,k) .LE. 0.0) reffclwc(i,k)=0.0 IF (lcc3d(i,k) .LE. 0.0) lcc3d(i,k)=0.0 IF (lcc3dcon(i,k) .LE. 0.0) lcc3dcon(i,k)=0.0 IF (lcc3dstra(i,k) .LE. 0.0) lcc3dstra(i,k)=0.0 ENDDO IF (reffclwtop(i) .LE. 0.0) reffclwtop(i)=0.0 IF (cldncl(i) .LE. 0.0) cldncl(i)=0.0 IF (cldnvi(i) .LE. 0.0) cldnvi(i)=0.0 IF (lcc(i) .LE. 0.0) lcc(i)=0.0 ENDDO !c ENDIF !ok_cdnc !c RETURN !c END SUBROUTINE newmicro