Changeset 2283 for trunk/LMDZ.GENERIC/libf/phystd/setspi.F90

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

File:

• trunk/LMDZ.GENERIC/libf/phystd/setspi.F90 (modified) (8 diffs)

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:'