Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90 (revision 5429) +++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90 (revision 5430) @@ -494,4 +494,5 @@ zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, & zqvapclr, zqupnew, & + cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), & zrfl, zrflclr, zrflcld, & zifl, ziflclr, ziflcld, & Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90 =================================================================== --- /LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90 (revision 5429) +++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90 (revision 5430) @@ -714,4 +714,5 @@ klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, & qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, & + cldfra, rvc_seri, qliq, qice, & rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub & ) @@ -720,5 +721,6 @@ USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG -USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub +USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub, ok_ice_supersat, ok_unadjusted_clouds +USE lmdz_lscp_ini, ONLY : eps, temp_nowater USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf @@ -747,10 +749,16 @@ REAL, INTENT(IN), DIMENSION(klon) :: qtotupnew !--total specific humidity IN THE LAYER ABOVE [kg/kg] +REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-] +REAL, INTENT(IN), DIMENSION(klon) :: rvc_seri !--cloud water vapor at the beginning of lscp (ratio wrt total water) - used only if the cloud properties are advected [kg/kg] +REAL, INTENT(IN), DIMENSION(klon) :: qliq !--liquid water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg] +REAL, INTENT(IN), DIMENSION(klon) :: qice !--ice water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg] + + REAL, INTENT(INOUT), DIMENSION(klon) :: rain !--flux of rain gridbox-mean coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2] -REAL, INTENT(IN), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2] +REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: snow !--flux of snow gridbox-mean coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2] -REAL, INTENT(IN), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2] +REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2] REAL, INTENT(OUT), DIMENSION(klon) :: dqreva !--rain tendency due to evaporation [kg/kg/s] @@ -762,13 +770,17 @@ !--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation REAL, DIMENSION(klon) :: dhum_to_dflux -!-- +!--Air density [kg/m3] and layer thickness [m] REAL, DIMENSION(klon) :: rho, dz +!--Temporary precip fractions and fluxes for the evaporation +REAL, DIMENSION(klon) :: precipfracclr_tmp, precipfraccld_tmp +REAL, DIMENSION(klon) :: rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp !--Saturation values REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati -!--Vapor in the clear sky -REAL :: qvapclr +!--Vapor in the clear sky and cloud +REAL :: qvapclr, qvapcld !--Fluxes tendencies because of evaporation and sublimation -REAL :: dprecip_evasub_max, draineva, dsnowsub, dprecip_evasub_tot +REAL :: dprecip_evasub_max, dprecip_evasub_tot +REAL :: drainclreva, draincldeva, dsnowclrsub, dsnowcldsub !--Specific humidity tendencies because of evaporation and sublimation REAL :: dqrevap, dqssubl @@ -828,26 +840,69 @@ ENDIF -! TODO Probleme : on utilise qvap total dans la maille pour l'evap / sub -! alors qu'on n'evap / sub que dans le ciel clair -! deux options pour cette routine : -! - soit on diagnostique le nuage AVANT l'evap / sub et on estime donc -! la fraction precipitante ciel clair dans la maille, ce qui permet de travailler -! avec des fractions, des fluxs et surtout un qvap dans le ciel clair -! - soit on pousse la param de Ludo au bout, et on prend un qvap de k+1 -! dans le ciel clair, avec un truc comme : -! qvapclr(k) = qvapclr(k+1)/qtot(k+1) * qtot(k) -! UPDATE : on code la seconde version. A voir si on veut mettre la premiere version. - - + +!--Initialise the precipitation fractions and fluxes DO i = 1, klon + precipfracclr_tmp(i) = precipfracclr(i) + precipfraccld_tmp(i) = precipfraccld(i) + rainclr_tmp(i) = rainclr(i) + snowclr_tmp(i) = snowclr(i) + raincld_tmp(i) = raincld(i) + snowcld_tmp(i) = snowcld(i) +ENDDO + +!--If we relax the LTP assumption that the cloud is the same than in the +!--gridbox above, we can change the tmp quantities +IF ( ok_ice_supersat ) THEN + !--Update the precipitation fraction using advected cloud fraction + CALL poprecip_fracupdate( & + klon, cldfra, precipfracclr_tmp, precipfraccld_tmp, & + rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp) + +! here we could put an ELSE statement and do one iteration of the condensation scheme +! (but it would only diagnose cloud from lognormal - link with thermals?) + + DO i = 1, klon + !--We cannot have a total fraction greater than the one in the gridbox above + !--because new cloud formation has not yet occured. + !--If this happens, we reduce the precip frac in the cloud, because it can + !--only mean that the cloud fraction at this level is greater than the total + !--precipitation fraction of the level above, and the precip frac in clear sky + !--is zero. + !--NB. this only works because we assume a max-random cloud and precip overlap + precipfraccld_tmp(i) = MIN( precipfraccld_tmp(i), precipfracclr(i) + precipfraccld(i) ) + ENDDO + +ENDIF +!--At this stage, we guarantee that +!-- precipfracclr + precipfraccld = precipfracclr_tmp + precipfraccld_tmp + + +DO i = 1, klon !--If there is precipitation from the layer above - ! NOTE TODO here we could replace the condition on precipfracclr(i) by a condition - ! such as eps or thresh_precip_frac, to remove the senseless barrier in the formulas - ! of evap / sublim - IF ( ( ( rain(i) + snow(i) ) .GT. 0. ) .AND. ( precipfracclr(i) .GT. 0. ) ) THEN - - IF ( ok_corr_vap_evasub ) THEN - !--Corrected version - we use the same water ratio between + IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN + + !--Init + drainclreva = 0. + draincldeva = 0. + dsnowclrsub = 0. + dsnowcldsub = 0. + + !--Init the properties of the air where reevaporation occurs + IF ( ok_ice_supersat ) THEN + !--We can diagnose the water vapor in the cloud using the advected + !--water vapor in the cloud and cloud fraction + IF ( ok_unadjusted_clouds .AND. ( cldfra(i) .GT. eps ) ) THEN + qvapcld = rvc_seri(i) * qvap(i) / cldfra(i) + ENDIF + !--We can diagnose completely the water vapor in clear sky, because all + !--the needed variables (ice, liq, vapor in cld, cloud fraction) are + !--advected + !--Note that qvap(i) is the total water in the gridbox + IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN + qvapclr = ( qvap(i) - qice(i) - qliq(i) - rvc_seri(i) * qvap(i) ) / ( 1. - cldfra(i) ) + ENDIF + ELSEIF ( ok_corr_vap_evasub ) THEN + !--Corrected version from Ludo - we use the same water ratio between !--the clear and the cloudy sky as in the layer above. This !--extends the assumption that the cloud fraction is the same @@ -856,5 +911,7 @@ !--Note that qvap(i) is the total water in the gridbox, and !--precipfraccld(i) is the cloud fraction in the layer above - qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) ) + IF ( ( 1. - precipfraccld(i) ) .GT. eps ) THEN + qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) ) + ENDIF ELSE !--Legacy version from Ludo - we use the total specific humidity @@ -863,56 +920,134 @@ ENDIF - !--Evaporation of liquid precipitation coming from above - !--in the clear sky only - !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva) - !--formula from Sundqvist 1988, Klemp & Wilhemson 1978 - !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly) - !--which does not need a barrier on rainclr, because included in the formula - draineva = precipfracclr(i) * ( MAX(0., & - - coef_eva * ( 1. - expo_eva ) * (1. - qvapclr / qsatl(i)) * dz(i) & - + ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_eva ) & - ) )**( 1. / ( 1. - expo_eva ) ) - rainclr(i) - - !--Evaporation is limited by 0 - draineva = MIN(0., draineva) - - - !--Sublimation of the solid precipitation coming from above - !--(same formula as for liquid precip) - !--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly) - !--which does not need a barrier on snowclr, because included in the formula - dsnowsub = precipfracclr(i) * ( MAX(0., & - - coef_sub * ( 1. - expo_sub ) * (1. - qvapclr / qsati(i)) * dz(i) & - + ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_sub ) & - ) )**( 1. / ( 1. - expo_sub ) ) - snowclr(i) - - !--Sublimation is limited by 0 - ! TODO: change max when we will allow for vapor deposition in supersaturated regions - dsnowsub = MIN(0., dsnowsub) - - !--Evaporation limit: we ensure that the layer's fraction below - !--the clear sky does not reach saturation. In this case, we - !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice - !--Max evaporation is computed not to saturate the clear sky precip fraction - !--(i.e., the fraction where evaporation occurs) - !--It is expressed as a max flux dprecip_evasub_max - - dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr(i)) & - * dhum_to_dflux(i) - dprecip_evasub_tot = draineva + dsnowsub - - !--Barriers - !--If activates if the total is LOWER than the max because - !--everything is negative - IF ( dprecip_evasub_tot .LT. dprecip_evasub_max ) THEN - draineva = dprecip_evasub_max * draineva / dprecip_evasub_tot - dsnowsub = dprecip_evasub_max * dsnowsub / dprecip_evasub_tot - ENDIF + !--------------------------- + !-- EVAP / SUBL IN CLEAR SKY + !--------------------------- + + IF ( precipfracclr_tmp(i) .GT. eps ) THEN + + !--Evaporation of liquid precipitation coming from above + !--in the clear sky only + !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva) + !--formula from Sundqvist 1988, Klemp & Wilhemson 1978 + !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly) + !--which does not need a barrier on rainclr, because included in the formula + drainclreva = precipfracclr_tmp(i) * MAX(0., & + - coef_eva * ( 1. - expo_eva ) * (1. - qvapclr / qsatl(i)) * dz(i) & + + ( rainclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_eva ) & + )**( 1. / ( 1. - expo_eva ) ) - rainclr_tmp(i) + + !--Evaporation is limited by 0 + !--NB. with ok_ice_supersat activated, this barrier should be useless + drainclreva = MIN(0., drainclreva) + + + !--Sublimation of the solid precipitation coming from above + !--(same formula as for liquid precip) + !--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly) + !--which does not need a barrier on snowclr, because included in the formula + dsnowclrsub = precipfracclr_tmp(i) * MAX(0., & + - coef_sub * ( 1. - expo_sub ) * (1. - qvapclr / qsati(i)) * dz(i) & + + ( snowclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_sub ) & + )**( 1. / ( 1. - expo_sub ) ) - snowclr_tmp(i) + + !--If ice supersaturation is activated, we allow vapor to condense on falling snow + !--i.e., having a positive dsnowclrsub + IF ( .NOT. ok_ice_supersat ) THEN + !--Sublimation is limited by 0 + dsnowclrsub = MIN(0., dsnowclrsub) + ENDIF + + + !--Evaporation limit: we ensure that the layer's fraction below + !--the clear sky does not reach saturation. In this case, we + !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice + !--Max evaporation is computed not to saturate the clear sky precip fraction + !--(i.e., the fraction where evaporation occurs) + !--It is expressed as a max flux dprecip_evasub_max + + IF ( .NOT. ok_ice_supersat ) THEN + dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr_tmp(i)) & + * dhum_to_dflux(i) + ELSE + dprecip_evasub_max = ( qvapclr - qsat(i) ) * precipfracclr_tmp(i) & + * dhum_to_dflux(i) + ENDIF + dprecip_evasub_tot = drainclreva + dsnowclrsub + + !--Barriers + !--If activates if the total is LOWER than the max because + !--everything is negative + IF ( (( dprecip_evasub_max .LT. 0. ) .AND. & + ( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. & + (( dprecip_evasub_max .GT. 0. ) .AND. & + ( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN + drainclreva = dprecip_evasub_max * drainclreva / dprecip_evasub_tot + dsnowclrsub = dprecip_evasub_max * dsnowclrsub / dprecip_evasub_tot + ENDIF + + ENDIF ! precipfracclr_tmp .GT. eps + + + !--------------------------- + !-- EVAP / SUBL IN THE CLOUD + !--------------------------- + + IF ( ok_unadjusted_clouds .AND. ( temp(i) .LE. temp_nowater ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN + !--Evaporation of liquid precipitation coming from above + !--in the cloud only + !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva) + !--formula from Sundqvist 1988, Klemp & Wilhemson 1978 + !--Exact explicit formulation (raincld is resolved exactly, qvap explicitly) + !--which does not need a barrier on raincld, because included in the formula + !draincldeva = precipfraccld_tmp(i) * MAX(0., & + ! - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) & + ! + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) & + ! )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i) + + !--Evaporation is limited by 0 + !--NB. with ok_ice_supersat activated, this barrier should be useless + !draincldeva = MIN(0., draincldeva) + draincldeva = 0. + + + !--Sublimation of the solid precipitation coming from above + !--(same formula as for liquid precip) + !--Exact explicit formulation (snowcld is resolved exactly, qvap explicitly) + !--which does not need a barrier on snowcld, because included in the formula + dsnowcldsub = precipfraccld_tmp(i) * MAX(0., & + - coef_sub * ( 1. - expo_sub ) * (1. - qvapcld / qsati(i)) * dz(i) & + + ( snowcld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_sub ) & + )**( 1. / ( 1. - expo_sub ) ) - snowcld_tmp(i) + + !--There is no barrier because we want deposition to occur if it is possible + + + !--Evaporation limit: we ensure that the layer's fraction below + !--the clear sky does not reach saturation. In this case, we + !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice + !--Max evaporation is computed not to saturate the clear sky precip fraction + !--(i.e., the fraction where evaporation occurs) + !--It is expressed as a max flux dprecip_evasub_max + + dprecip_evasub_max = ( qvapcld - qsat(i) ) * precipfraccld_tmp(i) * dhum_to_dflux(i) + dprecip_evasub_tot = draincldeva + dsnowcldsub + + !--Barriers + !--If activates if the total is LOWER than the max because + !--everything is negative + IF ( (( dprecip_evasub_max .LT. 0. ) .AND. & + ( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. & + (( dprecip_evasub_max .GT. 0. ) .AND. & + ( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN + draincldeva = dprecip_evasub_max * draincldeva / dprecip_evasub_tot + dsnowcldsub = dprecip_evasub_max * dsnowcldsub / dprecip_evasub_tot + ENDIF + + ENDIF ! ok_unadjusted_clouds !--New solid and liquid precipitation fluxes after evap and sublimation - dqrevap = draineva / dhum_to_dflux(i) - dqssubl = dsnowsub / dhum_to_dflux(i) - + dqrevap = ( drainclreva + draincldeva ) / dhum_to_dflux(i) + dqssubl = ( dsnowclrsub + dsnowcldsub ) / dhum_to_dflux(i) !--Vapor is updated after evaporation/sublimation (it is increased) @@ -929,11 +1064,45 @@ !--Add tendencies - !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) - 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 - !--sky is set to 0 + !--The MAX is needed because in some cases, the flux can be slightly + !--negative (numerical precision) + rainclr_tmp(i) = MAX(0., rainclr_tmp(i) + drainclreva) + snowclr_tmp(i) = MAX(0., snowclr_tmp(i) + dsnowclrsub) + raincld_tmp(i) = MAX(0., raincld_tmp(i) + draincldeva) + snowcld_tmp(i) = MAX(0., snowcld_tmp(i) + dsnowcldsub) + + + !--We reattribute fluxes according to initial fractions, using proportionnality + !--Note that trough this process, we conserve total precip fluxes + !-- rainclr_tmp + raincld_tmp = rainclr + raincld (same for snow) + !--If the cloud has increased, a part of the precip in clear sky falls in clear sky, + !--the rest falls in the cloud + IF ( ( precipfraccld(i) .GT. eps ) .AND. ( precipfraccld_tmp(i) .GT. precipfraccld(i) ) ) THEN + !--All the precip from cloud falls in the cloud + raincld(i) = MAX(0., raincld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i)) + rainclr(i) = MAX(0., rainclr_tmp(i) + raincld_tmp(i) - raincld(i)) + snowcld(i) = MAX(0., snowcld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i)) + snowclr(i) = MAX(0., snowclr_tmp(i) + snowcld_tmp(i) - snowcld(i)) + + !--If the cloud has narrowed, a part of the precip in cloud falls in clear sky, + !--the rest falls in the cloud + ELSEIF ( ( precipfracclr(i) .GT. eps ) .AND. ( precipfracclr_tmp(i) .GT. precipfracclr(i) ) ) THEN + !--All the precip from clear sky falls in clear sky + rainclr(i) = MAX(0., rainclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i)) + raincld(i) = MAX(0., raincld_tmp(i) + rainclr_tmp(i) - rainclr(i)) + snowclr(i) = MAX(0., snowclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i)) + snowcld(i) = MAX(0., snowcld_tmp(i) + snowclr_tmp(i) - snowclr(i)) + + ELSE + !--If the cloud stays the same or if there is no cloud above and + !--in the current layer, there is no reattribution + rainclr(i) = rainclr_tmp(i) + raincld(i) = raincld_tmp(i) + snowclr(i) = snowclr_tmp(i) + snowcld(i) = snowcld_tmp(i) + ENDIF + + !--If there is no more precip fluxes, the precipitation fraction is set to 0 IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0. + IF ( ( raincld(i) + snowcld(i) ) .LE. 0. ) precipfraccld(i) = 0. !--Calculation of the total fluxes @@ -955,4 +1124,133 @@ END SUBROUTINE poprecip_precld + + +!---------------------------------------------------------------- +! Computes the precipitation fraction update +! The goal of this routine is to reattribute precipitation fractions +! and fluxes to clear or cloudy air, depending on the variation of +! the cloud fraction on the vertical dimension. We assume a +! maximum-random overlap of the cloud cover (see Jakob and Klein, 2000, +! and LTP thesis, 2021) +! NB. in fact, we assume a maximum-random overlap of the total precip. frac +! +SUBROUTINE poprecip_fracupdate( & + klon, cldfra, precipfracclr, precipfraccld, & + rainclr, raincld, snowclr, snowcld) + +USE lmdz_lscp_ini, ONLY : eps + +IMPLICIT NONE + +INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-] + +REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction [-] +REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-] +REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-] + !--NB. at the end of the routine, becomes the fraction of precip + !--in the current layer + +REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2] +REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2] +REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2] +REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2] + +!--Local variables +INTEGER :: i +REAL :: dcldfra +REAL :: precipfractot +REAL :: dprecipfracclr, dprecipfraccld +REAL :: drainclr, dsnowclr +REAL :: draincld, dsnowcld + + +DO i = 1, klon + + !--Initialisation + precipfractot = precipfracclr(i) + precipfraccld(i) + + !--Instead of using the cloud cover which was use in LTP thesis, we use the + !--total precip. fraction to compute the maximum-random overlap. This is + !--because all the information of the cloud cover is embedded into + !--precipfractot, and this allows for taking into account the potential + !--reduction of the precipitation fraction because either the flux is too + !--small (see barrier at the end of poprecip_postcld) or the flux is completely + !--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 ) * & + ! ( 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. - precipfraccld(i) ) + ENDIF + + !--precipfraccld(i) is here the cloud fraction of the layer above + dcldfra = cldfra(i) - precipfraccld(i) + !--Tendency of the clear-sky precipitation fraction. We add a MAX on the + !--calculation of the current CS precip. frac. + !dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i) + !--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated + !--if precipfractot < cldfra) + dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i) + !--Tendency of the cloudy precipitation fraction. We add a MAX on the + !--calculation of the current CS precip. frac. + !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) ) + !--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated + !--if cldfra < 0) + dprecipfraccld = dcldfra + + + !--If the cloud extends + IF ( dprecipfraccld .GT. 0. ) THEN + !--If there is no CS precip, nothing happens. + !--If there is, we reattribute some of the CS precip flux + !--to the cloud precip flux, proportionnally to the + !--decrease of the CS precip fraction + IF ( precipfracclr(i) .LT. eps ) THEN + drainclr = 0. + dsnowclr = 0. + ELSE + drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i) + dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i) + ENDIF + !--If the cloud narrows + ELSEIF ( dprecipfraccld .LT. 0. ) THEN + !--We reattribute some of the cloudy precip flux + !--to the CS precip flux, proportionnally to the + !--decrease of the cloud precip fraction + IF ( precipfraccld(i) .LT. eps ) THEN + draincld = 0. + dsnowcld = 0. + ELSE + draincld = dprecipfraccld / precipfraccld(i) * raincld(i) + dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i) + ENDIF + drainclr = - draincld + dsnowclr = - dsnowcld + !--If the cloud stays the same or if there is no cloud above and + !--in the current layer, nothing happens + ELSE + drainclr = 0. + dsnowclr = 0. + ENDIF + + !--We add the tendencies + !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) + precipfraccld(i) = precipfraccld(i) + dprecipfraccld + precipfracclr(i) = precipfracclr(i) + dprecipfracclr + 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) + +ENDDO + +END SUBROUTINE poprecip_fracupdate @@ -1043,11 +1341,4 @@ REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip -!--Partition of the fluxes -REAL :: dcldfra -REAL :: precipfractot -REAL :: dprecipfracclr, dprecipfraccld -REAL :: drainclr, dsnowclr -REAL :: draincld, dsnowcld - !--Collection, aggregation and riming REAL :: eff_cldfra @@ -1098,4 +1389,10 @@ +!--Update the precipitation fraction following cloud formation +CALL poprecip_fracupdate( & + klon, cldfra, precipfracclr, precipfraccld, & + rainclr, raincld, snowclr, snowcld) + + DO i = 1, klon @@ -1111,92 +1408,4 @@ qtot(i) = qvap(i) + qliq(i) + qice(i) & + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i) - - !------------------------------------------------------------ - !-- PRECIPITATION FRACTIONS UPDATE - !------------------------------------------------------------ - !--The goal of this routine is to reattribute precipitation fractions - !--and fluxes to clear or cloudy air, depending on the variation of - !--the cloud fraction on the vertical dimension. We assume a - !--maximum-random overlap of the cloud cover (see Jakob and Klein, 2000, - !--and LTP thesis, 2021) - !--NB. in fact, we assume a maximum-random overlap of the total precip. frac - - !--Initialisation - precipfractot = precipfracclr(i) + precipfraccld(i) - - !--Instead of using the cloud cover which was use in LTP thesis, we use the - !--total precip. fraction to compute the maximum-random overlap. This is - !--because all the information of the cloud cover is embedded into - !--precipfractot, and this allows for taking into account the potential - !--reduction of the precipitation fraction because either the flux is too - !--small (see barrier at the end of poprecip_postcld) or the flux is completely - !--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 ) * & - ! ( 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. - precipfraccld(i) ) - ENDIF - - !--precipfraccld(i) is here the cloud fraction of the layer above - dcldfra = cldfra(i) - precipfraccld(i) - !--Tendency of the clear-sky precipitation fraction. We add a MAX on the - !--calculation of the current CS precip. frac. - !dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i) - !--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated - !--if precipfractot < cldfra) - dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i) - !--Tendency of the cloudy precipitation fraction. We add a MAX on the - !--calculation of the current CS precip. frac. - !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) ) - !--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated - !--if cldfra < 0) - dprecipfraccld = dcldfra - - - !--If the cloud extends - IF ( dprecipfraccld .GT. 0. ) THEN - !--If there is no CS precip, nothing happens. - !--If there is, we reattribute some of the CS precip flux - !--to the cloud precip flux, proportionnally to the - !--decrease of the CS precip fraction - IF ( precipfracclr(i) .LE. 0. ) THEN - drainclr = 0. - dsnowclr = 0. - ELSE - drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i) - dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i) - ENDIF - !--If the cloud narrows - ELSEIF ( dprecipfraccld .LT. 0. ) THEN - !--We reattribute some of the cloudy precip flux - !--to the CS precip flux, proportionnally to the - !--decrease of the cloud precip fraction - draincld = dprecipfraccld / precipfraccld(i) * raincld(i) - dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i) - drainclr = - draincld - dsnowclr = - dsnowcld - !--If the cloud stays the same or if there is no cloud above and - !--in the current layer, nothing happens - ELSE - drainclr = 0. - dsnowclr = 0. - ENDIF - - !--We add the tendencies - !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) - precipfraccld(i) = precipfraccld(i) + dprecipfraccld - precipfracclr(i) = precipfracclr(i) + dprecipfracclr - 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