source: LMDZ6/branches/IPSLCM6.0.14/libf/phymar/rrtm_ecrt_140gp.F90 @ 5420

Last change on this file since 5420 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR

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