Changeset 538 for trunk/LMDZ.GENERIC/libf

Timestamp:

Feb 17, 2012, 11:54:16 PM (13 years ago)

Author:

rwordsworth

Message:

calc_cpp_mugaz changed to a check only in 3D

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• calc_cpp_mugaz.F90 (modified) (4 diffs)
• callcorrk.F90 (modified) (11 diffs)
• callkeys.h (modified) (2 diffs)
• inifis.F (modified) (4 diffs)
• physiq.F90 (modified) (7 diffs)

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

@@ -4 +4 @@
 !     Purpose
 !     -------
-!     Compute the atmospheric specific heat capacity.
-!     Compute the atmospheric mean molar mass.
+!     Check to see if the atmospheric specific heat capacity and
+!     mean molar mass for the gas mixture defined in gases.def
+!     corresponds to what we're using. If it doesn't, abort run
+!     unless option 'check_cpp_match' is set to false in 
+!     callphys.def.
+!     NOTE: for now, in 1D we do as before. Jeremy, if you're 
+!     re-writing rcm1d you may want to alter this.
 !
 !     Authors
@@ -16 +21 @@
       implicit none
 
+#include "dimensions.h"
+#include "dimphys.h"
 #include "comcstfi.h"
-!#include "callkeys.h"
+#include "callkeys.h"
+
+      real cpp_c
+      real mugaz_c
 
       integer igas
 
-      cpp=0.0
-      mugaz=0.0
+      cpp_c   = 0.0
+      mugaz_c = 0.0
 
       do igas=1,ngasmx
@@ -31 +41 @@
             ! all values at 300 K from Engineering Toolbox
             if(gnom(igas).eq.'CO2')then
-               cpp   = cpp   + 0.846*gfrac(igas)
-               mugaz = mugaz + 44.01*gfrac(igas)
+               cpp_c   = cpp_c   + 0.846*gfrac(igas)
+               mugaz_c = mugaz_c + 44.01*gfrac(igas)
             elseif(gnom(igas).eq.'N2_')then
-               cpp   = cpp   + 1.040*gfrac(igas)
-               mugaz = mugaz + 28.01*gfrac(igas)
+               cpp_c   = cpp_c   + 1.040*gfrac(igas)
+               mugaz_c = mugaz_c + 28.01*gfrac(igas)
             elseif(gnom(igas).eq.'H2_')then
-               cpp   = cpp   + 14.31*gfrac(igas)
-               mugaz = mugaz + 2.01*gfrac(igas)
+               cpp_c   = cpp_c   + 14.31*gfrac(igas)
+               mugaz_c = mugaz_c + 2.01*gfrac(igas)
             elseif(gnom(igas).eq.'H2O')then
-               cpp   = cpp   + 1.864*gfrac(igas)
-               mugaz = mugaz + 18.02*gfrac(igas)
+               cpp_c   = cpp_c   + 1.864*gfrac(igas)
+               mugaz_c = mugaz_c + 18.02*gfrac(igas)
             elseif(gnom(igas).eq.'CH4')then
-               cpp   = cpp   + 2.226*gfrac(igas)
-               mugaz = mugaz + 16.04*gfrac(igas)
+               cpp_c   = cpp_c   + 2.226*gfrac(igas)
+               mugaz_c = mugaz_c + 16.04*gfrac(igas)
             elseif(gnom(igas).eq.'NH3')then
-               cpp   = cpp   + 2.175*gfrac(igas)
-               mugaz = mugaz + 17.03*gfrac(igas)
+               cpp_c   = cpp_c   + 2.175*gfrac(igas)
+               mugaz_c = mugaz_c + 17.03*gfrac(igas)
                print*,'WARNING, cpp for NH3 may be for liquid' 
             else
@@ -57 +67 @@
       enddo
 
-      cpp=1000.0*cpp
+      cpp_c = 1000.0*cpp_c
 
-      print*,'Cp in calc_cpp_mugaz is ',cpp,'J kg^-1 K^-1'
-      print*,'Mg in calc_cpp_mugaz is ',mugaz,'amu'
+      print*,'Cp in calc_cpp_mugaz is ',cpp_c,'J kg^-1 K^-1'
+      print*,'Mg in calc_cpp_mugaz is ',mugaz_c,'amu'
+      print*,'Predefined Cp in physics is ',cpp,'J kg^-1 K^-1'
+      print*,'Predefined Mg in physics is ',mugaz,'amu'
 
