source: LMDZ6/trunk/libf/phylmd/ecrad/rrtm_rrtm_140gp_mcica.F90.or @ 4062

Last change on this file since 4062 was 3908, checked in by idelkadi, 3 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

  • Inclusion of the ecrad directory containing the sources of the ECRAD code
    • interface routine : radiation_scheme.F90
  • Adaptation of compilation scripts :
    • compilation under CPP key CPP_ECRAD
    • compilation with option "-rad ecard" or "-ecard true"
    • The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
  • Runing LMDZ simulations with ecrad, you need :
    • logical key iflag_rrtm = 2 in physiq.def
    • namelist_ecrad (DefLists?)
    • the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
  • Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...
File size: 12.3 KB
Line 
1SUBROUTINE RRTM_RRTM_140GP_MCICA &
2 &( KIDIA , KFDIA , KLON , KLEV, KCOLS, KCLDCOL,&
3 &  PAER  , PAPH  , PAP  , PAERTAUL, PAERASYL, PAEROMGL, &
4 &  PTS   , PTH   , PT   , &
5 &  PEMIS , PEMIW ,&
6 &  PQ    , PCO2  , PCH4 , PN2O, PNO2 , PC11, PC12, PC22, PCL4, POZN   ,&
7 &  PCLDF , PTAUCLD, PCLFR,&
8 &  PEMIT , PFLUX , PFLUC, &
9 &  PLwDerivative) 
10
11! *** This program is the driver for the McICA version of RRTM_LW,
12!     the AER rapid model. 
13
14!     For each atmosphere the user wishes to analyze, this routine
15!     a) calls ECRTATM to read in the atmospheric profile
16!     b) calls SETCOEF to calculate various quantities needed for
17!        the radiative transfer algorithm
18!     c) calls RTRN to do the radiative transfer calculation for
19!        clear or cloudy sky
20!     d) writes out the upward, downward, and net flux for each
21!        level and the heating rate for each layer
22
23!     JJMorcrette 20050110 McICA version revisited (changes in RRTM_ECRT, RRTM_RTRN)
24!        NEC           25-Oct-2007 Optimisations
25!     JJMorcrette 20080424 3D fields of CO2, CH4, N2O, NO2, CFC11, 12, 22 and CCL4
26!     JJMorcrette 20110613 flexible number of g-points
27!     P Bechtold 14/05/2012 replace ZHEATF by core constants RG*RDAY/RCPD
28!                           and put arrays to scalars
29!     R Hogan    20/05/2014 pass partial derivatives back to calling function
30!-----------------------------------------------------------------------
31
32USE PARKIND1 , ONLY : JPIM, JPRB
33USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
34USE PARRRTM  , ONLY : JPBAND, JPXSEC, JPINPX
35USE YOERRTM  , ONLY : JPGPT 
36USE YOMCST   , ONLY : RG  ! , RDAYI, RCPD
37
38IMPLICIT NONE
39
40!------------------------------Arguments--------------------------------
41
42! Input arguments
43
44INTEGER(KIND=JPIM),INTENT(IN)    :: KLON! Number of atmospheres (longitudes)
45INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV! Number of atmospheric layers
46INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA ! First atmosphere index
47INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA ! Last atmosphere index
48INTEGER(KIND=JPIM),INTENT(IN)    :: KCOLS ! Number of columns on which to perform RT
49                                          ! should be the same as number of g-points, JPGPT
50INTEGER(KIND=JPIM),INTENT(IN)    :: KCLDCOL(KLON) ! cloudy column index: 1=cloud, 0: clear   
51
52REAL(KIND=JPRB)   ,INTENT(IN)    :: PAER(KLON,6,KLEV) ! Aerosol optical thickness
53REAL(KIND=JPRB)   ,INTENT(IN)    :: PAERTAUL(KLON,KLEV,16), PAERASYL(KLON,KLEV,16), PAEROMGL(KLON,KLEV,16)
54REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPH(KLON,KLEV+1) ! Interface pressures (Pa)
55REAL(KIND=JPRB)   ,INTENT(IN)    :: PAP(KLON,KLEV) ! Layer pressures (Pa)
56REAL(KIND=JPRB)   ,INTENT(IN)    :: PTS(KLON) ! Surface temperature (JK)
57REAL(KIND=JPRB)   ,INTENT(IN)    :: PTH(KLON,KLEV+1) ! Interface temperatures (JK)
58REAL(KIND=JPRB)   ,INTENT(IN)    :: PT(KLON,KLEV) ! Layer temperature (JK)
59REAL(KIND=JPRB)   ,INTENT(IN)    :: PEMIS(KLON) ! Non-window surface emissivity
60REAL(KIND=JPRB)   ,INTENT(IN)    :: PEMIW(KLON) ! Window surface emissivity
61REAL(KIND=JPRB)   ,INTENT(IN)    :: PQ(KLON,KLEV) ! H2O specific humidity (mmr)
62REAL(KIND=JPRB)   ,INTENT(IN)    :: PCO2(KLON,KLEV) ! CO2 mass mixing ratio
63REAL(KIND=JPRB)   ,INTENT(IN)    :: PCH4(KLON,KLEV) ! CH4 mass mixing ratio
64REAL(KIND=JPRB)   ,INTENT(IN)    :: PN2O(KLON,KLEV) ! N2O mass mixing ratio
65REAL(KIND=JPRB)   ,INTENT(IN)    :: PNO2(KLON,KLEV) ! NO2 mass mixing ratio
66REAL(KIND=JPRB)   ,INTENT(IN)    :: PC11(KLON,KLEV) ! CFC11 mass mixing ratio
67REAL(KIND=JPRB)   ,INTENT(IN)    :: PC12(KLON,KLEV) ! CFC12 mass mixing ratio
68REAL(KIND=JPRB)   ,INTENT(IN)    :: PC22(KLON,KLEV) ! CFC22 mass mixing ratio
69REAL(KIND=JPRB)   ,INTENT(IN)    :: PCL4(KLON,KLEV) ! CCL4  mass mixing ratio
70REAL(KIND=JPRB)   ,INTENT(IN)    :: POZN(KLON,KLEV) ! O3 mass mixing ratio
71REAL(KIND=JPRB)   ,INTENT(IN)    :: PCLFR(KLON,KLEV)
72
73REAL(KIND=JPRB)   ,INTENT(OUT)   :: PEMIT(KLON) ! Surface LW emissivity
74REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFLUX(KLON,2,KLEV+1) ! LW total sky flux (1=up, 2=down)
75REAL(KIND=JPRB)   ,INTENT(OUT)   :: PFLUC(KLON,2,KLEV+1) ! LW clear sky flux (1=up, 2=down)
76
77! Partial derivative of total upward flux at each level with respect
78! to upward flux at surface, used to correct heating rates at
79! gridpoints/timesteps between calls to the full radiation scheme:
80REAL(KIND=JPRB)   ,INTENT(OUT)   :: PLwDerivative(KLON,KLEV+1)
81
82!-- McICA ----------
83REAL(KIND=JPRB)   ,INTENT(IN)    :: PCLDF(KLON,KCOLS,KLEV)    ! Cloud fraction
84REAL(KIND=JPRB)   ,INTENT(IN)    :: PTAUCLD(KLON,KLEV,KCOLS)  ! Cloud optical depth
85
86REAL(KIND=JPRB) :: ZCLDFRAC(KIDIA:KFDIA,KCOLS,KLEV)   ! Cloud fraction
87REAL(KIND=JPRB) :: ZTAUCLD(KIDIA:KFDIA,KLEV,KCOLS)    ! Spectral optical thickness
88!-- McICA ----------
89
90REAL(KIND=JPRB) :: ZATR1(KIDIA:KFDIA,JPGPT,KLEV)
91
92REAL(KIND=JPRB) :: ZOD(KIDIA:KFDIA,JPGPT,KLEV)
93
94REAL(KIND=JPRB) :: ZTF1(KIDIA:KFDIA,JPGPT,KLEV)
95
96REAL(KIND=JPRB) :: ZCOLDRY(KIDIA:KFDIA,KLEV)
97REAL(KIND=JPRB) :: ZWBRODL(KIDIA:KFDIA,KLEV) !BROADENING GASES,column density (mol/cm2)
98REAL(KIND=JPRB) :: ZCOLBRD(KIDIA:KFDIA,KLEV) !BROADENING GASES, column amount
99REAL(KIND=JPRB) :: ZWKL(KIDIA:KFDIA,JPINPX,KLEV)
100
101REAL(KIND=JPRB) :: ZWX(KIDIA:KFDIA,JPXSEC,KLEV)         ! Amount of trace gases
102
103REAL(KIND=JPRB) :: ZTOTDFLUC(KIDIA:KFDIA,0:KLEV)
104REAL(KIND=JPRB) :: ZTOTDFLUX(KIDIA:KFDIA,0:KLEV)
105REAL(KIND=JPRB) :: ZTOTUFLUC(KIDIA:KFDIA,0:KLEV)
106REAL(KIND=JPRB) :: ZTOTUFLUX(KIDIA:KFDIA,0:KLEV)
107
108INTEGER(KIND=JPIM) :: JL, JK
109INTEGER(KIND=JPIM) :: ISTART
110INTEGER(KIND=JPIM) :: IEND
111
112REAL(KIND=JPRB) :: ZFLUXFAC, ZHEATFAC, ZPI
113REAL(KIND=JPRB) :: ZEPSEC
114
115!- from AER
116REAL(KIND=JPRB) :: ZTAUAERL(KIDIA:KFDIA,KLEV,JPBAND)
117
118!- from INTFAC     
119REAL(KIND=JPRB) :: ZFAC00(KIDIA:KFDIA,KLEV)
120REAL(KIND=JPRB) :: ZFAC01(KIDIA:KFDIA,KLEV)
121REAL(KIND=JPRB) :: ZFAC10(KIDIA:KFDIA,KLEV)
122REAL(KIND=JPRB) :: ZFAC11(KIDIA:KFDIA,KLEV)
123
124!- from FOR
125REAL(KIND=JPRB) :: ZFORFAC(KIDIA:KFDIA,KLEV)
126REAL(KIND=JPRB) :: ZFORFRAC(KIDIA:KFDIA,KLEV)
127INTEGER(KIND=JPIM) :: INDFOR(KIDIA:KFDIA,KLEV)
128
129!- from MINOR
130INTEGER(KIND=JPIM)  :: INDMINOR(KIDIA:KFDIA,KLEV)
131REAL(KIND=JPRB)     :: ZSCALEMINOR(KIDIA:KFDIA,KLEV)
132REAL(KIND=JPRB)     :: ZSCALEMINORN2(KIDIA:KFDIA,KLEV)
133REAL(KIND=JPRB)     :: ZMINORFRAC(KIDIA:KFDIA,KLEV)
134
135REAL(KIND=JPRB)     :: &                 
136                    &  ZRAT_H2OCO2(KIDIA:KFDIA,KLEV),ZRAT_H2OCO2_1(KIDIA:KFDIA,KLEV), &
137                    &  ZRAT_H2OO3(KIDIA:KFDIA,KLEV) ,ZRAT_H2OO3_1(KIDIA:KFDIA,KLEV), &
138                    &  ZRAT_H2ON2O(KIDIA:KFDIA,KLEV),ZRAT_H2ON2O_1(KIDIA:KFDIA,KLEV), &
139                    &  ZRAT_H2OCH4(KIDIA:KFDIA,KLEV),ZRAT_H2OCH4_1(KIDIA:KFDIA,KLEV), &
140                    &  ZRAT_N2OCO2(KIDIA:KFDIA,KLEV),ZRAT_N2OCO2_1(KIDIA:KFDIA,KLEV), &
141                    &  ZRAT_O3CO2(KIDIA:KFDIA,KLEV) ,ZRAT_O3CO2_1(KIDIA:KFDIA,KLEV)
142
143!- from INTIND
144INTEGER(KIND=JPIM) :: JP(KIDIA:KFDIA,KLEV)
145INTEGER(KIND=JPIM) :: JT(KIDIA:KFDIA,KLEV)
146INTEGER(KIND=JPIM) :: JT1(KIDIA:KFDIA,KLEV)
147
148!- from PRECISE             
149REAL(KIND=JPRB) :: ZONEMINUS
150
151!- from PROFDATA             
152REAL(KIND=JPRB) :: ZCOLH2O(KIDIA:KFDIA,KLEV)
153REAL(KIND=JPRB) :: ZCOLCO2(KIDIA:KFDIA,KLEV)
154REAL(KIND=JPRB) :: ZCOLO3(KIDIA:KFDIA,KLEV)
155REAL(KIND=JPRB) :: ZCOLN2O(KIDIA:KFDIA,KLEV)
156REAL(KIND=JPRB) :: ZCOLCH4(KIDIA:KFDIA,KLEV)
157REAL(KIND=JPRB) :: ZCOLO2(KIDIA:KFDIA,KLEV)
158REAL(KIND=JPRB) :: ZCO2MULT(KIDIA:KFDIA,KLEV)
159INTEGER(KIND=JPIM) :: ILAYTROP(KIDIA:KFDIA)
160INTEGER(KIND=JPIM) :: ILAYSWTCH(KIDIA:KFDIA)
161INTEGER(KIND=JPIM) :: ILAYLOW(KIDIA:KFDIA)
162
163!- from PROFILE             
164REAL(KIND=JPRB) :: ZPAVEL(KIDIA:KFDIA,KLEV)
165REAL(KIND=JPRB) :: ZTAVEL(KIDIA:KFDIA,KLEV)
166REAL(KIND=JPRB) :: ZPZ(KIDIA:KFDIA,0:KLEV)
167REAL(KIND=JPRB) :: ZTZ(KIDIA:KFDIA,0:KLEV)
168REAL(KIND=JPRB) :: ZTBOUND(KIDIA:KFDIA)
169
170!- from SELF             
171REAL(KIND=JPRB) :: ZSELFFAC(KIDIA:KFDIA,KLEV)
172REAL(KIND=JPRB) :: ZSELFFRAC(KIDIA:KFDIA,KLEV)
173INTEGER(KIND=JPIM) :: INDSELF(KIDIA:KFDIA,KLEV)
174
175!- from SP             
176REAL(KIND=JPRB) :: ZPFRAC(KIDIA:KFDIA,JPGPT,KLEV)
177
178!- from SURFACE             
179REAL(KIND=JPRB) :: ZSEMISS(KIDIA:KFDIA,JPBAND)
180REAL(KIND=JPRB) :: ZSEMISLW(KIDIA:KFDIA)
181INTEGER(KIND=JPIM) :: IREFLECT(KIDIA:KFDIA)
182
183! Local variable required in case KFDIA /= KLON
184REAL(KIND=JPRB) :: ZLwDerivative(KIDIA:KFDIA,KLEV+1)
185
186LOGICAL            :: LLPRINT
187
188REAL(KIND=JPRB) :: ZHOOK_HANDLE
189
190#include "rrtm_ecrt_140gp_mcica.intfb.h"
191#include "rrtm_gasabs1a_140gp.intfb.h"
192#include "rrtm_rtrn1a_140gp_mcica.intfb.h"
193#include "rrtm_setcoef_140gp.intfb.h"
194
195!     HEATFAC is the factor by which one must multiply delta-flux/
196!     delta-pressure, with flux in w/m-2 and pressure in mbar, to get
197!     the heating rate in units of degrees/day.  It is equal to
198!           (g)x(#sec/day)x(1e-5)/(specific heat of air at const. p)
199!        =  (9.8066)(86400)(1e-5)/(1.004)
200
201IF (LHOOK) CALL DR_HOOK('RRTM_RRTM_140GP_MCICA',0,ZHOOK_HANDLE)
202
203ASSOCIATE(NFLEVG=>KLEV)
204
205ZEPSEC = 1.E-06_JPRB
206ZONEMINUS = 1.0_JPRB - ZEPSEC
207ZPI = 2.0_JPRB*ASIN(1.0_JPRB)
208ZFLUXFAC = ZPI * 2.E+4
209!ZHEATFAC = 8.4391_JPRB
210!ZHEATFAC = RG*RDAYI/RCPD*1.E-2_JPRB
211
212! *** mji ***
213
214! For use with ECRT, this loop is over atmospheres (or longitudes)
215LLPRINT=.TRUE.
216
217!  do JK=1,KLEV
218!    print 9901,JK,PT(JL,JK),PQ(JL,JK),POZN(JL,JK),PCLDF(JL,JK,1),PTAUCLD(JL,JK,1)
219!  enddo
220
221! *** mji ***
222!- Prepare atmospheric profile from ECRT for use in RRTM, and define
223!  other RRTM input parameters.  Arrays are passed back through the
224!  existing RRTM commons and arrays.
225
226  CALL RRTM_ECRT_140GP_MCICA &
227   &( KIDIA, KFDIA, KLON , KLEV, KCOLS , &
228   &  PAER , PAPH , PAP , PAERTAUL, PAERASYL, PAEROMGL, &
229   &  pts  , PTH  , PT  , &
230   &  PEMIS, PEMIW, &
231   &  PQ   , PCO2 , PCH4, PN2O, PNO2, PC11, PC12, PC22, PCL4, POZN , PCLDF, PTAUCLD, &
232   &  ZCLDFRAC, ZTAUCLD, ZCOLDRY, ZWBRODL,ZWKL, ZWX, &
233   &  ZTAUAERL, ZPAVEL , ZTAVEL , ZPZ , ZTZ, ZTBOUND, ZSEMISS, IREFLECT) 
234
235  ISTART = 1
236  IEND   = 16
237
238!  Calculate information needed by the radiative transfer routine
239!  that is specific to this atmosphere, especially some of the
240!  coefficients and indices needed to compute the optical depths
241!  by interpolating data from stored reference atmospheres.
242
243  CALL RRTM_SETCOEF_140GP &
244   &( KIDIA  , KFDIA    , KLEV   , ZCOLDRY  , ZWBRODL , ZWKL , &
245   &  ZFAC00 , ZFAC01   , ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, JP, JT, JT1 , &
246   &  ZCOLH2O, ZCOLCO2  , ZCOLO3 , ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT , ZCOLBRD, &
247   &  ILAYTROP,ILAYSWTCH, ILAYLOW, ZPAVEL , ZTAVEL , ZSELFFAC, ZSELFFRAC, INDSELF, &
248   &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
249   &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
250   &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
251   &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)   
252
253  CALL RRTM_GASABS1A_140GP &
254   &( KIDIA   , KFDIA  , KLEV, ZATR1, ZOD, ZTF1, ZPAVEL, ZCOLDRY, ZCOLBRD, ZWX ,&
255   &  ZTAUAERL, ZFAC00 , ZFAC01, ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, JP, JT, JT1, ZONEMINUS ,&
256   &  ZCOLH2O , ZCOLCO2, ZCOLO3, ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT ,&
257   &  ILAYTROP, ILAYSWTCH,ILAYLOW, ZSELFFAC, ZSELFFRAC, INDSELF, ZPFRAC, &
258   &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
259   &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
260   &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
261   &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)     
262
263!- Call the radiative transfer routine.
264
265!  Clear and cloudy parts of column are treated together in RTRN.
266
267!  print 9901,JL,ZTBOUND
268
269  CALL RRTM_RTRN1A_140GP_MCICA &
270   &( KIDIA, KFDIA, KLEV, ISTART, IEND, KCOLS ,&
271   &  ZCLDFRAC, ZTAUCLD, ZATR1 ,&
272   &  ZOD , ZTF1 , &
273   &  ZTOTDFLUC, ZTOTDFLUX, ZTOTUFLUC, ZTOTUFLUX ,&
274   &  ZTAVEL, ZTZ, ZTBOUND, ZPFRAC, ZSEMISS, ZSEMISLW ,&
275   &  ZLwDerivative ) 
276
277! ***   Pass clear sky and total sky up and down flux profiles to ECRT
278!       output arrays (zflux, zfluc). Array indexing from bottom to top
279!       is preserved for ECRT.
280!       Invert down flux arrays for consistency with ECRT sign conventions.
281
282DO JL = KIDIA,KFDIA
283
284  PEMIT(JL) = ZSEMISLW(JL)
285  DO JK = 0, KLEV
286    PFLUC(JL,1,JK+1) =  ZTOTUFLUC(JL,JK)*ZFLUXFAC
287    PFLUC(JL,2,JK+1) = -ZTOTDFLUC(JL,JK)*ZFLUXFAC
288    PFLUX(JL,1,JK+1) =  ZTOTUFLUX(JL,JK)*ZFLUXFAC
289    PFLUX(JL,2,JK+1) = -ZTOTDFLUX(JL,JK)*ZFLUXFAC
290  ENDDO
291
292  ! Copy to output array, noting that they may be dimensioned
293  ! differently
294  PLwDerivative(JL,:) = ZLwDerivative(JL,:)
295
296ENDDO
297
2989901 FORMAT(1X,'rrtm:',I4,12E12.5)
299
300!------------------------------------------------------------------------
301END ASSOCIATE
302
303IF (LHOOK) CALL DR_HOOK('RRTM_RRTM_140GP_MCICA',1,ZHOOK_HANDLE)
304END SUBROUTINE RRTM_RRTM_140GP_MCICA
Note: See TracBrowser for help on using the repository browser.