Changeset 622

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

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/phystd/callcorrk.F90 (modified) (1 diff)
• libf/phystd/inifis.F (modified) (4 diffs)
• libf/phystd/physiq.F90 (modified) (9 diffs)
• libf/phystd/rain.F90 (modified) (5 diffs)

trunk/LMDZ.GENERIC/README

@@ -652 +652 @@
 corresponding variable, sigma) should be used. Now done in physics, vdifc and turbdiff.
   
+== 16/04/2012 == JL
+- Added consistency checks for calculations including water and global1d+diurnal.
+

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

@@ -274 +274 @@
            call getin("global1d",global1d)
            write(*,*) "global1d = ",global1d
+           ! Test of incompatibility:
+           ! if global1d is true, there should not be any diurnal cycle
+           if (global1d.and.diurnal) then
+            print*,'if global1d is true, diurnal must be set to false'
+            stop
+           endif
+
            if (global1d) then
              PRINT *,'Solar Zenith angle (deg.) ?'

trunk/LMDZ.GENERIC/libf/phystd/inifis.F

@@ -438 +438 @@
          endif
 
-! Test of incompatibility:
-
          write(*,*) "Gravitationnal sedimentation ?"
          sedimentation=.true. ! default value
          call getin("sedimentation",sedimentation)
          write(*,*) " sedimentation = ",sedimentation
-
-
-! Test of incompatibility:
 
          write(*,*) "Compute water cycle ?"
@@ -453 +448 @@
          write(*,*) " water = ",water
          
+! Test of incompatibility:
+! if water is true, there should be at least a tracer
+         if (water.and.(.not.tracer)) then
+           print*,'if water is ON, tracer must be ON too!'
+           stop
+         endif
+
          write(*,*) "Include water condensation ?"
          watercond=.false. ! default value
@@ -458 +460 @@
          write(*,*) " watercond = ",watercond
 
+! Test of incompatibility:
+! if watercond is used, then water should be used too
+         if (watercond.and.(.not.water)) then
+           print*,'if watercond is used, water should be used too'
+           stop
+         endif
+
          write(*,*) "Include water precipitation ?"
          waterrain=.false. ! default value
@@ -497 +506 @@
          call getin("Tsaldiff",Tsaldiff)
          write(*,*) " Tsaldiff = ",Tsaldiff
-
-! Test of incompatibility:
-! if watercond is used, then water should be used too
-
-         if (watercond.and.(.not.watercond)) then
-           print*,'if watercond is used, water should be used too'
-           stop
-         endif
 
          write(*,*) "Does user want to force cpp and mugaz?"

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

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

@@ -51 +51 @@
       PARAMETER (seuil_neb=0.001)
 
-!      REAL ct                           ! Inverse of cloud precipitation time
-!      PARAMETER (ct=1./1800.)
-!      PARAMETER (ct=1./1849.479)
-
-      REAL cl                           ! Precipitation threshold
-      PARAMETER (cl=2.0e-4)
+      REAL,PARAMETER :: cl=2.0e-4                          ! Precipitation threshold
+      REAL,PARAMETER :: ct=1./1800.                        ! Precipitation rate
 
       INTEGER ninter
@@ -62 +58 @@
 
       logical simple                    ! Use very simple Emanuel scheme?
-      parameter(simple=.false.)
+      parameter(simple=.true.)
 
       logical evap_prec                 ! Does the rain evaporate?
@@ -120 +116 @@
          PRINT*, 'in rain.F, ninter=', ninter
          PRINT*, 'in rain.F, evap_prec=', evap_prec
-
-         !print*,ptimestep
-         !print*,1./ct
-         !if(.not.simple)then
-         !   IF (ABS(ptimestep-1./ct).GT.0.001) THEN
-         !      PRINT*, 'Must talk to Laurent Li!!!'
-         !      call abort
-         !   ENDIF
-         !endif
 
          firstcall = .false.
@@ -245 +232 @@
                   zneb(i)  = MAX(rneb(i,k), seuil_neb)
                   IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
-                     d_ql(i,k) = -MAX((ql(i,k)-lconvert),0.0)
+                     d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
                      zrfl(i)   = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
                   ENDIF
@@ -268 +255 @@
                DO i = 1, ngridmx
                   IF (rneb(i,k).GT.0.0) THEN
-                     zchau(i) = (1./FLOAT(ninter)) * zoliq(i)      &
+                     zchau(i) = (ct*ptimestep/FLOAT(ninter)) * zoliq(i)      &
                           * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)**2)) * zfrac(i)
-                     ! warning: this may give dodgy results for physics calls .ne. 48 per day...
 
                      ! this is the ONLY place where zneb, ct and cl are used