source: LMDZ6/trunk/libf/phylmdiso/rrtm/rrtm_setcoef_140gp.F90 @ 3927

Last change on this file since 3927 was 3927, checked in by Laurent Fairhead, 3 years ago

Initial import of the physics wih isotopes from Camille Risi
CR
File size: 7.5 KB
Line 
1SUBROUTINE RRTM_SETCOEF_140GP (KLEV,P_COLDRY,P_WKL,&
2 & P_FAC00,P_FAC01,P_FAC10,P_FAC11,P_FORFAC,K_JP,K_JT,K_JT1,&
3 & P_COLH2O,P_COLCO2,P_COLO3,P_COLN2O,P_COLCH4,P_COLO2,P_CO2MULT,&
4 & K_LAYTROP,K_LAYSWTCH,K_LAYLOW,PAVEL,P_TAVEL,P_SELFFAC,P_SELFFRAC,K_INDSELF) 
5
6!     Reformatted for F90 by JJMorcrette, ECMWF, 980714
7
8!     Purpose:  For a given atmosphere, calculate the indices and
9!     fractions related to the pressure and temperature interpolations.
10!     Also calculate the values of the integrated Planck functions
11!     for each band at the level and layer temperatures.
12
13USE PARKIND1  ,ONLY : JPIM     ,JPRB
14USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
15
16USE PARRRTM  , ONLY : JPLAY     ,JPINPX
17USE YOERRTRF , ONLY :       PREFLOG   ,TREF
18
19IMPLICIT NONE
20
21INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV
22REAL(KIND=JPRB)   ,INTENT(IN)    :: P_COLDRY(JPLAY)
23REAL(KIND=JPRB)   ,INTENT(IN)    :: P_WKL(JPINPX,JPLAY)
24REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_FAC00(JPLAY)
25REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_FAC01(JPLAY)
26REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_FAC10(JPLAY)
27REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_FAC11(JPLAY)
28REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_FORFAC(JPLAY)
29INTEGER(KIND=JPIM),INTENT(OUT)   :: K_JP(JPLAY)
30INTEGER(KIND=JPIM),INTENT(OUT)   :: K_JT(JPLAY)
31INTEGER(KIND=JPIM),INTENT(OUT)   :: K_JT1(JPLAY)
32REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_COLH2O(JPLAY)
33REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_COLCO2(JPLAY)
34REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_COLO3(JPLAY)
35REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_COLN2O(JPLAY)
36REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_COLCH4(JPLAY)
37REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_COLO2(JPLAY)
38REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_CO2MULT(JPLAY)
39INTEGER(KIND=JPIM),INTENT(OUT)   :: K_LAYTROP
40INTEGER(KIND=JPIM),INTENT(OUT)   :: K_LAYSWTCH
41INTEGER(KIND=JPIM),INTENT(OUT)   :: K_LAYLOW
42REAL(KIND=JPRB)   ,INTENT(IN)    :: PAVEL(JPLAY)
43REAL(KIND=JPRB)   ,INTENT(IN)    :: P_TAVEL(JPLAY)
44REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_SELFFAC(JPLAY)
45REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_SELFFRAC(JPLAY)
46INTEGER(KIND=JPIM),INTENT(OUT)   :: K_INDSELF(JPLAY)
47!- from INTFAC     
48!- from INTIND
49!- from PROFDATA             
50!- from PROFILE             
51!- from SELF             
52INTEGER(KIND=JPIM) :: JP1, I_LAY
53
54REAL(KIND=JPRB) :: Z_CO2REG, Z_COMPFP, Z_FACTOR, Z_FP, Z_FT, Z_FT1, Z_PLOG, Z_SCALEFAC, Z_STPFAC, Z_WATER
55REAL(KIND=JPRB) :: ZHOOK_HANDLE
56
57!#include "yoeratm.h"   
58
59IF (LHOOK) CALL DR_HOOK('RRTM_SETCOEF_140GP',0,ZHOOK_HANDLE)
60Z_STPFAC = 296._JPRB/1013._JPRB
61
62K_LAYTROP  = 0
63K_LAYSWTCH = 0
64K_LAYLOW   = 0
65DO I_LAY = 1, KLEV
66!        Find the two reference pressures on either side of the
67!        layer pressure.  Store them in JP and JP1.  Store in FP the
68!        fraction of the difference (in ln(pressure)) between these
69!        two values that the layer pressure lies.
70  Z_PLOG = LOG(PAVEL(I_LAY))
71  K_JP(I_LAY) = INT(36._JPRB - 5*(Z_PLOG+0.04_JPRB))
72  IF (K_JP(I_LAY)  <  1) THEN
73    K_JP(I_LAY) = 1
74  ELSEIF (K_JP(I_LAY)  >  58) THEN
75    K_JP(I_LAY) = 58
76  ENDIF
77  JP1 = K_JP(I_LAY) + 1
78  Z_FP = 5._JPRB * (PREFLOG(K_JP(I_LAY)) - Z_PLOG)
79
80!        Determine, for each reference pressure (JP and JP1), which
81!        reference temperature (these are different for each 
82!        reference pressure) is nearest the layer temperature but does
83!        not exceed it.  Store these indices in JT and JT1, resp.
84!        Store in FT (resp. FT1) the fraction of the way between JT
85!        (JT1) and the next highest reference temperature that the
86!        layer temperature falls.
87
88  K_JT(I_LAY) = INT(3._JPRB + (P_TAVEL(I_LAY)-TREF(K_JP(I_LAY)))/15._JPRB)
89  IF (K_JT(I_LAY)  <  1) THEN
90    K_JT(I_LAY) = 1
91  ELSEIF (K_JT(I_LAY)  >  4) THEN
92    K_JT(I_LAY) = 4
93  ENDIF
94  Z_FT = ((P_TAVEL(I_LAY)-TREF(K_JP(I_LAY)))/15._JPRB) - REAL(K_JT(I_LAY)-3)
95  K_JT1(I_LAY) = INT(3._JPRB + (P_TAVEL(I_LAY)-TREF(JP1))/15._JPRB)
96  IF (K_JT1(I_LAY)  <  1) THEN
97    K_JT1(I_LAY) = 1
98  ELSEIF (K_JT1(I_LAY)  >  4) THEN
99    K_JT1(I_LAY) = 4
100  ENDIF
101  Z_FT1 = ((P_TAVEL(I_LAY)-TREF(JP1))/15._JPRB) - REAL(K_JT1(I_LAY)-3)
102
103  Z_WATER = P_WKL(1,I_LAY)/P_COLDRY(I_LAY)
104  Z_SCALEFAC = PAVEL(I_LAY) * Z_STPFAC / P_TAVEL(I_LAY)
105
106!        If the pressure is less than ~100mb, perform a different
107!        set of species interpolations.
108!         IF (PLOG .LE. 4.56) GO TO 5300
109!--------------------------------------         
110  IF (Z_PLOG  >  4.56_JPRB) THEN
111    K_LAYTROP =  K_LAYTROP + 1
112!        For one band, the "switch" occurs at ~300 mb.
113    IF (Z_PLOG  >=  5.76_JPRB) K_LAYSWTCH = K_LAYSWTCH + 1
114    IF (Z_PLOG  >=  6.62_JPRB) K_LAYLOW = K_LAYLOW + 1
115
116    P_FORFAC(I_LAY) = Z_SCALEFAC / (1.0_JPRB+Z_WATER)
117
118!        Set up factors needed to separately include the water vapor
119!        self-continuum in the calculation of absorption coefficient.
120!C           SELFFAC(LAY) = WATER * SCALEFAC / (1.+WATER)
121    P_SELFFAC(I_LAY) = Z_WATER * P_FORFAC(I_LAY)
122    Z_FACTOR = (P_TAVEL(I_LAY)-188.0_JPRB)/7.2_JPRB
123    K_INDSELF(I_LAY) = MIN(9, MAX(1, INT(Z_FACTOR)-7))
124    P_SELFFRAC(I_LAY) = Z_FACTOR - REAL(K_INDSELF(I_LAY) + 7)
125
126!        Calculate needed column amounts.
127    P_COLH2O(I_LAY) = 1.E-20_JPRB * P_WKL(1,I_LAY)
128    P_COLCO2(I_LAY) = 1.E-20_JPRB * P_WKL(2,I_LAY)
129    P_COLO3(I_LAY)  = 1.E-20_JPRB * P_WKL(3,I_LAY)
130    P_COLN2O(I_LAY) = 1.E-20_JPRB * P_WKL(4,I_LAY)
131    P_COLCH4(I_LAY) = 1.E-20_JPRB * P_WKL(6,I_LAY)
132    P_COLO2(I_LAY)  = 1.E-20_JPRB * P_WKL(7,I_LAY)
133    IF (P_COLCO2(I_LAY)  ==  0.0_JPRB) P_COLCO2(I_LAY) = 1.E-32_JPRB * P_COLDRY(I_LAY)
134    IF (P_COLN2O(I_LAY)  ==  0.0_JPRB) P_COLN2O(I_LAY) = 1.E-32_JPRB * P_COLDRY(I_LAY)
135    IF (P_COLCH4(I_LAY)  ==  0.0_JPRB) P_COLCH4(I_LAY) = 1.E-32_JPRB * P_COLDRY(I_LAY)
136!        Using E = 1334.2 cm-1.
137    Z_CO2REG = 3.55E-24_JPRB * P_COLDRY(I_LAY)
138    P_CO2MULT(I_LAY)= (P_COLCO2(I_LAY) - Z_CO2REG) *&
139     & 272.63_JPRB*EXP(-1919.4_JPRB/P_TAVEL(I_LAY))/(8.7604E-4_JPRB*P_TAVEL(I_LAY)) 
140!         GO TO 5400
141!------------------
142  ELSE
143!        Above LAYTROP.
144! 5300    CONTINUE
145
146!        Calculate needed column amounts.
147    P_FORFAC(I_LAY) = Z_SCALEFAC / (1.0_JPRB+Z_WATER)
148
149    P_COLH2O(I_LAY) = 1.E-20_JPRB * P_WKL(1,I_LAY)
150    P_COLCO2(I_LAY) = 1.E-20_JPRB * P_WKL(2,I_LAY)
151    P_COLO3(I_LAY)  = 1.E-20_JPRB * P_WKL(3,I_LAY)
152    P_COLN2O(I_LAY) = 1.E-20_JPRB * P_WKL(4,I_LAY)
153    P_COLCH4(I_LAY) = 1.E-20_JPRB * P_WKL(6,I_LAY)
154    P_COLO2(I_LAY)  = 1.E-20_JPRB * P_WKL(7,I_LAY)
155    IF (P_COLCO2(I_LAY)  ==  0.0_JPRB) P_COLCO2(I_LAY) = 1.E-32_JPRB * P_COLDRY(I_LAY)
156    IF (P_COLN2O(I_LAY)  ==  0.0_JPRB) P_COLN2O(I_LAY) = 1.E-32_JPRB * P_COLDRY(I_LAY)
157    IF (P_COLCH4(I_LAY)  ==  0.0_JPRB) P_COLCH4(I_LAY) = 1.E-32_JPRB * P_COLDRY(I_LAY)
158    Z_CO2REG = 3.55E-24_JPRB * P_COLDRY(I_LAY)
159    P_CO2MULT(I_LAY)= (P_COLCO2(I_LAY) - Z_CO2REG) *&
160     & 272.63_JPRB*EXP(-1919.4_JPRB/P_TAVEL(I_LAY))/(8.7604E-4_JPRB*P_TAVEL(I_LAY)) 
161!----------------     
162  ENDIF
163! 5400    CONTINUE
164
165!        We have now isolated the layer ln pressure and temperature,
166!        between two reference pressures and two reference temperatures
167!        (for each reference pressure).  We multiply the pressure
168!        fraction FP with the appropriate temperature fractions to get
169!        the factors that will be needed for the interpolation that yields
170!        the optical depths (performed in routines TAUGBn for band n).
171
172  Z_COMPFP = 1.0_JPRB - Z_FP
173  P_FAC10(I_LAY) = Z_COMPFP * Z_FT
174  P_FAC00(I_LAY) = Z_COMPFP * (1.0_JPRB - Z_FT)
175  P_FAC11(I_LAY) = Z_FP * Z_FT1
176  P_FAC01(I_LAY) = Z_FP * (1.0_JPRB - Z_FT1)
177
178ENDDO
179
180! MT 981104
181!-- Set LAYLOW for profiles with surface pressure less than 750 hPa.
182IF (K_LAYLOW == 0) K_LAYLOW=1
183
184IF (LHOOK) CALL DR_HOOK('RRTM_SETCOEF_140GP',1,ZHOOK_HANDLE)
185END SUBROUTINE RRTM_SETCOEF_140GP
Note: See TracBrowser for help on using the repository browser.