Changeset 2283 for trunk/LMDZ.GENERIC

Timestamp:

Apr 9, 2020, 6:34:53 PM (5 years ago)

Author:

jvatant

Message:

Set the temperature boundaries and step for Planck function integration as input in callphys.def, for more flexibility.

+ User can now set them by tplanckmin, tplanckmax and dtplanck
+ Default values are a wide range 30-1500K
+ Add a sanity check in callcorrk instead of leaving out-of-bounds planckir indexes.

--JVO

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/phystd/callcorrk.F90 (modified) (3 diffs)
• libf/phystd/callkeys_mod.F90 (modified) (1 diff)
• libf/phystd/gfluxi.F (modified) (2 diffs)
• libf/phystd/inifis_mod.F90 (modified) (2 diffs)
• libf/phystd/radcommon_h.F90 (modified) (2 diffs)
• libf/phystd/radinc_h.F90 (modified) (2 diffs)
• libf/phystd/setspi.F90 (modified) (8 diffs)
• libf/phystd/sfluxi.F (modified) (1 diff)

trunk/LMDZ.GENERIC/README

@@ -1525 +1525 @@
 - Get rid of "zerophys"; modern Fortran syntax is better to initialize an
   array.
+
+== 09/04/2020 (r2283) == JVO
++ Set the temperature boundaries and step for Planck function
+  integration as input in callphys.def, for more flexibility.
+  + User can now set them by tplanckmin, tplanckmax and dtplanck
+  + Default values are a wide range 30-1500K
+  + Add a sanity check in callcorrk instead of leaving out-of-bounds planckir indexes.
+
+

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

@@ -32 +32 @@
       use comcstfi_mod, only: pi, mugaz, cpp
       use callkeys_mod, only: varactive,diurnal,tracer,water,varfixed,satval,diagdtau,    &
-                              kastprof,strictboundcorrk,specOLR,CLFvarying
+                              kastprof,strictboundcorrk,specOLR,CLFvarying,               &
+                              tplanckmin,tplanckmax
       use optcv_mod, only: optcv
       use optci_mod, only: optci
@@ -739 +740 @@
 
       ! Test for out-of-bounds temperature.
+      ! -- JVO 20 : Also add a sanity test checking that tlevrad is
+      !             within Planck function temperature boundaries,
+      !             which would cause gfluxi/sfluxi to crash.
       do k=1,L_LEVELS
+
          if(tlevrad(k).lt.tgasmin)then
             print*,'Minimum temperature is outside the radiative'
@@ -773 +778 @@
             endif
          endif
+
+         if (tlevrad(k).lt.tplanckmin) then
+            print*,'Minimum temperature is outside the boundaries for'
+            print*,'Planck function integration set in callphys.def, aborting.'
+            print*,"k=",k," tlevrad(k)=",tlevrad(k)
+            print*,"tplanckmin=",tplanckmin
+            message="Minimum temperature outside Planck function bounds - Change tplanckmin in callphys.def"
+            call abort_physic(subname,message,1)
+          else if (tlevrad(k).gt.tplanckmax) then
+            print*,'Maximum temperature is outside the boundaries for'
+            print*,'Planck function integration set in callphys.def, aborting.'
+            print*,"k=",k," tlevrad(k)=",tlevrad(k)
+            print*,"tplanckmax=",tplanckmax
+            message="Maximum temperature outside Planck function bounds - Change tplanckmax in callphys.def"
+            call abort_physic(subname,message,1)
+          endif
+
       enddo
+
       do k=1,L_NLAYRAD+1
          if(tmid(k).lt.tgasmin)then

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

@@ -69 +69 @@
 !$OMP THREADPRIVATE(versH2H2cia,iddist,iaervar,iradia,startype)
 
+      real,save :: tplanckmin
+      real,save :: tplanckmax
+      real,save :: dtplanck
+!$OMP THREADPRIVATE(tplanckmin,tplanckmax,dtplanck)
       real,save :: topdustref
       real,save :: Nmix_co2

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

@@ -85 +85 @@
          LAMDA(L) = ALPHA(L)*(1.0D0-W0(L)*COSBAR(L))/UBARI
          
-         NT    = int(TLEV(2*L)*NTfac)   - NTstar+1
-         NT2   = int(TLEV(2*L+2)*NTfac) - NTstar+1
+         NT    = int(TLEV(2*L)*NTfac)   - NTstart+1
+         NT2   = int(TLEV(2*L+2)*NTfac) - NTstart+1
          
 ! AB : PLANCKIR(NW,NT) is replaced by P1, the linear interpolation result for a temperature NT
