- Timestamp:
- Jul 29, 2024, 11:01:04 PM (3 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_lscp_tools.F90
r5143 r5144 1 1 MODULE lmdz_lscp_tools 2 2 3 3 IMPLICIT NONE 4 4 5 5 CONTAINS 6 6 7 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++8 SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo)7 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 8 SUBROUTINE FALLICE_VELOCITY(klon, iwc, temp, rho, pres, ptconv, velo) 9 9 10 10 ! Ref: … … 15 15 ! Advances in Modeling Earth Systems, 11, 16 16 ! 3212–3234. https://doi.org/10.1029/2019MS001642 17 17 18 18 USE lmdz_lscp_ini, ONLY: iflag_vice, ffallv_con, ffallv_lsc 19 19 USE lmdz_lscp_ini, ONLY: cice_velo, dice_velo … … 30 30 REAL, INTENT(OUT), DIMENSION(klon) :: velo ! fallspeed velocity of crystals [m/s] 31 31 32 33 32 INTEGER i 34 REAL logvm, iwcg,tempc,phpa,fallv_tun33 REAL logvm, iwcg, tempc, phpa, fallv_tun 35 34 REAL m2ice, m2snow, vmice, vmsnow 36 35 REAL aice, bice, asnow, bsnow 37 38 39 DO i=1,klon 40 41 IF (ptconv(i)) THEN 42 fallv_tun=ffallv_con 43 ELSE 44 fallv_tun=ffallv_lsc 45 ENDIF 46 47 tempc=temp(i)-273.15 ! celcius temp 48 iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0 49 phpa=pres(i)/100. ! pressure in hPa 50 51 IF (iflag_vice == 1) THEN 36 37 DO i = 1, klon 38 39 IF (ptconv(i)) THEN 40 fallv_tun = ffallv_con 41 ELSE 42 fallv_tun = ffallv_lsc 43 ENDIF 44 45 tempc = temp(i) - 273.15 ! celcius temp 46 iwcg = MAX(iwc(i) * 1000., 1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0 47 phpa = pres(i) / 100. ! pressure in hPa 48 49 IF (iflag_vice == 1) THEN 52 50 ! so-called 'empirical parameterization' in Stubenrauch et al. 2019 53 51 IF (tempc >= -60.0) THEN 54 logvm= -0.0000414122*tempc*tempc*log(iwcg)-0.00538922*tempc*log(iwcg) &55 -0.0516344*log(iwcg)+0.00216078*tempc + 1.971456 velo(i)=exp(logvm)52 logvm = -0.0000414122 * tempc * tempc * log(iwcg) - 0.00538922 * tempc * log(iwcg) & 53 - 0.0516344 * log(iwcg) + 0.00216078 * tempc + 1.9714 54 velo(i) = exp(logvm) 57 55 else 58 velo(i)=65.0*(iwcg**0.2)*(150./phpa)**0.1556 velo(i) = 65.0 * (iwcg**0.2) * (150. / phpa)**0.15 59 57 endif 60 61 velo(i) =fallv_tun*velo(i)/100.0 ! from cm/s to m/s62 63 ELSE IF (iflag_vice == 2) THEN58 59 velo(i) = fallv_tun * velo(i) / 100.0 ! from cm/s to m/s 60 61 ELSE IF (iflag_vice == 2) THEN 64 62 ! so called PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019 65 aice=0.587 66 bice=2.45 67 asnow=0.0444 68 bsnow=2.1 69 70 m2ice=((iwcg*0.001/aice)/(exp(13.6-bice*7.76+0.479*bice**2)* & 71 exp((-0.0361+bice*0.0151+0.00149*bice**2)*tempc))) & 72 **(1./(0.807+bice*0.00581+0.0457*bice**2)) 73 74 vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+ & 75 (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc) & 76 *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice) 77 78 79 vmice=vmice*((1000./phpa)**0.2) 80 81 m2snow=((iwcg*0.001/asnow)/(exp(13.6-bsnow*7.76+0.479*bsnow**2)* & 82 exp((-0.0361+bsnow*0.0151+0.00149*bsnow**2)*tempc))) & 83 **(1./(0.807+bsnow*0.00581+0.0457*bsnow**2)) 84 85 86 vmsnow=100.*14.3*exp(13.6-(bsnow+.416)*7.76+0.479*(bsnow+.416)**2)& 87 *exp((-0.0361+(bsnow+.416)*0.0151+0.00149*(bsnow+.416)**2)*tempc)& 88 *(m2snow**(0.807+(bsnow+.416)*0.00581+0.0457*(bsnow+.416)**2))/(iwcg*0.001/asnow) 89 90 vmsnow=vmsnow*((1000./phpa)**0.35) 91 velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s 92 93 ELSE 63 aice = 0.587 64 bice = 2.45 65 asnow = 0.0444 66 bsnow = 2.1 67 68 m2ice = ((iwcg * 0.001 / aice) / (exp(13.6 - bice * 7.76 + 0.479 * bice**2) * & 69 exp((-0.0361 + bice * 0.0151 + 0.00149 * bice**2) * tempc))) & 70 **(1. / (0.807 + bice * 0.00581 + 0.0457 * bice**2)) 71 72 vmice = 100. * 1042.4 * exp(13.6 - (bice + 1) * 7.76 + 0.479 * (bice + 1.)**2) * exp((-0.0361 + & 73 (bice + 1.) * 0.0151 + 0.00149 * (bice + 1.)**2) * tempc) & 74 * (m2ice**(0.807 + (bice + 1.) * 0.00581 + 0.0457 * (bice + 1.)**2)) / (iwcg * 0.001 / aice) 75 76 vmice = vmice * ((1000. / phpa)**0.2) 77 78 m2snow = ((iwcg * 0.001 / asnow) / (exp(13.6 - bsnow * 7.76 + 0.479 * bsnow**2) * & 79 exp((-0.0361 + bsnow * 0.0151 + 0.00149 * bsnow**2) * tempc))) & 80 **(1. / (0.807 + bsnow * 0.00581 + 0.0457 * bsnow**2)) 81 82 vmsnow = 100. * 14.3 * exp(13.6 - (bsnow + .416) * 7.76 + 0.479 * (bsnow + .416)**2)& 83 * exp((-0.0361 + (bsnow + .416) * 0.0151 + 0.00149 * (bsnow + .416)**2) * tempc)& 84 * (m2snow**(0.807 + (bsnow + .416) * 0.00581 + 0.0457 * (bsnow + .416)**2)) / (iwcg * 0.001 / asnow) 85 86 vmsnow = vmsnow * ((1000. / phpa)**0.35) 87 velo(i) = fallv_tun * min(vmsnow, vmice) / 100. ! to m/s 88 89 ELSE 94 90 ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990 95 velo(i) = fallv_tun *cice_velo*((iwcg/1000.)**dice_velo)96 ENDIF91 velo(i) = fallv_tun * cice_velo * ((iwcg / 1000.)**dice_velo) 92 ENDIF 97 93 ENDDO 98 94 99 END SUBROUTINE FALLICE_VELOCITY 100 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 101 102 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 103 SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT) 104 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 105 106 ! Compute the ice fraction 1-xliq (see e.g. 107 ! Doutriaux-Boucher & Quaas 2004, section 2.2.) 108 ! as a function of temperature 109 ! see also Fig 3 of Madeleine et al. 2020, JAMES 110 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 111 95 END SUBROUTINE FALLICE_VELOCITY 96 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 97 98 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 99 SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT) 100 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 101 102 ! Compute the ice fraction 1-xliq (see e.g. 103 ! Doutriaux-Boucher & Quaas 2004, section 2.2.) 104 ! as a function of temperature 105 ! see also Fig 3 of Madeleine et al. 2020, JAMES 106 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 112 107 113 108 USE lmdz_print_control, ONLY: lunout, prt_level … … 118 113 IMPLICIT NONE 119 114 120 121 INTEGER, INTENT(IN) :: klon ! number of horizontal grid points 122 REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature 123 REAL, INTENT(IN), DIMENSION(klon) :: distcltop ! distance to cloud top 124 REAL, INTENT(IN), DIMENSION(klon) :: temp_cltop ! temperature of cloud top 125 INTEGER, INTENT(IN) :: iflag_ice_thermo 126 REAL, INTENT(OUT), DIMENSION(klon) :: icefrac 127 REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT 128 115 INTEGER, INTENT(IN) :: klon ! number of horizontal grid points 116 REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature 117 REAL, INTENT(IN), DIMENSION(klon) :: distcltop ! distance to cloud top 118 REAL, INTENT(IN), DIMENSION(klon) :: temp_cltop ! temperature of cloud top 119 INTEGER, INTENT(IN) :: iflag_ice_thermo 120 REAL, INTENT(OUT), DIMENSION(klon) :: icefrac 121 REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT 129 122 130 123 INTEGER i 131 124 REAL liqfrac_tmp, dicefrac_tmp 132 REAL Dv, denomdep, beta,qsi,dqsidt125 REAL Dv, denomdep, beta, qsi, dqsidt 133 126 LOGICAL ice_thermo 134 127 … … 137 130 138 131 IF ((iflag_t_glace<2)) THEN !.OR. (iflag_t_glace.GT.6)) THEN 139 140 CALL abort_physic(modname,abort_message,1)132 abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6' 133 CALL abort_physic(modname, abort_message, 1) 141 134 ENDIF 142 135 143 136 IF (.NOT.((iflag_ice_thermo == 1).OR.(iflag_ice_thermo >= 3))) THEN 144 145 CALL abort_physic(modname,abort_message,1)137 abort_message = 'lscp cannot be used without ice thermodynamics' 138 CALL abort_physic(modname, abort_message, 1) 146 139 ENDIF 147 140 148 149 DO i=1,klon 150 151 ! old function with sole dependence upon temperature 152 IF (iflag_t_glace == 2) THEN 153 liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) 154 liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) 155 icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace 156 IF (icefrac(i) >0.) THEN 157 dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) & 158 / (t_glace_min - t_glace_max) 141 DO i = 1, klon 142 143 ! old function with sole dependence upon temperature 144 IF (iflag_t_glace == 2) THEN 145 liqfrac_tmp = (temp(i) - t_glace_min) / (t_glace_max - t_glace_min) 146 liqfrac_tmp = MIN(MAX(liqfrac_tmp, 0.0), 1.0) 147 icefrac(i) = (1.0 - liqfrac_tmp)**exposant_glace 148 IF (icefrac(i) >0.) THEN 149 dicefracdT(i) = exposant_glace * (icefrac(i)**(exposant_glace - 1.)) & 150 / (t_glace_min - t_glace_max) 151 ENDIF 152 153 IF ((icefrac(i)==0).OR.(icefrac(i)==1)) THEN 154 dicefracdT(i) = 0. 155 ENDIF 156 157 ENDIF 158 159 ! function of temperature used in CMIP6 physics 160 IF (iflag_t_glace == 3) THEN 161 liqfrac_tmp = (temp(i) - t_glace_min) / (t_glace_max - t_glace_min) 162 liqfrac_tmp = MIN(MAX(liqfrac_tmp, 0.0), 1.0) 163 icefrac(i) = 1.0 - liqfrac_tmp**exposant_glace 164 IF ((icefrac(i) >0.) .AND. (liqfrac_tmp > 0.)) THEN 165 dicefracdT(i) = exposant_glace * ((liqfrac_tmp)**(exposant_glace - 1.)) & 166 / (t_glace_min - t_glace_max) 167 ELSE 168 dicefracdT(i) = 0. 169 ENDIF 170 ENDIF 171 172 ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top 173 ! and then decreases with decreasing height 174 175 !with linear function of temperature at cloud top 176 IF (iflag_t_glace == 4) THEN 177 liqfrac_tmp = (temp(i) - t_glace_min) / (t_glace_max - t_glace_min) 178 liqfrac_tmp = MIN(MAX(liqfrac_tmp, 0.0), 1.0) 179 icefrac(i) = MAX(MIN(1., 1.0 - liqfrac_tmp * exp(-distcltop(i) / dist_liq)), 0.) 180 dicefrac_tmp = - temp(i) / (t_glace_max - t_glace_min) 181 dicefracdT(i) = dicefrac_tmp * exp(-distcltop(i) / dist_liq) 182 IF ((liqfrac_tmp <=0) .OR. (liqfrac_tmp >= 1)) THEN 183 dicefracdT(i) = 0. 184 ENDIF 185 ENDIF 186 187 ! with CMIP6 function of temperature at cloud top 188 IF ((iflag_t_glace == 5) .OR. (iflag_t_glace == 7)) THEN 189 liqfrac_tmp = (temp(i) - t_glace_min) / (t_glace_max - t_glace_min) 190 liqfrac_tmp = MIN(MAX(liqfrac_tmp, 0.0), 1.0) 191 liqfrac_tmp = liqfrac_tmp**exposant_glace 192 icefrac(i) = MAX(MIN(1., 1.0 - liqfrac_tmp * exp(-distcltop(i) / dist_liq)), 0.) 193 IF ((liqfrac_tmp <=0) .OR. (liqfrac_tmp >= 1)) THEN 194 dicefracdT(i) = 0. 195 ELSE 196 dicefracdT(i) = exposant_glace * ((liqfrac_tmp)**(exposant_glace - 1.)) / (t_glace_min - t_glace_max) & 197 * exp(-distcltop(i) / dist_liq) 198 ENDIF 199 ENDIF 200 201 ! with modified function of temperature at cloud top 202 ! to get largere values around 260 K, works well with t_glace_min = 241K 203 IF (iflag_t_glace == 6) THEN 204 IF (temp(i) > t_glace_max) THEN 205 liqfrac_tmp = 1. 206 ELSE 207 liqfrac_tmp = -((temp(i) - t_glace_max) / (t_glace_max - t_glace_min))**2 + 1. 208 ENDIF 209 liqfrac_tmp = MIN(MAX(liqfrac_tmp, 0.0), 1.0) 210 icefrac(i) = MAX(MIN(1., 1.0 - liqfrac_tmp * exp(-distcltop(i) / dist_liq)), 0.) 211 IF ((liqfrac_tmp <=0) .OR. (liqfrac_tmp >= 1)) THEN 212 dicefracdT(i) = 0. 213 ELSE 214 dicefracdT(i) = 2 * ((temp(i) - t_glace_max) / (t_glace_max - t_glace_min)) / (t_glace_max - t_glace_min) & 215 * exp(-distcltop(i) / dist_liq) 216 ENDIF 217 ENDIF 218 219 ! if temperature of cloud top <-40°C, 220 IF (iflag_t_glace >= 4) THEN 221 IF ((temp_cltop(i) <= temp_nowater) .AND. (temp(i) <= t_glace_max)) THEN 222 icefrac(i) = 1. 223 dicefracdT(i) = 0. 224 ENDIF 225 ENDIF 226 227 ENDDO ! klon 228 229 END SUBROUTINE ICEFRAC_LSCP 230 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 231 232 SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, temp, pplay, paprsdn, paprsup, qice_ini, snowcld, qtot_incl, cldfra, tke, tke_dissip, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb) 233 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 234 ! Compute the liquid, ice and vapour content (+ice fraction) based 235 ! on turbulence (see Fields 2014, Furtado 2016) 236 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 237 238 USE lmdz_lscp_ini, ONLY: prt_level, lunout 239 USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI 240 USE lmdz_lscp_ini, ONLY: seuil_neb, temp_nowater 241 USE lmdz_lscp_ini, ONLY: tau_mixenv, lmix_mpc, naero5, gamma_snwretro, gamma_taud, capa_crystal 242 USE lmdz_lscp_ini, ONLY: eps 243 244 IMPLICIT NONE 245 246 INTEGER, INTENT(IN) :: klon !--number of horizontal grid points 247 REAL, INTENT(IN) :: dtime !--time step [s] 248 249 REAL, INTENT(IN), DIMENSION(klon) :: temp !--temperature 250 REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa] 251 REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa] 252 REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa] 253 REAL, INTENT(IN), DIMENSION(klon) :: qtot_incl !--specific total cloud water content, in-cloud content [kg/kg] 254 REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction in gridbox [-] 255 REAL, INTENT(IN), DIMENSION(klon) :: tke !--turbulent kinetic energy [m2/s2] 256 REAL, INTENT(IN), DIMENSION(klon) :: tke_dissip !--TKE dissipation [m2/s3] 257 258 REAL, INTENT(IN), DIMENSION(klon) :: qice_ini !--initial specific ice content gridbox-mean [kg/kg] 259 REAL, INTENT(IN), DIMENSION(klon) :: snowcld 260 REAL, INTENT(OUT), DIMENSION(klon) :: qliq !--specific liquid content gridbox-mean [kg/kg] 261 REAL, INTENT(OUT), DIMENSION(klon) :: qvap_cld !--specific cloud vapor content, gridbox-mean [kg/kg] 262 REAL, INTENT(OUT), DIMENSION(klon) :: qice !--specific ice content gridbox-mean [kg/kg] 263 REAL, INTENT(OUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-] 264 REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT 265 266 REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra which is liquid only 267 REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Temporary 268 REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Temporary 269 270 REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati !--specific humidity saturation values 271 INTEGER :: i 272 273 REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl !--In-cloud specific quantities [kg/kg] 274 REAL :: qsnowcld_incl 275 !REAL :: capa_crystal !--Capacitance of ice crystals [-] 276 REAL :: water_vapor_diff !--Water-vapour diffusion coefficient in air [m2/s] (function of T&P) 277 REAL :: air_thermal_conduct !--Thermal conductivity of air [J/m/K/s] (function of T) 278 REAL :: C0 !--Lagrangian structure function [-] 279 REAL :: tau_mixingenv 280 REAL :: tau_dissipturb 281 REAL :: invtau_phaserelax 282 REAL :: sigma2_pdf, mean_pdf 283 REAL :: ai, bi, B0 284 REAL :: sursat_iceliq 285 REAL :: sursat_env 286 REAL :: liqfra_max 287 REAL :: sursat_iceext 288 REAL :: nb_crystals !--number concentration of ice crystals [#/m3] 289 REAL :: moment1_PSD !--1st moment of ice PSD 290 REAL :: N0_PSD, lambda_PSD !--parameters of the exponential PSD 291 292 REAL :: rho_ice !--ice density [kg/m3] 293 REAL :: cldfra1D 294 REAL :: deltaz, rho_air 295 REAL :: psati !--saturation vapor pressure wrt i [Pa] 296 297 C0 = 10. !--value assumed in Field2014 298 rho_ice = 950. 299 sursat_iceext = -0.1 300 !capa_crystal = 1. !r_ice 301 qzero(:) = 0. 302 cldfraliq(:) = 0. 303 icefrac(:) = 0. 304 dicefracdT(:) = 0. 305 306 sigma2_icefracturb(:) = 0. 307 mean_icefracturb(:) = 0. 308 309 !--wrt liquid water 310 CALL calc_qsat_ecmwf(klon, temp(:), qzero(:), pplay(:), RTT, 1, .FALSE., qsatl(:), dqsatl(:)) 311 !--wrt ice 312 CALL calc_qsat_ecmwf(klon, temp(:), qzero(:), pplay(:), RTT, 2, .FALSE., qsati(:), dqsati(:)) 313 314 DO i = 1, klon 315 316 rho_air = pplay(i) / temp(i) / RD 317 !deltaz = ( paprsdn(i) - paprsup(i) ) / RG / rho_air(i) 318 ! because cldfra is intent in, but can be locally modified due to test 319 cldfra1D = cldfra(i) 320 IF (cldfra(i) <= 0.) THEN 321 qvap_cld(i) = 0. 322 qliq(i) = 0. 323 qice(i) = 0. 324 cldfraliq(i) = 0. 325 icefrac(i) = 0. 326 dicefracdT(i) = 0. 327 328 ! If there is a cloud 329 ELSE 330 IF (cldfra(i) >= 1.0) THEN 331 cldfra1D = 1.0 332 END IF 333 334 ! T>0°C, no ice allowed 335 IF (temp(i) >= RTT) THEN 336 qvap_cld(i) = qsatl(i) * cldfra1D 337 qliq(i) = MAX(0.0, qtot_incl(i) - qsatl(i)) * cldfra1D 338 qice(i) = 0. 339 cldfraliq(i) = 1. 340 icefrac(i) = 0. 341 dicefracdT(i) = 0. 342 343 ! T<-38°C, no liquid allowed 344 ELSE IF (temp(i) <= temp_nowater) THEN 345 qvap_cld(i) = qsati(i) * cldfra1D 346 qliq(i) = 0. 347 qice(i) = MAX(0.0, qtot_incl(i) - qsati(i)) * cldfra1D 348 cldfraliq(i) = 0. 349 icefrac(i) = 1. 350 dicefracdT(i) = 0. 351 352 ! MPC temperature 353 ELSE 354 ! Not enough TKE 355 IF (tke_dissip(i) <= eps) THEN 356 qvap_cld(i) = qsati(i) * cldfra1D 357 qliq(i) = 0. 358 qice(i) = MAX(0., qtot_incl(i) - qsati(i)) * cldfra1D 359 cldfraliq(i) = 0. 360 icefrac(i) = 1. 361 dicefracdT(i) = 0. 362 363 ! Enough TKE 364 ELSE 365 !--------------------------------------------------------- 366 !-- ICE SUPERSATURATION PDF 367 !--------------------------------------------------------- 368 !--If -38°C< T <0°C and there is enough turbulence, 369 !--we compute the cloud liquid properties with a Gaussian PDF 370 !--of ice supersaturation F(Si) (Field2014, Furtado2016). 371 !--Parameters of the PDF are function of turbulence and 372 !--microphysics/existing ice. 373 374 sursat_iceliq = qsatl(i) / qsati(i) - 1. 375 psati = qsati(i) * pplay(i) / (RD / RV) 376 377 !-------------- MICROPHYSICAL TERMS -------------- 378 !--We assume an exponential ice PSD whose parameters 379 !--are computed following Morrison&Gettelman 2008 380 !--Ice number density is assumed equals to INP density 381 !--which is a function of temperature (DeMott 2010) 382 !--bi and B0 are microphysical function characterizing 383 !--vapor/ice interactions 384 !--tau_phase_relax is the typical time of vapor deposition 385 !--onto ice crystals 386 387 qiceini_incl = qice_ini(i) / cldfra1D 388 qsnowcld_incl = snowcld(i) * RG * dtime / (paprsdn(i) - paprsup(i)) / cldfra1D 389 sursat_env = max(0., (qtot_incl(i) - qiceini_incl) / qsati(i) - 1.) 390 IF (qiceini_incl > eps) THEN 391 nb_crystals = 1.e3 * 5.94e-5 * (RTT - temp(i))**3.33 * naero5**(0.0264 * (RTT - temp(i)) + 0.0033) 392 lambda_PSD = ((RPI * rho_ice * nb_crystals * 24.) / (6. * (qiceini_incl + gamma_snwretro * qsnowcld_incl))) ** (1. / 3.) 393 N0_PSD = nb_crystals * lambda_PSD 394 moment1_PSD = N0_PSD / 2. / lambda_PSD**2 395 ELSE 396 moment1_PSD = 0. 159 397 ENDIF 160 398 161 IF ((icefrac(i)==0).OR.(icefrac(i)==1)) THEN 162 dicefracdT(i)=0. 399 !--Formulae for air thermal conductivity and water vapor diffusivity 400 !--comes respectively from Beard and Pruppacher (1971) 401 !--and Hall and Pruppacher (1976) 402 403 air_thermal_conduct = (5.69 + 0.017 * (temp(i) - RTT)) * 1.e-3 * 4.184 404 water_vapor_diff = 2.11 * 1e-5 * (temp(i) / RTT)**1.94 * (101325 / pplay(i)) 405 406 bi = 1. / ((qsati(i) + qsatl(i)) / 2.) + RLSTT**2 / RCPD / RV / temp(i)**2 407 B0 = 4. * RPI * capa_crystal * 1. / (RLSTT**2 / air_thermal_conduct / RV / temp(i)**2 & 408 + RV * temp(i) / psati / water_vapor_diff) 409 410 invtau_phaserelax = (bi * B0 * moment1_PSD) 411 412 ! Old way of estimating moment1 : spherical crystals + monodisperse PSD 413 ! nb_crystals = rho_air * qiceini_incl / ( 4. / 3. * RPI * r_ice**3. * rho_ice ) 414 ! moment1_PSD = nb_crystals * r_ice 415 416 !----------------- TURBULENT SOURCE/SINK TERMS ----------------- 417 !--Tau_mixingenv is the time needed to homogeneize the parcel 418 !--with its environment by turbulent diffusion over the parcel 419 !--length scale 420 !--if lmix_mpc <0, tau_mixigenv value is prescribed 421 !--else tau_mixigenv value is derived from tke_dissip and lmix_mpc 422 !--Tau_dissipturb is the time needed turbulence to decay due to 423 !--viscosity 424 425 ai = RG / RD / temp(i) * (RD * RLSTT / RCPD / RV / temp(i) - 1.) 426 IF (lmix_mpc > 0) THEN 427 tau_mixingenv = (lmix_mpc**2. / tke_dissip(i))**(1. / 3.) 428 ELSE 429 tau_mixingenv = tau_mixenv 163 430 ENDIF 164 431 165 ENDIF 166 167 ! function of temperature used in CMIP6 physics 168 IF (iflag_t_glace == 3) THEN 169 liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) 170 liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) 171 icefrac(i) = 1.0-liqfrac_tmp**exposant_glace 172 IF ((icefrac(i) >0.) .AND. (liqfrac_tmp > 0.)) THEN 173 dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) & 174 / (t_glace_min - t_glace_max) 432 tau_dissipturb = gamma_taud * 2. * 2. / 3. * tke(i) / tke_dissip(i) / C0 433 434 !--------------------- PDF COMPUTATIONS --------------------- 435 !--Formulae for sigma2_pdf (variance), mean of PDF in Furtado2016 436 !--cloud liquid fraction and in-cloud liquid content are given 437 !--by integrating resp. F(Si) and Si*F(Si) 438 !--Liquid is limited by the available water vapor trough a 439 !--maximal liquid fraction 440 441 liqfra_max = MAX(0., (MIN (1., (qtot_incl(i) - qiceini_incl - qsati(i) * (1 + sursat_iceext)) / (qsatl(i) - qsati(i))))) 442 sigma2_pdf = 1. / 2. * (ai**2) * 2. / 3. * tke(i) * tau_dissipturb / (invtau_phaserelax + 1. / tau_mixingenv) 443 mean_pdf = sursat_env * 1. / tau_mixingenv / (invtau_phaserelax + 1. / tau_mixingenv) 444 cldfraliq(i) = 0.5 * (1. - erf((sursat_iceliq - mean_pdf) / (SQRT(2. * sigma2_pdf)))) 445 !IF (cldfraliq(i) .GT. liqfra_max) THEN 446 ! cldfraliq(i) = liqfra_max 447 !ENDIF 448 449 qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2. * RPI) * EXP(-1. * (sursat_iceliq - mean_pdf)**2. / (2. * sigma2_pdf)) & 450 - qsati(i) * cldfraliq(i) * (sursat_iceliq - mean_pdf) 451 452 sigma2_icefracturb(i) = sigma2_pdf 453 mean_icefracturb(i) = mean_pdf 454 !------------ ICE AMOUNT AND WATER CONSERVATION ------------ 455 456 IF ((qliq_incl <= eps) .OR. (cldfraliq(i) <= eps)) THEN 457 qliq_incl = 0. 458 cldfraliq(i) = 0. 459 END IF 460 461 !--Choice for in-cloud vapor : 462 !--1.Weighted mean between qvap in MPC parts and in ice-only parts 463 !--2.Always at ice saturation 464 qvap_incl = MAX(qsati(i), (1. - cldfraliq(i)) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i)) 465 466 IF (qvap_incl >= qtot_incl(i)) THEN 467 qvap_incl = qsati(i) 468 qliq_incl = qtot_incl(i) - qvap_incl 469 qice_incl = 0. 470 471 ELSEIF ((qvap_incl + qliq_incl) >= qtot_incl(i)) THEN 472 qliq_incl = MAX(0.0, qtot_incl(i) - qvap_incl) 473 qice_incl = 0. 175 474 ELSE 176 dicefracdT(i)=0. 177 ENDIF 178 ENDIF 179 180 ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top 181 ! and then decreases with decreasing height 182 183 !with linear function of temperature at cloud top 184 IF (iflag_t_glace == 4) THEN 185 liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) 186 liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) 187 icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) 188 dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min) 189 dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq) 190 IF ((liqfrac_tmp <=0) .OR. (liqfrac_tmp >= 1)) THEN 191 dicefracdT(i) = 0. 192 ENDIF 193 ENDIF 194 195 ! with CMIP6 function of temperature at cloud top 196 IF ((iflag_t_glace == 5) .OR. (iflag_t_glace == 7)) THEN 197 liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) 198 liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) 199 liqfrac_tmp = liqfrac_tmp**exposant_glace 200 icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) 201 IF ((liqfrac_tmp <=0) .OR. (liqfrac_tmp >= 1)) THEN 202 dicefracdT(i) = 0. 203 ELSE 204 dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace-1.))/(t_glace_min- t_glace_max) & 205 *exp(-distcltop(i)/dist_liq) 206 ENDIF 207 ENDIF 208 209 ! with modified function of temperature at cloud top 210 ! to get largere values around 260 K, works well with t_glace_min = 241K 211 IF (iflag_t_glace == 6) THEN 212 IF (temp(i) > t_glace_max) THEN 213 liqfrac_tmp = 1. 214 ELSE 215 liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1. 216 ENDIF 217 liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) 218 icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) 219 IF ((liqfrac_tmp <=0) .OR. (liqfrac_tmp >= 1)) THEN 220 dicefracdT(i) = 0. 221 ELSE 222 dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) & 223 *exp(-distcltop(i)/dist_liq) 224 ENDIF 225 ENDIF 226 227 ! if temperature of cloud top <-40°C, 228 IF (iflag_t_glace >= 4) THEN 229 IF ((temp_cltop(i) <= temp_nowater) .AND. (temp(i) <= t_glace_max)) THEN 230 icefrac(i) = 1. 231 dicefracdT(i) = 0. 232 ENDIF 233 ENDIF 234 235 236 ENDDO ! klon 237 238 239 END SUBROUTINE ICEFRAC_LSCP 240 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 241 242 SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, temp, pplay, paprsdn, paprsup, qice_ini, snowcld, qtot_incl, cldfra, tke, tke_dissip, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb) 243 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 244 ! Compute the liquid, ice and vapour content (+ice fraction) based 245 ! on turbulence (see Fields 2014, Furtado 2016) 246 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 247 248 249 USE lmdz_lscp_ini, ONLY: prt_level, lunout 250 USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI 251 USE lmdz_lscp_ini, ONLY: seuil_neb, temp_nowater 252 USE lmdz_lscp_ini, ONLY: tau_mixenv, lmix_mpc, naero5, gamma_snwretro, gamma_taud, capa_crystal 253 USE lmdz_lscp_ini, ONLY: eps 254 255 IMPLICIT NONE 256 257 INTEGER, INTENT(IN) :: klon !--number of horizontal grid points 258 REAL, INTENT(IN) :: dtime !--time step [s] 259 260 REAL, INTENT(IN), DIMENSION(klon) :: temp !--temperature 261 REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa] 262 REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa] 263 REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa] 264 REAL, INTENT(IN), DIMENSION(klon) :: qtot_incl !--specific total cloud water content, in-cloud content [kg/kg] 265 REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction in gridbox [-] 266 REAL, INTENT(IN), DIMENSION(klon) :: tke !--turbulent kinetic energy [m2/s2] 267 REAL, INTENT(IN), DIMENSION(klon) :: tke_dissip !--TKE dissipation [m2/s3] 268 269 REAL, INTENT(IN), DIMENSION(klon) :: qice_ini !--initial specific ice content gridbox-mean [kg/kg] 270 REAL, INTENT(IN), DIMENSION(klon) :: snowcld 271 REAL, INTENT(OUT), DIMENSION(klon) :: qliq !--specific liquid content gridbox-mean [kg/kg] 272 REAL, INTENT(OUT), DIMENSION(klon) :: qvap_cld !--specific cloud vapor content, gridbox-mean [kg/kg] 273 REAL, INTENT(OUT), DIMENSION(klon) :: qice !--specific ice content gridbox-mean [kg/kg] 274 REAL, INTENT(OUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-] 275 REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT 276 277 REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra which is liquid only 278 REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Temporary 279 REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Temporary 280 281 REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati !--specific humidity saturation values 282 INTEGER :: i 283 284 REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl !--In-cloud specific quantities [kg/kg] 285 REAL :: qsnowcld_incl 286 !REAL :: capa_crystal !--Capacitance of ice crystals [-] 287 REAL :: water_vapor_diff !--Water-vapour diffusion coefficient in air [m2/s] (function of T&P) 288 REAL :: air_thermal_conduct !--Thermal conductivity of air [J/m/K/s] (function of T) 289 REAL :: C0 !--Lagrangian structure function [-] 290 REAL :: tau_mixingenv 291 REAL :: tau_dissipturb 292 REAL :: invtau_phaserelax 293 REAL :: sigma2_pdf, mean_pdf 294 REAL :: ai, bi, B0 295 REAL :: sursat_iceliq 296 REAL :: sursat_env 297 REAL :: liqfra_max 298 REAL :: sursat_iceext 299 REAL :: nb_crystals !--number concentration of ice crystals [#/m3] 300 REAL :: moment1_PSD !--1st moment of ice PSD 301 REAL :: N0_PSD, lambda_PSD !--parameters of the exponential PSD 302 303 REAL :: rho_ice !--ice density [kg/m3] 304 REAL :: cldfra1D 305 REAL :: deltaz, rho_air 306 REAL :: psati !--saturation vapor pressure wrt i [Pa] 307 308 C0 = 10. !--value assumed in Field2014 309 rho_ice = 950. 310 sursat_iceext = -0.1 311 !capa_crystal = 1. !r_ice 312 qzero(:) = 0. 313 cldfraliq(:) = 0. 314 icefrac(:) = 0. 315 dicefracdT(:) = 0. 316 317 sigma2_icefracturb(:) = 0. 318 mean_icefracturb(:) = 0. 319 320 !--wrt liquid water 321 CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.FALSE.,qsatl(:),dqsatl(:)) 322 !--wrt ice 323 CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.FALSE.,qsati(:),dqsati(:)) 324 325 326 DO i=1,klon 327 328 329 rho_air = pplay(i) / temp(i) / RD 330 !deltaz = ( paprsdn(i) - paprsup(i) ) / RG / rho_air(i) 331 ! because cldfra is intent in, but can be locally modified due to test 332 cldfra1D = cldfra(i) 333 IF (cldfra(i) <= 0.) THEN 334 qvap_cld(i) = 0. 335 qliq(i) = 0. 336 qice(i) = 0. 337 cldfraliq(i) = 0. 338 icefrac(i) = 0. 339 dicefracdT(i) = 0. 340 341 ! If there is a cloud 342 ELSE 343 IF (cldfra(i) >= 1.0) THEN 344 cldfra1D = 1.0 345 END IF 346 347 ! T>0°C, no ice allowed 348 IF ( temp(i) >= RTT ) THEN 349 qvap_cld(i) = qsatl(i) * cldfra1D 350 qliq(i) = MAX(0.0,qtot_incl(i)-qsatl(i)) * cldfra1D 351 qice(i) = 0. 352 cldfraliq(i) = 1. 353 icefrac(i) = 0. 354 dicefracdT(i) = 0. 355 356 ! T<-38°C, no liquid allowed 357 ELSE IF ( temp(i) <= temp_nowater) THEN 358 qvap_cld(i) = qsati(i) * cldfra1D 359 qliq(i) = 0. 360 qice(i) = MAX(0.0,qtot_incl(i)-qsati(i)) * cldfra1D 361 cldfraliq(i) = 0. 362 icefrac(i) = 1. 363 dicefracdT(i) = 0. 364 365 ! MPC temperature 366 ELSE 367 ! Not enough TKE 368 IF ( tke_dissip(i) <= eps ) THEN 369 qvap_cld(i) = qsati(i) * cldfra1D 370 qliq(i) = 0. 371 qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D 372 cldfraliq(i) = 0. 373 icefrac(i) = 1. 374 dicefracdT(i) = 0. 375 376 ! Enough TKE 377 ELSE 378 !--------------------------------------------------------- 379 !-- ICE SUPERSATURATION PDF 380 !--------------------------------------------------------- 381 !--If -38°C< T <0°C and there is enough turbulence, 382 !--we compute the cloud liquid properties with a Gaussian PDF 383 !--of ice supersaturation F(Si) (Field2014, Furtado2016). 384 !--Parameters of the PDF are function of turbulence and 385 !--microphysics/existing ice. 386 387 sursat_iceliq = qsatl(i)/qsati(i) - 1. 388 psati = qsati(i) * pplay(i) / (RD/RV) 389 390 !-------------- MICROPHYSICAL TERMS -------------- 391 !--We assume an exponential ice PSD whose parameters 392 !--are computed following Morrison&Gettelman 2008 393 !--Ice number density is assumed equals to INP density 394 !--which is a function of temperature (DeMott 2010) 395 !--bi and B0 are microphysical function characterizing 396 !--vapor/ice interactions 397 !--tau_phase_relax is the typical time of vapor deposition 398 !--onto ice crystals 399 400 qiceini_incl = qice_ini(i) / cldfra1D 401 qsnowcld_incl = snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) ) / cldfra1D 402 sursat_env = max(0., (qtot_incl(i) - qiceini_incl)/qsati(i) - 1.) 403 IF ( qiceini_incl > eps ) THEN 404 nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )**3.33 * naero5**(0.0264*(RTT-temp(i))+0.0033) 405 lambda_PSD = ( (RPI*rho_ice*nb_crystals*24.) / (6.*(qiceini_incl + gamma_snwretro * qsnowcld_incl)) ) ** (1./3.) 406 N0_PSD = nb_crystals * lambda_PSD 407 moment1_PSD = N0_PSD/2./lambda_PSD**2 408 ELSE 409 moment1_PSD = 0. 410 ENDIF 411 412 !--Formulae for air thermal conductivity and water vapor diffusivity 413 !--comes respectively from Beard and Pruppacher (1971) 414 !--and Hall and Pruppacher (1976) 415 416 air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184 417 water_vapor_diff = 2.11*1e-5 * ( temp(i) / RTT )**1.94 * ( 101325 / pplay(i) ) 418 419 bi = 1./((qsati(i)+qsatl(i))/2.) + RLSTT**2 / RCPD / RV / temp(i)**2 420 B0 = 4. * RPI * capa_crystal * 1. / ( RLSTT**2 / air_thermal_conduct / RV / temp(i)**2 & 421 + RV * temp(i) / psati / water_vapor_diff ) 422 423 invtau_phaserelax = (bi * B0 * moment1_PSD ) 424 425 ! Old way of estimating moment1 : spherical crystals + monodisperse PSD 426 ! nb_crystals = rho_air * qiceini_incl / ( 4. / 3. * RPI * r_ice**3. * rho_ice ) 427 ! moment1_PSD = nb_crystals * r_ice 428 429 !----------------- TURBULENT SOURCE/SINK TERMS ----------------- 430 !--Tau_mixingenv is the time needed to homogeneize the parcel 431 !--with its environment by turbulent diffusion over the parcel 432 !--length scale 433 !--if lmix_mpc <0, tau_mixigenv value is prescribed 434 !--else tau_mixigenv value is derived from tke_dissip and lmix_mpc 435 !--Tau_dissipturb is the time needed turbulence to decay due to 436 !--viscosity 437 438 ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. ) 439 IF ( lmix_mpc > 0 ) THEN 440 tau_mixingenv = ( lmix_mpc**2. / tke_dissip(i) )**(1./3.) 441 ELSE 442 tau_mixingenv = tau_mixenv 443 ENDIF 444 445 tau_dissipturb = gamma_taud * 2. * 2./3. * tke(i) / tke_dissip(i) / C0 446 447 !--------------------- PDF COMPUTATIONS --------------------- 448 !--Formulae for sigma2_pdf (variance), mean of PDF in Furtado2016 449 !--cloud liquid fraction and in-cloud liquid content are given 450 !--by integrating resp. F(Si) and Si*F(Si) 451 !--Liquid is limited by the available water vapor trough a 452 !--maximal liquid fraction 453 454 liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - qiceini_incl - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) ) 455 sigma2_pdf = 1./2. * ( ai**2 ) * 2./3. * tke(i) * tau_dissipturb / ( invtau_phaserelax + 1./tau_mixingenv ) 456 mean_pdf = sursat_env * 1./tau_mixingenv / ( invtau_phaserelax + 1./tau_mixingenv ) 457 cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - mean_pdf) / (SQRT(2.* sigma2_pdf) ) ) ) 458 !IF (cldfraliq(i) .GT. liqfra_max) THEN 459 ! cldfraliq(i) = liqfra_max 460 !ENDIF 461 462 qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.*RPI) * EXP( -1.*(sursat_iceliq - mean_pdf)**2. / (2.*sigma2_pdf) ) & 463 - qsati(i) * cldfraliq(i) * (sursat_iceliq - mean_pdf ) 464 465 sigma2_icefracturb(i)= sigma2_pdf 466 mean_icefracturb(i) = mean_pdf 467 !------------ ICE AMOUNT AND WATER CONSERVATION ------------ 468 469 IF ( (qliq_incl <= eps) .OR. (cldfraliq(i) <= eps) ) THEN 470 qliq_incl = 0. 471 cldfraliq(i) = 0. 472 END IF 473 474 !--Choice for in-cloud vapor : 475 !--1.Weighted mean between qvap in MPC parts and in ice-only parts 476 !--2.Always at ice saturation 477 qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) ) 478 479 IF ( qvap_incl >= qtot_incl(i) ) THEN 480 qvap_incl = qsati(i) 481 qliq_incl = qtot_incl(i) - qvap_incl 482 qice_incl = 0. 483 484 ELSEIF ( (qvap_incl + qliq_incl) >= qtot_incl(i) ) THEN 485 qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl) 486 qice_incl = 0. 487 ELSE 488 qice_incl = qtot_incl(i) - qvap_incl - qliq_incl 489 END IF 490 491 qvap_cld(i) = qvap_incl * cldfra1D 492 qliq(i) = qliq_incl * cldfra1D 493 qice(i) = qice_incl * cldfra1D 494 icefrac(i) = qice(i) / ( qice(i) + qliq(i) ) 495 dicefracdT(i) = 0. 496 !PRINT*,'MPC turb' 497 498 END IF ! Enough TKE 475 qice_incl = qtot_incl(i) - qvap_incl - qliq_incl 476 END IF 477 478 qvap_cld(i) = qvap_incl * cldfra1D 479 qliq(i) = qliq_incl * cldfra1D 480 qice(i) = qice_incl * cldfra1D 481 icefrac(i) = qice(i) / (qice(i) + qliq(i)) 482 dicefracdT(i) = 0. 483 !PRINT*,'MPC turb' 484 485 END IF ! Enough TKE 499 486 500 487 END IF ! MPC temperature 501 488 502 END IF ! cldfra503 504 ENDDO ! klon505 END SUBROUTINE ICEFRAC_LSCP_TURB506 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++507 508 509 SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs)510 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++489 END IF ! cldfra 490 491 ENDDO ! klon 492 END SUBROUTINE ICEFRAC_LSCP_TURB 493 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 494 495 496 SUBROUTINE CALC_QSAT_ECMWF(klon, temp, qtot, pressure, tref, phase, flagth, qs, dqs) 497 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 511 498 ! Calculate qsat following ECMWF method 512 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 513 USE lmdz_YOETHF 514 USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep 499 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 500 USE lmdz_yoethf 501 USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep 502 USE lmdz_yomcst 515 503 516 504 IMPLICIT NONE 517 518 include "YOMCST.h"519 505 520 506 INTEGER, INTENT(IN) :: klon ! number of horizontal grid points … … 522 508 REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg 523 509 REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa 524 REAL, INTENT(IN) 510 REAL, INTENT(IN) :: tref ! reference temperature in K 525 511 LOGICAL, INTENT(IN) :: flagth ! flag for qsat calculation for thermals 526 INTEGER, INTENT(IN) :: phase 512 INTEGER, INTENT(IN) :: phase 527 513 ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid) 528 514 ! 1=liquid … … 535 521 INTEGER i 536 522 537 DO i =1,klon538 539 IF (phase == 1) THEN540 delta =0.541 ELSEIF (phase == 2) THEN542 delta =1.543 ELSE544 delta =MAX(0.,SIGN(1.,tref-temp(i)))545 ENDIF546 547 IF (flagth) THEN548 cvm5=R5LES*(1.-delta) + R5IES*delta549 ELSE550 cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta551 cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i)))552 ENDIF553 554 qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i)555 qs(i)=MIN(0.5,qs(i))556 cor=1./(1.-RETV*qs(i))557 qs(i)=qs(i)*cor558 dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor)523 DO i = 1, klon 524 525 IF (phase == 1) THEN 526 delta = 0. 527 ELSEIF (phase == 2) THEN 528 delta = 1. 529 ELSE 530 delta = MAX(0., SIGN(1., tref - temp(i))) 531 ENDIF 532 533 IF (flagth) THEN 534 cvm5 = R5LES * (1. - delta) + R5IES * delta 535 ELSE 536 cvm5 = R5LES * RLVTT * (1. - delta) + R5IES * RLSTT * delta 537 cvm5 = cvm5 / RCPD / (1.0 + RVTMP2 * (qtot(i))) 538 ENDIF 539 540 qs(i) = R2ES * FOEEW(temp(i), delta) / pressure(i) 541 qs(i) = MIN(0.5, qs(i)) 542 cor = 1. / (1. - RETV * qs(i)) 543 qs(i) = qs(i) * cor 544 dqs(i) = FOEDE(temp(i), delta, cvm5, qs(i), cor) 559 545 560 546 END DO 561 547 562 END SUBROUTINE CALC_QSAT_ECMWF563 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++564 565 566 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++567 SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt)568 569 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++548 END SUBROUTINE CALC_QSAT_ECMWF 549 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 550 551 552 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 553 SUBROUTINE CALC_GAMMASAT(klon, temp, qtot, pressure, gammasat, dgammasatdt) 554 555 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 570 556 ! programme that calculates the gammasat parameter that determines the 571 557 ! homogeneous condensation thresholds for cold (<0oC) clouds 572 558 ! condensation at q>gammasat*qsat 573 559 ! Etienne Vignon, March 2021 574 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++560 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 575 561 576 562 USE lmdz_lscp_ini, ONLY: iflag_gammasat, t_glace_min, RTT 577 563 578 564 IMPLICIT NONE 579 580 565 581 566 INTEGER, INTENT(IN) :: klon ! number of horizontal grid points … … 588 573 REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature 589 574 590 REAL, DIMENSION(klon) :: qsi,qsl,dqsl,dqsi575 REAL, DIMENSION(klon) :: qsi, qsl, dqsl, dqsi 591 576 REAL fcirrus, fac 592 REAL, PARAMETER :: acirrus =2.349593 REAL, PARAMETER :: bcirrus =259.0577 REAL, PARAMETER :: acirrus = 2.349 578 REAL, PARAMETER :: bcirrus = 259.0 594 579 595 580 INTEGER i 596 597 CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.FALSE.,qsl,dqsl)598 CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.FALSE.,qsi,dqsi)599 600 DO i =1,klon601 602 603 604 gammasat(i)=1.605 dgammasatdt(i)=0.606 607 608 609 610 gammasat(i)=qsl(i)/qsi(i)611 dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)612 613 gammasat(i)=1.614 dgammasatdt(i)=0.615 616 617 618 619 620 621 622 623 fcirrus=acirrus-temp(i)/bcirrus624 IF (fcirrus > qsl(i)/qsi(i)) THEN625 gammasat(i)=qsl(i)/qsi(i)626 dgammasatdt(i)=(dqsl(i)*qsi(i)-dqsi(i)*qsl(i))/qsi(i)/qsi(i)627 628 gammasat(i)=fcirrus629 dgammasatdt(i)=-1.0/bcirrus630 631 632 633 634 gammasat(i)=1.635 dgammasatdt(i)=0.636 637 638 639 640 581 582 CALL CALC_QSAT_ECMWF(klon, temp, qtot, pressure, RTT, 1, .FALSE., qsl, dqsl) 583 CALL CALC_QSAT_ECMWF(klon, temp, qtot, pressure, RTT, 2, .FALSE., qsi, dqsi) 584 585 DO i = 1, klon 586 587 IF (temp(i) >= RTT) THEN 588 ! warm clouds: condensation at saturation wrt liquid 589 gammasat(i) = 1. 590 dgammasatdt(i) = 0. 591 592 ELSEIF ((temp(i) < RTT) .AND. (temp(i) > t_glace_min)) THEN 593 594 IF (iflag_gammasat >= 2) THEN 595 gammasat(i) = qsl(i) / qsi(i) 596 dgammasatdt(i) = (dqsl(i) * qsi(i) - dqsi(i) * qsl(i)) / qsi(i) / qsi(i) 597 ELSE 598 gammasat(i) = 1. 599 dgammasatdt(i) = 0. 600 ENDIF 601 602 ELSE 603 604 IF (iflag_gammasat >=1) THEN 605 ! homogeneous freezing of aerosols, according to 606 ! Koop, 2000 and Karcher 2008, QJRMS 607 ! 'Cirrus regime' 608 fcirrus = acirrus - temp(i) / bcirrus 609 IF (fcirrus > qsl(i) / qsi(i)) THEN 610 gammasat(i) = qsl(i) / qsi(i) 611 dgammasatdt(i) = (dqsl(i) * qsi(i) - dqsi(i) * qsl(i)) / qsi(i) / qsi(i) 612 ELSE 613 gammasat(i) = fcirrus 614 dgammasatdt(i) = -1.0 / bcirrus 615 ENDIF 616 617 ELSE 618 619 gammasat(i) = 1. 620 dgammasatdt(i) = 0. 621 622 ENDIF 623 624 ENDIF 625 641 626 END DO 642 627 643 644 END SUBROUTINE CALC_GAMMASAT 645 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 646 647 648 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 649 SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop) 650 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 651 652 USE lmdz_lscp_ini, ONLY: rd,rg,tresh_cl 653 654 IMPLICIT NONE 655 656 INTEGER, INTENT(IN) :: klon,klev !number of horizontal and vertical grid points 657 INTEGER, INTENT(IN) :: k ! vertical index 658 REAL, INTENT(IN), DIMENSION(klon,klev) :: temp ! temperature in K 659 REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa 660 REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa 661 REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb ! cloud fraction 662 663 REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D ! distance from cloud top 664 REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop ! temperature of cloud top 665 666 REAL dzlay(klon,klev) 667 REAL zlay(klon,klev) 668 REAL dzinterf 669 INTEGER i,k_top, kvert 670 LOGICAL bool_cl 671 672 673 DO i=1,klon 674 ! Initialization height middle of first layer 675 dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2)) 676 zlay(i,1) = dzlay(i,1)/2 677 678 DO kvert=2,klev 679 IF (kvert==klev) THEN 680 dzlay(i,kvert) = 2*(rd * temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert))) 681 ELSE 682 dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1)) 683 ENDIF 684 dzinterf = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert)) 685 zlay(i,kvert) = zlay(i,kvert-1) + dzinterf 686 ENDDO 687 ENDDO 688 689 DO i=1,klon 690 k_top = k 691 IF (rneb(i,k) <= tresh_cl) THEN 692 bool_cl = .FALSE. 693 ELSE 694 bool_cl = .TRUE. 695 ENDIF 696 697 DO WHILE ((bool_cl) .AND. (k_top <= klev)) 698 ! find cloud top 699 IF (rneb(i,k_top) > tresh_cl) THEN 700 k_top = k_top + 1 701 ELSE 702 bool_cl = .FALSE. 703 k_top = k_top - 1 704 ENDIF 705 ENDDO 706 k_top=min(k_top,klev) 707 708 !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to 709 !interf for layer of cloud top 710 distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2 711 temp_cltop(i) = temp(i,k_top) 712 ENDDO ! klon 713 714 END SUBROUTINE DISTANCE_TO_CLOUD_TOP 715 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 628 END SUBROUTINE CALC_GAMMASAT 629 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 630 631 632 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 633 SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon, klev, k, temp, pplay, paprs, rneb, distcltop1D, temp_cltop) 634 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 635 636 USE lmdz_lscp_ini, ONLY: rd, rg, tresh_cl 637 638 IMPLICIT NONE 639 640 INTEGER, INTENT(IN) :: klon, klev !number of horizontal and vertical grid points 641 INTEGER, INTENT(IN) :: k ! vertical index 642 REAL, INTENT(IN), DIMENSION(klon, klev) :: temp ! temperature in K 643 REAL, INTENT(IN), DIMENSION(klon, klev) :: pplay ! pressure middle layer in Pa 644 REAL, INTENT(IN), DIMENSION(klon, klev + 1) :: paprs ! pressure interfaces in Pa 645 REAL, INTENT(IN), DIMENSION(klon, klev) :: rneb ! cloud fraction 646 647 REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D ! distance from cloud top 648 REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop ! temperature of cloud top 649 650 REAL dzlay(klon, klev) 651 REAL zlay(klon, klev) 652 REAL dzinterf 653 INTEGER i, k_top, kvert 654 LOGICAL bool_cl 655 656 DO i = 1, klon 657 ! Initialization height middle of first layer 658 dzlay(i, 1) = Rd * temp(i, 1) / rg * log(paprs(i, 1) / paprs(i, 2)) 659 zlay(i, 1) = dzlay(i, 1) / 2 660 661 DO kvert = 2, klev 662 IF (kvert==klev) THEN 663 dzlay(i, kvert) = 2 * (rd * temp(i, kvert) / rg * log(paprs(i, kvert) / pplay(i, kvert))) 664 ELSE 665 dzlay(i, kvert) = rd * temp(i, kvert) / rg * log(paprs(i, kvert) / paprs(i, kvert + 1)) 666 ENDIF 667 dzinterf = rd * temp(i, kvert) / rg * log(pplay(i, kvert - 1) / pplay(i, kvert)) 668 zlay(i, kvert) = zlay(i, kvert - 1) + dzinterf 669 ENDDO 670 ENDDO 671 672 DO i = 1, klon 673 k_top = k 674 IF (rneb(i, k) <= tresh_cl) THEN 675 bool_cl = .FALSE. 676 ELSE 677 bool_cl = .TRUE. 678 ENDIF 679 680 DO WHILE ((bool_cl) .AND. (k_top <= klev)) 681 ! find cloud top 682 IF (rneb(i, k_top) > tresh_cl) THEN 683 k_top = k_top + 1 684 ELSE 685 bool_cl = .FALSE. 686 k_top = k_top - 1 687 ENDIF 688 ENDDO 689 k_top = min(k_top, klev) 690 691 !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to 692 !interf for layer of cloud top 693 distcltop1D(i) = zlay(i, k_top) - zlay(i, k) + dzlay(i, k_top) / 2 694 temp_cltop(i) = temp(i, k_top) 695 ENDDO ! klon 696 697 END SUBROUTINE DISTANCE_TO_CLOUD_TOP 698 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 716 699 717 700 END MODULE lmdz_lscp_tools
Note: See TracChangeset
for help on using the changeset viewer.
