Index: LMDZ6/branches/contrails/DefLists/file_def_histday_lmdz.xml =================================================================== --- LMDZ6/branches/contrails/DefLists/file_def_histday_lmdz.xml (revision 5612) +++ LMDZ6/branches/contrails/DefLists/file_def_histday_lmdz.xml (revision 5613) @@ -504,13 +504,4 @@ - - - - - - - - - Index: LMDZ6/branches/contrails/DefLists/file_def_histhf_lmdz.xml =================================================================== --- LMDZ6/branches/contrails/DefLists/file_def_histhf_lmdz.xml (revision 5612) +++ LMDZ6/branches/contrails/DefLists/file_def_histhf_lmdz.xml (revision 5613) @@ -479,5 +479,4 @@ - Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90 =================================================================== --- LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90 (revision 5612) +++ LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90 (revision 5613) @@ -174,18 +174,36 @@ pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qcritcont(i) / qsatl(i) * 100. - pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) - pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) - pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) - pdf_fra_above_qcritcont = EXP( - pdf_y ) * clrfra(i) - pdf_q_above_qcritcont = ( pdf_e3 * clrfra(i) & - + pdf_loc * pdf_fra_above_qcritcont ) * qsatl(i) / 100. + pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) + IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows + pdf_fra_above_qcritcont = 0. + pdf_q_above_qcritcont = 0. + ELSEIF ( pdf_y .LT. -10. ) THEN + pdf_fra_above_qcritcont = 1. + pdf_q_above_qcritcont = qclr(i) / clrfra(i) + ELSE + pdf_y = EXP( pdf_y ) + pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) + pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) + pdf_fra_above_qcritcont = EXP( - pdf_y ) * clrfra(i) + pdf_q_above_qcritcont = ( pdf_e3 * clrfra(i) & + + pdf_loc * pdf_fra_above_qcritcont ) * qsatl(i) / 100. + ENDIF pdf_x = qsat(i) / qsatl(i) * 100. - pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) - pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) - pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) - pdf_fra_above_qsat = EXP( - pdf_y ) * clrfra(i) - pdf_q_above_qsat = ( pdf_e3 * clrfra(i) & - + pdf_loc * pdf_fra_above_qsat ) * qsatl(i) / 100. + pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) + IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows + pdf_fra_above_qsat = 0. + pdf_q_above_qsat = 0. + ELSEIF ( pdf_y .LT. -10. ) THEN + pdf_fra_above_qsat = 1. + pdf_q_above_qsat = qclr(i) / clrfra(i) + ELSE + pdf_y = EXP( pdf_y ) + pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) + pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) + pdf_fra_above_qsat = EXP( - pdf_y ) * clrfra(i) + pdf_q_above_qsat = ( pdf_e3 * clrfra(i) & + + pdf_loc * pdf_fra_above_qsat ) * qsatl(i) / 100. + ENDIF potcontfraNP(i) = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat) Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90 =================================================================== --- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90 (revision 5612) +++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90 (revision 5613) @@ -239,5 +239,5 @@ ! ! for deposition / sublimation -REAL :: pres_sat, kappa_depsub, tau_depsub +REAL :: pres_sat, kappa_depsub, tauinv_depsub REAL :: air_thermal_conduct, water_vapor_diff REAL :: N_ice @@ -264,5 +264,5 @@ REAL :: pdf_fra_above_nuc, pdf_q_above_nuc REAL :: pdf_fra_above_lim, pdf_q_above_lim -REAL :: pdf_fra_below_lim, pdf_q_below_lim +REAL :: pdf_fra_below_lim REAL :: mixed_fraction ! @@ -545,12 +545,12 @@ !--If the cloud is initially supersaturated !--Exact explicit integration (qvc exact, qice explicit) - tau_depsub = 1. / ( depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub ) + tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub ELSE !--If the cloud is initially subsaturated !--Exact explicit integration (qvc exact, qice explicit) !--Same but depo_coef_cirrus = 1 - tau_depsub = 1. / ( qiceincld**(1./3.) * kappa_depsub ) + tauinv_depsub = qiceincld**(1./3.) * kappa_depsub ENDIF ! qvapincld .GT. qsat - qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime / tau_depsub ) + qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub ) !--If all the ice is sublimated IF ( qvapincld_new .GE. ( qvapincld + qiceincld ) ) qvapincld_new = 0. @@ -641,5 +641,5 @@ !--New PDF rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. - rhl_clr = MAX(0., MIN(1000., rhl_clr)) + rhl_clr = MAX(0., MIN(150., rhl_clr)) !--Calculation of the properties of the PDF @@ -662,4 +662,5 @@ pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. ) pdf_alpha(i) = EXP( rhl_clr / 100. ) * pdf_e3 + pdf_alpha(i) = MIN(10., pdf_alpha(i)) !--Avoid overflows !IF ( ok_warm_cloud ) THEN @@ -671,5 +672,6 @@ pdf_gamma(i) = GAMMA(1. + 1. / pdf_alpha(i)) - pdf_scale(i) = MAX(eps, pdf_std / SQRT( & + !--Barrier to avoid overflows + pdf_scale(i) = MAX(0.01, pdf_std / SQRT( & GAMMA(1. + 2. / pdf_alpha(i)) - pdf_gamma(i)**2 )) pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) @@ -680,9 +682,18 @@ pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100. - pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) - pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) - pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) - pdf_fra_above_nuc = EXP( - pdf_y ) - pdf_q_above_nuc = ( pdf_e3 + pdf_loc * pdf_fra_above_nuc ) * qsatl(i) / 100. + pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) + IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows + pdf_fra_above_nuc = 0. + pdf_q_above_nuc = 0. + ELSEIF ( pdf_y .LT. -10. ) THEN + pdf_fra_above_nuc = 1. + pdf_q_above_nuc = qclr(i) / clrfra(i) + ELSE + pdf_y = EXP( pdf_y ) + pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) + pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) + pdf_fra_above_nuc = EXP( - pdf_y ) + pdf_q_above_nuc = ( pdf_e3 + pdf_loc * pdf_fra_above_nuc ) * qsatl(i) / 100. + ENDIF IF ( pdf_fra_above_nuc .GT. eps ) THEN @@ -701,6 +712,7 @@ qvapincld = gamma_cond(i) * qsat(i) qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i) - tau_depsub = 1. / ( depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub ) - qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime / 2. / tau_depsub ) + tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub + qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) & + * EXP( - dtime / 2. * tauinv_depsub ) ELSE !--We keep the saturation adjustment hypothesis, and the vapor in the @@ -720,6 +732,11 @@ qclr(i) = qclr(i) - dqt_con - ENDIF ! ok_warm_cloud, cf_cond .GT. eps - ENDIF ! ( 1. - cldfra(i) ) .GT. eps + ENDIF ! pdf_fra_above_nuc .GT. eps + ELSE + !--Default value for the clear sky distribution: homogeneous distribution + pdf_alpha(i) = 1. + pdf_gamma(i) = 1. + pdf_scale(i) = 0.01 + ENDIF ! clrfra(i) .GT. eps @@ -826,10 +843,10 @@ IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN qiceinmix = ( qcld(i) - qvc(i) ) / cldfra(i) / ( 1. + clrfra_mix / cldfra_mix ) - tau_depsub = 1. / ( qiceinmix**(1./3.) * kappa_depsub ) - qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime / tau_depsub ) ) + tauinv_depsub = qiceinmix**(1./3.) * kappa_depsub + qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime * tauinv_depsub ) ) qvapinclr_lim = qvapinmix_lim * ( 1. + cldfra_mix / clrfra_mix ) & - qvc(i) / cldfra(i) * cldfra_mix / clrfra_mix ELSE - !--NB. if tau_depsub = 0, we get the same result as above + !--NB. if tau_depsub = 0 (ie tauinv_depsub = inf), we get the same result as above qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) & - qcld(i) / cldfra(i) * cldfra_mix / clrfra_mix @@ -851,13 +868,21 @@ pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qvapinclr_lim / qsatl(i) * 100. - pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) - pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) - pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) - pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) - pdf_q_above_lim = ( pdf_e3 * clrfra(i) & - + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. + pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) + IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows + pdf_fra_above_lim = 0. + pdf_q_above_lim = 0. + ELSEIF ( pdf_y .LT. -10. ) THEN + pdf_fra_above_lim = 1. + pdf_q_above_lim = qclr(i) / clrfra(i) + ELSE + pdf_y = EXP( pdf_y ) + pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) + pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) + pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) + pdf_q_above_lim = ( pdf_e3 * clrfra(i) & + + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. + ENDIF pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim - pdf_q_below_lim = qclr(i) - pdf_q_above_lim !--sigma_mix is the ratio of the surroundings of the clouds where mixing @@ -1000,10 +1025,19 @@ pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qvapinclr_lim / qsatl(i) * 100. - pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) - pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) - pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) - pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) - pdf_q_above_lim = ( pdf_e3 * clrfra(i) & - + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. + pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) + IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows + pdf_fra_above_lim = 0. + pdf_q_above_lim = 0. + ELSEIF ( pdf_y .LT. -10. ) THEN + pdf_fra_above_lim = 1. + pdf_q_above_lim = qclr(i) / clrfra(i) + ELSE + pdf_y = EXP( pdf_y ) + pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) + pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) + pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) + pdf_q_above_lim = ( pdf_e3 * clrfra(i) & + + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. + ENDIF pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim @@ -1086,6 +1120,6 @@ qiceincld = qice_sedim(i) / cldfra_above(i) IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN - tau_depsub = 1. / ( qiceincld**(1./3.) * kappa_depsub ) - qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime / tau_depsub ) ) + tauinv_depsub = qiceincld**(1./3.) * kappa_depsub + qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) ) ELSE qvapinclr_lim = qsat(i) - qiceincld @@ -1095,10 +1129,19 @@ pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qvapinclr_lim / qsatl(i) * 100. - pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) - pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) - pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) - pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) - pdf_q_above_lim = ( pdf_e3 * clrfra(i) & - + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. + pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) + IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows + pdf_fra_above_lim = 0. + pdf_q_above_lim = 0. + ELSEIF ( pdf_y .LT. -10. ) THEN + pdf_fra_above_lim = 1. + pdf_q_above_lim = qclr(i) / clrfra(i) + ELSE + pdf_y = EXP( pdf_y ) + pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) + pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) + pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) + pdf_q_above_lim = ( pdf_e3 * clrfra(i) & + + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. + ENDIF IF ( pdf_fra_above_lim .GT. eps ) THEN @@ -1140,9 +1183,18 @@ pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) pdf_x = qsat(i) / qsatl(i) * 100. - pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i) - pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) - pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) - issrfra(i) = EXP( - pdf_y ) * clrfra(i) - qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc * issrfra(i) ) * qsatl(i) / 100. + pdf_y = LOG( MAX( pdf_x - pdf_loc, eps ) / pdf_scale(i) ) * pdf_alpha(i) + IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows + issrfra(i) = 0. + qissr(i) = 0. + ELSEIF ( pdf_y .LT. -10. ) THEN + issrfra(i) = 1. + qissr(i) = qclr(i) / clrfra(i) + ELSE + pdf_y = EXP( pdf_y ) + pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) + pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) + issrfra(i) = EXP( - pdf_y ) * clrfra(i) + qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc * issrfra(i) ) * qsatl(i) / 100. + ENDIF ELSE issrfra(i) = 0. Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_precip.f90 =================================================================== --- LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_precip.f90 (revision 5612) +++ LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_precip.f90 (revision 5613) @@ -697,18 +697,20 @@ iwcg = MAX( zrho(i) * zoliqi(i) / zneb(i) * 1000., eps ) ! iwc in g/m3 - tempc = zt(i) - RTT ! celcius temp - ! so-called 'empirical parameterization' in Stubenrauch et al. 2019 - IF ( tempc .GE. -60. ) THEN - icevelo = EXP( 1.9714 + 0.00216078 * tempc & - - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) ) - ELSE - icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15 - ENDIF - icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s + !tempc = zt(i) - RTT ! celcius temp + !! so-called 'empirical parameterization' in Stubenrauch et al. 2019 + !IF ( tempc .GE. -60. ) THEN + ! icevelo = EXP( 1.9714 + 0.00216078 * tempc & + ! - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) ) + !ELSE + ! icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15 + !ENDIF + !icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s + + icevelo = fallice_sedim * 1.645 * ( iwcg / 1000. )**0.16 qice_sedim = zoliqi(i) * MIN(1., icevelo * dtime / zdz(i)) !--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two !--times in the same timestep) - !qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime) + qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime) !--Add tendencies zoliqi(i) = zoliqi(i) - qice_sedim @@ -1964,16 +1966,21 @@ rho = pplay(i) / temp(i) / RD iwcg = MAX( rho * qice(i) / cldfra(i) * 1000., eps ) ! iwc in g/m3 - tempc = temp(i) - RTT ! celcius temp - ! so-called 'empirical parameterization' in Stubenrauch et al. 2019 - IF ( tempc .GE. -60. ) THEN - icevelo = EXP( 1.9714 + 0.00216078 * tempc & - - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) ) - ELSE - icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15 - ENDIF - icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s + !tempc = temp(i) - RTT ! celcius temp + !! so-called 'empirical parameterization' in Stubenrauch et al. 2019 + !IF ( tempc .GE. -60. ) THEN + ! icevelo = EXP( 1.9714 + 0.00216078 * tempc & + ! - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) ) + !ELSE + ! icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15 + !ENDIF + !icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s + + icevelo = fallice_sedim * 1.645 * ( iwcg / 1000. )**0.16 dz = ( paprsdn(i) - paprsup(i) ) / RG / rho qice_sedim = qice(i) * MIN(1., icevelo * dtime / dz) + !--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two + !--times in the same timestep) + qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime) !--Add tendencies qice(i) = qice(i) - qice_sedim