source: LMDZ6/branches/contrails/libf/phylmd/ecrad/ifsrrtm/srtm_setcoef.F90 @ 5440

Last change on this file since 5440 was 4773, checked in by idelkadi, 12 months ago
  • Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (​https://github.com/lguez/ecrad.git)
  • Implementation of the double call of Ecrad in LMDZ


File size: 8.5 KB
Line 
1SUBROUTINE SRTM_SETCOEF &
2 & ( KIDIA   , KFDIA    , KLEV    ,&
3 &   PAVEL   , PTAVEL   ,&
4 &   PCOLDRY , PWKL     ,&
5 &   KLAYTROP,&
6 &   PCOLCH4  , PCOLCO2 , PCOLH2O , PCOLMOL  , PCOLO2 , PCOLO3 ,&
7 &   PFORFAC , PFORFRAC , KINDFOR , PSELFFAC, PSELFFRAC, KINDSELF ,&
8 &   PFAC00  , PFAC01   , PFAC10  , PFAC11  ,&
9 &   KJP     , KJT      , KJT1    , PRMU0    &
10 & ) 
11
12!     J. Delamere, AER, Inc. (version 2.5, 02/04/01)
13
14!     Modifications:
15!     JJMorcrette 030224   rewritten / adapted to ECMWF F90 system
16!        M.Hamrud      01-Oct-2003 CY28 Cleaning
17!        D.Salmond  31-Oct-2007 Vector version in the style of RRTM from Meteo France & NEC
18
19!     Purpose:  For a given atmosphere, calculate the indices and
20!     fractions related to the pressure and temperature interpolations.
21
22USE PARKIND1 , ONLY : JPIM, JPRB
23USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK
24USE YOESRTWN , ONLY : PREFLOG, TREF
25!!  USE YOESWN  , ONLY : NDBUG
26
27IMPLICIT NONE
28
29!-- Input arguments
30
31INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA, KFDIA
32INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV
33REAL(KIND=JPRB)   ,INTENT(IN)    :: PAVEL(KIDIA:KFDIA,KLEV)
34REAL(KIND=JPRB)   ,INTENT(IN)    :: PTAVEL(KIDIA:KFDIA,KLEV)
35REAL(KIND=JPRB)   ,INTENT(IN)    :: PCOLDRY(KIDIA:KFDIA,KLEV)
36REAL(KIND=JPRB)   ,INTENT(IN)    :: PWKL(KIDIA:KFDIA,35,KLEV)
37INTEGER(KIND=JPIM),INTENT(OUT)   :: KLAYTROP(KIDIA:KFDIA)
38REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLCH4(KIDIA:KFDIA,KLEV)
39REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLCO2(KIDIA:KFDIA,KLEV)
40REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLH2O(KIDIA:KFDIA,KLEV)
41REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLMOL(KIDIA:KFDIA,KLEV)
42REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLO2(KIDIA:KFDIA,KLEV)
43REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCOLO3(KIDIA:KFDIA,KLEV)
44REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFORFAC(KIDIA:KFDIA,KLEV)
45REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFORFRAC(KIDIA:KFDIA,KLEV)
46INTEGER(KIND=JPIM),INTENT(OUT)   :: KINDFOR(KIDIA:KFDIA,KLEV)
47REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSELFFAC(KIDIA:KFDIA,KLEV)
48REAL(KIND=JPRB)   ,INTENT(OUT)   :: PSELFFRAC(KIDIA:KFDIA,KLEV)
49INTEGER(KIND=JPIM),INTENT(OUT)   :: KINDSELF(KIDIA:KFDIA,KLEV)
50REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFAC00(KIDIA:KFDIA,KLEV)
51REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFAC01(KIDIA:KFDIA,KLEV)
52REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFAC10(KIDIA:KFDIA,KLEV)
53REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFAC11(KIDIA:KFDIA,KLEV)
54INTEGER(KIND=JPIM),INTENT(OUT)   :: KJP(KIDIA:KFDIA,KLEV)
55INTEGER(KIND=JPIM),INTENT(OUT)   :: KJT(KIDIA:KFDIA,KLEV)
56INTEGER(KIND=JPIM),INTENT(OUT)   :: KJT1(KIDIA:KFDIA,KLEV)
57REAL(KIND=JPRB)   ,INTENT(IN)    :: PRMU0(KIDIA:KFDIA)
58!-- Output arguments
59
60!-- local integers
61
62INTEGER(KIND=JPIM) :: I_NLAYERS, JK, JL, JP1
63
64!-- local reals
65
66REAL(KIND=JPRB) :: Z_STPFAC, Z_PLOG
67REAL(KIND=JPRB) :: Z_FP, Z_FT, Z_FT1, Z_WATER, Z_SCALEFAC
68REAL(KIND=JPRB) :: Z_FACTOR, Z_CO2REG, Z_COMPFP
69!REAL(KIND=JPRB) :: Z_TBNDFRAC, Z_T0FRAC
70REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
71
72
73IF (LHOOK) CALL DR_HOOK('SRTM_SETCOEF',0,ZHOOK_HANDLE)
74
75Z_STPFAC = 296._JPRB/1013._JPRB
76I_NLAYERS = KLEV
77
78DO JL = KIDIA, KFDIA
79  IF (PRMU0(JL) > 0.0_JPRB) THEN
80    KLAYTROP(JL)  = 0
81  ENDIF
82ENDDO
83
84DO JK = 1, I_NLAYERS
85  DO JL = KIDIA, KFDIA
86    IF (PRMU0(JL) > 0.0_JPRB) THEN
87      !        Find the two reference pressures on either side of the
88      !        layer pressure.  Store them in JP and JP1.  Store in FP the
89      !        fraction of the difference (in ln(pressure)) between these
90      !        two values that the layer pressure lies.
91
92      Z_PLOG = LOG(PAVEL(JL,JK))
93      KJP(JL,JK) = INT(36._JPRB - 5._JPRB*(Z_PLOG+0.04_JPRB))
94      IF (KJP(JL,JK) < 1) THEN
95        KJP(JL,JK) = 1
96      ELSEIF (KJP(JL,JK) > 58) THEN
97        KJP(JL,JK) = 58
98      ENDIF
99      JP1 = KJP(JL,JK) + 1
100      Z_FP = 5. * (PREFLOG(KJP(JL,JK)) - Z_PLOG)
101
102      !        Determine, for each reference pressure (JP and JP1), which
103      !        reference temperature (these are different for each 
104      !        reference pressure) is nearest the layer temperature but does
105      !        not exceed it.  Store these indices in JT and JT1, resp.
106      !        Store in FT (resp. FT1) the fraction of the way between JT
107      !        (JT1) and the next highest reference temperature that the
108      !        layer temperature falls.
109
110      KJT(JL,JK) = INT(3. + (PTAVEL(JL,JK)-TREF(KJP(JL,JK)))/15.)
111      IF (KJT(JL,JK) < 1) THEN
112        KJT(JL,JK) = 1
113      ELSEIF (KJT(JL,JK) > 4) THEN
114        KJT(JL,JK) = 4
115      ENDIF
116      Z_FT = ((PTAVEL(JL,JK)-TREF(KJP(JL,JK)))/15.) - REAL(KJT(JL,JK)-3)
117      KJT1(JL,JK) = INT(3. + (PTAVEL(JL,JK)-TREF(JP1))/15.)
118      IF (KJT1(JL,JK) < 1) THEN
119        KJT1(JL,JK) = 1
120      ELSEIF (KJT1(JL,JK) > 4) THEN
121        KJT1(JL,JK) = 4
122      ENDIF
123      Z_FT1 = ((PTAVEL(JL,JK)-TREF(JP1))/15.) - REAL(KJT1(JL,JK)-3)
124
125      Z_WATER = PWKL(JL,1,JK)/PCOLDRY(JL,JK)
126      Z_SCALEFAC = PAVEL(JL,JK) * Z_STPFAC / PTAVEL(JL,JK)
127
128      !        If the pressure is less than ~100mb, perform a different
129      !        set of species interpolations.
130
131      ! Olivier Marsden (15 June 2017)
132      !IF (Z_PLOG <= 4.56_JPRB) GO TO 5300
133      IF (KJP(JL,JK) >= 13) GO TO 5300
134     
135      KLAYTROP(JL) =  KLAYTROP(JL) + 1
136
137      !        Set up factors needed to separately include the water vapor
138      !        foreign-continuum in the calculation of absorption coefficient.
139
140      PFORFAC(JL,JK) = Z_SCALEFAC / (1.+Z_WATER)
141      Z_FACTOR = (332.0-PTAVEL(JL,JK))/36.0
142      KINDFOR(JL,JK) = MIN(2, MAX(1, INT(Z_FACTOR)))
143      PFORFRAC(JL,JK) = Z_FACTOR - REAL(KINDFOR(JL,JK))
144
145      !        Set up factors needed to separately include the water vapor
146      !        self-continuum in the calculation of absorption coefficient.
147
148      PSELFFAC(JL,JK) = Z_WATER * PFORFAC(JL,JK)
149      Z_FACTOR = (PTAVEL(JL,JK)-188.0)/7.2
150      KINDSELF(JL,JK) = MIN(9, MAX(1, INT(Z_FACTOR)-7))
151      PSELFFRAC(JL,JK) = Z_FACTOR - REAL(KINDSELF(JL,JK) + 7)
152
153      !        Calculate needed column amounts.
154
155      PCOLH2O(JL,JK) = 1.E-20 * PWKL(JL,1,JK)
156      PCOLCO2(JL,JK) = 1.E-20 * PWKL(JL,2,JK)
157      PCOLO3(JL,JK) = 1.E-20 * PWKL(JL,3,JK)
158      !         COLO3(LAY) = 0.
159      !         COLO3(LAY) = colo3(lay)/1.16
160      PCOLCH4(JL,JK) = 1.E-20 * PWKL(JL,6,JK)
161      PCOLO2(JL,JK) = 1.E-20 * PWKL(JL,7,JK)
162      PCOLMOL(JL,JK) = 1.E-20 * PCOLDRY(JL,JK) + PCOLH2O(JL,JK)
163      !         colco2(lay) = 0.
164      !         colo3(lay) = 0.
165      !         colch4(lay) = 0.
166      !         colo2(lay) = 0.
167      !         colmol(lay) = 0.
168      IF (PCOLCO2(JL,JK) == 0.) PCOLCO2(JL,JK) = 1.E-32 * PCOLDRY(JL,JK)
169      IF (PCOLCH4(JL,JK) == 0.) PCOLCH4(JL,JK) = 1.E-32 * PCOLDRY(JL,JK)
170      IF (PCOLO2(JL,JK) == 0.) PCOLO2(JL,JK) = 1.E-32 * PCOLDRY(JL,JK)
171      !        Using E = 1334.2 cm-1.
172      Z_CO2REG = 3.55E-24 * PCOLDRY(JL,JK)
173      GO TO 5400
174
175      !        Above LAYTROP.
1765300  CONTINUE
177
178      !        Set up factors needed to separately include the water vapor
179      !        foreign-continuum in the calculation of absorption coefficient.
180
181      PFORFAC(JL,JK) = Z_SCALEFAC / (1.+Z_WATER)
182      Z_FACTOR = (PTAVEL(JL,JK)-188.0)/36.0
183      KINDFOR(JL,JK) = 3
184      PFORFRAC(JL,JK) = Z_FACTOR - 1.0
185
186      !        Calculate needed column amounts.
187
188      PCOLH2O(JL,JK) = 1.E-20 * PWKL(JL,1,JK)
189      PCOLCO2(JL,JK) = 1.E-20 * PWKL(JL,2,JK)
190      PCOLO3(JL,JK)  = 1.E-20 * PWKL(JL,3,JK)
191      PCOLCH4(JL,JK) = 1.E-20 * PWKL(JL,6,JK)
192      PCOLO2(JL,JK)  = 1.E-20 * PWKL(JL,7,JK)
193      PCOLMOL(JL,JK) = 1.E-20 * PCOLDRY(JL,JK) + PCOLH2O(JL,JK)
194      IF (PCOLCO2(JL,JK) == 0.) PCOLCO2(JL,JK) = 1.E-32 * PCOLDRY(JL,JK)
195      IF (PCOLCH4(JL,JK) == 0.) PCOLCH4(JL,JK) = 1.E-32 * PCOLDRY(JL,JK)
196      IF (PCOLO2(JL,JK) == 0.) PCOLO2(JL,JK)  = 1.E-32 * PCOLDRY(JL,JK)
197      Z_CO2REG = 3.55E-24 * PCOLDRY(JL,JK)
198
199      PSELFFAC(JL,JK) =0.0_JPRB
200      PSELFFRAC(JL,JK)=0.0_JPRB
201      KINDSELF(JL,JK) = 0
202
2035400  CONTINUE
204
205      !        We have now isolated the layer ln pressure and temperature,
206      !        between two reference pressures and two reference temperatures
207      !        (for each reference pressure).  We multiply the pressure
208      !        fraction FP with the appropriate temperature fractions to get
209      !        the factors that will be needed for the interpolation that yields
210      !        the optical depths (performed in routines TAUGBn for band n).
211
212      Z_COMPFP = 1. - Z_FP
213      PFAC10(JL,JK) = Z_COMPFP * Z_FT
214      PFAC00(JL,JK) = Z_COMPFP * (1. - Z_FT)
215      PFAC11(JL,JK) = Z_FP * Z_FT1
216      PFAC01(JL,JK) = Z_FP * (1. - Z_FT1)
217
218    ENDIF
219  ENDDO
220ENDDO
221
222!-----------------------------------------------------------------------
223IF (LHOOK) CALL DR_HOOK('SRTM_SETCOEF',1,ZHOOK_HANDLE)
224
225END SUBROUTINE SRTM_SETCOEF
Note: See TracBrowser for help on using the repository browser.