Changeset 596 for trunk/LMDZ.GENERIC/libf/phystd

Timestamp:

Mar 26, 2012, 5:44:40 PM (13 years ago)

Author:

jleconte

Message:

• Added double gray case (if graybody=true in callphys.def):
  • opacities are set to a constant value in sugas_corrk.
  • the values are kappa_IR m^{2/kg in the infrared (to be read in callphys.def)}

kappa_VI m^{2/kg in the visible (to be read in callphys.def)}

• Cleaned continuum part in optc*
• Added .def files for a typical 1d earth case in deftank (dry case for the moment)

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• callcorrk.F90 (modified) (1 diff)
• callkeys.h (modified) (4 diffs)
• inifis.F (modified) (2 diffs)
• optci.F90 (modified) (2 diffs)
• optcv.F90 (modified) (2 diffs)
• physiq.F90 (modified) (6 diffs)
• sugas_corrk.F90 (modified) (8 diffs)

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

@@ -341 +341 @@
 
             do ig=1,ngrid
-               if(tracer)then
-               !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

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

@@ -3 +3 @@
 ! symbols '&' in columns 73 and 6
 !
-      COMMON/callkeys/callrad,corrk,calldifv,calladj                    &
+      COMMON/callkeys/callrad,corrk,calldifv,UseTurbDiff,calladj        &
      &   , co2cond,callsoil                                             &
      &   , season,diurnal,tlocked,iradia,lwrite                         &
      &   , iaervar,iddist,topdustref,callstats,calleofdump              &
      &   , enertest                                                     &
-     &   , callgasvis                                                   &
-     &   , Continuum                                                    &
+     &   , callgasvis,Continuum, graybody                               &
      &   , Nmix_co2, Nmix_h2o                                           &
      &   , dusttau                                                      &
@@ -17 +16 @@
      &   , kastprof                                                     &
      &   , noradsurf                                                    &
-     &   , graybody                                                     &
      &   , Tstrat                                                       &
      &   , newtonian                                                    &
@@ -51 +49 @@
      &   , pceil                                                    
 
-      logical callrad,corrk,calldifv,calladj,co2cond,callsoil           &
+      logical callrad,corrk,calldifv,UseTurbDiff                        &
+     &   , calladj,co2cond,callsoil                                     &
      &   , season,diurnal,tlocked,lwrite                                &
      &   , callstats,calleofdump                                        &
-     &   , callgasvis,Continuum
+     &   , callgasvis,Continuum,graybody
 
       logical enertest
@@ -79 +78 @@
       logical CLFvarying
       logical noradsurf
-      logical graybody
 
       integer iddist

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

@@ -210 +210 @@
          write(*,*) " calldifv = ",calldifv
 
+         write(*,*) "use turbdiff instead of vdifc ?"
+         UseTurbDiff=.true. ! default value
+         call getin("UseTurbDiff",UseTurbDiff)
+         write(*,*) " UseTurbDiff = ",UseTurbDiff
+
          write(*,*) "call convective adjustment ?"
          calladj=.true. ! default value
@@ -280 +285 @@
          write(*,*)" kastprof = ",kastprof
 
-         write(*,*)"Uniform absorption coefficient in IR?"
+         write(*,*)"Uniform absorption in radiative transfer?"
          graybody=.false.
          call getin("graybody",graybody)

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

@@ -137 +137 @@
             DCONT = 0.0 ! continuum absorption
 
+            if(Continuum.and.(.not.graybody))then
             ! include continua if necessary
-            wn_cont = dble(wnoi(nw))
-            T_cont  = dble(TMID(k))
-            do igas=1,ngasmx
-
-               if(gfrac(igas).eq.-1)then ! variable
-                  p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
-               else
-                  p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
-               endif
-
-               dtemp=0.0
-               if(gnom(igas).eq.'H2_')then
-                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
-               elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
-                  p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
-                  call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
-
-               endif
-
-               DCONT = DCONT + dtemp
-
-            enddo
-
-
-            DCONT = DCONT*dz(k) 
-            
-            if(.not.Continuum)then
-               DCONT=0.0
+               wn_cont = dble(wnoi(nw))
+               T_cont  = dble(TMID(k))
+               do igas=1,ngasmx
+
+                  if(gfrac(igas).eq.-1)then ! variable
+                     p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
+                  else
+                     p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
+                  endif
+
+                  dtemp=0.0
+                  if(gnom(igas).eq.'H2_')then
+                     call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
+                  elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
+                     p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
+                     call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
+                  endif
+
+                  DCONT = DCONT + dtemp
+               enddo
+
+               DCONT = DCONT*dz(k) 
             endif
 
