source: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad/ifsrrtm/rrtm_taumol16.F90 @ 5097

Last change on this file since 5097 was 4773, checked in by idelkadi, 11 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: 10.7 KB
Line 
1!----------------------------------------------------------------------------
2SUBROUTINE RRTM_TAUMOL16 (KIDIA,KFDIA,KLEV,P_TAU,&
3 & P_TAUAERL,P_FAC00,P_FAC01,P_FAC10,P_FAC11,P_FORFAC,P_FORFRAC,K_INDFOR,K_JP,K_JT,K_JT1,P_ONEMINUS,&
4 & P_COLH2O,P_COLCH4,K_LAYTROP,P_SELFFAC,P_SELFFRAC,K_INDSELF,PFRAC, &
5 & P_RAT_H2OCH4,P_RAT_H2OCH4_1) 
6
7!     BAND 16:  2600-3000 cm-1 (low - H2O,CH4; high - nothing)
8
9!     AUTHOR.
10!     -------
11!      JJMorcrette, ECMWF
12
13!     MODIFICATIONS.
14!     --------------
15!      M.Hamrud      01-Oct-2003 CY28 Cleaning
16!      NEC           25-Oct-2007 Optimisations
17!      JJMorcrette 20110613 flexible number of g-points
18!      ABozzo 201306 updated to rrtmg v4.85
19!      band 16:  2600-3250 cm-1 (low key- h2o,ch4; high key - ch4)
20! ---------------------------------------------------------------------------
21
22USE PARKIND1  ,ONLY : JPIM     ,JPRB
23USE YOMHOOK   ,ONLY : LHOOK, DR_HOOK, JPHOOK
24
25USE PARRRTM  , ONLY : JPBAND
26USE YOERRTM  , ONLY : JPGPT  ,NGS15  ,NG16
27USE YOERRTWN , ONLY : NSPA,NSPB   
28USE YOERRTA16, ONLY : ABSA,ABSB,FRACREFA,FRACREFB,SELFREF,FORREF
29USE YOERRTRF, ONLY : CHI_MLS
30
31IMPLICIT NONE
32
33INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA
34INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA
35INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV
36REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TAU(KIDIA:KFDIA,JPGPT,KLEV)
37REAL(KIND=JPRB)   ,INTENT(IN)    :: P_TAUAERL(KIDIA:KFDIA,KLEV,JPBAND)
38REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC00(KIDIA:KFDIA,KLEV)
39REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC01(KIDIA:KFDIA,KLEV)
40REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC10(KIDIA:KFDIA,KLEV)
41REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC11(KIDIA:KFDIA,KLEV)
42INTEGER(KIND=JPIM),INTENT(IN)    :: K_JP(KIDIA:KFDIA,KLEV)
43INTEGER(KIND=JPIM),INTENT(IN)    :: K_JT(KIDIA:KFDIA,KLEV)
44INTEGER(KIND=JPIM),INTENT(IN)    :: K_JT1(KIDIA:KFDIA,KLEV)
45REAL(KIND=JPRB)   ,INTENT(IN)    :: P_ONEMINUS
46REAL(KIND=JPRB)   ,INTENT(IN)    :: P_COLH2O(KIDIA:KFDIA,KLEV)
47REAL(KIND=JPRB)   ,INTENT(IN)    :: P_COLCH4(KIDIA:KFDIA,KLEV)
48INTEGER(KIND=JPIM),INTENT(IN)    :: K_LAYTROP(KIDIA:KFDIA)
49REAL(KIND=JPRB)   ,INTENT(IN)    :: P_SELFFAC(KIDIA:KFDIA,KLEV)
50REAL(KIND=JPRB)   ,INTENT(IN)    :: P_SELFFRAC(KIDIA:KFDIA,KLEV)
51INTEGER(KIND=JPIM),INTENT(IN)    :: K_INDSELF(KIDIA:KFDIA,KLEV)
52REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFRAC(KIDIA:KFDIA,JPGPT,KLEV)
53
54REAL(KIND=JPRB)   ,INTENT(IN)   :: P_RAT_H2OCH4(KIDIA:KFDIA,KLEV)
55REAL(KIND=JPRB)   ,INTENT(IN)   :: P_RAT_H2OCH4_1(KIDIA:KFDIA,KLEV)
56INTEGER(KIND=JPIM),INTENT(IN)   :: K_INDFOR(KIDIA:KFDIA,KLEV)
57REAL(KIND=JPRB)   ,INTENT(IN)   :: P_FORFAC(KIDIA:KFDIA,KLEV)
58REAL(KIND=JPRB)   ,INTENT(IN)   :: P_FORFRAC(KIDIA:KFDIA,KLEV)
59
60! ---------------------------------------------------------------------------
61
62INTEGER(KIND=JPIM) :: IG, IND0, IND1, INDS,INDF, JS,JS1,JPL,JLAY
63INTEGER(KIND=JPIM) :: JLON
64
65REAL(KIND=JPRB) ::  Z_FAC000, Z_FAC100, Z_FAC200,&
66 & Z_FAC010, Z_FAC110, Z_FAC210, &
67 & Z_FAC001, Z_FAC101, Z_FAC201, &
68 & Z_FAC011, Z_FAC111, Z_FAC211
69REAL(KIND=JPRB) :: ZP, ZP4, ZFK0, ZFK1, ZFK2
70
71REAL(KIND=JPRB) :: ZREFRAT_PLANCK_A
72REAL(KIND=JPRB) :: ZTAUFOR,ZTAUSELF,ZTAU_MAJOR,ZTAU_MAJOR1
73REAL(KIND=JPRB) :: Z_FS, Z_SPECMULT, Z_SPECPARM,Z_SPECCOMB,  &
74& Z_FS1, Z_SPECMULT1, Z_SPECPARM1,Z_SPECCOMB1, &
75& Z_FPL, Z_SPECMULT_PLANCK, Z_SPECPARM_PLANCK,Z_SPECCOMB_PLANCK
76REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
77
78IF (LHOOK) CALL DR_HOOK('RRTM_TAUMOL16',0,ZHOOK_HANDLE)
79
80! Calculate reference ratio to be used in calculation of Planck
81! fraction in lower atmosphere.
82
83! P = 387. mb (Level 6)
84      ZREFRAT_PLANCK_A = CHI_MLS(1,6)/CHI_MLS(6,6)
85
86! Compute the optical depth by interpolating in ln(pressure),
87! temperature,and appropriate species.  Below laytrop, the water
88! vapor self-continuum and foreign continuum is interpolated
89! (in temperature) separately. 
90
91DO JLAY = 1, KLEV
92  DO JLON = KIDIA, KFDIA
93    IF (JLAY <= K_LAYTROP(JLON)) THEN
94      Z_SPECCOMB = P_COLH2O(JLON,JLAY) + P_RAT_H2OCH4(JLON,JLAY)*P_COLCH4(JLON,JLAY)
95      Z_SPECPARM = P_COLH2O(JLON,JLAY)/Z_SPECCOMB
96      Z_SPECPARM = MIN(Z_SPECPARM,P_ONEMINUS)
97      Z_SPECMULT = 8._JPRB*(Z_SPECPARM)
98      JS = 1 + INT(Z_SPECMULT)
99      Z_FS = MOD(Z_SPECMULT,1.0_JPRB)
100
101      Z_SPECCOMB1 = P_COLH2O(JLON,JLAY) + P_RAT_H2OCH4_1(JLON,JLAY)*P_COLCH4(JLON,JLAY)
102      Z_SPECPARM1 = P_COLH2O(JLON,JLAY)/Z_SPECCOMB1
103      IF (Z_SPECPARM1 >= P_ONEMINUS) Z_SPECPARM1 = P_ONEMINUS
104      Z_SPECMULT1 = 8._JPRB*(Z_SPECPARM1)
105      JS1 = 1 + INT(Z_SPECMULT1)
106      Z_FS1 = MOD(Z_SPECMULT1,1.0_JPRB)
107
108      Z_SPECCOMB_PLANCK = P_COLH2O(JLON,JLAY)+ZREFRAT_PLANCK_A*P_COLCH4(JLON,JLAY)
109      Z_SPECPARM_PLANCK = P_COLH2O(JLON,JLAY)/Z_SPECCOMB_PLANCK
110      IF (Z_SPECPARM_PLANCK >= P_ONEMINUS) Z_SPECPARM_PLANCK=P_ONEMINUS
111      Z_SPECMULT_PLANCK = 8._JPRB*Z_SPECPARM_PLANCK
112      JPL= 1 + INT(Z_SPECMULT_PLANCK)
113      Z_FPL = MOD(Z_SPECMULT_PLANCK,1.0_JPRB)
114
115      IND0 = ((K_JP(JLON,JLAY)-1)*5+(K_JT(JLON,JLAY)-1))*NSPA(16) + JS
116      IND1 = (K_JP(JLON,JLAY)*5+(K_JT1(JLON,JLAY)-1))*NSPA(16) + JS1
117      INDS = K_INDSELF(JLON,JLAY)
118      INDF = K_INDFOR(JLON,JLAY)
119
120IF (Z_SPECPARM < 0.125_JPRB) THEN
121            ZP = Z_FS - 1
122            ZP4 = ZP**4
123            ZFK0 = ZP4
124            ZFK1 = 1 - ZP - 2.0_JPRB*ZP4
125            ZFK2 = ZP + ZP4
126            Z_FAC000 = ZFK0*P_FAC00(JLON,JLAY)
127            Z_FAC100 = ZFK1*P_FAC00(JLON,JLAY)
128            Z_FAC200 = ZFK2*P_FAC00(JLON,JLAY)
129            Z_FAC010 = ZFK0*P_FAC10(JLON,JLAY)
130            Z_FAC110 = ZFK1*P_FAC10(JLON,JLAY)
131            Z_FAC210 = ZFK2*P_FAC10(JLON,JLAY)
132      ELSEIF (Z_SPECPARM > 0.875_JPRB) THEN
133            ZP = -Z_FS
134            ZP4 = ZP**4
135            ZFK0 = ZP4
136            ZFK1 = 1 - ZP - 2.0_JPRB*ZP4
137            ZFK2 = ZP + ZP4
138            Z_FAC000 = ZFK0*P_FAC00(JLON,JLAY)
139            Z_FAC100 = ZFK1*P_FAC00(JLON,JLAY)
140            Z_FAC200 = ZFK2*P_FAC00(JLON,JLAY)
141            Z_FAC010 = ZFK0*P_FAC10(JLON,JLAY)
142            Z_FAC110 = ZFK1*P_FAC10(JLON,JLAY)
143            Z_FAC210 = ZFK2*P_FAC10(JLON,JLAY)
144      ELSE
145            Z_FAC000 = (1._JPRB - Z_FS) * P_FAC00(JLON,JLAY)
146            Z_FAC010 = (1._JPRB - Z_FS) * P_FAC10(JLON,JLAY)
147            Z_FAC100 = Z_FS * P_FAC00(JLON,JLAY)
148            Z_FAC110 = Z_FS * P_FAC10(JLON,JLAY)
149      ENDIF
150      IF (Z_SPECPARM1 < 0.125_JPRB) THEN
151            ZP = Z_FS1 - 1
152            ZP4 = ZP**4
153            ZFK0 = ZP4
154            ZFK1 = 1 - ZP - 2.0_JPRB*ZP4
155            ZFK2 = ZP + ZP4
156            Z_FAC001 = ZFK0*P_FAC01(JLON,JLAY)
157            Z_FAC101 = ZFK1*P_FAC01(JLON,JLAY)
158            Z_FAC201 = ZFK2*P_FAC01(JLON,JLAY)
159            Z_FAC011 = ZFK0*P_FAC11(JLON,JLAY)
160            Z_FAC111 = ZFK1*P_FAC11(JLON,JLAY)
161            Z_FAC211 = ZFK2*P_FAC11(JLON,JLAY)
162      ELSEIF (Z_SPECPARM1 > 0.875_JPRB) THEN
163            ZP = -Z_FS1
164            ZP4 = ZP**4
165            ZFK0 = ZP4
166            ZFK1 = 1 - ZP - 2.0_JPRB*ZP4
167            ZFK2 = ZP + ZP4
168            Z_FAC001 = ZFK0*P_FAC01(JLON,JLAY)
169            Z_FAC101 = ZFK1*P_FAC01(JLON,JLAY)
170            Z_FAC201 = ZFK2*P_FAC01(JLON,JLAY)
171            Z_FAC011 = ZFK0*P_FAC11(JLON,JLAY)
172            Z_FAC111 = ZFK1*P_FAC11(JLON,JLAY)
173            Z_FAC211 = ZFK2*P_FAC11(JLON,JLAY)
174      ELSE
175            Z_FAC001 = (1._JPRB - Z_FS1) * P_FAC01(JLON,JLAY)
176            Z_FAC011 = (1._JPRB - Z_FS1) * P_FAC11(JLON,JLAY)
177            Z_FAC101 = Z_FS1 * P_FAC01(JLON,JLAY)
178            Z_FAC111 = Z_FS1 * P_FAC11(JLON,JLAY)
179      ENDIF
180
181
182      DO IG = 1, NG16
183!- - DS_990714
184        ZTAUSELF = P_SELFFAC(JLON,JLAY)* (SELFREF(INDS,IG) + P_SELFFRAC(JLON,JLAY) * &
185          &       (SELFREF(INDS+1,IG) - SELFREF(INDS,IG)))
186        ZTAUFOR = P_FORFAC(JLON,JLAY) * (FORREF(INDF,IG) + P_FORFRAC(JLON,JLAY) * &
187          &       (FORREF(INDF+1,IG) - FORREF(INDF,IG)))
188
189     IF (Z_SPECPARM < 0.125_JPRB) THEN
190               ZTAU_MAJOR = Z_SPECCOMB * &
191                 &   (Z_FAC000 * ABSA(IND0,IG) + &
192                 &   Z_FAC100 * ABSA(IND0+1,IG) + &
193                 &   Z_FAC200 * ABSA(IND0+2,IG) + &
194                 &   Z_FAC010 * ABSA(IND0+9,IG) + &
195                 &   Z_FAC110 * ABSA(IND0+10,IG) + &
196                 &   Z_FAC210 * ABSA(IND0+11,IG))
197            ELSEIF (Z_SPECPARM > 0.875_JPRB) THEN
198               ZTAU_MAJOR = Z_SPECCOMB * &
199                 &   (Z_FAC200 * ABSA(IND0-1,IG) + &
200                 &   Z_FAC100 * ABSA(IND0,IG) + &
201                 &   Z_FAC000 * ABSA(IND0+1,IG) + &
202                 &   Z_FAC210 * ABSA(IND0+8,IG) + &
203                 &   Z_FAC110 * ABSA(IND0+9,IG) + &
204                 &   Z_FAC010 * ABSA(IND0+10,IG))
205            ELSE
206               ZTAU_MAJOR = Z_SPECCOMB * &
207                 &   (Z_FAC000 * ABSA(IND0,IG) + &
208                 &   Z_FAC100 * ABSA(IND0+1,IG) + &
209                 &   Z_FAC010 * ABSA(IND0+9,IG) + &
210                 &   Z_FAC110 * ABSA(IND0+10,IG))
211            ENDIF
212
213            IF (Z_SPECPARM1 < 0.125_JPRB) THEN
214               ZTAU_MAJOR1 = Z_SPECCOMB1 * &
215                &    (Z_FAC001 * ABSA(IND1,IG) + &
216                &    Z_FAC101 * ABSA(IND1+1,IG) + &
217                &    Z_FAC201 * ABSA(IND1+2,IG) + &
218                &    Z_FAC011 * ABSA(IND1+9,IG) + &
219                &    Z_FAC111 * ABSA(IND1+10,IG) + &
220                &    Z_FAC211 * ABSA(IND1+11,IG))
221            ELSEIF (Z_SPECPARM1 > 0.875_JPRB) THEN
222               ZTAU_MAJOR1 = Z_SPECCOMB1 * &
223                &    (Z_FAC201 * ABSA(IND1-1,IG) + &
224                &    Z_FAC101 * ABSA(IND1,IG) + &
225                &    Z_FAC001 * ABSA(IND1+1,IG) + &
226                &    Z_FAC211 * ABSA(IND1+8,IG) + &
227                &    Z_FAC111 * ABSA(IND1+9,IG) + &
228                &    Z_FAC011 * ABSA(IND1+10,IG))
229            ELSE
230               ZTAU_MAJOR1 = Z_SPECCOMB1 * &
231                &    (Z_FAC001 * ABSA(IND1,IG) +  &
232                &    Z_FAC101 * ABSA(IND1+1,IG) + &
233                &    Z_FAC011 * ABSA(IND1+9,IG) + &
234                &    Z_FAC111 * ABSA(IND1+10,IG))
235            ENDIF
236
237
238        P_TAU(JLON,NGS15+IG,JLAY) = ZTAU_MAJOR + ZTAU_MAJOR1 &
239                & + ZTAUSELF + ZTAUFOR &
240                & + P_TAUAERL(JLON,JLAY,16) 
241        PFRAC(JLON,NGS15+IG,JLAY) = FRACREFA(IG,JPL) + Z_FPL * &
242         & (FRACREFA(IG,JPL+1) - FRACREFA(IG,JPL))
243
244        ENDDO
245
246!-- DS_990714
247    ENDIF
248
249    IF (JLAY > K_LAYTROP(JLON)) THEN
250      IND0 = ((K_JP(JLON,JLAY)-13)*5+(K_JT(JLON,JLAY)-1))*NSPB(16) + 1
251      IND1 = ((K_JP(JLON,JLAY)-12)*5+(K_JT1(JLON,JLAY)-1))*NSPB(16) + 1
252      DO IG = 1, NG16
253
254        P_TAU(JLON,NGS15+IG,JLAY) =  P_COLCH4(JLON,JLAY) * &
255               &  (P_FAC00(JLON,JLAY) * ABSB(IND0,IG) + &
256               &  P_FAC10(JLON,JLAY) * ABSB(IND0+1,IG) + &
257               &  P_FAC01(JLON,JLAY) * ABSB(IND1,IG) + &
258               &  P_FAC11(JLON,JLAY) * ABSB(IND1+1,IG)) + &
259               &  P_TAUAERL(JLON,JLAY,16)
260        PFRAC(JLON,NGS15+IG,JLAY) = FRACREFB(IG)
261
262      ENDDO
263    ENDIF
264  ENDDO
265ENDDO
266
267IF (LHOOK) CALL DR_HOOK('RRTM_TAUMOL16',1,ZHOOK_HANDLE)
268
269END SUBROUTINE RRTM_TAUMOL16
Note: See TracBrowser for help on using the repository browser.