Changeset 731

Timestamp:

Jul 19, 2012, 12:04:02 PM (12 years ago)

Author:

jleconte

Message:

• Corrected precipitation evaporation scheme + snow fall (JL+BC)
• Commented some unecessary writediagfi calls that where used for test.

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/phystd/physiq.F90 (modified) (4 diffs)
• libf/phystd/rain.F90 (modified) (9 diffs)
• libf/phystd/turbdiff.F90 (modified) (2 diffs)

trunk/LMDZ.GENERIC/README

@@ -738 +738 @@
 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.

trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

@@ -1644 +1644 @@
         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 +1659 @@
           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 +1672 @@
              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 +1709 @@
 
           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

trunk/LMDZ.GENERIC/libf/phystd/rain.F90

@@ -102 +102 @@
 
 !     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 +200 @@
       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 +215 @@
 
                   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 +271 @@
          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 +281 @@
                   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 +299 @@
                      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 +321 @@
                      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 +334 @@
  !precip scheme modified from Boucher 95
            CASE(4)
-
-            !recalculate liquid water particle radii
-            call h2o_cloudrad(ql,reffh2oliq,reffh2oice)
 
             DO n = 1, ninter
@@ -347 +344 @@
                      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)
 

trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90

@@ -517 +517 @@
                 ! 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 +719 @@
          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