@@ -213 +208 @@
                
                TAUGAS  = U(k)*ANS
-
-               if(graybody)then
-                  TAUGAS = 0.0
-                  DCONT  = U(k)*3.3e-26
-               endif
 
                TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT

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

@@ -132 +132 @@
             DCONT = 0.0 ! continuum absorption
 
+            if(callgasvis.and.Continuum.and.(.not.graybody))then
             ! include continua if necessary
-            wn_cont = dble(wnov(nw))
-            T_cont  = dble(TMID(k))
-            do igas=1,ngasmx
-
-               if(gfrac(igas).eq.-1)then ! variable
-                  p_cont  = dble(PMID(k)*scalep*QVAR(K)) ! qvar = mol/mol
-               else
-                  p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
-               endif
-
-               dtemp=0.0
-               if(gnom(igas).eq.'H2_')then 
-                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
-               elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
-                  p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
-                  call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
-               endif
-
-               DCONT = DCONT + dtemp
-
-            enddo
-
-            DCONT = DCONT*dz(k) 
-
-            if(.not.(callgasvis.and.Continuum))then
-               DCONT=0.0
+               wn_cont = dble(wnov(nw))
+               T_cont  = dble(TMID(k))
+               do igas=1,ngasmx
+   
+                  if(gfrac(igas).eq.-1)then ! variable
+                     p_cont  = dble(PMID(k)*scalep*QVAR(K)) ! qvar = mol/mol
+                  else
+                     p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
+                  endif
+
+                  dtemp=0.0
+                  if(gnom(igas).eq.'H2_')then 
+                     call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
+                  elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
+                     p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
+                     call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
+                  endif
+
+                  DCONT = DCONT + dtemp
+
+               enddo
+
+               DCONT = DCONT*dz(k) 
             endif
 
@@ -198 +196 @@
 
                TAUGAS          = U(k)*ANS
-
-               if(graybody)then
-                  TAUGAS = 0.0
-                  DCONT  = 0.0
-               endif
-
 
 

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

@@ -315 +315 @@
       real dEtot, dEtots, masse, AtmToSurf_TurbFlux
       real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
-      real dEzRad(ngridmx,nlayermx),dEzdif(ngridmx,nlayermx)
+      real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdif(ngridmx,nlayermx)
       real madjdE(ngridmx), lscaledE(ngridmx)
 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily
@@ -400 +400 @@
       logical ice_update
       save ice_update
-
-!JL12 temporary to test difference in diffusion schemes
-      logical UseTurbDiff
 
 !=======================================================================
@@ -822 +819 @@
             dEtotsSW = 0.0
             dEtotsLW = 0.0
+            dEzRadsw(:,:)=0.
+            dEzRadlw(:,:)=0.
             do ig = 1, ngrid
                do l = 1, nlayer
@@ -827 +826 @@
                   dEtotSW = dEtotSW + cpp*masse*zdtsw(ig,l)*area(ig)
                   dEtotLW = dEtotLW + cpp*masse*zdtlw(ig,l)*area(ig)
+                  dEzRadsw(ig,l)=cpp*masse*zdtsw(ig,l)
+                  dEzRadlw(ig,l)=cpp*masse*zdtlw(ig,l)
                enddo
                dEtotsSW = dEtotsSW + fluxsurf_sw(ig)*(1.-albedo(ig))*area(ig)
@@ -862 +863 @@
 
 !JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff      
-         UseTurbDiff=.true.
          if (UseTurbDiff) then
          
@@ -2008 +2008 @@
         
         if(enertest) then
-          call writediagfi(ngridmx,"dEzdif","atm vdifc energy change","w.m^-2",3,dEzdif)
+          if (calldifv) call writediagfi(ngridmx,"dEzdif","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdif)
+          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)
+          endif 
           if(watercond) then
              call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 

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

@@ -12 +12 @@
 !     -------
 !     Adapted and generalised from the NASA Ames code by Robin Wordsworth (2009)
