Changeset 622 for trunk/LMDZ.GENERIC/libf/phystd/physiq.F90

Timestamp:

Apr 16, 2012, 12:22:40 PM (13 years ago)

Author:

jleconte

Message:

• Added consistency checks for calculations including water and global1d+diurnal.
• Corrected small bugs in precipitation scheme

File:

• trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 (modified) (9 diffs)

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

@@ -316 +316 @@
       real dEtot, dEtots, masse, AtmToSurf_TurbFlux
       real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
-      real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdif(ngridmx,nlayermx)
+      real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdiff(ngridmx,nlayermx)
+      real dEdiffs(ngridmx),dEdiff(ngridmx)
       real madjdE(ngridmx), lscaledE(ngridmx)
 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily
@@ -559 +560 @@
 
             do ig=1,ngrid
-               qsurf_hist(ig,igcm_h2o_vap) = qsurf(ig,igcm_h2o_vap)
+!already done above (JL12)
+!               qsurf_hist(ig,igcm_h2o_vap) = qsurf(ig,igcm_h2o_vap)
                if(ice_update)then
                   ice_initial(ig)=qsurf(ig,igcm_h2o_ice)
@@ -877 +879 @@
            do l=1,nlayer
              do ig=1,ngrid 
-               zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
                zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
              enddo
@@ -931 +932 @@
             dEtot=0.0
             dEtots=0.0
-            dEzdif(:,:)=0.
+            dEzdiff(:,:)=0.
             AtmToSurf_TurbFlux=0.
-            
+            dEdiffs(:)=0.
+            dEdiff(:)=0.
+                    
             do ig = 1, ngrid
                do l = 1, nlayer
                   masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                  dEtot = dEtot + cpp*masse*zdtdif(ig,l)*area(ig)
-                  dEzdif (ig,l)= cpp*masse*(zdtdif(ig,l))
-               enddo
-               dEzdif(ig,1)= dEzdif(ig,1)+ sensibFlux(ig)! subtract flux to the ground
-               dEtot = dEtot + sensibFlux(ig)*area(ig)! subtract flux to the ground
-               dEtots = dEtots + capcal(ig)*zdtsdif(ig)*area(ig)-zflubid(ig)*area(ig)-sensibFlux(ig)*area(ig)
+                  dEzdiff (ig,l)= cpp*masse*(zdtdif(ig,l))
+                  dEdiff(ig)=dEdiff(ig)+dEzdiff (ig,l)
+               enddo
+               dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
+               dEdiff(ig)=dEdiff(ig)+ sensibFlux(ig)! subtract flux to the ground
+               dEtot = dEtot + dEdiff(ig)*area(ig)
+               dEdiffs(ig)=capcal(ig)*zdtsdif(ig)-zflubid(ig)-sensibFlux(ig)
+               dEtots = dEtots + dEdiffs(ig)*area(ig)
                AtmToSurf_TurbFlux=AtmToSurf_TurbFlux+sensibFlux(ig)*area(ig)
             enddo
@@ -1342 +1347 @@
                  rainout(ig)             = zdqsrain(ig)+zdqssnow(ig) ! diagnostic   
               enddo
-
-              !-------------------------
-              ! test energy conservation
-              if(enertest)then
-                 dEtot=0.0
-                 do ig = 1, ngrid
-                    do l = 1, nlayer
-                       masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                       dEtot = dEtot + cpp*masse*zdtrain(ig,l)*area(ig)
-                    enddo
-                 enddo
-                 dEtot=dEtot/totarea
-                 print*,'In rain atmospheric energy change       =',dEtot,' W m-2'
-              endif
-              !-------------------------
 
 
@@ -1572 +1562 @@
          ! content of ocean gridpoints back to zero, in order
          ! to avoid rounding errors in vdifc, rain
-         do ig = 1, ngrid
-            do iq = 1, nq
-               if(iq.eq.igcm_h2o_vap .and. rnat(ig).eq.0)then ! if liquid water and terrain = ocean 
-                  qsurf_hist(ig,iq) = qsurf(ig,iq)
-                  !qsurf(ig,iq)      = qcol(ig,iq) 
-                  ! the value of qsurf we choose here makes NO DIFFERENCE TO ANYTHING AT ALL
-               else
-                  qsurf_hist(ig,iq) = qsurf(ig,iq)
-               endif
-            enddo
-         enddo
+         qsurf_hist(:,:) = qsurf(:,:)
 
          if(ice_update)then
@@ -2016 +1996 @@
         
         if(enertest) then
-          if (calldifv) call writediagfi(ngridmx,"dEzdif","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdif)
+          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,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
+          endif
           if (corrk) then
              call writediagfi(ngridmx,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
@@ -2023 +2009 @@
           if(watercond) then
              call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
-             call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE) 
+             call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
+             call writediagfi(ngridmx,"RH","relative humidity"," ",3,RH)
+             call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
           endif 
         endif
@@ -2054 +2042 @@
                           'kg m^-2',2,qcol(1,iq) )
 
-          if(water)then
-             call writediagfi(ngridmx,"H2Omaxcol","max. poss. H2O column","kg m^-2",2,H2Omaxcol)
-          endif
-
           if(watercond.or.CLFvarying)then
-             !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