source: LMDZ6/trunk/libf/phylmd/ecrad/ifsrrtm/rrtm_taumol15.F90 @ 4999

Last change on this file since 4999 was 4773, checked in by idelkadi, 7 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: 12.0 KB
Line 
1!----------------------------------------------------------------------------
2SUBROUTINE RRTM_TAUMOL15 (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_COLCO2,P_COLN2O,K_LAYTROP,P_SELFFAC,P_SELFFRAC,K_INDSELF,PFRAC, &
5 & PRAT_N2OCO2, PRAT_N2OCO2_1,PMINORFRAC,KINDMINOR,PSCALEMINOR,PCOLBRD) 
6
7!     BAND 15:  2380-2600 cm-1 (low - N2O,CO2; 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 2001306 updated to rrtmg v4.85
19!     band 15:  2380-2600 cm-1 (low - n2o,co2; low minor - n2)
20!                              (high - nothing)
21! ---------------------------------------------------------------------------
22
23USE PARKIND1  ,ONLY : JPIM     ,JPRB
24USE YOMHOOK   ,ONLY : LHOOK, DR_HOOK, JPHOOK
25
26USE PARRRTM  , ONLY : JPBAND
27USE YOERRTM  , ONLY : JPGPT  ,NGS14  ,NG15
28USE YOERRTWN , ONLY : NSPA   
29USE YOERRTA15, ONLY : ABSA   ,KA_MN2,FRACREFA,SELFREF,FORREF
30USE YOERRTRF, ONLY : CHI_MLS
31
32IMPLICIT NONE
33
34INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA
35INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA
36INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV
37REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TAU(KIDIA:KFDIA,JPGPT,KLEV)
38REAL(KIND=JPRB)   ,INTENT(IN)    :: P_TAUAERL(KIDIA:KFDIA,KLEV,JPBAND)
39REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC00(KIDIA:KFDIA,KLEV)
40REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC01(KIDIA:KFDIA,KLEV)
41REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC10(KIDIA:KFDIA,KLEV)
42REAL(KIND=JPRB)   ,INTENT(IN)    :: P_FAC11(KIDIA:KFDIA,KLEV)
43INTEGER(KIND=JPIM),INTENT(IN)    :: K_JP(KIDIA:KFDIA,KLEV)
44INTEGER(KIND=JPIM),INTENT(IN)    :: K_JT(KIDIA:KFDIA,KLEV)
45INTEGER(KIND=JPIM),INTENT(IN)    :: K_JT1(KIDIA:KFDIA,KLEV)
46REAL(KIND=JPRB)   ,INTENT(IN)    :: P_ONEMINUS
47REAL(KIND=JPRB)   ,INTENT(IN)    :: P_COLH2O(KIDIA:KFDIA,KLEV)
48REAL(KIND=JPRB)   ,INTENT(IN)    :: P_COLCO2(KIDIA:KFDIA,KLEV)
49REAL(KIND=JPRB)   ,INTENT(IN)    :: P_COLN2O(KIDIA:KFDIA,KLEV)
50INTEGER(KIND=JPIM),INTENT(IN)    :: K_LAYTROP(KIDIA:KFDIA)
51REAL(KIND=JPRB)   ,INTENT(IN)    :: P_SELFFAC(KIDIA:KFDIA,KLEV)
52REAL(KIND=JPRB)   ,INTENT(IN)    :: P_SELFFRAC(KIDIA:KFDIA,KLEV)
53INTEGER(KIND=JPIM),INTENT(IN)    :: K_INDSELF(KIDIA:KFDIA,KLEV)
54REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFRAC(KIDIA:KFDIA,JPGPT,KLEV)
55
56REAL(KIND=JPRB)   ,INTENT(IN)   :: PRAT_N2OCO2(KIDIA:KFDIA,KLEV)
57REAL(KIND=JPRB)   ,INTENT(IN)   :: PRAT_N2OCO2_1(KIDIA:KFDIA,KLEV)
58INTEGER(KIND=JPIM),INTENT(IN)   :: K_INDFOR(KIDIA:KFDIA,KLEV)
59REAL(KIND=JPRB)   ,INTENT(IN)   :: P_FORFAC(KIDIA:KFDIA,KLEV)
60REAL(KIND=JPRB)   ,INTENT(IN)   :: P_FORFRAC(KIDIA:KFDIA,KLEV)
61REAL(KIND=JPRB)   ,INTENT(IN)   :: PMINORFRAC(KIDIA:KFDIA,KLEV)
62INTEGER(KIND=JPIM),INTENT(IN)   :: KINDMINOR(KIDIA:KFDIA,KLEV)
63REAL(KIND=JPRB)   ,INTENT(IN)    :: PSCALEMINOR(KIDIA:KFDIA,KLEV)
64REAL(KIND=JPRB)   ,INTENT(IN)    :: PCOLBRD(KIDIA:KFDIA,KLEV)         
65! ---------------------------------------------------------------------------
66
67INTEGER(KIND=JPIM) :: IG, IND0, IND1, INDS,INDF,INDM, JS,JS1,JPL,JMN2, JLAY
68INTEGER(KIND=JPIM) :: JLON
69REAL(KIND=JPRB) :: ZREFRAT_PLANCK_A, ZREFRAT_M_A
70REAL(KIND=JPRB) :: ZTAUFOR,ZTAUSELF,ZTAU_MAJOR,ZTAU_MAJOR1, ZN2M1, ZN2M2, ZTAUN2,ZSCALEN2
71REAL(KIND=JPRB) ::  Z_FAC000, Z_FAC100, Z_FAC200,&
72 & Z_FAC010, Z_FAC110, Z_FAC210, &
73 & Z_FAC001, Z_FAC101, Z_FAC201, &
74 & Z_FAC011, Z_FAC111, Z_FAC211
75REAL(KIND=JPRB) :: ZP, ZP4, ZFK0, ZFK1, ZFK2
76
77REAL(KIND=JPRB) :: Z_FS, Z_SPECMULT, Z_SPECPARM,Z_SPECCOMB,  &
78& Z_FS1, Z_SPECMULT1, Z_SPECPARM1,Z_SPECCOMB1, &
79& Z_FMN2, Z_SPECMULT_MN2, Z_SPECPARM_MN2,Z_SPECCOMB_MN2, &
80& Z_FPL, Z_SPECMULT_PLANCK, Z_SPECPARM_PLANCK,Z_SPECCOMB_PLANCK
81REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
82
83IF (LHOOK) CALL DR_HOOK('RRTM_TAUMOL15',0,ZHOOK_HANDLE)
84
85! ---------------------------------------------------------------------------
86! Minor gas mapping level :
87!     Lower - Nitrogen Continuum, P = 1053., T = 294.
88
89! Calculate reference ratio to be used in calculation of Planck
90! fraction in lower atmosphere.
91! P = 1053. mb (Level 1)
92      ZREFRAT_PLANCK_A = CHI_MLS(4,1)/CHI_MLS(2,1)
93
94! P = 1053.
95      ZREFRAT_M_A = CHI_MLS(4,1)/CHI_MLS(2,1)
96
97! Compute the optical depth by interpolating in ln(pressure),
98! temperature, and appropriate species.  Below laytrop, the water
99! vapor self-continuum and foreign continuum is interpolated
100! (in temperature) separately. 
101 
102DO JLAY = 1, KLEV
103  DO JLON = KIDIA, KFDIA
104    IF (JLAY <= K_LAYTROP(JLON)) THEN
105
106      Z_SPECCOMB = P_COLN2O(JLON,JLAY) + PRAT_N2OCO2(JLON,JLAY)*P_COLCO2(JLON,JLAY)
107      Z_SPECPARM = P_COLN2O(JLON,JLAY)/Z_SPECCOMB
108      Z_SPECPARM = MIN(Z_SPECPARM,P_ONEMINUS)
109      Z_SPECMULT = 8._JPRB*(Z_SPECPARM)
110      JS = 1 + INT(Z_SPECMULT)
111      Z_FS = MOD(Z_SPECMULT,1.0_JPRB)
112
113      Z_SPECCOMB1 = P_COLN2O(JLON,JLAY) + PRAT_N2OCO2_1(JLON,JLAY)*P_COLCO2(JLON,JLAY)
114      Z_SPECPARM1 = P_COLN2O(JLON,JLAY)/Z_SPECCOMB1
115      IF (Z_SPECPARM1 >= P_ONEMINUS) Z_SPECPARM1 = P_ONEMINUS
116      Z_SPECMULT1 = 8._JPRB*(Z_SPECPARM1)
117      JS1 = 1 + INT(Z_SPECMULT1)
118      Z_FS1 = MOD(Z_SPECMULT1,1.0_JPRB)
119
120      Z_SPECCOMB_MN2 = P_COLN2O(JLON,JLAY) + ZREFRAT_M_A*P_COLCO2(JLON,JLAY)
121      Z_SPECPARM_MN2 = P_COLN2O(JLON,JLAY)/Z_SPECCOMB_MN2
122      IF (Z_SPECPARM_MN2 >= P_ONEMINUS) Z_SPECPARM_MN2 = P_ONEMINUS
123      Z_SPECMULT_MN2 = 8._JPRB*Z_SPECPARM_MN2
124      JMN2 = 1 + INT(Z_SPECMULT_MN2)
125      Z_FMN2 = MOD(Z_SPECMULT_MN2,1.0_JPRB)
126     
127      Z_SPECCOMB_PLANCK = P_COLN2O(JLON,JLAY)+ZREFRAT_PLANCK_A*P_COLCO2(JLON,JLAY)
128      Z_SPECPARM_PLANCK = P_COLN2O(JLON,JLAY)/Z_SPECCOMB_PLANCK
129      IF (Z_SPECPARM_PLANCK >= P_ONEMINUS) Z_SPECPARM_PLANCK=P_ONEMINUS
130      Z_SPECMULT_PLANCK = 8._JPRB*Z_SPECPARM_PLANCK
131      JPL= 1 + INT(Z_SPECMULT_PLANCK)
132      Z_FPL = MOD(Z_SPECMULT_PLANCK,1.0_JPRB)
133
134      IND0 = ((K_JP(JLON,JLAY)-1)*5+(K_JT(JLON,JLAY)-1))*NSPA(15) + JS
135      IND1 = (K_JP(JLON,JLAY)*5+(K_JT1(JLON,JLAY)-1))*NSPA(15) + JS1
136      INDS = K_INDSELF(JLON,JLAY)
137      INDF = K_INDFOR(JLON,JLAY)
138      INDM = KINDMINOR(JLON,JLAY)
139 
140      ZSCALEN2 = PCOLBRD(JLON,JLAY)*PSCALEMINOR(JLON,JLAY)
141
142IF (Z_SPECPARM < 0.125_JPRB) THEN
143            ZP = Z_FS - 1
144            ZP4 = ZP**4
145            ZFK0 = ZP4
146            ZFK1 = 1 - ZP - 2.0_JPRB*ZP4
147            ZFK2 = ZP + ZP4
148            Z_FAC000 = ZFK0*P_FAC00(JLON,JLAY)
149            Z_FAC100 = ZFK1*P_FAC00(JLON,JLAY)
150            Z_FAC200 = ZFK2*P_FAC00(JLON,JLAY)
151            Z_FAC010 = ZFK0*P_FAC10(JLON,JLAY)
152            Z_FAC110 = ZFK1*P_FAC10(JLON,JLAY)
153            Z_FAC210 = ZFK2*P_FAC10(JLON,JLAY)
154      ELSEIF (Z_SPECPARM > 0.875_JPRB) THEN
155            ZP = -Z_FS
156            ZP4 = ZP**4
157            ZFK0 = ZP4
158            ZFK1 = 1 - ZP - 2.0_JPRB*ZP4
159            ZFK2 = ZP + ZP4
160            Z_FAC000 = ZFK0*P_FAC00(JLON,JLAY)
161            Z_FAC100 = ZFK1*P_FAC00(JLON,JLAY)
162            Z_FAC200 = ZFK2*P_FAC00(JLON,JLAY)
163            Z_FAC010 = ZFK0*P_FAC10(JLON,JLAY)
164            Z_FAC110 = ZFK1*P_FAC10(JLON,JLAY)
165            Z_FAC210 = ZFK2*P_FAC10(JLON,JLAY)
166      ELSE
167            Z_FAC000 = (1._JPRB - Z_FS) * P_FAC00(JLON,JLAY)
168            Z_FAC010 = (1._JPRB - Z_FS) * P_FAC10(JLON,JLAY)
169            Z_FAC100 = Z_FS * P_FAC00(JLON,JLAY)
170            Z_FAC110 = Z_FS * P_FAC10(JLON,JLAY)
171      ENDIF
172      IF (Z_SPECPARM1 < 0.125_JPRB) THEN
173            ZP = Z_FS1 - 1
174            ZP4 = ZP**4
175            ZFK0 = ZP4
176            ZFK1 = 1 - ZP - 2.0_JPRB*ZP4
177            ZFK2 = ZP + ZP4
178            Z_FAC001 = ZFK0*P_FAC01(JLON,JLAY)
179            Z_FAC101 = ZFK1*P_FAC01(JLON,JLAY)
180            Z_FAC201 = ZFK2*P_FAC01(JLON,JLAY)
181            Z_FAC011 = ZFK0*P_FAC11(JLON,JLAY)
182            Z_FAC111 = ZFK1*P_FAC11(JLON,JLAY)
183            Z_FAC211 = ZFK2*P_FAC11(JLON,JLAY)
184      ELSEIF (Z_SPECPARM1 > 0.875_JPRB) THEN
185            ZP = -Z_FS1
186            ZP4 = ZP**4
187            ZFK0 = ZP4
188            ZFK1 = 1 - ZP - 2.0_JPRB*ZP4
189            ZFK2 = ZP + ZP4
190            Z_FAC001 = ZFK0*P_FAC01(JLON,JLAY)
191            Z_FAC101 = ZFK1*P_FAC01(JLON,JLAY)
192            Z_FAC201 = ZFK2*P_FAC01(JLON,JLAY)
193            Z_FAC011 = ZFK0*P_FAC11(JLON,JLAY)
194            Z_FAC111 = ZFK1*P_FAC11(JLON,JLAY)
195            Z_FAC211 = ZFK2*P_FAC11(JLON,JLAY)
196      ELSE
197            Z_FAC001 = (1._JPRB - Z_FS1) * P_FAC01(JLON,JLAY)
198            Z_FAC011 = (1._JPRB - Z_FS1) * P_FAC11(JLON,JLAY)
199            Z_FAC101 = Z_FS1 * P_FAC01(JLON,JLAY)
200            Z_FAC111 = Z_FS1 * P_FAC11(JLON,JLAY)
201      ENDIF
202
203
204!-- DS_990714 
205!-- jjm20110728 re-establishing the loop instead of specified IG to allow a flexible number of NG15
206      DO IG = 1, NG15
207        ZTAUSELF = P_SELFFAC(JLON,JLAY)* (SELFREF(INDS,IG) + P_SELFFRAC(JLON,JLAY) * &
208          &       (SELFREF(INDS+1,IG) - SELFREF(INDS,IG)))
209        ZTAUFOR = P_FORFAC(JLON,JLAY) * (FORREF(INDF,IG) + P_FORFRAC(JLON,JLAY) * &
210          &       (FORREF(INDF+1,IG) - FORREF(INDF,IG)))
211        ZN2M1 = KA_MN2(JMN2,INDM,IG) + Z_FMN2 * &
212          &       (KA_MN2(JMN2+1,INDM,IG) - KA_MN2(JMN2,INDM,IG))
213        ZN2M2 = KA_MN2(JMN2,INDM+1,IG) + Z_FMN2 * &
214          &       (KA_MN2(JMN2+1,INDM+1,IG) - KA_MN2(JMN2,INDM+1,IG))
215        ZTAUN2 = ZSCALEN2 * (ZN2M1 + PMINORFRAC(JLON,JLAY) * (ZN2M2 - ZN2M1))
216
217     IF (Z_SPECPARM < 0.125_JPRB) THEN
218               ZTAU_MAJOR = Z_SPECCOMB * &
219                 &   (Z_FAC000 * ABSA(IND0,IG) + &
220                 &   Z_FAC100 * ABSA(IND0+1,IG) + &
221                 &   Z_FAC200 * ABSA(IND0+2,IG) + &
222                 &   Z_FAC010 * ABSA(IND0+9,IG) + &
223                 &   Z_FAC110 * ABSA(IND0+10,IG) + &
224                 &   Z_FAC210 * ABSA(IND0+11,IG))
225            ELSEIF (Z_SPECPARM > 0.875_JPRB) THEN
226               ZTAU_MAJOR = Z_SPECCOMB * &
227                 &   (Z_FAC200 * ABSA(IND0-1,IG) + &
228                 &   Z_FAC100 * ABSA(IND0,IG) + &
229                 &   Z_FAC000 * ABSA(IND0+1,IG) + &
230                 &   Z_FAC210 * ABSA(IND0+8,IG) + &
231                 &   Z_FAC110 * ABSA(IND0+9,IG) + &
232                 &   Z_FAC010 * ABSA(IND0+10,IG))
233            ELSE
234               ZTAU_MAJOR = Z_SPECCOMB * &
235                 &   (Z_FAC000 * ABSA(IND0,IG) + &
236                 &   Z_FAC100 * ABSA(IND0+1,IG) + &
237                 &   Z_FAC010 * ABSA(IND0+9,IG) + &
238                 &   Z_FAC110 * ABSA(IND0+10,IG))
239            ENDIF
240
241            IF (Z_SPECPARM1 < 0.125_JPRB) THEN
242               ZTAU_MAJOR1 = Z_SPECCOMB1 * &
243                &    (Z_FAC001 * ABSA(IND1,IG) + &
244                &    Z_FAC101 * ABSA(IND1+1,IG) + &
245                &    Z_FAC201 * ABSA(IND1+2,IG) + &
246                &    Z_FAC011 * ABSA(IND1+9,IG) + &
247                &    Z_FAC111 * ABSA(IND1+10,IG) + &
248                &    Z_FAC211 * ABSA(IND1+11,IG))
249            ELSEIF (Z_SPECPARM1 > 0.875_JPRB) THEN
250               ZTAU_MAJOR1 = Z_SPECCOMB1 * &
251                &    (Z_FAC201 * ABSA(IND1-1,IG) + &
252                &    Z_FAC101 * ABSA(IND1,IG) + &
253                &    Z_FAC001 * ABSA(IND1+1,IG) + &
254                &    Z_FAC211 * ABSA(IND1+8,IG) + &
255                &    Z_FAC111 * ABSA(IND1+9,IG) + &
256                &    Z_FAC011 * ABSA(IND1+10,IG))
257            ELSE
258               ZTAU_MAJOR1 = Z_SPECCOMB1 * &
259                &    (Z_FAC001 * ABSA(IND1,IG) +  &
260                &    Z_FAC101 * ABSA(IND1+1,IG) + &
261                &    Z_FAC011 * ABSA(IND1+9,IG) + &
262                &    Z_FAC111 * ABSA(IND1+10,IG))
263            ENDIF
264
265
266        P_TAU(JLON,NGS14+IG,JLAY) = ZTAU_MAJOR + ZTAU_MAJOR1 &
267               & + ZTAUSELF + ZTAUFOR + ZTAUN2 &
268               & + P_TAUAERL(JLON,JLAY,15) 
269       PFRAC(JLON,NGS14+IG,JLAY) = FRACREFA(IG,JPL) + Z_FPL * &
270         & (FRACREFA(IG,JPL+1) - FRACREFA(IG,JPL)) 
271      ENDDO
272!-- jjm20110728
273!-- DS_990714 
274    ENDIF
275
276    IF (JLAY > K_LAYTROP(JLON)) THEN
277      DO IG = 1, NG15
278!-- jjm20110728 re-establishing the loop instead of specified IG to allow a flexible number of NG15
279!-- DS_990714 
280        P_TAU(JLON,NGS14+IG,JLAY) = P_TAUAERL(JLON,JLAY,15)
281        PFRAC(JLON,NGS14+IG,JLAY) = 0.0_JPRB
282      ENDDO
283!-- jjm20110728
284!-- DS_990714 
285    ENDIF
286  ENDDO
287ENDDO
288
289IF (LHOOK) CALL DR_HOOK('RRTM_TAUMOL15',1,ZHOOK_HANDLE)
290
291END SUBROUTINE RRTM_TAUMOL15
Note: See TracBrowser for help on using the repository browser.