Changeset 2582 for trunk/LMDZ.GENERIC/libf/phystd/tpindex.F

Timestamp:

Nov 9, 2021, 12:56:40 PM (3 years ago)

Author:

emillour

Message:

Generic GCM:
Fixed bug in tpindex (for low temperatures, between first and second
reference temperatures, temperature was wrongly set to tref(1)).
The input temperature was also allowed to be modified by the routine,
which is probably not a good idea and no longer the case.
Took this opportunity to turn tpindex into a module.
GM+EM

File:

• trunk/LMDZ.GENERIC/libf/phystd/tpindex.F (modified) (4 diffs)

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

@@ +1 @@
+      module tpindex_mod
+      
+      implicit none
+      
+      contains
+      
       subroutine tpindex(pw,tw,qvar,pref,tref,wrefvar,LCOEF,MT,MP,
      &     NVAR,wratio)
@@ -34 +40 @@
 !     PW                 - The pressure to interpolate to
 !     TW                 - The temperature to interpolate to
-!     Pref(NP)           - The pressure grid array
-!     Tref(NT)           - The temperature grid array
+!     Qvar               - The abundance to interpolate to
+!     Pref               - The tabulated (log) pressure grid
+!     Tref               - The tabulated temperature grid
+!     WREFVAR            - The tabulated abundance grid
 !   
 !     Outputs
 !     -------
-!     TI                 - Interpolation term (pressure)
-!     UI                 - Interpolation term (temperature)
 !     MT                 - Temperature index (bottom left temperature)
 !                          of bounding box
 !     MP                 - Pressure index (bottom left pressure)
 !                          of bounding box
+!     NVAR               -  index (bottom left )
+!                          of bounding box
+!     LCOEF(4)           - Interpolation coefficients
+!     WRATIO             - weight for the interpolation of abundance
 !
 !     Authors
@@ -52 +62 @@
 !==================================================================
 
-      use radinc_h
+      use radinc_h, only: L_PINT, L_NTREF, L_REFVAR
 
       implicit none
 
-      real*8 Tref(L_NTREF)
-      real*8 pref(L_PINT)
-      real*8 wrefvar(L_REFVAR)
+      REAL,INTENT(IN) :: Tref(L_NTREF)
+      REAL,INTENT(IN) :: pref(L_PINT)
+      REAL,INTENT(IN) :: wrefvar(L_REFVAR)
+      REAL,INTENT(IN) :: PW, TW, Qvar
 
-      integer MT, MP, N, M, NP, NVAR
-      real*8  PW, TW, Qvar, wratio
-      real*8  PWL, LCOEF(4), T, U
+      INTEGER,INTENT(OUT) :: MT, MP, NVAR
+      REAL,INTENT(OUT) :: LCOEF(4), wratio
+      
+      INTEGER :: N
+      REAL :: PWL ! local value for log10(PW)
+      REAL :: TWL ! local value for TW
+      REAL :: T ! linear interpolation weight along log(pressure)
+      REAL :: U ! linear interpolation weight along temperature
 
-C======================================================================C
- 
+!     1. For the interpolation along temperature 
 !     Get the upper and lower temperature grid indicies that bound the
-!     requested temperature. If the requested temperature is outside
-!     the T-grid, set up to extrapolate from the appropriate end.
-!     TW : temperature to be interpolated
-!     TREF : grid array
-!     MT : index of TREF for bounding new temperature
-!     U : new index (real) for temperature interpolated
+!     requested temperature. If the sought temperature is below
+!     the T-grid minimum, assume that lower value. If the sought teperature
+!     is above the T-grid maximum, extrapolate.
+!     TW : temperature to interpolate to
+!     TREF : reference temperature array
+!     MT : index of TREF for (lower) bounding temperature
+!     U : weight for interpolation in temperature
 
-      IF(TW.LE.TREF(1)) THEN
+      TWL = TW
+
+      IF(TWL.LT.TREF(1)) THEN
+        ! Below lowest tabulated temperature
         MT = 1
-        IF (TW.LT.TREF(1)) THEN
-         write(*,*) 'tpindex: Caution! Temperature of upper levels lower 
+        TWL=TREF(1)*1.00
+        write(*,*) 'tpindex: Caution! Temperature of upper levels lower 
      $ than ref temperature for k-coef: k-coeff fixed for upper levels'
-         write(*,*) "         TW=",TW
-         write(*,*) "         TREF(1)=",TREF(1)
-        ENDIF
-      ELSE
+        write(*,*) "         TW=",TWL
+        write(*,*) "         TREF(1)=",TREF(1)
+      ELSEIF(TWL.GE.TREF(L_NTREF)) THEN
+       ! Above highest tabulated temperature
+       MT = L_NTREF-1
+      ELSE 
+        ! Regular case, find encompassing tabulated temperatures
         do n=1,L_NTREF-1
-          if(tw.gt.Tref(n) .and. TW.LE.TREF(N+1)) then
+          IF (TWL.GE.TREF(N) .and. TWL.LT.TREF(N+1)) THEN
             MT = n