-      R = 8.314511E+0 *1000.E+0/mugaz
-      rcp = R/cpp
+      if(ngridmx.eq.1)then
+         print*,'Automatically setting cpp & mugaz to calculated values in calc_cpp_mugaz' 
+         cpp   = cpp_c
+         mugaz = mugaz_c
+         R     = 8.314511E+0 *1000.E+0/mugaz
+         rcp   = R/cpp
+      elseif((cpp.ne.cpp_c) .or. (mugaz.ne.mugaz_c))then
+         if(check_cpp_match)then
+            print*,'Values do not match!'
+            print*,'Either adjust cpp / mugaz via newstart to calculated values,'
+            print*,'or set check_cpp_match to .false. in callphys.def.'
+            stop
+         endif
+      endif
 
       return

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

@@ -4 +4 @@
           dtlw,dtsw,fluxsurf_lw,                               &
           fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,        &
-          OLR_nu,OSR_nu,                                       &
+          OLR_nu,OSR_nu,                                       &
           reffrad,tau_col,cloudfrac,totcloudfrac,              &
           clearsky,firstcall,lastcall)
@@ -195 +195 @@
       qsiaer(:,:,:)=0.0
       giaer(:,:,:) =0.0
+
 
       if(firstcall) then
@@ -264 +265 @@
             stop
          endif
-         
-         OLR_nu=0.
-         OSR_nu=0.
+
+         OLR_nu=0.
+         OSR_nu=0.
 
          firstcall=.false.   
@@ -326 +327 @@
 
             do ig=1,ngrid
-               if(tracer.and.igcm_co2_ice.gt.0)then
+               if(tracer)then
+               !if(tracer.and.igcm_co2_ice.gt.0)then
+
+                  if(igcm_co2_ice.lt.1)then
+                     print*,'Tracers but no CO2 ice still seems to be a problem...'
+                     print*,'Aborting in callcorrk.'
+                     stop
+                  endif
+
                   reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
                        (4*Nmix_co2*pi*rho_co2) )
@@ -556 +565 @@
          end do
          qvar(1)=qvar(2)
-!         qvar(2*nlayermx+1)=qsurf(ig,i_var) !JL12 not very good to compare kg/kg and kg/m2???
 
       elseif(varfixed)then
@@ -585 +593 @@
          Ttemp = tsurf(ig)
          call watersat(Ttemp,ptemp,qsat)
-         
 
          !qvar(2*nlayermx+1)=qsat      ! fully saturated everywhere
          qvar(2*nlayermx+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!  JL: for completely constant water profile uncoment the following line
-!        qvar=0.005
+         !qvar=0.005                   ! completely constant profile (JL)
 
       else
@@ -600 +604 @@
       end if
 
+      if(.not.kastprof)then
       ! IMPORTANT: Now convert from kg/kg to mol/mol
       do k=1,L_LEVELS
          qvar(k) = qvar(k)/epsi
       end do
+      end if
 
 !-----------------------------------------------------------------------
@@ -634 +640 @@
 
          print*,'Recalculating qvar with VARIABLE epsi for kastprof'
-         i_var=igcm_h2o_vap
+         print*,'Assumes that the variable gas is H2O!!!'
+         print*,'Assumes that there is only one tracer'
+         !i_var=igcm_h2o_vap
+         i_var=1
+         if(nqmx.gt.1)then
+            print*,'Need 1 tracer only to run kcm1d.e' 
+            stop
+         endif
          do l=1,nlayer
             vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O
@@ -644 +657 @@
          end do
          qvar(1)=qvar(2)
-         qvar(2*nlayermx+1)=qsurf(ig,i_var)*muvar(ig,1)/mH2O
+
+         print*,'Warning: reducing qvar in callcorrk.'
+         print*,'Temperature profile no longer consistent ', &
+                            'with saturated H2O.'
+         do k=1,L_LEVELS
+            qvar(k) = qvar(k)*satval
+         end do
 
       endif
@@ -766 +785 @@
               qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
               taugsurfi,qvar,muvarrad)
-              
+
          call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
               wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, & 
@@ -859 +878 @@
            close(116)
 
