source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/rrtm/rrtm_ecrt_140gp.F90 @ 4249

Last change on this file since 4249 was 3990, checked in by millour, 8 years ago

An updated LMDZ5 (based on LMDZ rev 2786) to play with.
EM
File size: 13.4 KB
Line 
1!
2! $Id: rrtm_ecrt_140gp.F90 2626 2016-09-15 14:20:56Z musat $
3!
4!****************** SUBROUTINE RRTM_ECRT_140GP **************************
5
6SUBROUTINE RRTM_ECRT_140GP &
7 & ( K_IPLON, klon , klev, kcld,&
8 & paer , paph , pap,&
9 & pts  , pth  , pt,&
10 & P_ZEMIS, P_ZEMIW,&
11 & pq   , pcco2, pozn, pcldf, ptaucld, ptclear,&
12 & P_CLDFRAC,P_TAUCLD,&
13 & PTAU_LW,&
14 & P_COLDRY,P_WKL,P_WX,&
15 & P_TAUAERL,PAVEL,P_TAVEL,PZ,P_TZ,P_TBOUND,K_NLAYERS,P_SEMISS,K_IREFLECT ) 
16
17!     Reformatted for F90 by JJMorcrette, ECMWF, 980714
18
19!     Read in atmospheric profile from ECMWF radiation code, and prepare it
20!     for use in RRTM.  Set other RRTM input parameters.  Values are passed
21!     back through existing RRTM arrays and commons.
22
23!- Modifications
24
25!     2000-05-15 Deborah Salmond  Speed-up
26
27USE PARKIND1  ,ONLY : JPIM     ,JPRB
28USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
29
30USE PARRRTM  , ONLY : JPBAND   ,JPXSEC   ,JPLAY   ,&
31 & JPINPX 
32USE YOERAD   , ONLY : NLW      ,NOVLP
33!MPL/IM 20160915 on prend GES de phylmd USE YOERDI   , ONLY :    RCH4     ,RN2O    ,RCFC11  ,RCFC12
34USE YOESW    , ONLY : RAER
35
36!------------------------------Arguments--------------------------------
37
38IMPLICIT NONE
39
40
41INTEGER(KIND=JPIM),INTENT(IN)    :: KLON! Number of atmospheres (longitudes)
42INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV! Number of atmospheric layers
43INTEGER(KIND=JPIM),INTENT(IN)    :: K_IPLON
44INTEGER(KIND=JPIM),INTENT(OUT)   :: KCLD
45REAL(KIND=JPRB)   ,INTENT(IN)    :: PAER(KLON,6,KLEV) ! Aerosol optical thickness
46REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPH(KLON,KLEV+1) ! Interface pressures (Pa)
47REAL(KIND=JPRB)   ,INTENT(IN)    :: PAP(KLON,KLEV) ! Layer pressures (Pa)
48REAL(KIND=JPRB)   ,INTENT(IN)    :: PTS(KLON) ! Surface temperature (K)
49REAL(KIND=JPRB)   ,INTENT(IN)    :: PTH(KLON,KLEV+1) ! Interface temperatures (K)
50REAL(KIND=JPRB)   ,INTENT(IN)    :: PT(KLON,KLEV) ! Layer temperature (K)
51REAL(KIND=JPRB)   ,INTENT(IN)    :: P_ZEMIS(KLON) ! Non-window surface emissivity
52REAL(KIND=JPRB)   ,INTENT(IN)    :: P_ZEMIW(KLON) ! Window surface emissivity
53REAL(KIND=JPRB)   ,INTENT(IN)    :: PQ(KLON,KLEV) ! H2O specific humidity (mmr)
54REAL(KIND=JPRB)   ,INTENT(IN)    :: PCCO2 ! CO2 mass mixing ratio
55REAL(KIND=JPRB)   ,INTENT(IN)    :: POZN(KLON,KLEV) ! O3 mass mixing ratio
56REAL(KIND=JPRB)   ,INTENT(IN)    :: PCLDF(KLON,KLEV) ! Cloud fraction
57REAL(KIND=JPRB)   ,INTENT(IN)    :: PTAUCLD(KLON,KLEV,JPBAND) ! Cloud optical depth
58!--C.Kleinschmitt
59REAL(KIND=JPRB)   ,INTENT(IN)    :: PTAU_LW(KLON,KLEV,NLW) ! LW Optical depth of aerosols 
60!--end
61REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTCLEAR
62REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_CLDFRAC(JPLAY) ! Cloud fraction
63REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TAUCLD(JPLAY,JPBAND) ! Spectral optical thickness
64REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_COLDRY(JPLAY)
65REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_WKL(JPINPX,JPLAY)
66REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_WX(JPXSEC,JPLAY) ! Amount of trace gases
67REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TAUAERL(JPLAY,JPBAND)
68REAL(KIND=JPRB)   ,INTENT(OUT)   :: PAVEL(JPLAY)
69REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TAVEL(JPLAY)
70REAL(KIND=JPRB)   ,INTENT(OUT)   :: PZ(0:JPLAY)
71REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TZ(0:JPLAY)
72REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TBOUND
73INTEGER(KIND=JPIM),INTENT(OUT)   :: K_NLAYERS
74REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_SEMISS(JPBAND)
75INTEGER(KIND=JPIM),INTENT(OUT)   :: K_IREFLECT
76!      real rch4                       ! CH4 mass mixing ratio
77!      real rn2o                       ! N2O mass mixing ratio
78!      real rcfc11                     ! CFC11 mass mixing ratio
79!      real rcfc12                     ! CFC12 mass mixing ratio
80!- from AER
81!- from PROFILE             
82!- from SURFACE             
83REAL(KIND=JPRB) :: ztauaer(5)
84REAL(KIND=JPRB) :: zc1j(0:klev)               ! total cloud from top and level k
85REAL(KIND=JPRB) :: Z_AMD                  ! Effective molecular weight of dry air (g/mol)
86REAL(KIND=JPRB) :: Z_AMW                  ! Molecular weight of water vapor (g/mol)
87REAL(KIND=JPRB) :: Z_AMCO2                ! Molecular weight of carbon dioxide (g/mol)
88REAL(KIND=JPRB) :: Z_AMO                  ! Molecular weight of ozone (g/mol)
89REAL(KIND=JPRB) :: Z_AMCH4                ! Molecular weight of methane (g/mol)
90REAL(KIND=JPRB) :: Z_AMN2O                ! Molecular weight of nitrous oxide (g/mol)
91REAL(KIND=JPRB) :: Z_AMC11                ! Molecular weight of CFC11 (g/mol) - CFCL3
92REAL(KIND=JPRB) :: Z_AMC12                ! Molecular weight of CFC12 (g/mol) - CF2CL2
93REAL(KIND=JPRB) :: Z_AVGDRO               ! Avogadro's number (molecules/mole)
94REAL(KIND=JPRB) :: Z_GRAVIT               ! Gravitational acceleration (cm/sec2)
95
96! Atomic weights for conversion from mass to volume mixing ratios; these
97!  are the same values used in ECRT to assure accurate conversion to vmr
98data Z_AMD   /  28.970_JPRB    /
99data Z_AMW   /  18.0154_JPRB   /
100data Z_AMCO2 /  44.011_JPRB    /
101data Z_AMO   /  47.9982_JPRB   /
102data Z_AMCH4 /  16.043_JPRB    /
103data Z_AMN2O /  44.013_JPRB    /
104data Z_AMC11 / 137.3686_JPRB   /
105data Z_AMC12 / 120.9140_JPRB   /
106data Z_AVGDRO/ 6.02214E23_JPRB /
107data Z_GRAVIT/ 9.80665E02_JPRB /
108
109INTEGER(KIND=JPIM) :: IATM, IMOL, IXMAX, J1, J2, JAE, JB, JK, JL, I_L
110INTEGER(KIND=JPIM) :: I_NMOL, I_NXMOL
111
112REAL(KIND=JPRB) :: Z_AMM, ZCLDLY, ZCLEAR, ZCLOUD, ZEPSEC
113REAL(KIND=JPRB) :: ZHOOK_HANDLE
114
115!MPL/IM 20160915 on prend GES de phylmd
116#include "clesphys.h"
117! ***
118
119! *** mji
120! Initialize all molecular amounts and aerosol optical depths to zero here,
121! then pass ECRT amounts into RRTM arrays below.
122
123!      DATA ZWKL /MAXPRDW*0.0/
124!      DATA ZWX  /MAXPROD*0.0/
125!      DATA KREFLECT /0/
126
127! Activate cross section molecules:
128!     NXMOL     - number of cross-sections input by user
129!     IXINDX(I) - index of cross-section molecule corresponding to Ith
130!                 cross-section specified by user
131!                 = 0 -- not allowed in RRTM
132!                 = 1 -- CCL4
133!                 = 2 -- CFC11
134!                 = 3 -- CFC12
135!                 = 4 -- CFC22
136!      DATA KXMOL  /2/
137!      DATA KXINDX /0,2,3,0,31*0/
138
139!      IREFLECT=KREFLECT
140!      NXMOL=KXMOL
141
142IF (LHOOK) CALL DR_HOOK('RRTM_ECRT_140GP',0,ZHOOK_HANDLE)
143K_IREFLECT=0
144I_NXMOL=2
145
146DO J1=1,35
147! IXINDX(J1)=0
148  DO J2=1,KLEV
149    P_WKL(J1,J2)=0.0_JPRB
150  ENDDO
151ENDDO
152!IXINDX(2)=2
153!IXINDX(3)=3
154
155!     Set parameters needed for RRTM execution:
156IATM    = 0
157!      IXSECT  = 1
158!      NUMANGS = 0
159!      IOUT    = -1
160IXMAX   = 4
161
162!     Bands 6,7,8 are considered the 'window' and allowed to have a
163!     different surface emissivity (as in ECMWF).  Eli wrote this part....
164P_SEMISS(1)  = P_ZEMIS(K_IPLON)
165P_SEMISS(2)  = P_ZEMIS(K_IPLON)
166P_SEMISS(3)  = P_ZEMIS(K_IPLON)
167P_SEMISS(4)  = P_ZEMIS(K_IPLON)
168P_SEMISS(5)  = P_ZEMIS(K_IPLON)
169P_SEMISS(6)  = P_ZEMIW(K_IPLON)
170P_SEMISS(7)  = P_ZEMIW(K_IPLON)
171P_SEMISS(8)  = P_ZEMIW(K_IPLON)
172P_SEMISS(9)  = P_ZEMIS(K_IPLON)
173P_SEMISS(10) = P_ZEMIS(K_IPLON)
174P_SEMISS(11) = P_ZEMIS(K_IPLON)
175P_SEMISS(12) = P_ZEMIS(K_IPLON)
176P_SEMISS(13) = P_ZEMIS(K_IPLON)
177P_SEMISS(14) = P_ZEMIS(K_IPLON)
178P_SEMISS(15) = P_ZEMIS(K_IPLON)
179P_SEMISS(16) = P_ZEMIS(K_IPLON)
180
181!     Set surface temperature. 
182
183P_TBOUND = pts(K_IPLON)
184
185!     Install ECRT arrays into RRTM arrays for pressure, temperature,
186!     and molecular amounts.  Pressures are converted from Pascals
187!     (ECRT) to mb (RRTM).  H2O, CO2, O3 and trace gas amounts are
188!     converted from mass mixing ratio to volume mixing ratio.  CO2
189!     converted with same dry air and CO2 molecular weights used in
190!     ECRT to assure correct conversion back to the proper CO2 vmr.
191!     The dry air column COLDRY (in molec/cm2) is calculated from
192!     the level pressures PZ (in mb) based on the hydrostatic equation
193!     and includes a correction to account for H2O in the layer.  The
194!     molecular weight of moist air (amm) is calculated for each layer.
195!     Note: RRTM levels count from bottom to top, while the ECRT input
196!     variables count from the top down and must be reversed here.
197
198K_NLAYERS = klev
199I_NMOL = 6
200PZ(0) = paph(K_IPLON,klev+1)/100._JPRB
201P_TZ(0) = pth(K_IPLON,klev+1)
202DO I_L = 1, KLEV
203  PAVEL(I_L) = pap(K_IPLON,KLEV-I_L+1)/100._JPRB
204  P_TAVEL(I_L) = pt(K_IPLON,KLEV-I_L+1)
205  PZ(I_L) = paph(K_IPLON,KLEV-I_L+1)/100._JPRB
206  P_TZ(I_L) = pth(K_IPLON,KLEV-I_L+1)
207  P_WKL(1,I_L) = pq(K_IPLON,KLEV-I_L+1)*Z_AMD/Z_AMW
208  P_WKL(2,I_L) = pcco2*Z_AMD/Z_AMCO2
209  P_WKL(3,I_L) = pozn(K_IPLON,KLEV-I_L+1)*Z_AMD/Z_AMO
210  P_WKL(4,I_L) = rn2o*Z_AMD/Z_AMN2O
211  P_WKL(6,I_L) = rch4*Z_AMD/Z_AMCH4
212  Z_AMM = (1-P_WKL(1,I_L))*Z_AMD + P_WKL(1,I_L)*Z_AMW
213  P_COLDRY(I_L) = (PZ(I_L-1)-PZ(I_L))*1.E3_JPRB*Z_AVGDRO/(Z_GRAVIT*Z_AMM*(1+P_WKL(1,I_L)))
214ENDDO
215
216!- Fill RRTM aerosol arrays with operational ECMWF aerosols,
217!  do the mixing and distribute over the 16 spectral intervals
218
219DO I_L=1,KLEV
220  JK=KLEV-I_L+1
221!       DO JAE=1,5
222  JAE=1
223  ZTAUAER(JAE) =&
224   & RAER(JAE,1)*PAER(K_IPLON,1,JK)+RAER(JAE,2)*PAER(K_IPLON,2,JK)&
225   & +RAER(JAE,3)*PAER(K_IPLON,3,JK)+RAER(JAE,4)*PAER(K_IPLON,4,JK)&
226   & +RAER(JAE,5)*PAER(K_IPLON,5,JK)+RAER(JAE,6)*PAER(K_IPLON,6,JK) 
227  P_TAUAERL(I_L, 1)=ZTAUAER(1)
228  P_TAUAERL(I_L, 2)=ZTAUAER(1)
229  JAE=2
230  ZTAUAER(JAE) =&
231   & RAER(JAE,1)*PAER(K_IPLON,1,JK)+RAER(JAE,2)*PAER(K_IPLON,2,JK)&
232   & +RAER(JAE,3)*PAER(K_IPLON,3,JK)+RAER(JAE,4)*PAER(K_IPLON,4,JK)&
233   & +RAER(JAE,5)*PAER(K_IPLON,5,JK)+RAER(JAE,6)*PAER(K_IPLON,6,JK) 
234  P_TAUAERL(I_L, 3)=ZTAUAER(2)
235  P_TAUAERL(I_L, 4)=ZTAUAER(2)
236  P_TAUAERL(I_L, 5)=ZTAUAER(2)
237  JAE=3
238  ZTAUAER(JAE) =&
239   & RAER(JAE,1)*PAER(K_IPLON,1,JK)+RAER(JAE,2)*PAER(K_IPLON,2,JK)&
240   & +RAER(JAE,3)*PAER(K_IPLON,3,JK)+RAER(JAE,4)*PAER(K_IPLON,4,JK)&
241   & +RAER(JAE,5)*PAER(K_IPLON,5,JK)+RAER(JAE,6)*PAER(K_IPLON,6,JK) 
242  P_TAUAERL(I_L, 6)=ZTAUAER(3)
243  P_TAUAERL(I_L, 8)=ZTAUAER(3)
244  P_TAUAERL(I_L, 9)=ZTAUAER(3)
245  JAE=4
246  ZTAUAER(JAE) =&
247   & RAER(JAE,1)*PAER(K_IPLON,1,JK)+RAER(JAE,2)*PAER(K_IPLON,2,JK)&
248   & +RAER(JAE,3)*PAER(K_IPLON,3,JK)+RAER(JAE,4)*PAER(K_IPLON,4,JK)&
249   & +RAER(JAE,5)*PAER(K_IPLON,5,JK)+RAER(JAE,6)*PAER(K_IPLON,6,JK) 
250  P_TAUAERL(I_L, 7)=ZTAUAER(4)
251  JAE=5
252  ZTAUAER(JAE) =&
253   & RAER(JAE,1)*PAER(K_IPLON,1,JK)+RAER(JAE,2)*PAER(K_IPLON,2,JK)&
254   & +RAER(JAE,3)*PAER(K_IPLON,3,JK)+RAER(JAE,4)*PAER(K_IPLON,4,JK)&
255   & +RAER(JAE,5)*PAER(K_IPLON,5,JK)+RAER(JAE,6)*PAER(K_IPLON,6,JK) 
256!       END DO
257  P_TAUAERL(I_L,10)=ZTAUAER(5)
258  P_TAUAERL(I_L,11)=ZTAUAER(5)
259  P_TAUAERL(I_L,12)=ZTAUAER(5)
260  P_TAUAERL(I_L,13)=ZTAUAER(5)
261  P_TAUAERL(I_L,14)=ZTAUAER(5)
262  P_TAUAERL(I_L,15)=ZTAUAER(5)
263  P_TAUAERL(I_L,16)=ZTAUAER(5)
264ENDDO
265!--Use LW AOD from own Mie calculations (C. Kleinschmitt)
266DO I_L=1,KLEV
267  JK=KLEV-I_L+1
268  DO JAE=1, NLW
269    P_TAUAERL(I_L,JAE) = MAX( PTAU_LW(K_IPLON, JK, JAE), 1e-30 )
270  ENDDO
271ENDDO
272!--end C. Kleinschmitt
273
274DO J2=1,KLEV
275  DO J1=1,JPXSEC
276    P_WX(J1,J2)=0.0_JPRB
277  ENDDO
278ENDDO
279
280DO I_L = 1, KLEV
281!- Set cross section molecule amounts from ECRT; convert to vmr
282  P_WX(2,I_L) = rcfc11*Z_AMD/Z_AMC11
283  P_WX(3,I_L) = rcfc12*Z_AMD/Z_AMC12
284  P_WX(2,I_L) = P_COLDRY(I_L) * P_WX(2,I_L) * 1.E-20_JPRB
285  P_WX(3,I_L) = P_COLDRY(I_L) * P_WX(3,I_L) * 1.E-20_JPRB
286
287!- Here, all molecules in WKL and WX are in volume mixing ratio; convert to
288!  molec/cm2 based on COLDRY for use in RRTM
289
290  DO IMOL = 1, I_NMOL
291    P_WKL(IMOL,I_L) = P_COLDRY(I_L) * P_WKL(IMOL,I_L)
292  ENDDO 
293 
294! DO IX = 1,JPXSEC
295! IF (IXINDX(IX)  /=  0) THEN
296!     WX(IXINDX(IX),L) = COLDRY(L) * WX(IX,L) * 1.E-20_JPRB
297! ENDIF
298! END DO 
299
300ENDDO
301
302!- Approximate treatment for various cloud overlaps
303ZCLEAR=1.0_JPRB
304ZCLOUD=0.0_JPRB
305ZC1J(0)=0.0_JPRB
306ZEPSEC=1.E-03_JPRB
307JL=K_IPLON
308
309!++MODIFCODE
310IF ((NOVLP == 1).OR.(NOVLP ==6).OR.(NOVLP ==8)) THEN
311!--MODIFCODE
312
313  DO JK=1,KLEV
314    IF (pcldf(JL,JK) > ZEPSEC) THEN
315      ZCLDLY=pcldf(JL,JK)
316      ZCLEAR=ZCLEAR &
317       & *(1.0_JPRB-MAX( ZCLDLY , ZCLOUD ))&
318       & /(1.0_JPRB-MIN( ZCLOUD , 1.0_JPRB-ZEPSEC )) 
319      ZCLOUD = ZCLDLY
320      ZC1J(JK)= 1.0_JPRB - ZCLEAR
321    ELSE
322      ZCLDLY=0.0_JPRB
323      ZCLEAR=ZCLEAR &
324       & *(1.0_JPRB-MAX( ZCLDLY , ZCLOUD ))&
325       & /(1.0_JPRB-MIN( ZCLOUD , 1.0_JPRB-ZEPSEC )) 
326      ZCLOUD = ZCLDLY
327      ZC1J(JK)= 1.0_JPRB - ZCLEAR
328    ENDIF
329  ENDDO
330
331!++MODIFCODE
332ELSEIF ((NOVLP == 2).OR.(NOVLP ==7)) THEN
333!--MODIFCODE
334
335  DO JK=1,KLEV
336    IF (pcldf(JL,JK) > ZEPSEC) THEN
337      ZCLDLY=pcldf(JL,JK)
338      ZCLOUD = MAX( ZCLDLY , ZCLOUD )
339      ZC1J(JK) = ZCLOUD
340    ELSE
341      ZCLDLY=0.0_JPRB
342      ZCLOUD = MAX( ZCLDLY , ZCLOUD )
343      ZC1J(JK) = ZCLOUD
344    ENDIF
345  ENDDO
346
347!++MODIFCODE
348ELSEIF ((NOVLP == 3).OR.(NOVLP ==5)) THEN
349!--MODIFCODE
350
351  DO JK=1,KLEV
352    IF (pcldf(JL,JK) > ZEPSEC) THEN
353      ZCLDLY=pcldf(JL,JK)
354      ZCLEAR = ZCLEAR * (1.0_JPRB-ZCLDLY)
355      ZCLOUD = 1.0_JPRB - ZCLEAR
356      ZC1J(JK) = ZCLOUD
357    ELSE
358      ZCLDLY=0.0_JPRB
359      ZCLEAR = ZCLEAR * (1.0_JPRB-ZCLDLY)
360      ZCLOUD = 1.0_JPRB - ZCLEAR
361      ZC1J(JK) = ZCLOUD
362    ENDIF
363  ENDDO
364
365ELSEIF (NOVLP == 4) THEN
366
367ENDIF
368PTCLEAR=1.0_JPRB-ZC1J(KLEV)
369
370! Transfer cloud fraction and cloud optical depth to RRTM arrays;
371! invert array index for pcldf to go from bottom to top for RRTM
372
373!- clear-sky column
374IF (PTCLEAR  >  1.0_JPRB-ZEPSEC) THEN
375  KCLD=0
376  DO I_L = 1, KLEV
377    P_CLDFRAC(I_L) = 0.0_JPRB
378  ENDDO
379  DO JB=1,JPBAND
380    DO I_L=1,KLEV
381      P_TAUCLD(I_L,JB) = 0.0_JPRB
382    ENDDO
383  ENDDO
384
385ELSE
386
387!- cloudy column
388!   The diffusivity factor (Savijarvi, 1997) on the cloud optical
389!   thickness TAUCLD has already been applied in RADLSW
390
391  KCLD=1
392  DO I_L=1,KLEV
393    P_CLDFRAC(I_L) = pcldf(K_IPLON,I_L)
394  ENDDO
395  DO JB=1,JPBAND
396    DO I_L=1,KLEV
397      P_TAUCLD(I_L,JB) = ptaucld(K_IPLON,I_L,JB)
398    ENDDO
399  ENDDO
400
401ENDIF
402
403!     ------------------------------------------------------------------
404 
405IF (LHOOK) CALL DR_HOOK('RRTM_ECRT_140GP',1,ZHOOK_HANDLE)
406END SUBROUTINE RRTM_ECRT_140GP
Note: See TracBrowser for help on using the repository browser.