Changeset 1712 for LMDZ5/trunk
- Timestamp:
- Jan 18, 2013, 3:07:29 PM (12 years ago)
- Location:
- LMDZ5/trunk/libf/phylmd
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/phylmd/conf_phys_m.F90
r1687 r1712 18 18 iflag_cldcon, & 19 19 iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, & 20 ok_ade, ok_aie, aerosol_couple, &20 ok_ade, ok_aie, ok_cdnc, aerosol_couple, & 21 21 flag_aerosol, new_aod, & 22 22 bl95_b0, bl95_b1,& … … 60 60 ! ok_instan: sorties instantanees 61 61 ! ok_ade, ok_aie: apply or not aerosol direct and indirect effects 62 ! ok_cdnc, ok cloud droplet number concentration 62 63 ! bl95_b*: parameters in the formula to link CDNC to aerosol mass conc 63 64 ! … … 70 71 logical :: ok_LES 71 72 LOGICAL :: callstats 72 LOGICAL :: ok_ade, ok_aie, aerosol_couple73 LOGICAL :: ok_ade, ok_aie, ok_cdnc, aerosol_couple 73 74 INTEGER :: flag_aerosol 74 75 LOGICAL :: new_aod … … 84 85 logical,SAVE :: ok_LES_omp 85 86 LOGICAL,SAVE :: callstats_omp 86 LOGICAL,SAVE :: ok_ade_omp, ok_aie_omp, aerosol_couple_omp87 LOGICAL,SAVE :: ok_ade_omp, ok_aie_omp, ok_cdnc_omp, aerosol_couple_omp 87 88 INTEGER, SAVE :: flag_aerosol_omp 88 89 LOGICAL, SAVE :: new_aod_omp … … 272 273 call getin('ok_aie', ok_aie_omp) 273 274 275 ! 276 !Config Key = ok_cdnc 277 !Config Desc = ok cloud droplet number concentration 278 !Config Def = .false. 279 !Config Help = Used in newmicro.F 280 ! 281 ok_cdnc_omp = .false. 282 call getin('ok_cdnc', ok_cdnc_omp) 274 283 ! 275 284 !Config Key = aerosol_couple … … 1677 1686 ok_ade = ok_ade_omp 1678 1687 ok_aie = ok_aie_omp 1688 ok_cdnc = ok_cdnc_omp 1679 1689 aerosol_couple = aerosol_couple_omp 1680 1690 flag_aerosol=flag_aerosol_omp … … 1774 1784 END IF 1775 1785 END IF 1786 1787 ! ok_cdnc must be set to y if ok_aie is activated 1788 IF (ok_aie .AND. .NOT. ok_cdnc) THEN 1789 CALL abort_gcm('conf_phys', 'ok_cdnc must be set to y if ok_aie is activated',1) 1790 ENDIF 1776 1791 1777 1792 !$OMP MASTER -
LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90
r1687 r1712 99 99 REAL :: dummy 100 100 LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf 101 LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats101 LOGICAL :: ok_LES, ok_ade, ok_aie, ok_cdnc, aerosol_couple, new_aod, callstats 102 102 INTEGER :: iflag_radia, flag_aerosol 103 103 REAL :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut … … 136 136 iflag_cldcon, & 137 137 iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, & 138 ok_ade, ok_aie, aerosol_couple,&138 ok_ade, ok_aie, ok_cdnc, aerosol_couple, & 139 139 flag_aerosol, new_aod, & 140 140 bl95_b0, bl95_b1, & -
LMDZ5/trunk/libf/phylmd/newmicro.F
r1525 r1712 2 2 3 3 4 5 4 ! 6 SUBROUTINE newmicro (paprs, pplay,ok_newmicro, 5 SUBROUTINE newmicro (ok_cdnc, bl95_b0, bl95_b1, 6 . paprs, pplay, 7 7 . t, pqlwp, pclc, pcltau, pclemi, 8 8 . pch, pcl, pcm, pct, pctlwp, 9 s xflwp, xfiwp, xflwc, xfiwc, 10 e ok_aie, 11 e mass_solu_aero, mass_solu_aero_pi, 12 e bl95_b0, bl95_b1, 13 s cldtaupi, re, fl, reliq, reice) 14 9 . xflwp, xfiwp, xflwc, xfiwc, 10 . mass_solu_aero, mass_solu_aero_pi, 11 . pcldtaupi, re, fl, reliq, reice) 12 c 15 13 USE dimphy 16 14 USE phys_local_var_mod, only: scdnc,cldncl,reffclwtop,lcc, … … 21 19 c====================================================================== 22 20 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910 21 c O. Boucher (LMD/CNRS) mise a jour en 201212 23 22 c Objet: Calculer epaisseur optique et emmissivite des nuages 24 23 c====================================================================== 25 24 c Arguments: 25 c ok_cdnc-input-L-flag pour calculer les rayons a partir des aerosols 26 c 26 27 c t-------input-R-temperature 27 c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)28 c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere dans la partie nuageuse (kg/kg) 28 29 c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1) 29 c30 c ok_aie--input-L-apply aerosol indirect effect or not31 30 c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3] 32 c mass_solu_aero_pi--input-R-dito, pre-industrial value 33 c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land) 34 c bl95_b1-input-R-a parameter, may be varied for tests ( -"- ) 31 c mass_solu_aero_pi--input-R-ditto, pre-industrial value 32 c 33 c bl95_b0-input-R-a PARAMETER, may be varied for tests (s-sea, l-land) 34 c bl95_b1-input-R-a PARAMETER, may be varied for tests ( -"- ) 35 35 c 36 c cldtaupi-output-R-pre-industrial value of cloud optical thickness,37 c needed for the diagnostics of the aerosol indirect38 c radiative forcing (see radlwsw)39 36 c re------output-R-Cloud droplet effective radius multiplied by fl [um] 40 37 c fl------output-R-Denominator to re, introduced to avoid problems in 41 38 c the averaging of the output. fl is the fraction of liquid 42 39 c water clouds within a grid cell 40 c 43 41 c pcltau--output-R-epaisseur optique des nuages 44 42 c pclemi--output-R-emissivite des nuages (0 a 1) 43 c pcldtaupi-output-R-pre-industrial value of cloud optical thickness, 44 c 45 c pcl-output-R-2D low-level cloud cover 46 c pcm-output-R-2D mid-level cloud cover 47 c pch-output-R-2D high-level cloud cover 48 c pct-output-R-2D total cloud cover 45 49 c====================================================================== 46 50 C 47 51 #include "YOMCST.h" 48 c49 cym#include "dimensions.h"50 cym#include "dimphy.h"51 52 #include "nuage.h" 52 cIM cf. CR: include pour NOVLP et ZEPSEC53 53 #include "radepsi.h" 54 54 #include "radopt.h" 55 55 56 c choix de l'hypothese de recouvrememnt nuageuse 56 LOGICAL RANDOM,MAXIMUM_RANDOM,MAXIMUM 57 parameter (RANDOM=.FALSE., MAXIMUM_RANDOM=.TRUE., MAXIMUM=.FALSE.) 57 LOGICAL RANDOM, MAXIMUM_RANDOM, MAXIMUM 58 PARAMETER (RANDOM=.FALSE., MAXIMUM_RANDOM=.TRUE., MAXIMUM=.FALSE.) 59 c 58 60 LOGICAL, SAVE :: FIRST=.TRUE. 59 61 !$OMP THREADPRIVATE(FIRST) 60 c Hypoyhese de recouvrement : MAXIMUM_RANDOM61 62 INTEGER flag_max 62 REAL phase3d(klon, klev),dh(klon, klev),pdel(klon, klev), 63 . zrho(klon, klev) 64 REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon) 63 c 64 c threshold PARAMETERs 65 65 REAL thres_tau,thres_neb 66 66 PARAMETER (thres_tau=0.3, thres_neb=0.001) 67 REAL t_tmp 68 REAL gravit69 PARAMETER (gravit=9.80616) !m/s270 REAL pqlwpcon(klon, klev), pqlwpstra(klon, klev) 71 c 72 REAL p aprs(klon,klev+1), pplay(klon,klev)67 c 68 REAL phase3d(klon, klev) 69 REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon) 70 c 71 REAL paprs(klon,klev+1) 72 REAL pplay(klon,klev) 73 73 REAL t(klon,klev) 74 c75 74 REAL pclc(klon,klev) 76 75 REAL pqlwp(klon,klev) 77 REAL pcltau(klon,klev), pclemi(klon,klev) 78 c 79 REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon) 76 REAL pcltau(klon,klev) 77 REAL pclemi(klon,klev) 78 REAL pcldtaupi(klon, klev) 79 c 80 REAL pct(klon) 81 REAL pcl(klon) 82 REAL pcm(klon) 83 REAL pch(klon) 84 REAL pctlwp(klon) 80 85 c 81 86 LOGICAL lo … … 85 90 ! PARAMETER (cetahb = 0.45, cetamb = 0.80) 86 91 ! Remplacer 87 ! cetahb*paprs(i,1) par prmhc88 ! cetamb*paprs(i,1) par prlmc89 REAL prmhc ! Pressure between medium and high level cloud 90 REAL prlmc ! Pressure between low and medium level cloud 92 ! cetahb*paprs(i,1) par prmhc 93 ! cetamb*paprs(i,1) par prlmc 94 REAL prmhc ! Pressure between medium and high level cloud in Pa 95 REAL prlmc ! Pressure between low and medium level cloud in Pa 91 96 PARAMETER (prmhc = 440.*100., prlmc = 680.*100.) 92 93 97 C 94 98 INTEGER i, k 95 cIM: 091003 REAL zflwp, zradef, zfice, zmsac96 REAL zflwp(klon), zradef, zfice, zmsac97 cIM: 091003 rajout98 99 REAL xflwp(klon), xfiwp(klon) 99 100 REAL xflwc(klon,klev), xfiwc(klon,klev) 100 101 c 101 REAL radius, rad_chaud 102 cc PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0) 103 ccc PARAMETER (rad_chaud=15.0, rad_froid=35.0) 104 c sintex initial PARAMETER (rad_chaud=10.0, rad_froid=30.0) 105 REAL coef, coef_froi, coef_chau 102 REAL radius 103 c 104 REAL coef_froi, coef_chau 106 105 PARAMETER (coef_chau=0.13, coef_froi=0.09) 106 c 107 107 REAL seuil_neb 108 108 PARAMETER (seuil_neb=0.001) 109 c 109 110 INTEGER nexpo ! exponentiel pour glace/eau 110 111 PARAMETER (nexpo=6) 111 ccc PARAMETER (nexpo=1) 112 113 c -- sb: 114 logical ok_newmicro 115 c parameter (ok_newmicro=.FALSE.) 116 cIM: 091003 real rel, tc, rei, zfiwp 117 real rel, tc, rei, zfiwp(klon) 118 real k_liq, k_ice0, k_ice, DF 119 parameter (k_liq=0.0903, k_ice0=0.005) ! units=m2/g 120 parameter (DF=1.66) ! diffusivity factor 121 c sb -- 112 c PARAMETER (nexpo=1) 113 114 REAL rel, tc, rei 115 REAL k_ice0, k_ice, DF 116 PARAMETER (k_ice0=0.005) ! units=m2/g 117 PARAMETER (DF=1.66) ! diffusivity factor 118 c 122 119 cjq for the aerosol indirect effect 123 120 cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003 124 121 cjq 125 LOGICAL ok_aie ! Apply AIE or not?126 LOGICAL ok_a1lwpdep ! a1 LWP dependent?127 128 122 REAL mass_solu_aero(klon, klev) ! total mass concentration for all soluble aerosols [ug m-3] 129 123 REAL mass_solu_aero_pi(klon, klev) ! - " - (pre-industrial value) … … 135 129 REAL fl(klon, klev) ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell) 136 130 131 LOGICAL ok_cdnc 137 132 REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula 138 133 139 REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag140 134 cjq-end 141 135 cIM cf. CR:parametres supplementaires … … 145 139 REAL zcloudm(klon) 146 140 REAL zcloudl(klon) 147 148 149 c ************************** 150 c * * 151 c * DEBUT PARTIE OPTIMISEE * 152 c * * 153 c ************************** 154 155 REAL diff_paprs(klon, klev), zfice1, zfice2(klon, klev) 156 REAL rad_chaud_tab(klon, klev), zflwp_var, zfiwp_var 141 REAL rhodz(klon, klev) !--rho*dz pour la couche 142 REAL zrho(klon, klev) !--rho pour la couche 143 REAL dh(klon, klev) !--dz pour la couche 144 REAL zfice(klon, klev) 145 REAL rad_chaud(klon, klev) !--rayon pour les nuages chauds 146 REAL zflwp_var, zfiwp_var 157 147 REAL d_rei_dt 158 148 … … 171 161 ! Pour retrouver les résultats numériques de la version d'origine, 172 162 ! on impose 0.71 quand on est proche de 0.71 173 163 c 174 164 d_rei_dt=(rei_max-rei_min)/81.4 175 165 if (abs(d_rei_dt-0.71)<1.e-4) d_rei_dt=0.71 176 ! print*,'d_rei_dT ',d_rei_dt,rei_min,rei_max177 166 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 178 167 c 179 168 c Calculer l'epaisseur optique et l'emmissivite des nuages 180 c181 169 c IM inversion des DO 170 c 182 171 xflwp = 0.d0 183 172 xfiwp = 0.d0 184 173 xflwc = 0.d0 185 174 xfiwc = 0.d0 186 187 ! Initialisation 175 c 188 176 reliq=0. 189 177 reice=0. 190 178 c 191 179 DO k = 1, klev 192 DO i = 1, klon 193 diff_paprs(i,k) = (paprs(i,k)-paprs(i,k+1))/RG 180 DO i = 1, klon 181 c-layer calculation 182 rhodz(i,k) = (paprs(i,k)-paprs(i,k+1))/RG ! kg/m2 183 zrho(i,k)=pplay(i,k)/t(i,k)/RD ! kg/m3 184 dh(i,k)=rhodz(i,k)/zrho(i,k) ! m 185 c-Fraction of ice in cloud using a linear transition 186 zfice(i,k) = 1.0 - (t(i,k)-t_glace_min) / 187 & (t_glace_max-t_glace_min) 188 zfice(i,k) = MIN(MAX(zfice(i,k),0.0),1.0) 189 c-IM Total Liquid/Ice water content 190 xflwc(i,k) = (1.-zfice(i,k))*pqlwp(i,k) 191 xfiwc(i,k) = zfice(i,k)*pqlwp(i,k) 194 192 ENDDO 195 193 ENDDO 196 194 197 IF (ok_newmicro) THEN 198 199 200 DO k = 1, klev 201 DO i = 1, klon 202 c zfice2(i,k) = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace) 203 zfice2(i,k) = 1.0 - (t(i,k)-t_glace_min) / 204 & (t_glace_max-t_glace_min) 205 zfice2(i,k) = MIN(MAX(zfice2(i,k),0.0),1.0) 206 c IM Total Liquid/Ice water content 207 xflwc(i,k) = (1.-zfice2(i,k))*pqlwp(i,k) 208 xfiwc(i,k) = zfice2(i,k)*pqlwp(i,k) 209 c IM In-Cloud Liquid/Ice water content 210 c xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k) 211 c xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k) 212 ENDDO 213 ENDDO 214 215 IF (ok_aie) THEN 216 DO k = 1, klev 217 DO i = 1, klon 218 ! Formula "D" of Boucher and Lohmann, Tellus, 1995 219 ! 220 cdnc(i,k) = 10.**(bl95_b0+bl95_b1* 195 IF (ok_cdnc) THEN 196 c 197 c--we compute cloud properties as a function of the aerosol load 198 c 199 DO k = 1, klev 200 DO i = 1, klon 201 c 202 c Formula "D" of Boucher and Lohmann, Tellus, 1995 203 c Cloud droplet number concentration (CDNC) is restricted 204 c to be within [20, 1000 cm^3] 205 c 206 c--present-day case 207 cdnc(i,k) = 10.**(bl95_b0+bl95_b1* 221 208 & log(MAX(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3 222 ! Cloud droplet number concentration (CDNC) is restricted 223 ! to be within [20, 1000 cm^3] 224 ! 225 cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k))) 226 ! 227 ! 228 cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1* 209 cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k))) 210 c 211 c--pre-industrial case 212 cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1* 229 213 & log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.)) 230 214 & *1.e6 !-m-3 231 cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k))) 232 ENDDO 233 ENDDO 234 DO k = 1, klev 235 DO i = 1, klon 236 ! rad_chaud_tab(i,k) = 237 ! & MAX(1.1e6 238 ! & *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 239 ! & /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.),5.) 240 rad_chaud_tab(i,k) = 241 & 1.1 242 & *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 243 & /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.) 244 rad_chaud_tab(i,k) = MAX(rad_chaud_tab(i,k) * 1e6, 5.) 245 ENDDO 246 ENDDO 247 ELSE 248 DO k = 1, MIN(3,klev) 249 DO i = 1, klon 250 rad_chaud_tab(i,k) = rad_chau2 251 ENDDO 252 ENDDO 253 DO k = MIN(3,klev)+1, klev 254 DO i = 1, klon 255 rad_chaud_tab(i,k) = rad_chau1 256 ENDDO 257 ENDDO 258 259 ENDIF 260 261 DO k = 1, klev 262 ! IF(.not.ok_aie) THEN 263 rad_chaud = rad_chau1 264 IF (k.LE.3) rad_chaud = rad_chau2 265 ! ENDIF 266 DO i = 1, klon 267 IF (pclc(i,k) .LE. seuil_neb) THEN 268 269 c -- effective cloud droplet radius (microns): 270 271 c for liquid water clouds: 272 ! For output diagnostics 273 ! 274 ! Cloud droplet effective radius [um] 275 ! 276 ! we multiply here with f * xl (fraction of liquid water 277 ! clouds in the grid cell) to avoid problems in the 278 ! averaging of the output. 279 ! In the output of IOIPSL, derive the real cloud droplet 280 ! effective radius as re/fl 281 ! 282 283 fl(i,k) = seuil_neb*(1.-zfice2(i,k)) 284 re(i,k) = rad_chaud_tab(i,k)*fl(i,k) 285 286 rel = 0. 287 rei = 0. 288 pclc(i,k) = 0.0 289 pcltau(i,k) = 0.0 290 pclemi(i,k) = 0.0 291 cldtaupi(i,k) = 0.0 292 ELSE 293 215 cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k))) 216 c 217 c--present-day case 218 rad_chaud(i,k) = 219 & 1.1*((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 220 & /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.) 221 rad_chaud(i,k) = MAX(rad_chaud(i,k) * 1.e6, 5.) 222 c 223 c--pre-industrial case 224 radius = 225 & 1.1*((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 226 & /(4./3.*RPI*1000.*cdnc_pi(i,k)))**(1./3.) 227 radius = MAX(radius*1.e6, 5.) 228 c 229 c--pre-industrial case 230 c--liquid/ice cloud water paths: 231 IF (pclc(i,k) .LE. seuil_neb) THEN 232 c 233 pcldtaupi(i,k) = 0.0 234 c 235 ELSE 236 c 237 zflwp_var= 1000.*(1.-zfice(i,k))*pqlwp(i,k)/pclc(i,k) 238 & *rhodz(i,k) 239 zfiwp_var= 1000.*zfice(i,k)*pqlwp(i,k)/pclc(i,k) 240 & *rhodz(i,k) 241 tc = t(i,k)-273.15 242 rei = d_rei_dt*tc + rei_max 243 if (tc.le.-81.4) rei = rei_min 244 c 245 c-- cloud optical thickness : 246 c [for liquid clouds, traditional formula, 247 c for ice clouds, Ebert & Curry (1992)] 248 c 249 if (zflwp_var.eq.0.) radius = 1. 250 if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1. 251 pcldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius 252 & + zfiwp_var * (3.448e-03 + 2.431/rei) 253 c 254 ENDIF 255 c 256 ENDDO 257 ENDDO 258 c 259 ELSE !--not ok_cdnc 260 c 261 c-prescribed cloud droplet radius 262 c 263 DO k = 1, MIN(3,klev) 264 DO i = 1, klon 265 rad_chaud(i,k) = rad_chau2 266 ENDDO 267 ENDDO 268 DO k = MIN(3,klev)+1, klev 269 DO i = 1, klon 270 rad_chaud(i,k) = rad_chau1 271 ENDDO 272 ENDDO 273 274 ENDIF !--ok_cdnc 275 c 276 c--computation of cloud optical depth and emissivity 277 c--in the general case 278 c 279 DO k = 1, klev 280 DO i = 1, klon 281 c 282 IF (pclc(i,k) .LE. seuil_neb) THEN 283 c 284 c effective cloud droplet radius (microns) for liquid water clouds: 285 c For output diagnostics cloud droplet effective radius [um] 286 c we multiply here with f * xl (fraction of liquid water 287 c clouds in the grid cell) to avoid problems in the averaging of the output. 288 c In the output of IOIPSL, derive the REAL cloud droplet 289 c effective radius as re/fl 290 c 291 fl(i,k) = seuil_neb*(1.-zfice(i,k)) 292 re(i,k) = rad_chaud(i,k)*fl(i,k) 293 pclc(i,k) = 0.0 294 pcltau(i,k) = 0.0 295 pclemi(i,k) = 0.0 296 c 297 ELSE 298 c 294 299 c -- liquid/ice cloud water paths: 295 300 296 zflwp_var= 1000.*(1.-zfice2(i,k))*pqlwp(i,k)/pclc(i,k) 297 & *diff_paprs(i,k) 298 zfiwp_var= 1000.*zfice2(i,k)*pqlwp(i,k)/pclc(i,k) 299 & *diff_paprs(i,k) 300 301 c -- effective cloud droplet radius (microns): 302 303 c for liquid water clouds: 304 305 IF (ok_aie) THEN 306 radius = 307 & 1.1 308 & *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 309 & /(4./3.*RPI*1000.*cdnc_pi(i,k)))**(1./3.) 310 radius = MAX(radius*1e6, 5.) 311 312 tc = t(i,k)-273.15 313 rei = d_rei_dt*tc + rei_max 314 if (tc.le.-81.4) rei = rei_min 315 if (zflwp_var.eq.0.) radius = 1. 316 if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1. 317 cldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius 318 & + zfiwp_var * (3.448e-03 + 2.431/rei) 319 320 ENDIF ! ok_aie 321 ! For output diagnostics 322 ! 323 ! Cloud droplet effective radius [um] 324 ! 325 ! we multiply here with f * xl (fraction of liquid water 326 ! clouds in the grid cell) to avoid problems in the 327 ! averaging of the output. 328 ! In the output of IOIPSL, derive the real cloud droplet 329 ! effective radius as re/fl 330 ! 331 332 fl(i,k) = pclc(i,k)*(1.-zfice2(i,k)) 333 re(i,k) = rad_chaud_tab(i,k)*fl(i,k) 334 335 rel = rad_chaud_tab(i,k) 336 c for ice clouds: as a function of the ambiant temperature 337 c [formula used by Iacobellis and Somerville (2000), with an 338 c asymptotical value of 3.5 microns at T<-81.4 C added to be 339 c consistent with observations of Heymsfield et al. 1986]: 340 c 2011/05/24 : rei_min = 3.5 becomes a free parameter as well as rei_max=61.29 341 tc = t(i,k)-273.15 342 rei = d_rei_dt*tc + rei_max 343 if (tc.le.-81.4) rei = rei_min 344 c -- cloud optical thickness : 345 346 c [for liquid clouds, traditional formula, 347 c for ice clouds, Ebert & Curry (1992)] 348 301 zflwp_var= 1000.*(1.-zfice(i,k))*pqlwp(i,k)/pclc(i,k) 302 & *rhodz(i,k) 303 zfiwp_var= 1000.*zfice(i,k)*pqlwp(i,k)/pclc(i,k) 304 & *rhodz(i,k) 305 c 306 c effective cloud droplet radius (microns) for liquid water clouds: 307 c For output diagnostics cloud droplet effective radius [um] 308 c we multiply here with f * xl (fraction of liquid water 309 c clouds in the grid cell) to avoid problems in the averaging of the output. 310 c In the output of IOIPSL, derive the REAL cloud droplet 311 c effective radius as re/fl 312 c 313 fl(i,k) = pclc(i,k)*(1.-zfice(i,k)) 314 re(i,k) = rad_chaud(i,k)*fl(i,k) 315 c 316 rel = rad_chaud(i,k) 317 c 318 c for ice clouds: as a function of the ambiant temperature 319 c [formula used by Iacobellis and Somerville (2000), with an 320 c asymptotical value of 3.5 microns at T<-81.4 C added to be 321 c consistent with observations of Heymsfield et al. 1986]: 322 c 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as rei_max=61.29 323 c 324 tc = t(i,k)-273.15 325 rei = d_rei_dt*tc + rei_max 326 if (tc.le.-81.4) rei = rei_min 327 c 328 c-- cloud optical thickness : 329 c [for liquid clouds, traditional formula, 330 c for ice clouds, Ebert & Curry (1992)] 331 c 349 332 if (zflwp_var.eq.0.) rel = 1. 350 333 if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1. 351 334 pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel ) 352 335 & + zfiwp_var * (3.448e-03 + 2.431/rei) 336 c 353 337 c -- cloud infrared emissivity: 354 355 c [the broadband infrared absorption coefficient is parameterized 338 c [the broadband infrared absorption coefficient is PARAMETERized 356 339 c as a function of the effective cld droplet radius] 357 358 340 c Ebert and Curry (1992) formula as used by Kiehl & Zender (1995): 341 c 359 342 k_ice = k_ice0 + 1.0/rei 360 343 c 361 344 pclemi(i,k) = 1.0 362 345 & - EXP( -coef_chau*zflwp_var - DF*k_ice*zfiwp_var) 363 364 ENDIF 365 reliq(i,k)=rel 366 reice(i,k)=rei 367 ! if (i.eq.1) then 368 ! print*,'Dans newmicro rel, rei :',rel, rei 369 ! print*,'Dans newmicro reliq, reice :', 370 ! $ reliq(i,k),reice(i,k) 371 ! endif 372 373 ENDDO 374 ENDDO 375 346 c 347 ENDIF 348 c 349 reliq(i,k)=rel 350 reice(i,k)=rei 351 c 352 xflwp(i) = xflwp(i)+ xflwc(i,k) * rhodz(i,k) 353 xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * rhodz(i,k) 354 c 355 ENDDO 356 ENDDO 357 c 358 c--if cloud droplet radius is fixed, then pcldtaupi=pcltau 359 c 360 IF (.NOT.ok_cdnc) THEN 376 361 DO k = 1, klev 377 362 DO i = 1, klon 378 xflwp(i) = xflwp(i)+ xflwc(i,k) * diff_paprs(i,k) 379 xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * diff_paprs(i,k) 380 ENDDO 381 ENDDO 382 383 ELSE 384 DO k = 1, klev 385 rad_chaud = rad_chau1 386 IF (k.LE.3) rad_chaud = rad_chau2 387 DO i = 1, klon 388 389 IF (pclc(i,k) .LE. seuil_neb) THEN 390 391 pclc(i,k) = 0.0 392 pcltau(i,k) = 0.0 393 pclemi(i,k) = 0.0 394 cldtaupi(i,k) = 0.0 395 396 ELSE 397 398 zflwp_var = 1000.*pqlwp(i,k)*diff_paprs(i,k) 399 & /pclc(i,k) 400 401 zfice1 = MIN( 402 & MAX( 1.0 - (t(i,k)-t_glace_min) / 403 & (t_glace_max-t_glace_min),0.0),1.0)**nexpo 404 405 radius = rad_chaud * (1.-zfice1) + rad_froid * zfice1 406 coef = coef_chau * (1.-zfice1) + coef_froi * zfice1 407 408 pcltau(i,k) = 3.0 * zflwp_var / (2.0 * radius) 409 pclemi(i,k) = 1.0 - EXP( - coef * zflwp_var) 410 411 ENDIF 412 413 ENDDO 414 ENDDO 415 ENDIF 416 417 IF (.NOT.ok_aie) THEN 418 DO k = 1, klev 419 DO i = 1, klon 420 cldtaupi(i,k)=pcltau(i,k) 363 pcldtaupi(i,k)=pcltau(i,k) 421 364 ENDDO 422 365 ENDDO 423 366 ENDIF 424 425 ccc DO k = 1, klev426 ccc DO i = 1, klon427 ccc t(i,k) = t(i,k)428 ccc pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )429 ccc lo = pclc(i,k) .GT. (2.*1.e-5)430 ccc zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))431 ccc . /(rg*pclc(i,k))432 ccc zradef = 10.0 + (1.-sigs(k))*45.0433 ccc pcltau(i,k) = 1.5 * zflwp / zradef434 ccc zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)435 ccc zmsac = 0.13*(1.0-zfice) + 0.08*zfice436 ccc pclemi(i,k) = 1.-EXP(-zmsac*zflwp)437 ccc if (.NOT.lo) pclc(i,k) = 0.0438 ccc if (.NOT.lo) pcltau(i,k) = 0.0439 ccc if (.NOT.lo) pclemi(i,k) = 0.0440 ccc ENDDO441 ccc ENDDO442 ccccc print*, 'pas de nuage dans le rayonnement'443 ccccc DO k = 1, klev444 ccccc DO i = 1, klon445 ccccc pclc(i,k) = 0.0446 ccccc pcltau(i,k) = 0.0447 ccccc pclemi(i,k) = 0.0448 ccccc ENDDO449 ccccc ENDDO450 367 C 451 368 C COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS 452 C453 369 c IM cf. CR:test: calcul prenant ou non en compte le recouvrement 454 370 c initialisations 371 c 455 372 DO i=1,klon 456 373 zclear(i)=1. … … 465 382 ENDDO 466 383 C 467 cIM cf CR DO k=1,klev 384 c--calculation of liquid water path 385 c 468 386 DO k = klev, 1, -1 469 387 DO i = 1, klon 470 pctlwp(i) = pctlwp(i) 471 & + pqlwp(i,k)*diff_paprs(i,k) 388 pctlwp(i) = pctlwp(i)+ pqlwp(i,k)*rhodz(i,k) 472 389 ENDDO 473 390 ENDDO 474 c IM cf. CR 391 c 392 c--calculation of cloud properties with cloud overlap 393 c 475 394 IF (NOVLP.EQ.1) THEN 476 395 DO k = klev, 1, -1 477 396 DO i = 1, klon 478 397 zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i))) 479 & /(1.-MIN( real(zcloud(i), kind=8),1.-ZEPSEC))398 & /(1.-MIN(REAL(zcloud(i), kind=8),1.-ZEPSEC)) 480 399 pct(i)=1.-zclear(i) 481 400 IF (paprs(i,k).LT.prmhc) THEN 482 401 pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloudh(i))) 483 & /(1.-MIN( real(zcloudh(i), kind=8),1.-ZEPSEC))402 & /(1.-MIN(REAL(zcloudh(i), kind=8),1.-ZEPSEC)) 484 403 zcloudh(i)=pclc(i,k) 485 404 ELSE IF (paprs(i,k).GE.prmhc .AND. 486 405 & paprs(i,k).LT.prlmc) THEN 487 406 pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloudm(i))) 488 & /(1.-MIN( real(zcloudm(i), kind=8),1.-ZEPSEC))407 & /(1.-MIN(REAL(zcloudm(i), kind=8),1.-ZEPSEC)) 489 408 zcloudm(i)=pclc(i,k) 490 409 ELSE IF (paprs(i,k).GE.prlmc) THEN 491 410 pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloudl(i))) 492 & /(1.-MIN( real(zcloudl(i), kind=8),1.-ZEPSEC))411 & /(1.-MIN(REAL(zcloudl(i), kind=8),1.-ZEPSEC)) 493 412 zcloudl(i)=pclc(i,k) 494 413 endif … … 527 446 ENDDO 528 447 ENDIF 529 530 448 C 531 449 DO i = 1, klon 532 c IM cf. CR pct(i)=1.-pct(i)533 450 pch(i)=1.-pch(i) 534 451 pcm(i)=1.-pcm(i) 535 452 pcl(i)=1.-pcl(i) 536 453 ENDDO 537 454 c 538 455 c ======================================================== 539 !DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL456 c DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL 540 457 c ======================================================== 541 !! change by Nicolas Yan (LSCE) 542 !! Cloud Droplet Number Concentration (CDNC) : 3D variable 543 !! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable 544 !! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable 545 !! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable 546 !! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable 547 IF (ok_newmicro) THEN 548 IF (ok_aie) THEN 458 c change by Nicolas Yan (LSCE) 459 c Cloud Droplet Number Concentration (CDNC) : 3D variable 460 c Fractionnal cover by liquid water cloud (LCC3D) : 3D variable 461 c Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable 462 c Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable 463 c Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable 464 c 465 IF (ok_cdnc) THEN 466 c 549 467 DO k = 1, klev 550 468 DO i = 1, klon 551 phase3d(i,k)=1-zfice 2(i,k)469 phase3d(i,k)=1-zfice(i,k) 552 470 IF (pclc(i,k) .LE. seuil_neb) THEN 553 471 lcc3d(i,k)=seuil_neb*phase3d(i,k) … … 558 476 ENDDO 559 477 ENDDO 560 478 c 561 479 DO i=1,klon 562 480 lcc(i)=0. … … 566 484 IF(MAXIMUM) tcc(i) = 0. 567 485 ENDDO 568 569 486 c 570 487 DO i=1,klon 571 488 DO k=klev-1,1,-1 !From TOA down 572 573 489 c 574 490 ! Test, if the cloud optical depth exceeds the necessary 575 491 ! threshold: 576 492 577 IF (pcltau(i,k).GT.thres_tau .AND. pclc(i,k).GT.thres_neb) 578 . THEN 579 ! To calculate the right Temperature at cloud top, 580 ! interpolate it between layers: 581 t_tmp = t(i,k) + 582 . (paprs(i,k+1)-pplay(i,k))/(pplay(i,k+1)-pplay(i,k)) 583 . * ( t(i,k+1) - t(i,k) ) 584 585 IF(MAXIMUM) THEN 586 IF(FIRST) THEN 493 IF (pcltau(i,k).GT.thres_tau 494 . .AND. pclc(i,k).GT.thres_neb) THEN 495 496 IF (MAXIMUM) THEN 497 IF (FIRST) THEN 587 498 write(*,*)'Hypothese de recouvrement: MAXIMUM' 588 499 FIRST=.FALSE. … … 592 503 ENDIF 593 504 594 IF (RANDOM) THEN595 IF (FIRST) THEN505 IF (RANDOM) THEN 506 IF (FIRST) THEN 596 507 write(*,*)'Hypothese de recouvrement: RANDOM' 597 508 FIRST=.FALSE. … … 601 512 ENDIF 602 513 603 IF (MAXIMUM_RANDOM) THEN604 IF (FIRST) THEN514 IF (MAXIMUM_RANDOM) THEN 515 IF (FIRST) THEN 605 516 write(*,*)'Hypothese de recouvrement: MAXIMUM_ 606 517 . RANDOM' … … 613 524 ENDIF 614 525 c Effective radius of cloud droplet at top of cloud (m) 615 reffclwtop(i) = reffclwtop(i) + rad_chaud _tab(i,k) *526 reffclwtop(i) = reffclwtop(i) + rad_chaud(i,k) * 616 527 . 1.0E-06 * phase3d(i,k) * ( tcc(i) - ftmp(i))*flag_max 617 528 c CDNC at top of cloud (m-3) … … 626 537 ENDIF ! is there a visible, not-too-small cloud? 627 538 ENDDO ! loop over k 628 629 IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i) 539 c 540 IF (RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i) 541 c 630 542 ENDDO ! loop over i 631 543 … … 633 545 DO i = 1, klon 634 546 DO k = 1, klev 635 pqlwpcon(i,k)=rnebcon(i,k)*clwcon(i,k) ! fraction eau liquide convective 636 pqlwpstra(i,k)=pclc(i,k)*phase3d(i,k)-pqlwpcon(i,k) ! fraction eau liquide stratiforme 637 IF (pqlwpstra(i,k) .LE. 0.0) pqlwpstra(i,k)=0.0 547 ! Weight to be used for outputs: eau_liquide*couverture nuageuse 548 lcc3dcon(i,k) =rnebcon(i,k)*phase3d(i,k)*clwcon(i,k) ! eau liquide convective 549 lcc3dstra(i,k)=pclc(i,k)*pqlwp(i,k)*phase3d(i,k) 550 lcc3dstra(i,k)=lcc3dstra(i,k)-lcc3dcon(i,k) ! eau liquide stratiforme 551 lcc3dstra(i,k)=MAX(lcc3dstra(i,k),0.0) 552 ! Compute cloud droplet radius as above in meter 553 radius=1.1*((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k))) 554 & /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.) 555 radius=MAX(radius, 5.e-6) 638 556 ! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D 639 reffclwc(i,k)=1.1 640 & *((pqlwpcon(i,k)*pplay(i,k)/(RD * T(i,k))) 641 & /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.) 642 reffclwc(i,k) = MAX(reffclwc(i,k) * 1e6, 5.) 643 557 reffclwc(i,k)=radius 558 reffclwc(i,k)=reffclwc(i,k)*lcc3dcon(i,k) 644 559 ! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D 645 IF ((pclc(i,k)-rnebcon(i,k)) .LE. seuil_neb) THEN ! tout sous la forme convective 646 reffclws(i,k)=0.0 647 lcc3dstra(i,k)= 0.0 648 ELSE 649 reffclws(i,k) = (pclc(i,k)*phase3d(i,k)* 650 & rad_chaud_tab(i,k)- 651 & pqlwpcon(i,k)*reffclwc(i,k)) 652 IF(reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0 653 lcc3dstra(i,k)=pqlwpstra(i,k) 654 ENDIF 655 !Convertion from um to m 656 IF(rnebcon(i,k). LE. seuil_neb) THEN 657 reffclwc(i,k) = reffclwc(i,k)*seuil_neb*clwcon(i,k) 658 & *1.0E-06 659 lcc3dcon(i,k)= seuil_neb*clwcon(i,k) 660 ELSE 661 reffclwc(i,k) = reffclwc(i,k)*pqlwpcon(i,k) 662 & *1.0E-06 663 lcc3dcon(i,k) = pqlwpcon(i,k) 664 ENDIF 665 666 reffclws(i,k) = reffclws(i,k)*1.0E-06 667 560 reffclws(i,k)=radius 561 reffclws(i,k)=reffclws(i,k)*lcc3dstra(i,k) 668 562 ENDDO !klev 669 563 ENDDO !klon 670 671 !! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D 672 DO k = 1, klev 673 DO i = 1, klon 674 pdel(i,k) = paprs(i,k)-paprs(i,k+1) 675 zrho(i,k)=pplay(i,k)/t(i,k)/RD ! kg/m3 676 dh(i,k)=pdel(i,k)/(gravit*zrho(i,k)) ! hauteur de chaque boite (m) 677 ENDDO 678 ENDDO 564 c 565 c Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D 679 566 c 680 567 DO i = 1, klon … … 697 584 DO i = 1, klon 698 585 DO k = 1, klev 699 IF (scdnc(i,k) .LE. 0.0) scdnc(i,k)=0.0700 IF (reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0701 IF (reffclwc(i,k) .LE. 0.0) reffclwc(i,k)=0.0702 IF (lcc3d(i,k) .LE. 0.0) lcc3d(i,k)=0.0703 IF (lcc3dcon(i,k) .LE. 0.0) lcc3dcon(i,k)=0.0586 IF (scdnc(i,k) .LE. 0.0) scdnc(i,k)=0.0 587 IF (reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0 588 IF (reffclwc(i,k) .LE. 0.0) reffclwc(i,k)=0.0 589 IF (lcc3d(i,k) .LE. 0.0) lcc3d(i,k)=0.0 590 IF (lcc3dcon(i,k) .LE. 0.0) lcc3dcon(i,k)=0.0 704 591 IF (lcc3dstra(i,k) .LE. 0.0) lcc3dstra(i,k)=0.0 705 592 ENDDO 706 IF (reffclwtop(i) .LE. 0.0) reffclwtop(i)=0.0707 IF (cldncl(i) .LE. 0.0) cldncl(i)=0.0708 IF (cldnvi(i) .LE. 0.0) cldnvi(i)=0.0709 IF (lcc(i) .LE. 0.0) lcc(i)=0.0593 IF (reffclwtop(i) .LE. 0.0) reffclwtop(i)=0.0 594 IF (cldncl(i) .LE. 0.0) cldncl(i)=0.0 595 IF (cldnvi(i) .LE. 0.0) cldnvi(i)=0.0 596 IF (lcc(i) .LE. 0.0) lcc(i)=0.0 710 597 ENDDO 711 712 ENDIF !ok_aie 713 ENDIF !ok newmicro 714 c 715 C 598 c 599 ENDIF !ok_cdnc 600 c 716 601 RETURN 602 c 717 603 END -
LMDZ5/trunk/libf/phylmd/physiq.F
r1694 r1712 1094 1094 ! Parameters 1095 1095 LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not 1096 LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013) 1096 1097 REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995) 1097 SAVE ok_ade, ok_aie, bl95_b0, bl95_b11098 c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)1098 SAVE ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1 1099 c$OMP THREADPRIVATE(ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1) 1099 1100 LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA 1100 1101 ! false : lecture des aerosol dans un fichier … … 1251 1252 . fact_cldcon, facttemps,ok_newmicro,iflag_radia, 1252 1253 . iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, 1253 . ok_ade, ok_aie, aerosol_couple,1254 . ok_ade, ok_aie, ok_cdnc, aerosol_couple, 1254 1255 . flag_aerosol, new_aod, 1255 1256 . bl95_b0, bl95_b1, … … 3145 3146 3146 3147 if (ok_newmicro) then 3147 CALL newmicro (paprs, pplay,ok_newmicro, 3148 . t_seri, cldliq, cldfra, cldtau, cldemi, 3149 . cldh, cldl, cldm, cldt, cldq, 3150 . flwp, fiwp, flwc, fiwc, 3151 e ok_aie, 3152 e mass_solu_aero, mass_solu_aero_pi, 3153 e bl95_b0, bl95_b1, 3154 s cldtaupi, re, fl, ref_liq, ref_ice) 3148 CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, 3149 . paprs, pplay, t_seri, cldliq, cldfra, 3150 . cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, 3151 e flwp, fiwp, flwc, fiwc, 3152 e mass_solu_aero, mass_solu_aero_pi, 3153 s cldtaupi, re, fl, ref_liq, ref_ice) 3155 3154 else 3156 3155 CALL nuage (paprs, pplay, … … 3161 3160 e bl95_b0, bl95_b1, 3162 3161 s cldtaupi, re, fl) 3163 3164 3162 endif 3165 3163 c
Note: See TracChangeset
for help on using the changeset viewer.