Changeset 600 for trunk

Timestamp:

Mar 27, 2012, 5:02:10 PM (13 years ago)

Author:

jleconte

Message:

• Corrects the computation of planck function at the surface in sfluxi

so that its integral is equal to sigma Tsurf^4.

• This ensure that no flux is lost due to:

-truncation of the planck function at high/low wavenumber
-numerical error during first spectral computation of the planck function
-discrepancy between Tsurf and NTS/NTfac in sfluxi

• OLR now equal to LW net heating/cooling at equilibrium!
• As much as possible, only the value of the stephan boltzmann constant defined in racommon_h (and the

corresponding variable, sigma) should be used. Now done in physics, vdifc and turbdiff.

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/phystd/callcorrk.F90 (modified) (2 diffs)
• libf/phystd/comcstfi.h (modified) (1 diff)
• libf/phystd/physiq.F90 (modified) (5 diffs)
• libf/phystd/sfluxi.F (modified) (3 diffs)
• libf/phystd/turbdiff.F90 (modified) (2 diffs)
• libf/phystd/vdifc.F (modified) (2 diffs)

trunk/LMDZ.GENERIC/README

@@ -640 +640 @@
 - Added .def files for a typical 1d earth case in deftank (dry case for the moment)
 - Corrects a bug on potential temperature calculation in physic 
+
+== 27/03/2012 == JL
+- Corrects the computation of planck function at the surface in sfluxi
+     so that its integral is equal to sigma Tsurf^4.
+- This ensure that no flux is lost due to:
+     -truncation of the planck function at high/low wavenumber
+     -numerical error during first spectral computation of the planck function
+     -discrepancy between Tsurf and NTS/NTfac in sfluxi
+- OLR now equal to LW net heating/cooling at equilibrium!
+- As much as possible, only the value of the stephan boltzmann constant defined in racommon_h (and the
+corresponding variable, sigma) should be used. Now done in physics, vdifc and turbdiff.
+  

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

@@ -702 +702 @@
       END DO
       
-      plevrad(1) = 0
+      plevrad(1) = 0.
       plevrad(2) = max(pgasmin,0.0001*plevrad(3))
 
@@ -816 +816 @@
          fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
 
-
          if(fluxtop_dn(ig).lt.0.0)then
             print*,'Achtung! fluxtop_dn has lost the plot!'

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

@@ -4 +4 @@
       common/comcstfi/pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
       common/comcstfi/avocado!,molrad,visc
-
+      
       real pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
       real avocado!,molrad,visc

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

@@ -10 +10 @@
       use watercommon_h, only : mx_eau_sol, To, RLVTT, mH2O
       use gases_h
+      use radcommon_h, only: sigma
       implicit none
 
@@ -180 +181 @@
       save capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
 
-      real, parameter :: stephan = 5.67e-08 ! Stephan Boltzman constant (should really be elsewhere)
-
 ! Local variables :
 ! -----------------
@@ -551 +550 @@
 
             if(ngrid.eq.1)then
-               qzero1D               = 10000.0
+               qzero1D               = 0.0
                qsurf(1,igcm_h2o_vap) = qzero1D
                do l=1, nlayer
@@ -784 +783 @@
                  fluxsurf_sw(ig) = fluxtop_dn(ig)
                  fluxrad_sky(ig) = fluxtop_dn(ig)*(1.-albedo(ig))
-                 fluxtop_lw(ig)  = emis(ig)*stephan*tsurf(ig)**4
+                 fluxtop_lw(ig)  = emis(ig)*sigma*tsurf(ig)**4
               enddo ! radiation skips the atmosphere entirely
 
@@ -803 +802 @@
         do ig=1,ngrid
            zplanck(ig)=tsurf(ig)*tsurf(ig)
-           zplanck(ig)=emis(ig)*stephan*zplanck(ig)*zplanck(ig)
+           zplanck(ig)=emis(ig)*sigma*zplanck(ig)*zplanck(ig)
            fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
         enddo

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

@@ -5 +5 @@
 
       use radinc_h
-      use radcommon_h, only: planckir, tlimit
+      use radcommon_h, only: planckir, tlimit,sigma
 
       implicit none
@@ -33 +33 @@
 
       real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI)
+      real*8 PLANCKSUM,PLANCKREF
 
 
@@ -72 +73 @@
 !      NTT   = TTOP *10.0D0-499
 
+!JL12 corrects the surface planck function so that its integral is equal to sigma Tsurf^4
+!JL12   this ensure that no flux is lost due to:
+!JL12          -truncation of the planck function at high/low wavenumber
+!JL12          -numerical error during first spectral integration
+!JL12          -discrepancy between Tsurf and NTS/NTfac
+      PLANCKSUM=0.d0
+      PLANCKREF=TSURF*TSURF
+      PLANCKREF=sigma*PLANCKREF*PLANCKREF
+      DO NW=1,L_NSPECTI
+         PLANCKSUM=PLANCKSUM+PLANCKIR(NW,NTS)*DWNI(NW)
+      ENDDO
+      PLANCKSUM=PLANCKREF/(PLANCKSUM*Pi)
+!JL12
+
       DO 501 NW=1,L_NSPECTI
 
 C       SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS
-        BSURF = (1.-RSFI)*PLANCKIR(NW,NTS) ! interpolate for accuracy??
+        BSURF = (1.-RSFI)*PLANCKIR(NW,NTS)*PLANCKSUM !JL12 plancksum see above
         PLTOP = PLANCKIR(NW,NTT)
 

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

@@ -8 +8 @@
 
       use watercommon_h, only : RLVTT, To, RCPD, mx_eau_sol
+      use radcommon_h, only : sigma
 
       implicit none
@@ -123 +124 @@
       save ivap, iliq, iliq_surf,iice_surf
 
-      real, parameter :: sigma=5.67e-8
       real, parameter :: karman=0.4
       real cd0, roughratio

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

@@ -8 +8 @@
 
       use watercommon_h, only : RLVTT, To, RCPD, mx_eau_sol
+      use radcommon_h, only : sigma
 
       implicit none
@@ -134 +135 @@
       save ivap, iice
 
-      real, parameter :: sigma=5.67e-8
       real, parameter :: karman=0.4
       real cd0, roughratio