Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90 (revision 4831) +++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90 (revision 4832) @@ -161,4 +161,7 @@ !$OMP THREADPRIVATE(rho_rain) + REAL, SAVE, PROTECTED :: rho_ice=920. ! A COMMENTER TODO [kg/m3] + !$OMP THREADPRIVATE(rho_ice) + REAL, SAVE, PROTECTED :: rho_snow ! A COMMENTER TODO [kg/m3] !$OMP THREADPRIVATE(rho_snow) @@ -169,4 +172,7 @@ REAL, SAVE, PROTECTED :: r_snow=1.E-3 ! A COMMENTER TODO [m] !$OMP THREADPRIVATE(r_snow) + + REAL, SAVE, PROTECTED :: r_ice=50.E-6 ! A COMMENTER TODO [m] + !$OMP THREADPRIVATE(r_ice) REAL, SAVE, PROTECTED :: tau_auto_snow_min=100. ! A COMMENTER TODO [s] Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90 (revision 4831) +++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90 (revision 4832) @@ -136,5 +136,4 @@ !--multiplying by dz = - dP / g / rho (line 3-4) !--formula from Sundqvist 1988, Klemp & Wilhemson 1978 - draineva = - precipfracclr(i) * coef_eva * (1. - qvap(i) / qsatl(i)) & * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva & @@ -197,6 +196,7 @@ !--Add tendencies - rainclr(i) = rainclr(i) + draineva - snowclr(i) = snowclr(i) + dsnowsub + + rainclr(i) = MAX(0., rainclr(i) + draineva) + snowclr(i) = MAX(0., snowclr(i) + dsnowsub) !--If there is no more precip fluxes, the precipitation fraction in clear @@ -248,5 +248,5 @@ cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, & thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, & - rho_rain, rho_snow, r_rain, r_snow, & + rho_rain, rho_snow, r_rain, r_snow, rho_ice, r_ice, & tau_auto_snow_min, tau_auto_snow_max, & thresh_precip_frac, eps, air_thermal_conduct, & @@ -340,5 +340,5 @@ REAL :: dqrfreez_max REAL :: tau_freez -REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez +REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2 REAL :: coef_freez REAL :: dqifreez !--loss of ice cloud content due to collection of ice from rain [kg/kg/s] @@ -371,5 +371,5 @@ dqifreez= 0. - !--hum_to_flux: coef to convert a specific quantity to a flux + !--hum_to_flux: coef to convert a delta in specific quantity to a flux !-- hum_to_flux = rho * dz/dt = 1 / g * dP/dt hum_to_flux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime @@ -398,8 +398,16 @@ !--evaporated (see barrier at the end of poprecip_precld) !--NB. precipfraccld(i) is here the cloud fraction of the layer above - precipfractot = 1. - ( 1. - precipfractot ) * & + !precipfractot = 1. - ( 1. - precipfractot ) * & + ! ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) & + ! / ( 1. - MIN( precipfraccld(i), 1. - eps ) ) + + + IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN + precipfractot = 1. + ELSE + precipfractot = 1. - ( 1. - precipfractot ) * & ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) & - / ( 1. - MIN( precipfraccld(i), 1. - eps ) ) - + / ( 1. - precipfraccld(i) ) + ENDIF !--precipfraccld(i) is here the cloud fraction of the layer above @@ -451,10 +459,10 @@ precipfraccld(i) = precipfraccld(i) + dprecipfraccld precipfracclr(i) = precipfracclr(i) + dprecipfracclr - rainclr(i) = rainclr(i) + drainclr - snowclr(i) = snowclr(i) + dsnowclr - raincld(i) = raincld(i) - drainclr - snowcld(i) = snowcld(i) - dsnowclr + + rainclr(i) = MAX( 0. , rainclr(i) + drainclr ) + snowclr(i) = MAX( 0. , snowclr(i) + dsnowclr ) + raincld(i) = MAX( 0. , raincld(i) - drainclr ) + snowcld(i) = MAX( 0. , snowcld(i) - dsnowclr ) - !--If vertical heterogeneity is taken into account, we use !--the "true" volume fraction instead of a modified @@ -513,5 +521,5 @@ !--Add tendencies qliq(i) = qliq(i) + dqlcol - raincld(i) = raincld(i) - dqlcol * hum_to_flux(i) + raincld(i) = MAX( 0. , raincld(i) - dqlcol * hum_to_flux(i) ) !--Diagnostic tendencies @@ -535,5 +543,5 @@ !--Add tendencies qice(i) = qice(i) + dqiagg - snowcld(i) = snowcld(i) - dqiagg * hum_to_flux(i) + snowcld(i) = MAX( 0. , snowcld(i) - dqiagg * hum_to_flux(i) ) !--Diagnostic tendencies @@ -599,6 +607,6 @@ qliq(i) = qliq(i) + dqlauto qice(i) = qice(i) + dqiauto - raincld(i) = raincld(i) - dqlauto * hum_to_flux(i) - snowcld(i) = snowcld(i) - dqiauto * hum_to_flux(i) + raincld(i) = MAX( 0. , raincld(i) - dqlauto * hum_to_flux(i) ) + snowcld(i) = MAX( 0. , snowcld(i) - dqiauto * hum_to_flux(i) ) !--Diagnostic tendencies @@ -631,5 +639,5 @@ !--Add tendencies qliq(i) = qliq(i) + dqlrim - snowcld(i) = snowcld(i) - dqlrim * hum_to_flux(i) + snowcld(i) = MAX( 0. , snowcld(i) - dqlrim * hum_to_flux(i) ) !--Temperature adjustment with the release of latent @@ -699,5 +707,5 @@ !--Barrier to limit the melting flux to the clr snow flux in the mesh dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / hum_to_flux(i)) - ENDIF + ENDIF !--In cloudy air @@ -726,11 +734,14 @@ dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt dqstotmelt = dqsmelt_max + ENDIF !--Add tendencies - rainclr(i) = rainclr(i) - dqsclrmelt * hum_to_flux(i) - raincld(i) = raincld(i) - dqscldmelt * hum_to_flux(i) - snowclr(i) = snowclr(i) + dqsclrmelt * hum_to_flux(i) - snowcld(i) = snowcld(i) + dqscldmelt * hum_to_flux(i) + + rainclr(i) = MAX( 0. , rainclr(i) - dqsclrmelt * hum_to_flux(i) ) + raincld(i) = MAX( 0. , raincld(i) - dqscldmelt * hum_to_flux(i) ) + snowclr(i) = MAX( 0. , snowclr(i) + dqsclrmelt * hum_to_flux(i) ) + snowcld(i) = MAX( 0. , snowcld(i) + dqscldmelt * hum_to_flux(i) ) + !--Temperature adjustment with the release of latent @@ -749,8 +760,10 @@ !-- FREEZING !--------------------------------------------------------- - !--Process through which rain freezes into snow. This is - !--parameterized as an exponential decrease of the rain - !--water content. - !--The formula is homemade. + !--Process through which rain freezes into snow. + !-- We parameterize it as a 2 step process: + !-- first: freezing following collision with ice crystals + !-- second: immersion freezing following (inspired by Bigg 1953) + !-- the latter is parameterized as an exponential decrease of the rain + !-- water content with a homemade formulya !--This is based on a caracteritic time of freezing, which !--exponentially depends on temperature so that it is @@ -758,5 +771,4 @@ !--0 for RTT (=273.15 K). !--NB.: this process needs a temperature adjustment - !--dqrfreez_max: maximum rain freezing so that temperature !-- stays lower than 273 K [kg/kg] @@ -768,30 +780,25 @@ IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN - !--Computed according to - !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q - dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD & + + !--1st step: freezing following collision with ice crystals + !--Sub-process of freezing which quantifies the collision between + !--ice crystals in suspension and falling rain droplets. + !--The rain droplets freeze, becoming graupel, and carrying + !--the ice crystal (which acted as an ice nucleating particle). + !--The formula is adapted from the riming formula. + !-- it works only in the cloudy part + + dqifreez = 0. + dqrtotfreez_step1 = 0. + + IF ((qice(i) .GT. 0.) .AND. (raincld(i) .GT. 0.)) THEN + dqrclrfreez = 0. + dqrcldfreez = 0. + + !--Computed according to + !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q + dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD & * ( 1. + RVTMP2 * qtot(i) )) - tau_freez = 1. / ( beta_freez & - * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) ) - - !--Initialisation - dqrclrfreez = 0. - dqrcldfreez = 0. - - !--In clear air - IF ( rainclr(i) .GT. 0. ) THEN - !--Exact solution of dqrain/dt = -qrain/tau_freez - dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. ) - ENDIF - - !--In cloudy air - IF ( raincld(i) .GT. 0. ) THEN - !--Exact solution of dqrain/dt = -qrain/tau_freez - dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. ) - - !--------------------------------------------------------- - !-- FREEZING OF RAIN WITH CONTACT OF ICE CRYSTALS - !--------------------------------------------------------- !--Sub-process of freezing which quantifies the collision between !--ice crystals in suspension and falling rain droplets. @@ -800,43 +807,86 @@ !--The formula is adapted from the riming formula. - !--The sticking efficacy is perfect. + !--The collision efficiency is assumed unity Eff_rain_ice = 1. coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice - !-- ATTENTION Double implicit version? + barriers if needed !--Exact version, which does not need a barrier because of - !--the exponential decrease + !--the exponential decrease. dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. ) - !--We add the two freezing processes in cloud - dqrcldfreez = dqrcldfreez - dqifreez - ENDIF - - !--Barriers + !--We add the part of rain water that freezes, limited by a temperature barrier + !-- this quantity is calculated assuming that the number of drop that freeze correspond to the number + !-- of crystals collected (and assuming uniform distributions of ice crystals and rain drops) + dqrcldfreez = MAX(dqifreez*rho_rain*r_rain/(rho_ice*r_ice),-raincld(i)/hum_to_flux(i)) + dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max) + dqrtotfreez_step1 = dqrcldfreez + + !--Add tendencies + qice(i) = qice(i) + dqifreez + raincld(i) = MAX( 0. , raincld(i) + dqrcldfreez * hum_to_flux(i) ) + snowcld(i) = MAX( 0. , snowcld(i) - dqrcldfreez * hum_to_flux(i) - dqifreez * hum_to_flux(i) ) + temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD & + / ( 1. + RVTMP2 * qtot(i) ) + + ENDIF + + !-- Second step immersion freezing of rain drops + !-- with a homemade timeconstant depending on temperature + + dqrclrfreez = 0. + dqrcldfreez = 0. + dqrtotfreez_step2 = 0. + !--Computed according to + !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q + + dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD & + * ( 1. + RVTMP2 * qtot(i) )) + + + tau_freez = 1. / ( beta_freez & + * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) ) + + + !--In clear air + IF ( rainclr(i) .GT. 0. ) THEN + !--Exact solution of dqrain/dt = -qrain/tau_freez + dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. ) + ENDIF + + !--In cloudy air + IF ( raincld(i) .GT. 0. ) THEN + !--Exact solution of dqrain/dt = -qrain/tau_freez + dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. ) + ENDIF + + !--temperature barrie step 2 !--It is activated if the total is LOWER than the max !--because everything is negative - dqrtotfreez = dqrclrfreez + dqrcldfreez - IF ( dqrtotfreez .LT. dqrfreez_max ) THEN + dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez + + IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN !--We redistribute the max freezed rain keeping !--the clear/cloud partition of the freezing rain - dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez - dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez - dqrtotfreez = dqrfreez_max + dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2 + dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2 + dqrtotfreez_step2 = dqrfreez_max ENDIF !--Add tendencies - rainclr(i) = rainclr(i) + dqrclrfreez * hum_to_flux(i) - raincld(i) = raincld(i) + dqrcldfreez * hum_to_flux(i) - snowclr(i) = snowclr(i) - dqrclrfreez * hum_to_flux(i) - snowcld(i) = snowcld(i) - dqrcldfreez * hum_to_flux(i) + rainclr(i) = MAX( 0. , rainclr(i) + dqrclrfreez * hum_to_flux(i) ) + raincld(i) = MAX( 0. , raincld(i) + dqrcldfreez * hum_to_flux(i) ) + snowclr(i) = MAX( 0. , snowclr(i) - dqrclrfreez * hum_to_flux(i) ) + snowcld(i) = MAX( 0. , snowcld(i) - dqrcldfreez * hum_to_flux(i) ) + !--Temperature adjustment with the uptake of latent !--heat because of freezing - temp(i) = temp(i) - dqrtotfreez * RLMLT / RCPD & + temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD & / ( 1. + RVTMP2 * qtot(i) ) + dqrtotfreez=dqrtotfreez_step1+dqrtotfreez_step2 !--Diagnostic tendencies dqrfreez(i) = dqrtotfreez / dtime - dqsfreez(i) = - dqrtotfreez / dtime + dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime ENDIF