! $Id: nuage.f90 5305 2024-10-30 18:29:21Z evignon $ SUBROUTINE nuage(paprs, pplay, t, pqlwp,picefra, pclc, pcltau, pclemi, pch, pcl, pcm, & pct, pctlwp, ok_aie, mass_solu_aero, mass_solu_aero_pi, bl95_b0, bl95_b1, distcltop, & temp_cltop, cldtaupi, re, fl) USE clesphys_mod_h USE yomcst_mod_h USE dimphy USE lmdz_lscp_tools, only: icefrac_lscp USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14) USE lmdz_lscp_ini, only : iflag_t_glace USE phys_local_var_mod, ONLY: ptconv USE nuage_params_mod_h IMPLICIT NONE ! ====================================================================== ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910 ! Objet: Calculer epaisseur optique et emmissivite des nuages ! ====================================================================== ! Arguments: ! t-------input-R-temperature ! pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg) ! picefra--inout-R-fraction de glace dans les nuages (-) ! pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1) ! ok_aie--input-L-apply aerosol indirect effect or not ! mass_solu_aero-----input-R-total mass concentration for all soluble ! aerosols[ug/m^3] ! mass_solu_aero_pi--input-R-dito, pre-industrial value ! bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land) ! bl95_b1-input-R-a parameter, may be varied for tests ( -"- ) ! cldtaupi-output-R-pre-industrial value of cloud optical thickness, ! needed for the diagnostics of the aerosol indirect ! radiative forcing (see radlwsw) ! re------output-R-Cloud droplet effective radius multiplied by fl [um] ! fl------output-R-Denominator to re, introduced to avoid problems in ! the averaging of the output. fl is the fraction of liquid ! water clouds within a grid cell ! pcltau--output-R-epaisseur optique des nuages ! pclemi--output-R-emissivite des nuages (0 a 1) ! ====================================================================== REAL paprs(klon, klev+1), pplay(klon, klev) REAL t(klon, klev) REAL pclc(klon, klev) REAL pqlwp(klon, klev), picefra(klon,klev) REAL pcltau(klon, klev), pclemi(klon, klev) REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon) REAL distcltop(klon,klev) REAL temp_cltop(klon,klev) LOGICAL lo REAL cetahb, cetamb PARAMETER (cetahb=0.45, cetamb=0.80) INTEGER i, k REAL zflwp, zradef, zfice(klon), zmsac REAL radius, rad_chaud ! JBM (3/14) parameters already defined in nuage.h: ! REAL rad_froid, rad_chau1, rad_chau2 ! PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0) ! cc PARAMETER (rad_chaud=15.0, rad_froid=35.0) ! sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0) REAL coef, coef_froi, coef_chau PARAMETER (coef_chau=0.13, coef_froi=0.09) REAL seuil_neb PARAMETER (seuil_neb=0.001) ! JBM (3/14) nexpo is replaced by exposant_glace ! REAL nexpo ! exponentiel pour glace/eau ! PARAMETER (nexpo=6.) REAL, PARAMETER :: t_glace_min_old = 258. INTEGER, PARAMETER :: exposant_glace_old = 6 ! jq for the aerosol indirect effect ! jq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003 ! jq LOGICAL ok_aie ! Apply AIE or not? 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) REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag REAl dzfice(klon) ! jq-end ! cc PARAMETER (nexpo=1) ! Calculer l'epaisseur optique et l'emmissivite des nuages DO k = 1, klev IF (iflag_t_glace.EQ.0) THEN DO i = 1, klon zfice(i) = 1.0 - (t(i,k)-t_glace_min_old)/(273.13-t_glace_min_old) zfice(i) = min(max(zfice(i),0.0), 1.0) zfice(i) = zfice(i)**exposant_glace_old ENDDO ELSE ! of IF (iflag_t_glace.EQ.0) ! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod ! zfice(i) = icefrac_lsc(t(i,k), t_glace_min, & ! t_glace_max, exposant_glace) IF (ok_new_lscp) THEN CALL icefrac_lscp(klon,t(:,k),iflag_ice_thermo,distcltop(:,k),temp_cltop(:,k),zfice(:),dzfice(:)) ELSE CALL icefrac_lsc(klon,t(:,k),pplay(:,k)/paprs(:,1),zfice(:)) ENDIF IF (ok_new_lscp .AND. ok_icefra_lscp) THEN ! EV: take the ice fraction directly from the lscp code ! consistent only for non convective grid points ! critical for mixed phase clouds DO i=1,klon IF (.NOT. ptconv(i,k)) THEN zfice(i)=picefra(i,k) ENDIF ENDDO ENDIF ENDIF DO i = 1, klon rad_chaud = rad_chau1 IF (k<=3) rad_chaud = rad_chau2 pclc(i, k) = max(pclc(i,k), seuil_neb) zflwp = 1000.*pqlwp(i, k)/rg/pclc(i, k)*(paprs(i,k)-paprs(i,k+1)) IF (ok_aie) THEN ! Formula "D" of Boucher and Lohmann, Tellus, 1995 ! cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), & 1.E-4))/log(10.))*1.E6 !-m-3 ! Cloud droplet number concentration (CDNC) is restricted ! to be within [20, 1000 cm^3] ! cdnc(i, k) = min(1000.E6, max(20.E6,cdnc(i,k))) 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))) ! ! ! air density: pplay(i,k) / (RD * zT(i,k)) ! factor 1.1: derive effective radius from volume-mean radius ! factor 1000 is the water density ! _chaud means that this is the CDR for liquid water clouds ! rad_chaud = 1.1*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi*1000. & *cdnc(i,k)))**(1./3.) ! ! Convert to um. CDR shall be at least 3 um. ! rad_chaud = max(rad_chaud*1.E6, 3.) ! For output diagnostics ! ! Cloud droplet effective radius [um] ! ! we multiply here with f * xl (fraction of liquid water ! clouds in the grid cell) to avoid problems in the ! averaging of the output. ! In the output of IOIPSL, derive the real cloud droplet ! effective radius as re/fl ! fl(i, k) = pclc(i, k)*(1.-zfice(i)) re(i, k) = rad_chaud*fl(i, k) ! Pre-industrial cloud opt thickness ! ! "radius" is calculated as rad_chaud above (plus the ! ice cloud contribution) but using cdnc_pi instead of ! cdnc. radius = max(1.1E6*((pqlwp(i,k)*pplay(i,k)/(rd*t(i,k)))/(4./3.*rpi* & 1000.*cdnc_pi(i,k)))**(1./3.), 3.)*(1.-zfice(i)) + rad_froid*zfice(i) cldtaupi(i, k) = 3.0/2.0*zflwp/radius END IF ! ok_aie radius = rad_chaud*(1.-zfice(i)) + rad_froid*zfice(i) coef = coef_chau*(1.-zfice(i)) + coef_froi*zfice(i) pcltau(i, k) = 3.0/2.0*zflwp/radius pclemi(i, k) = 1.0 - exp(-coef*zflwp) lo = (pclc(i,k)<=seuil_neb) IF (lo) pclc(i, k) = 0.0 IF (lo) pcltau(i, k) = 0.0 IF (lo) pclemi(i, k) = 0.0 IF (.NOT. ok_aie) cldtaupi(i, k) = pcltau(i, k) END DO END DO ! cc DO k = 1, klev ! cc DO i = 1, klon ! cc t(i,k) = t(i,k) ! cc pclc(i,k) = MAX( 1.e-5 , pclc(i,k) ) ! cc lo = pclc(i,k) .GT. (2.*1.e-5) ! cc zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1)) ! cc . /(rg*pclc(i,k)) ! cc zradef = 10.0 + (1.-sigs(k))*45.0 ! cc pcltau(i,k) = 1.5 * zflwp / zradef ! cc zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0) ! cc zmsac = 0.13*(1.0-zfice) + 0.08*zfice ! cc pclemi(i,k) = 1.-EXP(-zmsac*zflwp) ! cc if (.NOT.lo) pclc(i,k) = 0.0 ! cc if (.NOT.lo) pcltau(i,k) = 0.0 ! cc if (.NOT.lo) pclemi(i,k) = 0.0 ! cc ENDDO ! cc ENDDO ! ccccc print*, 'pas de nuage dans le rayonnement' ! ccccc DO k = 1, klev ! ccccc DO i = 1, klon ! ccccc pclc(i,k) = 0.0 ! ccccc pcltau(i,k) = 0.0 ! ccccc pclemi(i,k) = 0.0 ! ccccc ENDDO ! ccccc ENDDO ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS DO i = 1, klon pct(i) = 1.0 pch(i) = 1.0 pcm(i) = 1.0 pcl(i) = 1.0 pctlwp(i) = 0.0 END DO DO k = klev, 1, -1 DO i = 1, klon pctlwp(i) = pctlwp(i) + pqlwp(i, k)*(paprs(i,k)-paprs(i,k+1))/rg pct(i) = pct(i)*(1.0-pclc(i,k)) IF (pplay(i,k)<=cetahb*paprs(i,1)) pch(i) = pch(i)*(1.0-pclc(i,k)) IF (pplay(i,k)>cetahb*paprs(i,1) .AND. pplay(i,k)<=cetamb*paprs(i,1)) & pcm(i) = pcm(i)*(1.0-pclc(i,k)) IF (pplay(i,k)>cetamb*paprs(i,1)) pcl(i) = pcl(i)*(1.0-pclc(i,k)) END DO END DO DO i = 1, klon pct(i) = 1. - pct(i) pch(i) = 1. - pch(i) pcm(i) = 1. - pcm(i) pcl(i) = 1. - pcl(i) END DO RETURN END SUBROUTINE nuage SUBROUTINE diagcld1(paprs, pplay, rain, snow, kbot, ktop, diafra, dialiq) USE dimphy USE yomcst_mod_h IMPLICIT NONE ! Laurent Li (LMD/CNRS), le 12 octobre 1998 ! (adaptation du code ECMWF) ! Dans certains cas, le schema pronostique des nuages n'est ! pas suffisament performant. On a donc besoin de diagnostiquer ! ces nuages. Je dois avouer que c'est une frustration. ! Arguments d'entree: REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche REAL pplay(klon, klev) ! pression (Pa) au milieu de couche REAL t(klon, klev) ! temperature (K) REAL q(klon, klev) ! humidite specifique (Kg/Kg) REAL rain(klon) ! pluie convective (kg/m2/s) REAL snow(klon) ! neige convective (kg/m2/s) INTEGER ktop(klon) ! sommet de la convection INTEGER kbot(klon) ! bas de la convection ! Arguments de sortie: REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee REAL dialiq(klon, klev) ! eau liquide nuageuse ! Constantes ajustables: REAL canva, canvb, canvh PARAMETER (canva=2.0, canvb=0.3, canvh=0.4) REAL cca, ccb, ccc PARAMETER (cca=0.125, ccb=1.5, ccc=0.8) REAL ccfct, ccscal PARAMETER (ccfct=0.400) PARAMETER (ccscal=1.0E+11) REAL cetahb, cetamb PARAMETER (cetahb=0.45, cetamb=0.80) REAL cclwmr PARAMETER (cclwmr=1.E-04) REAL zepscr PARAMETER (zepscr=1.0E-10) ! Variables locales: INTEGER i, k REAL zcc(klon) ! Initialisation: DO k = 1, klev DO i = 1, klon diafra(i, k) = 0.0 dialiq(i, k) = 0.0 END DO END DO DO i = 1, klon ! Calculer la fraction nuageuse zcc(i) = 0.0 IF ((rain(i)+snow(i))>0.) THEN zcc(i) = cca*log(max(zepscr,(rain(i)+snow(i))*ccscal)) - ccb zcc(i) = min(ccc, max(0.0,zcc(i))) END IF END DO DO i = 1, klon ! pour traiter les enclumes diafra(i, ktop(i)) = max(diafra(i,ktop(i)), zcc(i)*ccfct) IF ((zcc(i)>=canvh) .AND. (pplay(i,ktop(i))<=cetahb*paprs(i, & 1))) diafra(i, ktop(i)) = max(diafra(i,ktop(i)), max(zcc( & i)*ccfct,canva*(zcc(i)-canvb))) dialiq(i, ktop(i)) = cclwmr*diafra(i, ktop(i)) END DO DO k = 1, klev ! nuages convectifs (sauf enclumes) DO i = 1, klon IF (k=kbot(i)) THEN diafra(i, k) = max(diafra(i,k), zcc(i)*ccfct) dialiq(i, k) = cclwmr*diafra(i, k) END IF END DO END DO RETURN END SUBROUTINE diagcld1 SUBROUTINE diagcld2(paprs, pplay, t, q, diafra, dialiq) USE dimphy USE yomcst_mod_h USE yoethf_mod_h IMPLICIT NONE ! Arguments d'entree: REAL paprs(klon, klev+1) ! pression (Pa) a inter-couche REAL pplay(klon, klev) ! pression (Pa) au milieu de couche REAL t(klon, klev) ! temperature (K) REAL q(klon, klev) ! humidite specifique (Kg/Kg) ! Arguments de sortie: REAL diafra(klon, klev) ! fraction nuageuse diagnostiquee REAL dialiq(klon, klev) ! eau liquide nuageuse REAL cetamb PARAMETER (cetamb=0.80) REAL cloia, cloib, cloic, cloid PARAMETER (cloia=1.0E+02, cloib=-10.00, cloic=-0.6, cloid=5.0) ! cc PARAMETER (CLOIA=1.0E+02, CLOIB=-10.00, CLOIC=-0.9, CLOID=5.0) REAL rgammas PARAMETER (rgammas=0.05) REAL crhl PARAMETER (crhl=0.15) ! cc PARAMETER (CRHL=0.70) REAL t_coup PARAMETER (t_coup=234.0) ! Variables locales: INTEGER i, k, kb, invb(klon) REAL zqs, zrhb, zcll, zdthmin(klon), zdthdp REAL zdelta, zcor ! Fonctions thermodynamiques: include "FCTTRE.h" ! Initialisation: DO k = 1, klev DO i = 1, klon diafra(i, k) = 0.0 dialiq(i, k) = 0.0 END DO END DO DO i = 1, klon invb(i) = klev zdthmin(i) = 0.0 END DO DO k = 2, klev/2 - 1 DO i = 1, klon zdthdp = (t(i,k)-t(i,k+1))/(pplay(i,k)-pplay(i,k+1)) - & rd*0.5*(t(i,k)+t(i,k+1))/rcpd/paprs(i, k+1) zdthdp = zdthdp*cloia IF (pplay(i,k)>cetamb*paprs(i,1) .AND. zdthdp0.0 .AND. zrhb