Index: LMDZ5/trunk/libf/phylmd/cloudth_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/cloudth_mod.F90 (revision 2944) +++ LMDZ5/trunk/libf/phylmd/cloudth_mod.F90 (revision 2945) @@ -587,5 +587,5 @@ SUBROUTINE cloudth_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & - & qcloud,ctot,zpspsk,paprs,ztla,zthl, & + & qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, & & ratqs,zqs,t) @@ -624,14 +624,16 @@ REAL zqsatenv(ngrid,klev) - - REAL sigma1(ngrid,klev) + REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) REAL qlth(ngrid,klev) REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) - REAL cth(ngrid,klev) + REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) - REAL rneb(ngrid,klev) + REAL cth_vol(ngrid,klev) + REAL cenv_vol(ngrid,klev) + REAL ctot_vol(ngrid,klev) + REAL rneb(ngrid,klev) REAL t(ngrid,klev) REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi @@ -652,5 +654,5 @@ CALL cloudth_vert_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & - & qcloud,ctot,zpspsk,paprs,ztla,zthl, & + & qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, & & ratqs,zqs,t) RETURN @@ -672,4 +674,7 @@ cenv(:,:)=0. ctot(:,:)=0. + cth_vol(:,:)=0. + cenv_vol(:,:)=0. + ctot_vol(:,:)=0. qsatmmussig1=0. qsatmmussig2=0. @@ -747,4 +752,5 @@ cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) + ctot_vol(ind1,ind2)=ctot(ind1,ind2) qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2)) @@ -780,4 +786,5 @@ xenv=senv/(sqrt2*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) + ctot_vol(ind1,ind2)=ctot(ind1,ind2) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) @@ -801,5 +808,5 @@ SUBROUTINE cloudth_vert_v3(ngrid,klev,ind2, & & ztv,po,zqta,fraca, & - & qcloud,ctot,zpspsk,paprs,ztla,zthl, & + & qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, & & ratqs,zqs,t) @@ -838,5 +845,7 @@ REAL zqsatenv(ngrid,klev) - + + + REAL sigma1(ngrid,klev) REAL sigma2(ngrid,klev) @@ -844,7 +853,10 @@ REAL qlenv(ngrid,klev) REAL qltot(ngrid,klev) - REAL cth(ngrid,klev) + REAL cth(ngrid,klev) REAL cenv(ngrid,klev) REAL ctot(ngrid,klev) + REAL cth_vol(ngrid,klev) + REAL cenv_vol(ngrid,klev) + REAL ctot_vol(ngrid,klev) REAL rneb(ngrid,klev) REAL t(ngrid,klev) @@ -854,5 +866,5 @@ REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2 REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv - REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth + REAL IntJ,IntI1,IntI2,IntI3,IntJ_CF,IntI1_CF,IntI3_CF,coeffqlenv,coeffqlth REAL Tbef,zdelta,qsatbef,zcor REAL qlbef @@ -883,4 +895,7 @@ cenv(:,:)=0. ctot(:,:)=0. + cth_vol(:,:)=0. + cenv_vol(:,:)=0. + ctot_vol(:,:)=0. qsatmmussig1=0. qsatmmussig2=0. @@ -1032,39 +1047,53 @@ exp_xth2 = exp(-1.*xth2**2) + !CF_surfacique cth(ind1,ind2)=0.5*(1.-1.*erf(xth1)) cenv(ind1,ind2)=0.5*(1.-1.*erf(xenv1)) - ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) - - + ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2) + + + !CF_volumique & eau condense !environnement IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2 + IntJ_CF=0.5*(1.-1.*erf(xenv2)) if (deltasenv .lt. 1.e-10) then qlenv(ind1,ind2)=IntJ + cenv_vol(ind1,ind2)=IntJ_CF else IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1)) IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2) IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2) + IntI1_CF=((senv+deltasenv)*(erf(xenv2)-erf(xenv1)))/(4*deltasenv) + IntI3_CF=(sqrt2*sigma1s*(exp_xenv1-exp_xenv2))/(4*sqrtpi*deltasenv) qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 + cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF endif - !thermique IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2 - if (deltasth .lt. 1.e-10) then ! Seuil a choisir !!! + IntJ_CF=0.5*(1.-1.*erf(xth2)) + if (deltasth .lt. 1.e-10) then qlth(ind1,ind2)=IntJ + cth_vol(ind1,ind2)=IntJ_CF else IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1)) IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2) IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2) + IntI1_CF=((sth+deltasth)*(erf(xth2)-erf(xth1)))/(4*deltasth) + IntI3_CF=(sqrt2*sigma2s*(exp_xth1-exp_xth2))/(4*sqrtpi*deltasth) qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 + cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF endif + qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2) - + ctot_vol(ind1,ind2)=fraca(ind1,ind2)*cth_vol(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv_vol(ind1,ind2) ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4) + if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then ctot(ind1,ind2)=0. + ctot_vol(ind1,ind2)=0. qcloud(ind1)=zqsatenv(ind1,ind2) @@ -1100,4 +1129,5 @@ xenv=senv/(sqrt2*sigma1s) ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv)) + ctot_vol(ind1,ind2)=ctot(ind1,ind2) qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2)) Index: LMDZ5/trunk/libf/phylmd/conf_phys_m.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/conf_phys_m.F90 (revision 2944) +++ LMDZ5/trunk/libf/phylmd/conf_phys_m.F90 (revision 2945) @@ -169,4 +169,5 @@ INTEGER,SAVE :: iflag_t_glace_omp INTEGER,SAVE :: iflag_cloudth_vert_omp + INTEGER,SAVE :: iflag_rain_incloud_vol_omp REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp REAL,SAVE :: t_glace_min_omp, t_glace_max_omp @@ -1248,4 +1249,13 @@ ! + !Config Key = iflag_rain_incloud_vol + !Config Desc = + !Config Def = 0 + !Config Help = + ! + iflag_rain_incloud_vol_omp = 0 + CALL getin('iflag_rain_incloud_vol',iflag_rain_incloud_vol_omp) + + ! !Config Key = iflag_ice_thermo !Config Desc = @@ -2163,4 +2173,5 @@ iflag_t_glace = iflag_t_glace_omp iflag_cloudth_vert=iflag_cloudth_vert_omp + iflag_rain_incloud_vol=iflag_rain_incloud_vol_omp iflag_ice_thermo = iflag_ice_thermo_omp rei_min = rei_min_omp @@ -2511,4 +2522,5 @@ write(lunout,*)' iflag_t_glace = ',iflag_t_glace write(lunout,*)' iflag_cloudth_vert = ',iflag_cloudth_vert + write(lunout,*)' iflag_rain_incloud_vol = ',iflag_rain_incloud_vol write(lunout,*)' iflag_ice_thermo = ',iflag_ice_thermo write(lunout,*)' rei_min = ',rei_min Index: LMDZ5/trunk/libf/phylmd/fisrtilp.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/fisrtilp.F90 (revision 2944) +++ LMDZ5/trunk/libf/phylmd/fisrtilp.F90 (revision 2945) @@ -18,4 +18,5 @@ USE ioipsl_getin_p_mod, ONLY : getin_p USE phys_local_var_mod, ONLY: ql_seri,qs_seri + USE phys_local_var_mod, ONLY: rneblsvol ! flag to include modifications to ensure energy conservation (if flag >0) USE add_phys_tend_mod, only : fl_cor_ebil @@ -132,4 +133,5 @@ INTEGER ncoreczq REAL ctot(klon,klev) + REAL ctot_vol(klon,klev) REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 REAL zdqsdT_raw(klon) @@ -217,9 +219,11 @@ ! ice_thermo = iflag_ice_thermo .GE. 1 -itap=itap+1 -znebprecip(:)=0. - - ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3) + itap=itap+1 + znebprecip(:)=0. + + ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3) zdelq=0.0 + ctot_vol(1:klon,1:klev)=0.0 + rneblsvol(1:klon,1:klev)=0.0 if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM' @@ -727,9 +731,10 @@ call cloudth_v3(klon,klev,k,ztv, & zq,zqta,fraca, & - qcloud,ctot,zpspsk,paprs,ztla,zthl, & + qcloud,ctot,ctot_vol,zpspsk,paprs,ztla,zthl, & ratqs,zqs,t) endif do i=1,klon rneb(i,k)=ctot(i,k) + rneblsvol(i,k)=ctot_vol(i,k) zqn(i)=qcloud(i) enddo @@ -1148,6 +1153,16 @@ *fallvs(zrhol(i)) * zfice(i) endif - zchau = zct *dtime/REAL(ninter) * zoliq(i) & - *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )**2)) *(1.-zfice(i)) + + ! si l'heterogeneite verticale est active, on utilise + ! la fraction volumique "vraie" plutot que la fraction + ! surfacique modifiee, qui est plus grande et reduit + ! sinon l'eau in-cloud de facon artificielle + if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then + zchau = zct *dtime/REAL(ninter) * zoliq(i) & + *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl )**2)) *(1.-zfice(i)) + else + zchau = zct *dtime/REAL(ninter) * zoliq(i) & + *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )**2)) *(1.-zfice(i)) + endif !AJ< IF (.NOT. ice_thermo) THEN Index: LMDZ5/trunk/libf/phylmd/nuage.h =================================================================== --- LMDZ5/trunk/libf/phylmd/nuage.h (revision 2944) +++ LMDZ5/trunk/libf/phylmd/nuage.h (revision 2945) @@ -10,4 +10,5 @@ INTEGER iflag_t_glace, iflag_cloudth_vert, iflag_cld_cv + INTEGER iflag_rain_incloud_vol common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max, & @@ -15,4 +16,5 @@ & tau_cld_cv,coefw_cld_cv, & & tmax_fonte_cv, & - & iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv + & iflag_t_glace,iflag_cloudth_vert,iflag_cld_cv, & + & iflag_rain_incloud_vol !$OMP THREADPRIVATE(/nuagecom/) Index: LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 2944) +++ LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 2945) @@ -417,6 +417,6 @@ REAL,ALLOCATABLE,SAVE,DIMENSION(:,:) :: beta_prec !$OMP THREADPRIVATE(beta_prec) - REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rneb,rnebjn -!$OMP THREADPRIVATE(rneb,rnebjn) + REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rneb,rnebjn,rneblsvol +!$OMP THREADPRIVATE(rneb,rnebjn,rneblsvol) ! variables de sorties MM @@ -732,5 +732,5 @@ ALLOCATE(wdtrainA(klon,klev),wdtrainM(klon,klev)) ALLOCATE(beta_prec(klon,klev)) - ALLOCATE(rneb(klon,klev),rnebjn(klon,klev)) + ALLOCATE(rneb(klon,klev),rnebjn(klon,klev),rneblsvol(klon,klev)) Index: LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 2944) +++ LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 2945) @@ -1328,4 +1328,6 @@ TYPE(ctrl_out), SAVE :: o_rnebls = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), & 'rnebls', 'LS Cloud fraction', '-', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_rneblsvol = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), & + 'rneblsvol', 'LS Cloud fraction by volume', '-', (/ ('', i=1, 10) /)) TYPE(ctrl_out), SAVE :: o_rhum = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), & 'rhum', 'Relative humidity', '-', (/ ('', i=1, 10) /)) Index: LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 2944) +++ LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 2945) @@ -124,5 +124,5 @@ o_vitu, o_vitv, o_vitw, o_pres, o_paprs, & o_zfull, o_zhalf, o_rneb, o_rnebjn, o_rnebcon, & - o_rnebls, o_rhum, o_ozone, o_ozone_light, & + o_rnebls, o_rneblsvol, o_rhum, o_ozone, o_ozone_light, & o_duphy, o_dtphy, o_dqphy, o_dqphy2d, o_dqlphy, o_dqlphy2d, & o_dqsphy, o_dqsphy2d, o_albe_srf, o_z0m_srf, o_z0h_srf, & @@ -269,5 +269,5 @@ ql_seri, qs_seri, tr_seri, & zphi, u_seri, v_seri, omega, cldfra, & - rneb, rnebjn, zx_rh, d_t_dyn, & + rneb, rnebjn, rneblsvol, zx_rh, d_t_dyn, & d_q_dyn, d_ql_dyn, d_qs_dyn, & d_q_dyn2d, d_ql_dyn2d, d_qs_dyn2d, & @@ -1344,4 +1344,5 @@ CALL histwrite_phy(o_rnebcon, rnebcon) CALL histwrite_phy(o_rnebls, rneb) + CALL histwrite_phy(o_rneblsvol, rneblsvol) IF (vars_defined) THEN DO k=1, klev