-            goto 10
-          end if
-        end do
-
-        MT = L_NTREF-1
+            EXIT
+          END IF
+        END DO
       
-   10   continue
       END IF
 
-      !TB15 : case low temp : MT=1: fixed TW right above tref(1)
-      IF (MT.eq.1) THEN
-         TW=tref(1)*1.00
-!         write(*,*) 'tpindex: Caution! Temperature of upper levels lower 
-!     $than ref temperature for k-coef: k-coeff fixed for upper levels'
-!         write(*,*) "         TW=",TW
-!         write(*,*) "         TREF(1)=",TREF(1)
-      ENDIF
+!     weight for the temperature interpolation
+      U = (TWL-TREF(MT))/(TREF(MT+1)-TREF(MT))
 
-      U = (TW-TREF(MT))/(TREF(MT+1)-TREF(MT))
+!     2. For interpolation along pressure 
+!     Get the upper and lower pressure grid indicies that bound the
+!     requested pressure. If the sought pressure is below
+!     the P-grid minimum, assume that lower value. If the sought pressure
+!     is above the P-grid maximum, extrapolate.
 
-!     Get the upper and lower pressure grid indicies that bound the
-!     requested pressure. If the requested pressure is outside
-!     the P-grid, set up to extrapolate from the appropriate end.
+      PWL = log10(PW)
 
-      pwl = log10(pw)
+      IF(PWL.LT.PREF(1)) THEN
+        ! Below lowest tabulated pressure
+        MP = 1
+        PWL=Pref(1)*1.00
+        write(*,*) 'tpindex: Caution! Pressure of upper levels lower 
+     $than ref pressure for k-coef: k-coeff fixed for upper levels'
+        write(*,*) "         PWL=",PWL
+        write(*,*) "         PREF(1)=",PREF(1)
+      ELSEIF(PWL.GE.PREF(L_PINT)) THEN
+        ! Above highest tabulated pressure
+        MP = L_PINT-1
+      ELSE 
+        ! Regular case, find encompassing tabulated pressures
+        do n=1,L_PINT-1
+          IF (PWL.GE.PREF(N) .and. PWL.LT.PREF(N+1)) THEN
+            MP = n
+            EXIT
+          END IF
+        END DO
+      
+      END IF
 
-      do n=2,L_PINT-1
-        if(pwl.le.Pref(n)) then
-          MP = n-1
-          goto 20
-        end if
-      end do
-
-      MP = L_PINT-1
-
-   20 continue
-      
-      !TB15 : case low pressure : n=2 : fixed pwl, right above pref(1)
-      IF (MP.eq.1) THEN
-        IF (PWL.LT.PREF(1)) THEN
-         write(*,*) 'tpindex: Caution! Pressure of upper levels lower 
-     $than ref pressure for k-coef: k-coeff fixed for upper levels'
-         write(*,*) "         PWL=",PWL
-         write(*,*) "         PREF(1)=",PREF(1)
-        ENDIF
-        PWL=Pref(1)*1.00
-      ENDIF
-
-!     interpolated pressure
+!     weight for the pressure interpolation
       T = (PWL-PREF(MP))/(PREF(MP+1)-PREF(MP))
 
-!  Fill in the interpolation coefficients
+!  Build the interpolation coefficients (bilinear in temperature-log(pressure))
       LCOEF(1) = (1.0-T)*(1.0-U)
       LCOEF(2) = T*(1.0-U)
@@ -143 +156 @@
       LCOEF(4) = (1.0-T)*U
 
+!  3. For interpolation along abundances 
 !  Get the indicies for abundance of the varying species. There are 10 sets of 
 !  k-coefficients with differing amounts of variable vs. constant gas.
 
-      IF(QVAR.le.WREFVAR(1)) then
+      IF(QVAR.LE.WREFVAR(1)) THEN
+        ! Below lowest tabulated abundance
         NVAR   = 1
         WRATIO = 0.0D0      ! put all the weight on the first point
-      ELSEIF(QVAR.ge.WREFVAR(L_REFVAR)) then
+      ELSEIF(QVAR.GE.WREFVAR(L_REFVAR)) THEN
+        ! Above highest tabulated abundance
         NVAR   = L_REFVAR-1 ! TB16 in order to not oversize NVAr when doing
                                   !NVAR+1
         WRATIO = 1.00D0     ! put all the weight on the last point
       ELSE
+        ! Regular case, find encompassing tabulated abundances
+        ! and compute corresponding weight
         DO N=2,L_REFVAR
-          IF(QVAR.GE.WREFVAR(N-1) .and. QVAR.lt.WREFVAR(N)) then
+          IF(QVAR.GE.WREFVAR(N-1) .and. QVAR.LT.WREFVAR(N)) THEN
             NVAR   = N-1
             WRATIO = (QVAR - WREFVAR(N-1))/(WREFVAR(N) - WREFVAR(N-1))
-            GOTO 30
+            EXIT
           END IF
         END DO
       END IF
 
-   30 CONTINUE
-
-      return
-      end
+      end subroutine tpindex
+      
+      end module tpindex_mod
+