Index: /trunk/LMDZ.GENERIC/README =================================================================== --- /trunk/LMDZ.GENERIC/README (revision 730) +++ /trunk/LMDZ.GENERIC/README (revision 731) @@ -738,2 +738,5 @@ gcm.e runs in Earth and Early Mars case with CO2 and H2O cycle + dust. +== 19/07/2012 == JL +- Corrected precipitation evaporation scheme + snow fall (JL+BC) +- Commented some unecessary writediagfi calls that where used for test. Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 (revision 730) +++ /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 (revision 731) @@ -1644,9 +1644,9 @@ call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv) call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw) - call writediagfi(ngrid,'p','Pressure','Pa',3,pplay) + call writediagfi(ngrid,"p","Pressure","Pa",3,pplay) ! Total energy balance diagnostics if(callrad.and.(.not.newtonian))then - call writediagfi(ngrid,'ALB','Surface albedo',' ',2,albedo) + call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo) call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn) call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw) @@ -1659,12 +1659,12 @@ if (calldifv) then call writediagfi(ngridmx,"q2","turbulent kinetic energy","J.kg^-1",3,q2) - call writediagfi(ngridmx,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff) - call writediagfi(ngridmx,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff) - call writediagfi(ngridmx,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs) +! call writediagfi(ngridmx,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff) +! call writediagfi(ngridmx,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff) +! call writediagfi(ngridmx,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs) call writediagfi(ngridmx,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux) endif if (corrk) then - call writediagfi(ngridmx,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw) - call writediagfi(ngridmx,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw) +! call writediagfi(ngridmx,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw) +! call writediagfi(ngridmx,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw) endif if(watercond) then @@ -1672,6 +1672,6 @@ call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE) call writediagfi(ngridmx,"RH","relative humidity"," ",3,RH) - call writediagfi(ngridmx,"qsatatm","atm qsat"," ",3,qsat) - call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol) +! call writediagfi(ngridmx,"qsatatm","atm qsat"," ",3,qsat) +! call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol) endif endif @@ -1709,7 +1709,7 @@ if(watercond.or.CLFvarying)then - call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man) - call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc) - call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac) +! call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man) +! call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc) +! call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac) call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac) endif Index: /trunk/LMDZ.GENERIC/libf/phystd/rain.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/rain.F90 (revision 730) +++ /trunk/LMDZ.GENERIC/libf/phystd/rain.F90 (revision 731) @@ -102,6 +102,7 @@ ! Online functions - REAL fallv, zzz ! falling speed of ice crystals + REAL fallv, fall2v, zzz ! falling speed of ice crystals fallv (zzz) = 3.29 * ((zzz)**0.16) + fall2v (zzz) =10.6 * ((zzz)**0.31) !for use with radii DATA firstcall /.true./ @@ -199,5 +200,4 @@ DO k = 1, nlayermx DO i = 1, ngridmx - !l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/(g*ptimestep) l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g ENDDO @@ -215,12 +215,11 @@ if(zt(i,k).gt.Tsat(i,k))then -! print*,'in rain',i,k,zrfl(i),zt(i,k),Tsat(i,k) +! treat the case where all liquid water should boil zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*l2c(i,k)/RLVTT,zrfl(i)) zrfl(i)=MAX(zrfl(i)-zqev,0.) d_q(i,k)=zqev/l2c(i,k) d_t(i,k) = - d_q(i,k) * RLVTT/RCPD -! print*,zqev,zrfl(i),d_q(i,k),d_t(i,k) else - zqev = MAX (0.0, (zqs(i,k)-q(i,k)))/ptimestep! there is a bug here + zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*l2c(i,k)/ptimestep !there was a bug here zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5 *sqrt(zrfl(i))*l2c(i,k)/pplay(i,k)*zt(i,k)*R ! BC modif here @@ -272,5 +271,5 @@ elseif (precip_scheme.ge.2) then - DO i = 1, ngridmx + DO i = 1, ngridmx IF (rneb(i,k).GT.0.0) THEN zoliq(i) = ql(i,k) @@ -282,5 +281,8 @@ zneb(i) = MAX(rneb(i,k), seuil_neb) ENDIF - ENDDO + ENDDO + + !recalculate liquid water particle radii + call h2o_cloudrad(ql,reffh2oliq,reffh2oice) SELECT CASE(precip_scheme) @@ -297,6 +299,5 @@ zrhol(i) = zrho(i) * zoliq(i) / zneb(i) zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & - *fallv(zrhol(i)) * (1.0-zfrac(i)) ! zfroi behaves oddly... -! * 0.1 * (1.0-zfrac(i)) + *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... ztot(i) = zchau(i) + zfroi(i) @@ -320,6 +321,5 @@ zrhol(i) = zrho(i) * zoliq(i) / zneb(i) zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & - *fallv(zrhol(i)) * (1.0-zfrac(i)) ! zfroi behaves oddly... -! * 0.1 * (1.0-zfrac(i)) + *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... ztot(i) = zchau(i) + zfroi(i) @@ -334,7 +334,4 @@ !precip scheme modified from Boucher 95 CASE(4) - - !recalculate liquid water particle radii - call h2o_cloudrad(ql,reffh2oliq,reffh2oice) DO n = 1, ninter @@ -347,6 +344,5 @@ zrhol(i) = zrho(i) * zoliq(i) / zneb(i) zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & - *fallv(zrhol(i)) * (1.0-zfrac(i)) ! zfroi behaves oddly... -! * 0.1 * (1.0-zfrac(i)) + *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... ztot(i) = zchau(i) + zfroi(i) Index: /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90 (revision 730) +++ /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90 (revision 731) @@ -517,14 +517,4 @@ ! calculate dQsat / dT by finite differences ! we cannot use the updated temperature value yet... -! if(qsat2(ig).gt.1.) then -! qsat(ig)=1. -! dqsat(ig)=0. -! else -! qsat(ig)=qsat2(ig) -! endif - enddo - - call writediagfi(ngrid,'qsatused','saturation mixing ratio at surface','',2,qsat) - call writediagfi(ngrid,'psat','saturation pressure at surface','',2,psat) ! coefficients for q @@ -729,5 +719,4 @@ if (tracer) then call writediagfi(ngrid,'dqevap','evaporated water vapor specific concentration','s-1',3,pdqevap) - call writediagfi(ngrid,'h2oevap','evaporated water vapor mass','kg/m2',2,zdmassevap) endif endif