@@ -112 +112 @@
       ! -- same results for most thin atmospheres
       ! -- and stabilizes integrations
-      NT    = int(TLEV(2*L+1)*NTfac) - NTstar+1
+      NT    = int(TLEV(2*L+1)*NTfac) - NTstart+1
       !! For deep, opaque, thick first layers (e.g. Saturn)
       !! what is below works much better, not unstable, ...
       !! ... and actually fully accurate because 1st layer temp (JL) 
-      !NT    = int(TLEV(2*L)*NTfac) - NTstar+1
+      !NT    = int(TLEV(2*L)*NTfac) - NTstart+1
       !! (or this one yields same results
-      !NT    = int( (TLEV(2*L)+TLEV(2*L+1))*0.5*NTfac ) - NTstar+1
-      
-      NT2   = int(TLEV(2*L)*NTfac)   - NTstar+1
+      !NT    = int( (TLEV(2*L)+TLEV(2*L+1))*0.5*NTfac ) - NTstart+1
+      
+      NT2   = int(TLEV(2*L)*NTfac)   - NTstart+1
       
 ! AB : PLANCKIR(NW,NT) is replaced by P1, the linear interpolation result for a temperature NT

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

@@ -233 +233 @@
      write(*,*) "strictboundcorrk = ",strictboundcorrk
 
+     write(*,*) "Minimum atmospheric temperature for Planck function integration ?"
+     tplanckmin=30.0 ! default value
+     call getin_p("tplanckmin",tplanckmin)
+     write(*,*) " tplanckmin = ",tplanckmin
+ 
+     write(*,*) "Maximum atmospheric temperature for Planck function integration ?"
+     tplanckmax=1500.0 ! default value
+     call getin_p("tplanckmax",tplanckmax)
+     write(*,*) " tplanckmax = ",tplanckmax
+ 
+     write(*,*) "Temperature step for Planck function integration ?"
+     dtplanck=0.1 ! default value
+     call getin_p("dtplanck",dtplanck)
+     write(*,*) " dtplanck = ",dtplanck
+ 
      write(*,*) "call gaseous absorption in the visible bands?", &
                     "(matters only if callrad=T)"
@@ -809 +824 @@
 #endif
   ! initialize variables in radinc_h
-  call ini_radinc_h(nlayer)
+  call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
  
   ! allocate "comsoil_h" arrays

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

@@ -1 +1 @@
 module radcommon_h
-      use radinc_h, only: L_NSPECTI, L_NSPECTV, NTstar, NTstop, &
+      use radinc_h, only: L_NSPECTI, L_NSPECTV, NTstart, NTstop, &
                           naerkind, nsizemax
       implicit none
@@ -123 +123 @@
       REAL,SAVE :: tstellar ! Stellar brightness temperature (SW)
 
-      real*8,save :: planckir(L_NSPECTI,NTstop-NTstar+1)
+      REAL*8, DIMENSION(:,:), ALLOCATABLE, SAVE :: planckir
 
       real*8,save :: PTOP

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

@@ -72 +72 @@
 
       ! For Planck function integration: 
-      ! equivalent temperatures are 1/NTfac of these values
-      integer, parameter :: NTstar = 500
-      integer, parameter :: NTstop = 15000 ! new default for all non hot Jupiter runs
-      real*8, parameter :: NTfac = 1.0D+1  
-      !integer, parameter :: NTstar = 1000
-      !integer, parameter :: NTstop = 25000
-      !real*8,parameter :: NTfac = 5.0D+1    
-      !integer, parameter :: NTstar = 2000
-      !integer, parameter :: NTstop = 50000
-      !real*8,parameter :: NTfac = 1.0D+2    
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      ! Integration boundary temperatures are NTstart/NTfac and Ntstop/NTfac
+      ! -- JVO 20 : Now read boundary T and integration dT as inputs in callphys.def
+      !             NTstart, Nstop and NTfac then set by ini_radinc_h
+      !             Smart user can adjust values depending he's running hot or cold atm
+      !             Default is wide range : 30K-1500K, with 0.1K step
+      !             ->  NTstart=300, Nstop=15000, NTfac=10
+      integer :: NTstart, NTstop
+      real*8  :: NTfac
 
       ! Maximum number of grain size classes for aerosol convolution:
@@ -96 +95 @@
 contains
 
-  subroutine ini_radinc_h(nbp_lev)
+  subroutine ini_radinc_h(nbp_lev,tplanckmin,tplanckmax,dtplanck)
   ! Initialize module variables
   implicit none
   integer,intent(in) :: nbp_lev
+  real*8, intent(in) :: tplanckmin
+  real*8, intent(in) :: tplanckmax
+  real*8, intent(in) :: dtplanck
   
   L_NLAYRAD = nbp_lev
-  L_LEVELS = 2*(nbp_lev-1)+3
+  L_LEVELS  = 2*(nbp_lev-1)+3
   L_NLEVRAD = nbp_lev+1
-  
+
+  NTfac   = 1.D0 / dtplanck
+  NTstart = int(tplanckmin * NTfac)
+  NTstop  = int(tplanckmax * NTfac)
+ 
   end subroutine
 

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

@@ -22 +22 @@
 !==================================================================
 
-      use radinc_h,    only: L_NSPECTI,corrkdir,banddir,NTstar,NTstop,NTfac
+      use radinc_h,    only: L_NSPECTI,corrkdir,banddir,NTstart,NTstop,NTfac
       use radcommon_h, only: BWNI,BLAMI,WNOI,DWNI,WAVEI,planckir,sigma
       use datafile_mod, only: datadir
@@ -149 +149 @@
       print*,''
       print*,'setspi: Current Planck integration range:'
-      print*,'T = ',dble(NTstar)/NTfac, ' to ',dble(NTstop)/NTfac,' K.'
+      print*,'T = ',dble(NTstart)/NTfac, ' to ',dble(NTstop)/NTfac,' K.'
+
+      IF(.NOT.ALLOCATED(planckir)) ALLOCATE(planckir(L_NSPECTI,NTstop-NTstart+1))
 
       do NW=1,L_NSPECTI
@@ -156 +158 @@
          bpa = (b+a)/2.0D0
          bma = (b-a)/2.0D0
-         do nt=NTstar,NTstop
+         do nt=NTstart,NTstop
             T   = dble(NT)/NTfac
             ans = 0.0D0
@@ -165 +167 @@
             end do
 
-            planckir(NW,nt-NTstar+1) = ans*bma/(PI*DWNI(NW))
+            planckir(NW,nt-NTstart+1) = ans*bma/(PI*DWNI(NW))
          end do
       end do
@@ -172 +174 @@
       if(forceEC)then
          print*,'setspi: Force F=sigma*eps*T^4 for all values of T!'
-         do nt=NTstar,NTstop
+         do nt=NTstart,NTstop
             plancksum=0.0D0
             T=dble(NT)/NTfac
@@ -178 +180 @@
             do NW=1,L_NSPECTI
                plancksum=plancksum+  &
-                  planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
+                  planckir(NW,nt-NTstart+1)*DWNI(NW)*pi
             end do
 
             do NW=1,L_NSPECTI
-               planckir(NW,nt-NTstar+1)=     &
-                  planckir(NW,nt-NTstar+1)*  &
+               planckir(NW,nt-NTstart+1)=     &
+                  planckir(NW,nt-NTstart+1)*  &
                           sigma*(dble(nt)/NTfac)**4/plancksum
             end do
@@ -192 +194 @@
          ! check energy conservation at lower temperature boundary
          plancksum=0.0D0
-         nt=NTstar
+         nt=NTstart
          do NW=1,L_NSPECTI
-            plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
+            plancksum=plancksum+planckir(NW,nt-NTstart+1)*DWNI(NW)*pi
          end do
          print*,'setspi: At lower limit:'
@@ -204 +206 @@
          nt=NTstop
          do NW=1,L_NSPECTI
-            plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
+            plancksum=plancksum+planckir(NW,nt-NTstart+1)*DWNI(NW)*pi
          end do
          print*,'setspi: At upper limit:'

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

@@ -67 +67 @@
       TSURF = TLEV(L_LEVELS)
 
-      NTS   = int(TSURF*NTfac)-NTstar+1
-      NTT   = int(TTOP *NTfac)-NTstar+1
+      NTS   = int(TSURF*NTfac)-NTstart+1
+      NTT   = int(TTOP *NTfac)-NTstart+1
 
 !JL12 corrects the surface planck function so that its integral is equal to sigma Tsurf^4