source: LMDZ5/branches/IPSLCM5A2.1/libf/phylmd/rrtm/rrtm_ecrt_140gp.F90 @ 5421

Last change on this file since 5421 was 2152, checked in by fhourdin, 10 years ago

Corrections cosmétiques pour RRTM
Bug fixing for RRTM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:keywords set to Author Date Id Revi
File size: 13.2 KB
Line 
1!
2! $Id: rrtm_ecrt_140gp.F90 2152 2014-11-19 16:52:35Z fhourdin $
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
33USE 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! ***
116
117! *** mji
118! Initialize all molecular amounts and aerosol optical depths to zero here,
119! then pass ECRT amounts into RRTM arrays below.
120
121!      DATA ZWKL /MAXPRDW*0.0/
122!      DATA ZWX  /MAXPROD*0.0/
123!      DATA KREFLECT /0/
124
125! Activate cross section molecules:
126!     NXMOL     - number of cross-sections input by user
127!     IXINDX(I) - index of cross-section molecule corresponding to Ith
128!                 cross-section specified by user
129!                 = 0 -- not allowed in RRTM
130!                 = 1 -- CCL4
131!                 = 2 -- CFC11
132!                 = 3 -- CFC12
133!                 = 4 -- CFC22
134!      DATA KXMOL  /2/
135!      DATA KXINDX /0,2,3,0,31*0/
136
137!      IREFLECT=KREFLECT
138!      NXMOL=KXMOL
139
140IF (LHOOK) CALL DR_HOOK('RRTM_ECRT_140GP',0,ZHOOK_HANDLE)
141K_IREFLECT=0
142I_NXMOL=2
143
144DO J1=1,35
145! IXINDX(J1)=0
146  DO J2=1,KLEV
147    P_WKL(J1,J2)=0.0_JPRB
148  ENDDO
149ENDDO
150!IXINDX(2)=2
151!IXINDX(3)=3
152
153!     Set parameters needed for RRTM execution:
154IATM    = 0
155!      IXSECT  = 1
156!      NUMANGS = 0
157!      IOUT    = -1
158IXMAX   = 4
159
160!     Bands 6,7,8 are considered the 'window' and allowed to have a
161!     different surface emissivity (as in ECMWF).  Eli wrote this part....
162P_SEMISS(1)  = P_ZEMIS(K_IPLON)
163P_SEMISS(2)  = P_ZEMIS(K_IPLON)
164P_SEMISS(3)  = P_ZEMIS(K_IPLON)
165P_SEMISS(4)  = P_ZEMIS(K_IPLON)
166P_SEMISS(5)  = P_ZEMIS(K_IPLON)
167P_SEMISS(6)  = P_ZEMIW(K_IPLON)
168P_SEMISS(7)  = P_ZEMIW(K_IPLON)
169P_SEMISS(8)  = P_ZEMIW(K_IPLON)
170P_SEMISS(9)  = P_ZEMIS(K_IPLON)
171P_SEMISS(10) = P_ZEMIS(K_IPLON)
172P_SEMISS(11) = P_ZEMIS(K_IPLON)
173P_SEMISS(12) = P_ZEMIS(K_IPLON)
174P_SEMISS(13) = P_ZEMIS(K_IPLON)
175P_SEMISS(14) = P_ZEMIS(K_IPLON)
176P_SEMISS(15) = P_ZEMIS(K_IPLON)
177P_SEMISS(16) = P_ZEMIS(K_IPLON)
178
179!     Set surface temperature. 
180
181P_TBOUND = pts(K_IPLON)
182
183!     Install ECRT arrays into RRTM arrays for pressure, temperature,
184!     and molecular amounts.  Pressures are converted from Pascals
185!     (ECRT) to mb (RRTM).  H2O, CO2, O3 and trace gas amounts are
186!     converted from mass mixing ratio to volume mixing ratio.  CO2
187!     converted with same dry air and CO2 molecular weights used in
188!     ECRT to assure correct conversion back to the proper CO2 vmr.
189!     The dry air column COLDRY (in molec/cm2) is calculated from
190!     the level pressures PZ (in mb) based on the hydrostatic equation
191!     and includes a correction to account for H2O in the layer.  The
192!     molecular weight of moist air (amm) is calculated for each layer.
193!     Note: RRTM levels count from bottom to top, while the ECRT input
194!     variables count from the top down and must be reversed here.
195
196K_NLAYERS = klev
197I_NMOL = 6
198PZ(0) = paph(K_IPLON,klev+1)/100._JPRB
199P_TZ(0) = pth(K_IPLON,klev+1)
200DO I_L = 1, KLEV
201  PAVEL(I_L) = pap(K_IPLON,KLEV-I_L+1)/100._JPRB
202  P_TAVEL(I_L) = pt(K_IPLON,KLEV-I_L+1)
203  PZ(I_L) = paph(K_IPLON,KLEV-I_L+1)/100._JPRB
204  P_TZ(I_L) = pth(K_IPLON,KLEV-I_L+1)
205  P_WKL(1,I_L) = pq(K_IPLON,KLEV-I_L+1)*Z_AMD/Z_AMW
206  P_WKL(2,I_L) = pcco2*Z_AMD/Z_AMCO2
207  P_WKL(3,I_L) = pozn(K_IPLON,KLEV-I_L+1)*Z_AMD/Z_AMO
208  P_WKL(4,I_L) = rn2o*Z_AMD/Z_AMN2O
209  P_WKL(6,I_L) = rch4*Z_AMD/Z_AMCH4
210  Z_AMM = (1-P_WKL(1,I_L))*Z_AMD + P_WKL(1,I_L)*Z_AMW
211  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)))
212ENDDO
213
214!- Fill RRTM aerosol arrays with operational ECMWF aerosols,
215!  do the mixing and distribute over the 16 spectral intervals
216
217DO I_L=1,KLEV
218  JK=KLEV-I_L+1
219!       DO JAE=1,5
220  JAE=1
221  ZTAUAER(JAE) =&
222   & RAER(JAE,1)*PAER(K_IPLON,1,JK)+RAER(JAE,2)*PAER(K_IPLON,2,JK)&
223   & +RAER(JAE,3)*PAER(K_IPLON,3,JK)+RAER(JAE,4)*PAER(K_IPLON,4,JK)&
224   & +RAER(JAE,5)*PAER(K_IPLON,5,JK)+RAER(JAE,6)*PAER(K_IPLON,6,JK) 
225  P_TAUAERL(I_L, 1)=ZTAUAER(1)
226  P_TAUAERL(I_L, 2)=ZTAUAER(1)
227  JAE=2
228  ZTAUAER(JAE) =&
229   & RAER(JAE,1)*PAER(K_IPLON,1,JK)+RAER(JAE,2)*PAER(K_IPLON,2,JK)&
230   & +RAER(JAE,3)*PAER(K_IPLON,3,JK)+RAER(JAE,4)*PAER(K_IPLON,4,JK)&
231   & +RAER(JAE,5)*PAER(K_IPLON,5,JK)+RAER(JAE,6)*PAER(K_IPLON,6,JK) 
232  P_TAUAERL(I_L, 3)=ZTAUAER(2)
233  P_TAUAERL(I_L, 4)=ZTAUAER(2)
234  P_TAUAERL(I_L, 5)=ZTAUAER(2)
235  JAE=3
236  ZTAUAER(JAE) =&
237   & RAER(JAE,1)*PAER(K_IPLON,1,JK)+RAER(JAE,2)*PAER(K_IPLON,2,JK)&
238   & +RAER(JAE,3)*PAER(K_IPLON,3,JK)+RAER(JAE,4)*PAER(K_IPLON,4,JK)&
239   & +RAER(JAE,5)*PAER(K_IPLON,5,JK)+RAER(JAE,6)*PAER(K_IPLON,6,JK) 
240  P_TAUAERL(I_L, 6)=ZTAUAER(3)
241  P_TAUAERL(I_L, 8)=ZTAUAER(3)
242  P_TAUAERL(I_L, 9)=ZTAUAER(3)
243  JAE=4
244  ZTAUAER(JAE) =&
245   & RAER(JAE,1)*PAER(K_IPLON,1,JK)+RAER(JAE,2)*PAER(K_IPLON,2,JK)&
246   & +RAER(JAE,3)*PAER(K_IPLON,3,JK)+RAER(JAE,4)*PAER(K_IPLON,4,JK)&
247   & +RAER(JAE,5)*PAER(K_IPLON,5,JK)+RAER(JAE,6)*PAER(K_IPLON,6,JK) 
248  P_TAUAERL(I_L, 7)=ZTAUAER(4)
249  JAE=5
250  ZTAUAER(JAE) =&
251   & RAER(JAE,1)*PAER(K_IPLON,1,JK)+RAER(JAE,2)*PAER(K_IPLON,2,JK)&
252   & +RAER(JAE,3)*PAER(K_IPLON,3,JK)+RAER(JAE,4)*PAER(K_IPLON,4,JK)&
253   & +RAER(JAE,5)*PAER(K_IPLON,5,JK)+RAER(JAE,6)*PAER(K_IPLON,6,JK) 
254!       END DO
255  P_TAUAERL(I_L,10)=ZTAUAER(5)
256  P_TAUAERL(I_L,11)=ZTAUAER(5)
257  P_TAUAERL(I_L,12)=ZTAUAER(5)
258  P_TAUAERL(I_L,13)=ZTAUAER(5)
259  P_TAUAERL(I_L,14)=ZTAUAER(5)
260  P_TAUAERL(I_L,15)=ZTAUAER(5)
261  P_TAUAERL(I_L,16)=ZTAUAER(5)
262ENDDO
263!--Use LW AOD from own Mie calculations (C. Kleinschmitt)
264DO I_L=1,KLEV
265  JK=KLEV-I_L+1
266  DO JAE=1, NLW
267    P_TAUAERL(I_L,JAE) = MAX( PTAU_LW(K_IPLON, JK, JAE), 1e-30 )
268  ENDDO
269ENDDO
270!--end C. Kleinschmitt
271
272DO J2=1,KLEV
273  DO J1=1,JPXSEC
274    P_WX(J1,J2)=0.0_JPRB
275  ENDDO
276ENDDO
277
278DO I_L = 1, KLEV
279!- Set cross section molecule amounts from ECRT; convert to vmr
280  P_WX(2,I_L) = rcfc11*Z_AMD/Z_AMC11
281  P_WX(3,I_L) = rcfc12*Z_AMD/Z_AMC12
282  P_WX(2,I_L) = P_COLDRY(I_L) * P_WX(2,I_L) * 1.E-20_JPRB
283  P_WX(3,I_L) = P_COLDRY(I_L) * P_WX(3,I_L) * 1.E-20_JPRB
284
285!- Here, all molecules in WKL and WX are in volume mixing ratio; convert to
286!  molec/cm2 based on COLDRY for use in RRTM
287
288  DO IMOL = 1, I_NMOL
289    P_WKL(IMOL,I_L) = P_COLDRY(I_L) * P_WKL(IMOL,I_L)
290  ENDDO 
291 
292! DO IX = 1,JPXSEC
293! IF (IXINDX(IX)  /=  0) THEN
294!     WX(IXINDX(IX),L) = COLDRY(L) * WX(IX,L) * 1.E-20_JPRB
295! ENDIF
296! END DO 
297
298ENDDO
299
300!- Approximate treatment for various cloud overlaps
301ZCLEAR=1.0_JPRB
302ZCLOUD=0.0_JPRB
303ZC1J(0)=0.0_JPRB
304ZEPSEC=1.E-03_JPRB
305JL=K_IPLON
306
307!++MODIFCODE
308IF ((NOVLP == 1).OR.(NOVLP ==6).OR.(NOVLP ==8)) THEN
309!--MODIFCODE
310
311  DO JK=1,KLEV
312    IF (pcldf(JL,JK) > ZEPSEC) THEN
313      ZCLDLY=pcldf(JL,JK)
314      ZCLEAR=ZCLEAR &
315       & *(1.0_JPRB-MAX( ZCLDLY , ZCLOUD ))&
316       & /(1.0_JPRB-MIN( ZCLOUD , 1.0_JPRB-ZEPSEC )) 
317      ZCLOUD = ZCLDLY
318      ZC1J(JK)= 1.0_JPRB - ZCLEAR
319    ELSE
320      ZCLDLY=0.0_JPRB
321      ZCLEAR=ZCLEAR &
322       & *(1.0_JPRB-MAX( ZCLDLY , ZCLOUD ))&
323       & /(1.0_JPRB-MIN( ZCLOUD , 1.0_JPRB-ZEPSEC )) 
324      ZCLOUD = ZCLDLY
325      ZC1J(JK)= 1.0_JPRB - ZCLEAR
326    ENDIF
327  ENDDO
328
329!++MODIFCODE
330ELSEIF ((NOVLP == 2).OR.(NOVLP ==7)) THEN
331!--MODIFCODE
332
333  DO JK=1,KLEV
334    IF (pcldf(JL,JK) > ZEPSEC) THEN
335      ZCLDLY=pcldf(JL,JK)
336      ZCLOUD = MAX( ZCLDLY , ZCLOUD )
337      ZC1J(JK) = ZCLOUD
338    ELSE
339      ZCLDLY=0.0_JPRB
340      ZCLOUD = MAX( ZCLDLY , ZCLOUD )
341      ZC1J(JK) = ZCLOUD
342    ENDIF
343  ENDDO
344
345!++MODIFCODE
346ELSEIF ((NOVLP == 3).OR.(NOVLP ==5)) THEN
347!--MODIFCODE
348
349  DO JK=1,KLEV
350    IF (pcldf(JL,JK) > ZEPSEC) THEN
351      ZCLDLY=pcldf(JL,JK)
352      ZCLEAR = ZCLEAR * (1.0_JPRB-ZCLDLY)
353      ZCLOUD = 1.0_JPRB - ZCLEAR
354      ZC1J(JK) = ZCLOUD
355    ELSE
356      ZCLDLY=0.0_JPRB
357      ZCLEAR = ZCLEAR * (1.0_JPRB-ZCLDLY)
358      ZCLOUD = 1.0_JPRB - ZCLEAR
359      ZC1J(JK) = ZCLOUD
360    ENDIF
361  ENDDO
362
363ELSEIF (NOVLP == 4) THEN
364
365ENDIF
366PTCLEAR=1.0_JPRB-ZC1J(KLEV)
367
368! Transfer cloud fraction and cloud optical depth to RRTM arrays;
369! invert array index for pcldf to go from bottom to top for RRTM
370
371!- clear-sky column
372IF (PTCLEAR  >  1.0_JPRB-ZEPSEC) THEN
373  KCLD=0
374  DO I_L = 1, KLEV
375    P_CLDFRAC(I_L) = 0.0_JPRB
376  ENDDO
377  DO JB=1,JPBAND
378    DO I_L=1,KLEV
379      P_TAUCLD(I_L,JB) = 0.0_JPRB
380    ENDDO
381  ENDDO
382
383ELSE
384
385!- cloudy column
386!   The diffusivity factor (Savijarvi, 1997) on the cloud optical
387!   thickness TAUCLD has already been applied in RADLSW
388
389  KCLD=1
390  DO I_L=1,KLEV
391    P_CLDFRAC(I_L) = pcldf(K_IPLON,I_L)
392  ENDDO
393  DO JB=1,JPBAND
394    DO I_L=1,KLEV
395      P_TAUCLD(I_L,JB) = ptaucld(K_IPLON,I_L,JB)
396    ENDDO
397  ENDDO
398
399ENDIF
400
401!     ------------------------------------------------------------------
402
403IF (LHOOK) CALL DR_HOOK('RRTM_ECRT_140GP',1,ZHOOK_HANDLE)
404END SUBROUTINE RRTM_ECRT_140GP
Note: See TracBrowser for help on using the repository browser.