+!          I am useful, please don`t remove me!
 !           if(specOLR)then
 !               open(117,file='OLRnu.out')

trunk/LMDZ.GENERIC/libf/phystd/callkeys.h

@@ -20 +20 @@
      &   , Tstrat                                                       &
      &   , newtonian                                                    &
+     &   , check_cpp_match                                              &
      &   , tau_relax                                                    &
      &   , testradtimes                                                 &
@@ -61 +62 @@
       logical kastprof
       logical newtonian
+      logical check_cpp_match
       logical testradtimes
       logical rayleigh

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

@@ -177 +177 @@
          write(*,*) " enertest = ",enertest
 
+         write(*,*) "Check to see if cpp values used match gases.def ?"
+         check_cpp_match=.true. ! default value
+         call getin("check_cpp_match",check_cpp_match)
+         write(*,*) " check_cpp_match = ",check_cpp_match
+
+
          write(*,*) "call radiative transfer ?"
          callrad=.true. ! default value
@@ -192 +198 @@
          call getin("callgasvis",callgasvis)
          write(*,*) " callgasvis = ",callgasvis
-
+        
          write(*,*) "call continuum opacities in radiative transfer ?",
      &              "(matters only if callrad=T)"
-         Continuum=.false. ! default value
+         Continuum=.true. ! default value
          call getin("Continuum",Continuum)
          write(*,*) " Continuum = ",Continuum
-         
+
+ 
          write(*,*) "call turbulent vertical diffusion ?"
          calldifv=.true. ! default value
@@ -310 +317 @@
          endif
 
-
          write(*,*)"Use Newtonian cooling for radiative transfer?"
          newtonian=.false.
@@ -506 +512 @@
          endif
 
-         mugaz=0.
-         cpp=0.
+         mugaz=8.314*1000./pr
          call su_gases
          call calc_cpp_mugaz

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

@@ -440 +440 @@
 !        read startfi (initial parameters)
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-          print*,'in physiqu i am in phyetat0'
          call phyetat0("startfi.nc",0,0,nsoilmx,nq,           &
                day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
@@ -712 +711 @@
                  call abort
               endif
-              muvar(:,:)=0.0 ! only used for climate evolution for now
-              
-              
+              muvar(:,:)=0.0 ! only used for climate evolution studies (kcm1d) for now
+      
 !             standard callcorrk
               clearsky=.false.
-              call callcorrk(ngrid,nlayer,pq,nq,qsurf,         &
+              call callcorrk(ngrid,nlayer,pq,nq,qsurf,                   &
                       albedo,emis,mu0,pplev,pplay,pt,                    &
                       tsurf,fract,dist_star,aerosol,cpp3D,muvar,         &
                       zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
-                      fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
+                      fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
                       reffrad,tau_col,cloudfrac,totcloudfrac,            &
                       clearsky,firstcall,lastcall)      
-              
 
 !             Option to call scheme once more for clear regions
@@ -732 +729 @@
                  !!! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
                  clearsky=.true.
-                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
+                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,                  &
                       albedo,emis,mu0,pplev,pplay,pt,                      &
                       tsurf,fract,dist_star,aerosol,cpp3D,muvar,           &
@@ -738 +735 @@
                       fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
                       reffrad,tau_col1,cloudfrac,totcloudfrac,             &
-                      clearsky,firstcall,lastcall)
-                 clearsky = .false.  !! just in case.     
+                      clearsky,firstcall,lastcall)
+                 clearsky = .false.  !! just in case.     
 
                  ! Sum the fluxes and heating rates from cloudy/clear cases
@@ -757 +754 @@
                     enddo
 
-                    do nw=1,L_NSPECTV
+                    do nw=1,L_NSPECTV
                        OSR_nu(ig,nw) = ntf*OSR_nu1(ig,nw) + tf*OSR_nu(ig,nw)                   
-                    enddo
-                    do nw=1,L_NSPECTI
+                    enddo
+                    do nw=1,L_NSPECTI
                        OLR_nu(ig,nw) = ntf*OLR_nu1(ig,nw) + tf*OLR_nu(ig,nw)                   
-                    enddo
+                    enddo
 
                  enddo
@@ -1939 +1936 @@
                enddo
 
-
-
             endif
 
@@ -2097 +2092 @@
                 print*,'WARNING: tau_col=',tau_col(ig)
                 print*,'at ig=',ig,'in PHYSIQ'
-                !call abort
              endif
           end do