Index: LMDZ5/trunk/libf/phylmd/conf_phys_m.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/conf_phys_m.F90 (revision 1711) +++ LMDZ5/trunk/libf/phylmd/conf_phys_m.F90 (revision 1712) @@ -18,5 +18,5 @@ iflag_cldcon, & iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, & - ok_ade, ok_aie, aerosol_couple, & + ok_ade, ok_aie, ok_cdnc, aerosol_couple, & flag_aerosol, new_aod, & bl95_b0, bl95_b1,& @@ -60,4 +60,5 @@ ! ok_instan: sorties instantanees ! ok_ade, ok_aie: apply or not aerosol direct and indirect effects +! ok_cdnc, ok cloud droplet number concentration ! bl95_b*: parameters in the formula to link CDNC to aerosol mass conc ! @@ -70,5 +71,5 @@ logical :: ok_LES LOGICAL :: callstats - LOGICAL :: ok_ade, ok_aie, aerosol_couple + LOGICAL :: ok_ade, ok_aie, ok_cdnc, aerosol_couple INTEGER :: flag_aerosol LOGICAL :: new_aod @@ -84,5 +85,5 @@ logical,SAVE :: ok_LES_omp LOGICAL,SAVE :: callstats_omp - LOGICAL,SAVE :: ok_ade_omp, ok_aie_omp, aerosol_couple_omp + LOGICAL,SAVE :: ok_ade_omp, ok_aie_omp, ok_cdnc_omp, aerosol_couple_omp INTEGER, SAVE :: flag_aerosol_omp LOGICAL, SAVE :: new_aod_omp @@ -272,4 +273,12 @@ call getin('ok_aie', ok_aie_omp) +! +!Config Key = ok_cdnc +!Config Desc = ok cloud droplet number concentration +!Config Def = .false. +!Config Help = Used in newmicro.F +! + ok_cdnc_omp = .false. + call getin('ok_cdnc', ok_cdnc_omp) ! !Config Key = aerosol_couple @@ -1677,4 +1686,5 @@ ok_ade = ok_ade_omp ok_aie = ok_aie_omp + ok_cdnc = ok_cdnc_omp aerosol_couple = aerosol_couple_omp flag_aerosol=flag_aerosol_omp @@ -1774,4 +1784,9 @@ END IF END IF + +! ok_cdnc must be set to y if ok_aie is activated + IF (ok_aie .AND. .NOT. ok_cdnc) THEN + CALL abort_gcm('conf_phys', 'ok_cdnc must be set to y if ok_aie is activated',1) + ENDIF !$OMP MASTER Index: LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90 (revision 1711) +++ LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90 (revision 1712) @@ -99,5 +99,5 @@ REAL :: dummy LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf - LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats + LOGICAL :: ok_LES, ok_ade, ok_aie, ok_cdnc, aerosol_couple, new_aod, callstats INTEGER :: iflag_radia, flag_aerosol REAL :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut @@ -136,5 +136,5 @@ iflag_cldcon, & iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, & - ok_ade, ok_aie, aerosol_couple, & + ok_ade, ok_aie, ok_cdnc, aerosol_couple, & flag_aerosol, new_aod, & bl95_b0, bl95_b1, & Index: LMDZ5/trunk/libf/phylmd/newmicro.F =================================================================== --- LMDZ5/trunk/libf/phylmd/newmicro.F (revision 1711) +++ LMDZ5/trunk/libf/phylmd/newmicro.F (revision 1712) @@ -2,15 +2,13 @@ - ! - SUBROUTINE newmicro (paprs, pplay,ok_newmicro, + SUBROUTINE newmicro (ok_cdnc, bl95_b0, bl95_b1, + . paprs, pplay, . t, pqlwp, pclc, pcltau, pclemi, . pch, pcl, pcm, pct, pctlwp, - s xflwp, xfiwp, xflwc, xfiwc, - e ok_aie, - e mass_solu_aero, mass_solu_aero_pi, - e bl95_b0, bl95_b1, - s cldtaupi, re, fl, reliq, reice) - + . 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, @@ -21,61 +19,68 @@ 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 (kg/kg) +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 -c ok_aie--input-L-apply aerosol indirect effect or not c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3] -c mass_solu_aero_pi--input-R-dito, pre-industrial value -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 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 cldtaupi-output-R-pre-industrial value of cloud optical thickness, -c needed for the diagnostics of the aerosol indirect -c radiative forcing (see radlwsw) 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" -c -cym#include "dimensions.h" -cym#include "dimphy.h" #include "nuage.h" -cIM cf. CR: include pour NOVLP et ZEPSEC #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.) + LOGICAL RANDOM, MAXIMUM_RANDOM, MAXIMUM + PARAMETER (RANDOM=.FALSE., MAXIMUM_RANDOM=.TRUE., MAXIMUM=.FALSE.) +c LOGICAL, SAVE :: FIRST=.TRUE. !$OMP THREADPRIVATE(FIRST) -c Hypoyhese de recouvrement : MAXIMUM_RANDOM INTEGER flag_max - REAL phase3d(klon, klev),dh(klon, klev),pdel(klon, klev), - . zrho(klon, klev) - REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon) +c +c threshold PARAMETERs REAL thres_tau,thres_neb PARAMETER (thres_tau=0.3, thres_neb=0.001) - REAL t_tmp - REAL gravit - PARAMETER (gravit=9.80616) !m/s2 - REAL pqlwpcon(klon, klev), pqlwpstra(klon, klev) -c - REAL paprs(klon,klev+1), pplay(klon,klev) +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) -c REAL pclc(klon,klev) REAL pqlwp(klon,klev) - REAL pcltau(klon,klev), pclemi(klon,klev) -c - REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon) + 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 @@ -85,45 +90,34 @@ ! 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 - REAL prlmc ! Pressure between low and medium level cloud +! 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 -cIM: 091003 REAL zflwp, zradef, zfice, zmsac - REAL zflwp(klon), zradef, zfice, zmsac -cIM: 091003 rajout REAL xflwp(klon), xfiwp(klon) REAL xflwc(klon,klev), xfiwc(klon,klev) c - REAL radius, rad_chaud -cc PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0) -ccc PARAMETER (rad_chaud=15.0, rad_froid=35.0) -c sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0) - REAL coef, coef_froi, coef_chau + 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) -ccc PARAMETER (nexpo=1) - -c -- sb: - logical ok_newmicro -c parameter (ok_newmicro=.FALSE.) -cIM: 091003 real rel, tc, rei, zfiwp - real rel, tc, rei, zfiwp(klon) - real k_liq, k_ice0, k_ice, DF - parameter (k_liq=0.0903, k_ice0=0.005) ! units=m2/g - parameter (DF=1.66) ! diffusivity factor -c sb -- +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 - LOGICAL ok_aie ! Apply AIE or not? - LOGICAL ok_a1lwpdep ! a1 LWP dependent? - 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) @@ -135,7 +129,7 @@ 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 - REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag cjq-end cIM cf. CR:parametres supplementaires @@ -145,14 +139,10 @@ REAL zcloudm(klon) REAL zcloudl(klon) - - -c ************************** -c * * -c * DEBUT PARTIE OPTIMISEE * -c * * -c ************************** - - REAL diff_paprs(klon, klev), zfice1, zfice2(klon, klev) - REAL rad_chaud_tab(klon, klev), zflwp_var, zfiwp_var + 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 @@ -171,286 +161,213 @@ ! 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 -! print*,'d_rei_dT ',d_rei_dt,rei_min,rei_max !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c c Calculer l'epaisseur optique et l'emmissivite des nuages -c c IM inversion des DO +c xflwp = 0.d0 xfiwp = 0.d0 xflwc = 0.d0 xfiwc = 0.d0 - -! Initialisation +c reliq=0. reice=0. - +c DO k = 1, klev - DO i = 1, klon - diff_paprs(i,k) = (paprs(i,k)-paprs(i,k+1))/RG + 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_newmicro) THEN - - - DO k = 1, klev - DO i = 1, klon -c zfice2(i,k) = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace) - zfice2(i,k) = 1.0 - (t(i,k)-t_glace_min) / - & (t_glace_max-t_glace_min) - zfice2(i,k) = MIN(MAX(zfice2(i,k),0.0),1.0) -c IM Total Liquid/Ice water content - xflwc(i,k) = (1.-zfice2(i,k))*pqlwp(i,k) - xfiwc(i,k) = zfice2(i,k)*pqlwp(i,k) -c IM In-Cloud Liquid/Ice water content -c xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k) -c xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k) - ENDDO - ENDDO - - IF (ok_aie) THEN - DO k = 1, klev - DO i = 1, klon - ! Formula "D" of Boucher and Lohmann, Tellus, 1995 - ! - cdnc(i,k) = 10.**(bl95_b0+bl95_b1* + 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 - ! 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* + 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))) - ENDDO - ENDDO - DO k = 1, klev - DO i = 1, klon -! rad_chaud_tab(i,k) = -! & MAX(1.1e6 -! & *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) -! & /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.),5.) - rad_chaud_tab(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_tab(i,k) = MAX(rad_chaud_tab(i,k) * 1e6, 5.) - ENDDO - ENDDO - ELSE - DO k = 1, MIN(3,klev) - DO i = 1, klon - rad_chaud_tab(i,k) = rad_chau2 - ENDDO - ENDDO - DO k = MIN(3,klev)+1, klev - DO i = 1, klon - rad_chaud_tab(i,k) = rad_chau1 - ENDDO - ENDDO - - ENDIF - - DO k = 1, klev -! IF(.not.ok_aie) THEN - rad_chaud = rad_chau1 - IF (k.LE.3) rad_chaud = rad_chau2 -! ENDIF - DO i = 1, klon - IF (pclc(i,k) .LE. seuil_neb) THEN - -c -- effective cloud droplet radius (microns): - -c for liquid water clouds: - ! 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) = seuil_neb*(1.-zfice2(i,k)) - re(i,k) = rad_chaud_tab(i,k)*fl(i,k) - - rel = 0. - rei = 0. - pclc(i,k) = 0.0 - pcltau(i,k) = 0.0 - pclemi(i,k) = 0.0 - cldtaupi(i,k) = 0.0 - ELSE - + 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) + 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.-zfice2(i,k))*pqlwp(i,k)/pclc(i,k) - & *diff_paprs(i,k) - zfiwp_var= 1000.*zfice2(i,k)*pqlwp(i,k)/pclc(i,k) - & *diff_paprs(i,k) - -c -- effective cloud droplet radius (microns): - -c for liquid water clouds: - - IF (ok_aie) THEN - 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*1e6, 5.) - - tc = t(i,k)-273.15 - rei = d_rei_dt*tc + rei_max - if (tc.le.-81.4) rei = rei_min - if (zflwp_var.eq.0.) radius = 1. - if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1. - cldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius - & + zfiwp_var * (3.448e-03 + 2.431/rei) - - ENDIF ! ok_aie - ! 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.-zfice2(i,k)) - re(i,k) = rad_chaud_tab(i,k)*fl(i,k) - - rel = rad_chaud_tab(i,k) -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 - tc = t(i,k)-273.15 - rei = d_rei_dt*tc + rei_max - if (tc.le.-81.4) rei = rei_min -c -- cloud optical thickness : - -c [for liquid clouds, traditional formula, -c for ice clouds, Ebert & Curry (1992)] - + 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 [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) - - ENDIF - reliq(i,k)=rel - reice(i,k)=rei -! if (i.eq.1) then -! print*,'Dans newmicro rel, rei :',rel, rei -! print*,'Dans newmicro reliq, reice :', -! $ reliq(i,k),reice(i,k) -! endif - - ENDDO - ENDDO - +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 - xflwp(i) = xflwp(i)+ xflwc(i,k) * diff_paprs(i,k) - xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * diff_paprs(i,k) - ENDDO - ENDDO - - ELSE - DO k = 1, klev - rad_chaud = rad_chau1 - IF (k.LE.3) rad_chaud = rad_chau2 - DO i = 1, klon - - IF (pclc(i,k) .LE. seuil_neb) THEN - - pclc(i,k) = 0.0 - pcltau(i,k) = 0.0 - pclemi(i,k) = 0.0 - cldtaupi(i,k) = 0.0 - - ELSE - - zflwp_var = 1000.*pqlwp(i,k)*diff_paprs(i,k) - & /pclc(i,k) - - zfice1 = MIN( - & MAX( 1.0 - (t(i,k)-t_glace_min) / - & (t_glace_max-t_glace_min),0.0),1.0)**nexpo - - radius = rad_chaud * (1.-zfice1) + rad_froid * zfice1 - coef = coef_chau * (1.-zfice1) + coef_froi * zfice1 - - pcltau(i,k) = 3.0 * zflwp_var / (2.0 * radius) - pclemi(i,k) = 1.0 - EXP( - coef * zflwp_var) - - ENDIF - - ENDDO - ENDDO - ENDIF - - IF (.NOT.ok_aie) THEN - DO k = 1, klev - DO i = 1, klon - cldtaupi(i,k)=pcltau(i,k) + pcldtaupi(i,k)=pcltau(i,k) ENDDO ENDDO ENDIF - -ccc DO k = 1, klev -ccc DO i = 1, klon -ccc t(i,k) = t(i,k) -ccc pclc(i,k) = MAX( 1.e-5 , pclc(i,k) ) -ccc lo = pclc(i,k) .GT. (2.*1.e-5) -ccc zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1)) -ccc . /(rg*pclc(i,k)) -ccc zradef = 10.0 + (1.-sigs(k))*45.0 -ccc pcltau(i,k) = 1.5 * zflwp / zradef -ccc zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0) -ccc zmsac = 0.13*(1.0-zfice) + 0.08*zfice -ccc pclemi(i,k) = 1.-EXP(-zmsac*zflwp) -ccc if (.NOT.lo) pclc(i,k) = 0.0 -ccc if (.NOT.lo) pcltau(i,k) = 0.0 -ccc if (.NOT.lo) pclemi(i,k) = 0.0 -ccc ENDDO -ccc 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 C C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS -C c IM cf. CR:test: calcul prenant ou non en compte le recouvrement c initialisations +c DO i=1,klon zclear(i)=1. @@ -465,30 +382,32 @@ ENDDO C -cIM cf CR DO k=1,klev +c--calculation of liquid water path +c DO k = klev, 1, -1 DO i = 1, klon - pctlwp(i) = pctlwp(i) - & + pqlwp(i,k)*diff_paprs(i,k) + pctlwp(i) = pctlwp(i)+ pqlwp(i,k)*rhodz(i,k) ENDDO ENDDO -c IM cf. CR +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)) + & /(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)) + & /(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)) + & /(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)) + & /(1.-MIN(REAL(zcloudl(i), kind=8),1.-ZEPSEC)) zcloudl(i)=pclc(i,k) endif @@ -527,27 +446,26 @@ ENDDO ENDIF - C DO i = 1, klon -c IM cf. CR pct(i)=1.-pct(i) pch(i)=1.-pch(i) pcm(i)=1.-pcm(i) pcl(i)=1.-pcl(i) ENDDO - +c c ======================================================== -! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL +c DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL c ======================================================== -!! change by Nicolas Yan (LSCE) -!! Cloud Droplet Number Concentration (CDNC) : 3D variable -!! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable -!! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable -!! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable -!! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable - IF (ok_newmicro) THEN - IF (ok_aie) THEN +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-zfice2(i,k) + phase3d(i,k)=1-zfice(i,k) IF (pclc(i,k) .LE. seuil_neb) THEN lcc3d(i,k)=seuil_neb*phase3d(i,k) @@ -558,5 +476,5 @@ ENDDO ENDDO - +c DO i=1,klon lcc(i)=0. @@ -566,23 +484,16 @@ 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 - ! To calculate the right Temperature at cloud top, - ! interpolate it between layers: - t_tmp = t(i,k) + - . (paprs(i,k+1)-pplay(i,k))/(pplay(i,k+1)-pplay(i,k)) - . * ( t(i,k+1) - t(i,k) ) - - IF(MAXIMUM) THEN - IF(FIRST) THEN + 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. @@ -592,6 +503,6 @@ ENDIF - IF(RANDOM) THEN - IF(FIRST) THEN + IF (RANDOM) THEN + IF (FIRST) THEN write(*,*)'Hypothese de recouvrement: RANDOM' FIRST=.FALSE. @@ -601,6 +512,6 @@ ENDIF - IF(MAXIMUM_RANDOM) THEN - IF(FIRST) THEN + IF (MAXIMUM_RANDOM) THEN + IF (FIRST) THEN write(*,*)'Hypothese de recouvrement: MAXIMUM_ . RANDOM' @@ -613,5 +524,5 @@ ENDIF c Effective radius of cloud droplet at top of cloud (m) - reffclwtop(i) = reffclwtop(i) + rad_chaud_tab(i,k) * + 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) @@ -626,6 +537,7 @@ ENDIF ! is there a visible, not-too-small cloud? ENDDO ! loop over k - - IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i) +c + IF (RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i) +c ENDDO ! loop over i @@ -633,48 +545,23 @@ DO i = 1, klon DO k = 1, klev - pqlwpcon(i,k)=rnebcon(i,k)*clwcon(i,k) ! fraction eau liquide convective - pqlwpstra(i,k)=pclc(i,k)*phase3d(i,k)-pqlwpcon(i,k) ! fraction eau liquide stratiforme - IF (pqlwpstra(i,k) .LE. 0.0) pqlwpstra(i,k)=0.0 +! 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)=1.1 - & *((pqlwpcon(i,k)*pplay(i,k)/(RD * T(i,k))) - & /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.) - reffclwc(i,k) = MAX(reffclwc(i,k) * 1e6, 5.) - + reffclwc(i,k)=radius + reffclwc(i,k)=reffclwc(i,k)*lcc3dcon(i,k) ! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D - IF ((pclc(i,k)-rnebcon(i,k)) .LE. seuil_neb) THEN ! tout sous la forme convective - reffclws(i,k)=0.0 - lcc3dstra(i,k)= 0.0 - ELSE - reffclws(i,k) = (pclc(i,k)*phase3d(i,k)* - & rad_chaud_tab(i,k)- - & pqlwpcon(i,k)*reffclwc(i,k)) - IF(reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0 - lcc3dstra(i,k)=pqlwpstra(i,k) - ENDIF -!Convertion from um to m - IF(rnebcon(i,k). LE. seuil_neb) THEN - reffclwc(i,k) = reffclwc(i,k)*seuil_neb*clwcon(i,k) - & *1.0E-06 - lcc3dcon(i,k)= seuil_neb*clwcon(i,k) - ELSE - reffclwc(i,k) = reffclwc(i,k)*pqlwpcon(i,k) - & *1.0E-06 - lcc3dcon(i,k) = pqlwpcon(i,k) - ENDIF - - reffclws(i,k) = reffclws(i,k)*1.0E-06 - + reffclws(i,k)=radius + reffclws(i,k)=reffclws(i,k)*lcc3dstra(i,k) ENDDO !klev ENDDO !klon - -!! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D - DO k = 1, klev - DO i = 1, klon - pdel(i,k) = paprs(i,k)-paprs(i,k+1) - zrho(i,k)=pplay(i,k)/t(i,k)/RD ! kg/m3 - dh(i,k)=pdel(i,k)/(gravit*zrho(i,k)) ! hauteur de chaque boite (m) - ENDDO - ENDDO +c +c Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D c DO i = 1, klon @@ -697,21 +584,20 @@ 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 (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 + 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 - - ENDIF !ok_aie - ENDIF !ok newmicro -c -C +c + ENDIF !ok_cdnc +c RETURN +c END Index: LMDZ5/trunk/libf/phylmd/physiq.F =================================================================== --- LMDZ5/trunk/libf/phylmd/physiq.F (revision 1711) +++ LMDZ5/trunk/libf/phylmd/physiq.F (revision 1712) @@ -1094,7 +1094,8 @@ ! Parameters LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not + LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013) REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995) - SAVE ok_ade, ok_aie, bl95_b0, bl95_b1 -c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1) + SAVE ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1 +c$OMP THREADPRIVATE(ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1) LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA ! false : lecture des aerosol dans un fichier @@ -1251,5 +1252,5 @@ . fact_cldcon, facttemps,ok_newmicro,iflag_radia, . iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, - . ok_ade, ok_aie, aerosol_couple, + . ok_ade, ok_aie, ok_cdnc, aerosol_couple, . flag_aerosol, new_aod, . bl95_b0, bl95_b1, @@ -3145,12 +3146,10 @@ if (ok_newmicro) then - CALL newmicro (paprs, pplay,ok_newmicro, - . t_seri, cldliq, cldfra, cldtau, cldemi, - . cldh, cldl, cldm, cldt, cldq, - . flwp, fiwp, flwc, fiwc, - e ok_aie, - e mass_solu_aero, mass_solu_aero_pi, - e bl95_b0, bl95_b1, - s cldtaupi, re, fl, ref_liq, ref_ice) + CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, + . paprs, pplay, t_seri, cldliq, cldfra, + . cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, + e flwp, fiwp, flwc, fiwc, + e mass_solu_aero, mass_solu_aero_pi, + s cldtaupi, re, fl, ref_liq, ref_ice) else CALL nuage (paprs, pplay, @@ -3161,5 +3160,4 @@ e bl95_b0, bl95_b1, s cldtaupi, re, fl) - endif c