source: trunk/LMDZ.GENERIC/libf/phystd/tpindex.F @ 2537

Last change on this file since 2537 was 1616, checked in by emillour, 8 years ago

Generic physics: better control of warning messages in tpindex, only write these if strictly out of range.
EM

File size: 5.0 KB
RevLine 
[135]1      subroutine tpindex(pw,tw,qvar,pref,tref,wrefvar,LCOEF,MT,MP,
2     &     NVAR,wratio)
3
4!==================================================================
5!     
6!     Purpose
7!     -------
8!     Interpolate K-coefficients to the given P,T and Qvar values.
9!
10!     Notes
11!     -----
12!     The interpolation is the usual one in two dimensions given
13!     in "Numerical Recipes", where the "X" are P, the "Y" are
14!     T, and the F(X,Y) are the CO2 K-coefficients.
15!
16!     The interpolating box is:
17!
18!           (PL,TU)                        (PR,TU)
19!
20!                          (TW,PW)
21!
22!           
23!           (PL,TL)                        (PR,TL)
24!
25!      PL  - Pressure left
26!      PR  - Pressure right
27!      TL  - Temperature lower
28!      TU  - Temperature upper
29!      PW  - Pressure wanted
30!      TW  - Temperature wanted
31!
32!     Inputs
33!     ------
34!     PW                 - The pressure to interpolate to
35!     TW                 - The temperature to interpolate to
36!     Pref(NP)           - The pressure grid array
37!     Tref(NT)           - The temperature grid array
38!   
39!     Outputs
40!     -------
41!     TI                 - Interpolation term (pressure)
42!     UI                 - Interpolation term (temperature)
43!     MT                 - Temperature index (bottom left temperature)
44!                          of bounding box
45!     MP                 - Pressure index (bottom left pressure)
46!                          of bounding box
47!
48!     Authors
49!     -------
50!     Adapted from the NASA Ames code by R. Wordsworth (2009)
51!     
52!==================================================================
53
54      use radinc_h
55
56      implicit none
57
58      real*8 Tref(L_NTREF)
59      real*8 pref(L_PINT)
60      real*8 wrefvar(L_REFVAR)
61
62      integer MT, MP, N, M, NP, NVAR
63      real*8  PW, TW, Qvar, wratio
64      real*8  PWL, LCOEF(4), T, U
65
66C======================================================================C
67 
68!     Get the upper and lower temperature grid indicies that bound the
69!     requested temperature. If the requested temperature is outside
70!     the T-grid, set up to extrapolate from the appropriate end.
[1600]71!     TW : temperature to be interpolated
72!     TREF : grid array
73!     MT : index of TREF for bounding new temperature
74!     U : new index (real) for temperature interpolated
[135]75
76      IF(TW.LE.TREF(1)) THEN
77        MT = 1
[1616]78        IF (TW.LT.TREF(1)) THEN
79         write(*,*) 'tpindex: Caution! Temperature of upper levels lower
80     $ than ref temperature for k-coef: k-coeff fixed for upper levels'
81         write(*,*) "         TW=",TW
82         write(*,*) "         TREF(1)=",TREF(1)
83        ENDIF
[135]84      ELSE
85        do n=1,L_NTREF-1
86          if(tw.gt.Tref(n) .and. TW.LE.TREF(N+1)) then
87            MT = n
88            goto 10
89          end if
90        end do
91
92        MT = L_NTREF-1
93     
94   10   continue
95      END IF
96
[1600]97      !TB15 : case low temp : MT=1: fixed TW right above tref(1)
98      IF (MT.eq.1) THEN
99         TW=tref(1)*1.00
[1616]100!         write(*,*) 'tpindex: Caution! Temperature of upper levels lower
101!     $than ref temperature for k-coef: k-coeff fixed for upper levels'
102!         write(*,*) "         TW=",TW
103!         write(*,*) "         TREF(1)=",TREF(1)
[1600]104      ENDIF
105
[135]106      U = (TW-TREF(MT))/(TREF(MT+1)-TREF(MT))
107
108!     Get the upper and lower pressure grid indicies that bound the
109!     requested pressure. If the requested pressure is outside
110!     the P-grid, set up to extrapolate from the appropriate end.
111
112      pwl = log10(pw)
113
114      do n=2,L_PINT-1
115        if(pwl.le.Pref(n)) then
116          MP = n-1
117          goto 20
118        end if
119      end do
120
121      MP = L_PINT-1
122
123   20 continue
[1600]124     
125      !TB15 : case low pressure : n=2 : fixed pwl, right above pref(1)
126      IF (MP.eq.1) THEN
[1616]127        IF (PWL.LT.PREF(1)) THEN
[1600]128         write(*,*) 'tpindex: Caution! Pressure of upper levels lower
129     $than ref pressure for k-coef: k-coeff fixed for upper levels'
[1616]130         write(*,*) "         PWL=",PWL
131         write(*,*) "         PREF(1)=",PREF(1)
132        ENDIF
133        PWL=Pref(1)*1.00
[1600]134      ENDIF
[135]135
[1600]136!     interpolated pressure
[135]137      T = (PWL-PREF(MP))/(PREF(MP+1)-PREF(MP))
138
139!  Fill in the interpolation coefficients
140      LCOEF(1) = (1.0-T)*(1.0-U)
141      LCOEF(2) = T*(1.0-U)
142      LCOEF(3) = T*U
143      LCOEF(4) = (1.0-T)*U
144
145!  Get the indicies for abundance of the varying species. There are 10 sets of
146!  k-coefficients with differing amounts of variable vs. constant gas.
147
148      IF(QVAR.le.WREFVAR(1)) then
149        NVAR   = 1
[1600]150        WRATIO = 0.0D0      ! put all the weight on the first point
[135]151      ELSEIF(QVAR.ge.WREFVAR(L_REFVAR)) then
[1600]152        NVAR   = L_REFVAR-1 ! TB16 in order to not oversize NVAr when doing
153                                  !NVAR+1
154        WRATIO = 1.00D0     ! put all the weight on the last point
[135]155      ELSE
156        DO N=2,L_REFVAR
157          IF(QVAR.GE.WREFVAR(N-1) .and. QVAR.lt.WREFVAR(N)) then
158            NVAR   = N-1
159            WRATIO = (QVAR - WREFVAR(N-1))/(WREFVAR(N) - WREFVAR(N-1))
160            GOTO 30
161          END IF
162        END DO
163      END IF
164
165   30 CONTINUE
166
167      return
168      end
Note: See TracBrowser for help on using the repository browser.