+!     Added double gray case by Jeremy Leconte (2012)
 !
 !     Summary
@@ -25 +26 @@
       use datafile_mod, only: datadir
       use gases_h
+      use ioipsl_getincom 
       implicit none
 
 #include "callkeys.h"
-
+#include "comcstfi.h"
 !==================================================================
 
@@ -44 +46 @@
 
       real*8 x, xi(4), yi(4), ans, p
+      real kappa_IR, kappa_VI
 
       integer ngas, igas
@@ -53 +56 @@
       file_id='/corrk_data/' // TRIM(corrkdir) // '/Q.dat'
       file_path=TRIM(datadir)//TRIM(file_id)
-
 
       ! check that the file exists
@@ -266 +268 @@
 !-----------------------------------------------------------------------
 
-
-
 !=======================================================================
 !     Get gaseous k-coefficients and interpolate onto finer pressure grid
 
       ! VISIBLE
-      if (callgasvis.and..not.TRIM(corrkdir).eq.'null') then
+      if (callgasvis.and.(.not.TRIM(corrkdir).eq.'null').and.(.not.graybody)) then
          file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat' 
          file_path=TRIM(datadir)//TRIM(file_id)
@@ -288 +288 @@
          read(111,*) gasv8
          close(111)
-         
+         
+      else if (callgasvis.and.graybody) then
+!    constant absorption coefficient in visible
+         write(*,*)"constant absorption coefficient in visible:"
+         kappa_VI=-100000.
+         call getin("kappa_VI",kappa_VI)
+         write(*,*)" kappa_VI = ",kappa_VI
+         kappa_VI=kappa_VI*1.e4* mugaz * 1.672621e-27    ! conversion from m^2/kg to cm^2/molecule         
+         gasv8(:,:,:,:,:)=kappa_VI
       else
          print*,'Visible corrk gaseous absorption is set to zero.'
@@ -295 +303 @@
 
       ! INFRA-RED
-      if (.not.TRIM(corrkdir).eq.'null') then
+      if ((.not.TRIM(corrkdir).eq.'null').and.(.not.graybody)) then
          file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat' 
          file_path=TRIM(datadir)//TRIM(file_id)
@@ -311 +319 @@
          read(111,*) gasi8
          close(111)
-     
-         do nw=1,L_NSPECTI
-            fzeroi(nw) = 0.
-!            do nt=1,L_NTREF
-!               do np=1,L_NPREF
-!                  do nh=1,L_REFVAR
-!                     do ng = 1,L_NGAUSS
-!                        if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
-!                           fzeroi(nw)=fzeroi(nw)+1.
-!                        endif
-!                     end do
-!                  end do
-!               end do
-!            end do
-!            fzeroi(nw)=fzeroi(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
-         end do
-
-         do nw=1,L_NSPECTV
-            fzerov(nw) = 0.
-!            do nt=1,L_NTREF
-!               do np=1,L_NPREF
-!                  do nh=1,L_REFVAR
-!                     do ng = 1,L_NGAUSS
-!                        if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
-!                           fzerov(nw)=fzerov(nw)+1.
-!                        endif
-!                     end do
-!                  end do
-!               end do
-!            end do
-!            fzerov(nw)=fzerov(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
-         end do
-
-
-         do nw=1,L_NSPECTV
-            fzerov(nw) = 0.
-         end do
-
+
+      else if (graybody) then
+!     constant absorption coefficient in IR
+         write(*,*)"constant absorption coefficient in InfraRed:"
+         kappa_IR=-100000.
+         call getin("kappa_IR",kappa_IR)
+         write(*,*)" kappa_IR = ",kappa_IR       
+         kappa_IR=kappa_IR*1.e4* mugaz * 1.672621e-27    ! conversion from m^2/kg to cm^2/molecule       
+         gasi8(:,:,:,:,:)=kappa_IR
       else
          print*,'Infrared corrk gaseous absorption is set to zero.'
          gasi8(:,:,:,:,:)=0.0
       endif
+
+      do nw=1,L_NSPECTI
+         fzeroi(nw) = 0.
+      end do
+
+      do nw=1,L_NSPECTV
+         fzerov(nw) = 0.
+      end do
+
 
 !     Take log10 of the values - this is what we will interpolate.