source: LMDZ6/trunk/libf/phylmd/radlwsw_m.F90 @ 3908

Last change on this file since 3908 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, ...
  • 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: 61.7 KB
Line 
1!
2! $Id: radlwsw_m.F90 3908 2021-05-20 07:11:13Z idelkadi $
3!
4module radlwsw_m
5
6  IMPLICIT NONE
7
8contains
9
10SUBROUTINE radlwsw( &
11   dist, rmu0, fract, &
12!albedo SB >>>
13!  paprs, pplay,tsol,alb1, alb2, &
14   paprs, pplay,tsol,SFRWL,alb_dir, alb_dif, &
15!albedo SB <<<
16   t,q,wo,&
17   cldfra, cldemi, cldtaupd,&
18   ok_ade, ok_aie, ok_volcan, flag_aerosol,&
19   flag_aerosol_strat, flag_aer_feedback, &
20   tau_aero, piz_aero, cg_aero,&
21   tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! rajoute par OB RRTM
22   tau_aero_lw_rrtm, &              ! rajoute par C.Kleinschmitt pour RRTM
23   cldtaupi, &
24   qsat, flwc, fiwc, &
25   ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
26   heat,heat0,cool,cool0,albpla,&
27   heat_volc, cool_volc,&
28   topsw,toplw,solsw,solswfdiff,sollw,&
29   sollwdown,&
30   topsw0,toplw0,solsw0,sollw0,&
31   lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup,&
32   swdnc0, swdn0, swdn, swupc0, swup0, swup,&
33   topswad_aero, solswad_aero,&
34   topswai_aero, solswai_aero, &
35   topswad0_aero, solswad0_aero,&
36   topsw_aero, topsw0_aero,&
37   solsw_aero, solsw0_aero, &
38   topswcf_aero, solswcf_aero,&
39!-C. Kleinschmitt for LW diagnostics
40   toplwad_aero, sollwad_aero,&
41   toplwai_aero, sollwai_aero, &
42   toplwad0_aero, sollwad0_aero,&
43!-end
44   ZLWFT0_i, ZFLDN0, ZFLUP0,&
45   ZSWFT0_i, ZFSDN0, ZFSUP0)
46
47! Modules necessaires
48  USE DIMPHY
49  USE assert_m, ONLY : assert
50  USE infotrac_phy, ONLY : type_trac
51  USE write_field_phy
52
53#ifdef REPROBUS
54  USE CHEM_REP, ONLY : solaireTIME, ok_SUNTIME, ndimozon
55#endif
56
57#ifdef CPP_RRTM
58!    modules necessaires au rayonnement
59!    -----------------------------------------
60      USE YOERAD   , ONLY : NLW, LRRTM    ,LCCNL    ,LCCNO ,&
61          NRADIP   , NRADLP , NICEOPT, NLIQOPT ,RCCNLND  , RCCNSEA
62      USE YOELW    , ONLY : NSIL     ,NTRA     ,NUA      ,TSTAND   ,XP
63      USE YOESW    , ONLY : RYFWCA   ,RYFWCB   ,RYFWCC   ,RYFWCD,&   
64          RYFWCE   ,RYFWCF   ,REBCUA   ,REBCUB   ,REBCUC,&   
65          REBCUD   ,REBCUE   ,REBCUF   ,REBCUI   ,REBCUJ,& 
66          REBCUG   ,REBCUH   ,RHSAVI   ,RFULIO   ,RFLAA0,& 
67          RFLAA1   ,RFLBB0   ,RFLBB1   ,RFLBB2   ,RFLBB3,& 
68          RFLCC0   ,RFLCC1   ,RFLCC2   ,RFLCC3   ,RFLDD0,& 
69          RFLDD1   ,RFLDD2   ,RFLDD3   ,RFUETA   ,RASWCA,&
70          RASWCB   ,RASWCC   ,RASWCD   ,RASWCE   ,RASWCF
71      USE YOERDU   , ONLY : NUAER  ,NTRAER ,REPLOG ,REPSC  ,REPSCW ,DIFF
72      USE YOERRTWN , ONLY : DELWAVE   ,TOTPLNK     
73      USE YOMPHY3  , ONLY : RII0
74#endif
75      USE aero_mod
76
77! AI 02.2021
78! Besoin pour ECRAD de pctsrf, zmasq, longitude, altitude
79#ifdef CPP_ECRAD
80      USE geometry_mod, ONLY: latitude, longitude
81      USE phys_state_var_mod, ONLY: pctsrf
82      USE indice_sol_mod
83      USE time_phylmdz_mod, only: current_time
84      USE phys_cal_mod, only: day_cur
85#endif
86
87  !======================================================================
88  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
89  ! Objet: interface entre le modele et les rayonnements
90  ! Arguments:
91  !                  INPUTS
92  ! dist----- input-R- distance astronomique terre-soleil
93  ! rmu0----- input-R- cosinus de l'angle zenithal
94  ! fract---- input-R- duree d'ensoleillement normalisee
95  ! co2_ppm-- input-R- concentration du gaz carbonique (en ppm)
96  ! paprs---- input-R- pression a inter-couche (Pa)
97  ! pplay---- input-R- pression au milieu de couche (Pa)
98  ! tsol----- input-R- temperature du sol (en K)
99  ! alb1----- input-R- albedo du sol(entre 0 et 1) dans l'interval visible
100  ! alb2----- input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge   
101  ! t-------- input-R- temperature (K)
102  ! q-------- input-R- vapeur d'eau (en kg/kg)
103  ! cldfra--- input-R- fraction nuageuse (entre 0 et 1)
104  ! cldtaupd- input-R- epaisseur optique des nuages dans le visible (present-day value)
105  ! cldemi--- input-R- emissivite des nuages dans l'IR (entre 0 et 1)
106  ! ok_ade--- input-L- apply the Aerosol Direct Effect or not?
107  ! ok_aie--- input-L- apply the Aerosol Indirect Effect or not?
108  ! ok_volcan input-L- activate volcanic diags (SW heat & LW cool rate, SW & LW flux)
109  ! flag_aerosol input-I- aerosol flag from 0 to 6
110  ! flag_aerosol_strat input-I- use stratospheric aerosols flag (0, 1, 2)
111  ! flag_aer_feedback  input-I- activate aerosol radiative feedback (T, F)
112  ! tau_ae, piz_ae, cg_ae input-R- aerosol optical properties (calculated in aeropt.F)
113  ! cldtaupi  input-R- epaisseur optique des nuages dans le visible
114  !                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
115  !                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
116  !                   it is needed for the diagnostics of the aerosol indirect radiative forcing     
117  !
118  !                  OUTPUTS
119  ! heat-----output-R- echauffement atmospherique (visible) (K/jour)
120  ! cool-----output-R- refroidissement dans l'IR (K/jour)
121  ! albpla---output-R- albedo planetaire (entre 0 et 1)
122  ! topsw----output-R- flux solaire net au sommet de l'atm.
123  ! toplw----output-R- ray. IR montant au sommet de l'atmosphere
124  ! solsw----output-R- flux solaire net a la surface
125  ! solswfdiff----output-R- fraction de rayonnement diffus pour le flux solaire descendant a la surface
126  ! sollw----output-R- ray. IR montant a la surface
127  ! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
128  ! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
129  ! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
130  ! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
131  !
132  ! heat_volc-----output-R- echauffement atmospherique  du au forcage volcanique (visible) (K/s)
133  ! cool_volc-----output-R- refroidissement dans l'IR du au forcage volcanique (K/s)
134  !
135  ! ATTENTION: swai and swad have to be interpreted in the following manner:
136  ! ---------
137  ! ok_ade=F & ok_aie=F -both are zero
138  ! ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
139  !                        indirect is zero
140  ! ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
141  !                        direct is zero
142  ! ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
143  !                        aerosol direct forcing is F_{AD} = topswai-topswad
144  !
145  ! --------- RRTM: output RECMWFL
146  ! ZEMTD (KPROMA,KLEV+1)         ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
147  ! ZEMTU (KPROMA,KLEV+1)         ; TOTAL UPWARD   LONGWAVE EMISSIVITY
148  ! ZTRSO (KPROMA,KLEV+1)         ; TOTAL SHORTWAVE TRANSMISSIVITY
149  ! ZTH   (KPROMA,KLEV+1)         ; HALF LEVEL TEMPERATURE
150  ! ZCTRSO(KPROMA,2)              ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
151  ! ZCEMTR(KPROMA,2)              ; CLEAR-SKY NET LONGWAVE EMISSIVITY
152  ! ZTRSOD(KPROMA)                ; TOTAL-SKY SURFACE SW TRANSMISSITY
153  ! ZLWFC (KPROMA,2)              ; CLEAR-SKY LONGWAVE FLUXES
154  ! ZLWFT (KPROMA,KLEV+1)         ; TOTAL-SKY LONGWAVE FLUXES
155  ! ZLWFT0(KPROMA,KLEV+1)         ; CLEAR-SKY LONGWAVE FLUXES      ! added by MPL 090109
156  ! ZSWFC (KPROMA,2)              ; CLEAR-SKY SHORTWAVE FLUXES
157  ! ZSWFT (KPROMA,KLEV+1)         ; TOTAL-SKY SHORTWAVE FLUXES
158  ! ZSWFT0(KPROMA,KLEV+1)         ; CLEAR-SKY SHORTWAVE FLUXES     ! added by MPL 090109
159  ! ZFLUX (KLON,2,KLEV+1)         ; TOTAL LW FLUXES  1=up, 2=DWN   ! added by MPL 080411
160  ! ZFLUC (KLON,2,KLEV+1)         ; CLEAR SKY LW FLUXES            ! added by MPL 080411
161  ! ZFSDWN(klon,KLEV+1)           ; TOTAL SW  DWN FLUXES           ! added by MPL 080411
162  ! ZFCDWN(klon,KLEV+1)           ; CLEAR SKY SW  DWN FLUXES       ! added by MPL 080411
163  ! ZFCCDWN(klon,KLEV+1)          ; CLEAR SKY CLEAN (NO AEROSOL) SW  DWN FLUXES      ! added by OB 211117
164  ! ZFSUP (klon,KLEV+1)           ; TOTAL SW  UP  FLUXES           ! added by MPL 080411
165  ! ZFCUP (klon,KLEV+1)           ; CLEAR SKY SW  UP  FLUXES       ! added by MPL 080411
166  ! ZFCCUP (klon,KLEV+1)          ; CLEAR SKY CLEAN (NO AEROSOL) SW  UP  FLUXES      ! added by OB 211117
167  ! ZFLCCDWN(klon,KLEV+1)         ; CLEAR SKY CLEAN (NO AEROSOL) LW  DWN FLUXES      ! added by OB 211117
168  ! ZFLCCUP (klon,KLEV+1)         ; CLEAR SKY CLEAN (NO AEROSOL) LW  UP  FLUXES      ! added by OB 211117
169 
170  !======================================================================
171 
172  ! ====================================================================
173  ! Adapte au modele de chimie INCA par Celine Deandreis & Anne Cozic -- 2009
174  ! 1 = ZERO   
175  ! 2 = AER total   
176  ! 3 = NAT   
177  ! 4 = BC   
178  ! 5 = SO4   
179  ! 6 = POM   
180  ! 7 = DUST   
181  ! 8 = SS   
182  ! 9 = NO3   
183  !
184  ! ====================================================================
185
186! ==============
187! DECLARATIONS
188! ==============
189  include "YOETHF.h"
190  include "YOMCST.h"
191  include "clesphys.h"
192
193! Input arguments
194  REAL,    INTENT(in)  :: dist
195  REAL,    INTENT(in)  :: rmu0(KLON), fract(KLON)
196  REAL,    INTENT(in)  :: paprs(KLON,KLEV+1), pplay(KLON,KLEV)
197!albedo SB >>>
198! REAL,    INTENT(in)  :: alb1(KLON), alb2(KLON), tsol(KLON)
199  REAL,    INTENT(in)  :: tsol(KLON)
200  REAL,    INTENT(in) :: alb_dir(KLON,NSW),alb_dif(KLON,NSW)
201  REAL,    INTENT(in) :: SFRWL(6)
202!albedo SB <<<
203  REAL,    INTENT(in)  :: t(KLON,KLEV), q(KLON,KLEV)
204
205  REAL, INTENT(in):: wo(:, :, :) ! dimension(KLON,KLEV, 1 or 2)
206  ! column-density of ozone in a layer, in kilo-Dobsons
207  ! "wo(:, :, 1)" is for the average day-night field,
208  ! "wo(:, :, 2)" is for daylight time.
209
210  LOGICAL, INTENT(in)  :: ok_ade, ok_aie                                 ! switches whether to use aerosol direct (indirect) effects or not
211  LOGICAL, INTENT(in)  :: ok_volcan                                      ! produce volcanic diags (SW/LW heat flux and rate)
212  LOGICAL              :: lldebug
213  INTEGER, INTENT(in)  :: flag_aerosol                                   ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
214  INTEGER, INTENT(in)  :: flag_aerosol_strat                             ! use stratospheric aerosols
215  LOGICAL, INTENT(in)  :: flag_aer_feedback                              ! activate aerosol radiative feedback
216  REAL,    INTENT(in)  :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV)
217  REAL,    INTENT(in)  :: tau_aero(KLON,KLEV,naero_grp,2)                        ! aerosol optical properties (see aeropt.F)
218  REAL,    INTENT(in)  :: piz_aero(KLON,KLEV,naero_grp,2)                        ! aerosol optical properties (see aeropt.F)
219  REAL,    INTENT(in)  :: cg_aero(KLON,KLEV,naero_grp,2)                         ! aerosol optical properties (see aeropt.F)
220!--OB
221  REAL,    INTENT(in)  :: tau_aero_sw_rrtm(KLON,KLEV,2,NSW)                 ! aerosol optical properties RRTM
222  REAL,    INTENT(in)  :: piz_aero_sw_rrtm(KLON,KLEV,2,NSW)                 ! aerosol optical properties RRTM
223  REAL,    INTENT(in)  :: cg_aero_sw_rrtm(KLON,KLEV,2,NSW)                  ! aerosol optical properties RRTM
224!--OB fin
225
226!--C. Kleinschmitt
227#ifdef CPP_RRTM
228  REAL,    INTENT(in)  :: tau_aero_lw_rrtm(KLON,KLEV,2,NLW)                 ! LW aerosol optical properties RRTM
229#else
230  REAL,    INTENT(in)  :: tau_aero_lw_rrtm(KLON,KLEV,2,nbands_lw_rrtm)
231#endif
232!--C. Kleinschmitt end
233
234  REAL,    INTENT(in)  :: cldtaupi(KLON,KLEV)                            ! cloud optical thickness for pre-industrial aerosol concentrations
235  REAL,    INTENT(in)  :: qsat(klon,klev) ! Variable pour iflag_rrtm=1
236  REAL,    INTENT(in)  :: flwc(klon,klev) ! Variable pour iflag_rrtm=1
237  REAL,    INTENT(in)  :: fiwc(klon,klev) ! Variable pour iflag_rrtm=1
238  REAL,    INTENT(in)  :: ref_liq(klon,klev) ! cloud droplet radius present-day from newmicro
239  REAL,    INTENT(in)  :: ref_ice(klon,klev) ! ice crystal radius   present-day from newmicro
240  REAL,    INTENT(in)  :: ref_liq_pi(klon,klev) ! cloud droplet radius pre-industrial from newmicro
241  REAL,    INTENT(in)  :: ref_ice_pi(klon,klev) ! ice crystal radius   pre-industrial from newmicro
242
243! Output arguments
244  REAL,    INTENT(out) :: heat(KLON,KLEV), cool(KLON,KLEV)
245  REAL,    INTENT(out) :: heat0(KLON,KLEV), cool0(KLON,KLEV)
246  REAL,    INTENT(out) :: heat_volc(KLON,KLEV), cool_volc(KLON,KLEV) !NL
247  REAL,    INTENT(out) :: topsw(KLON), toplw(KLON)
248  REAL,    INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON), solswfdiff(KLON)
249  REAL,    INTENT(out) :: topsw0(KLON), toplw0(KLON), solsw0(KLON), sollw0(KLON)
250  REAL,    INTENT(out) :: sollwdown(KLON)
251  REAL,    INTENT(out) :: swdn(KLON,kflev+1),swdn0(KLON,kflev+1), swdnc0(KLON,kflev+1)
252  REAL,    INTENT(out) :: swup(KLON,kflev+1),swup0(KLON,kflev+1), swupc0(KLON,kflev+1)
253  REAL,    INTENT(out) :: lwdn(KLON,kflev+1),lwdn0(KLON,kflev+1), lwdnc0(KLON,kflev+1)
254  REAL,    INTENT(out) :: lwup(KLON,kflev+1),lwup0(KLON,kflev+1), lwupc0(KLON,kflev+1)
255  REAL,    INTENT(out) :: topswad_aero(KLON), solswad_aero(KLON)         ! output: aerosol direct forcing at TOA and surface
256  REAL,    INTENT(out) :: topswai_aero(KLON), solswai_aero(KLON)         ! output: aerosol indirect forcing atTOA and surface
257  REAL,    INTENT(out) :: toplwad_aero(KLON), sollwad_aero(KLON)         ! output: LW aerosol direct forcing at TOA and surface
258  REAL,    INTENT(out) :: toplwai_aero(KLON), sollwai_aero(KLON)         ! output: LW aerosol indirect forcing atTOA and surface
259  REAL, DIMENSION(klon), INTENT(out)    :: topswad0_aero
260  REAL, DIMENSION(klon), INTENT(out)    :: solswad0_aero
261  REAL, DIMENSION(klon), INTENT(out)    :: toplwad0_aero
262  REAL, DIMENSION(klon), INTENT(out)    :: sollwad0_aero
263  REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw_aero
264  REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw0_aero
265  REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw_aero
266  REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw0_aero
267  REAL, DIMENSION(kdlon,3), INTENT(out) :: topswcf_aero
268  REAL, DIMENSION(kdlon,3), INTENT(out) :: solswcf_aero
269  REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZSWFT0_i
270  REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZLWFT0_i
271
272! Local variables
273  REAL(KIND=8) ZFSUP(KDLON,KFLEV+1)
274  REAL(KIND=8) ZFSDN(KDLON,KFLEV+1)
275  REAL(KIND=8) ZFSUP0(KDLON,KFLEV+1)
276  REAL(KIND=8) ZFSDN0(KDLON,KFLEV+1)
277  REAL(KIND=8) ZFSUPC0(KDLON,KFLEV+1)
278  REAL(KIND=8) ZFSDNC0(KDLON,KFLEV+1)
279  REAL(KIND=8) ZFLUP(KDLON,KFLEV+1)
280  REAL(KIND=8) ZFLDN(KDLON,KFLEV+1)
281  REAL(KIND=8) ZFLUP0(KDLON,KFLEV+1)
282  REAL(KIND=8) ZFLDN0(KDLON,KFLEV+1)
283  REAL(KIND=8) ZFLUPC0(KDLON,KFLEV+1)
284  REAL(KIND=8) ZFLDNC0(KDLON,KFLEV+1)
285  REAL(KIND=8) zx_alpha1, zx_alpha2
286  INTEGER k, kk, i, j, iof, nb_gr
287  INTEGER ist,iend,ktdia,kmode
288  REAL(KIND=8) PSCT
289  REAL(KIND=8) PALBD(kdlon,2), PALBP(kdlon,2)
290!  MPL 06.01.09: pour RRTM, creation de PALBD_NEW et PALBP_NEW
291! avec NSW en deuxieme dimension       
292  REAL(KIND=8) PALBD_NEW(kdlon,NSW), PALBP_NEW(kdlon,NSW)
293  REAL(KIND=8) PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
294  REAL(KIND=8) PPSOL(kdlon), PDP(kdlon,KLEV)
295  REAL(KIND=8) PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
296  REAL(KIND=8) PTAVE(kdlon,kflev)
297  REAL(KIND=8) PWV(kdlon,kflev), PQS(kdlon,kflev)
298
299!!!!!!! Declarations specifiques pour ECRAD !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
300! AI 02.2021
301#ifdef CPP_ECRAD
302! ATTENTION les dimensions klon, kdlon ???
303! INPUTS
304  REAL(KIND=8) ZEMISW(klon), &              ! LW emissivity inside the window region
305               ZEMIS(klon)                  ! LW emissivity outside the window region
306  REAL(KIND=8) ZGELAM(klon), &              ! longitudes en rad
307               ZGEMU(klon)                  ! sin(latitude)
308  REAL(KIND=8) ZCO2(klon,klev), &           ! CO2 mass mixing ratios on full levels
309               ZCH4(klon,klev), &           ! CH4 mass mixing ratios on full levels
310               ZN2O(klon,klev), &           ! N2O mass mixing ratios on full levels
311               ZNO2(klon,klev), &           ! NO2 mass mixing ratios on full levels
312               ZCFC11(klon,klev), &         ! CFC11
313               ZCFC12(klon,klev), &         ! CFC12
314               ZHCFC22(klon,klev), &        ! HCFC22
315               ZCCL4(klon,klev), &          ! CCL4
316               ZO3_DP(klon,klev), ZO3_DP_i(klon,klev)            ! Ozone
317  REAL(KIND=8) ZQ_RAIN(klon,klev), &        ! Rain cloud mass mixing ratio (kg/kg) ?
318               ZQ_SNOW(klon,klev)           ! Snow cloud mass mixing ratio (kg/kg) ?
319  REAL(KIND=8) ZAEROSOL_OLD(KLON,6,KLEV), &  !
320               ZAEROSOL(KLON,KLEV,naero_tot) !
321! OUTPUTS
322  REAL(KIND=8) ZFLUX_DIR(klon), &           ! Direct compt of surf flux into horizontal plane
323               ZFLUX_DIR_CLEAR(klon), &     ! CS Direct
324               ZFLUX_DIR_INTO_SUN(klon), &  !
325               ZFLUX_UV(klon), &            ! UV flux
326               ZFLUX_PAR(klon), &           ! photosynthetically active radiation similarly
327               ZFLUX_PAR_CLEAR(klon), &     ! CS photosynthetically
328               ZFLUX_SW_DN_TOA(klon), &     ! DN SW flux at TOA
329               ZEMIS_OUT(klon)              ! effective broadband emissivity
330  REAL(KIND=8) ZLWDERIVATIVE(klon,klev+1)   ! LW derivatives
331  REAL(KIND=8) ZSWDIFFUSEBAND(klon,NSW), &  ! SW DN flux in diffuse albedo band
332               ZSWDIRECTBAND(klon,NSW)      ! SW DN flux in direct albedo band
333#endif
334!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
335
336  REAL(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
337  ! "POZON(:, :, 1)" is for the average day-night field,
338  ! "POZON(:, :, 2)" is for daylight time.
339!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6 
340  REAL(KIND=8) PAER(kdlon,kflev,6)
341  REAL(KIND=8) PCLDLD(kdlon,kflev)
342  REAL(KIND=8) PCLDLU(kdlon,kflev)
343  REAL(KIND=8) PCLDSW(kdlon,kflev)
344  REAL(KIND=8) PTAU(kdlon,2,kflev)
345  REAL(KIND=8) POMEGA(kdlon,2,kflev)
346  REAL(KIND=8) PCG(kdlon,2,kflev)
347  REAL(KIND=8) zfract(kdlon), zrmu0(kdlon), zdist
348  REAL(KIND=8) zheat(kdlon,kflev), zcool(kdlon,kflev)
349  REAL(KIND=8) zheat0(kdlon,kflev), zcool0(kdlon,kflev)
350  REAL(KIND=8) zheat_volc(kdlon,kflev), zcool_volc(kdlon,kflev) !NL
351  REAL(KIND=8) ztopsw(kdlon), ztoplw(kdlon)
352  REAL(KIND=8) zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon), zsolswfdiff(kdlon)
353  REAL(KIND=8) zsollwdown(kdlon)
354  REAL(KIND=8) ztopsw0(kdlon), ztoplw0(kdlon)
355  REAL(KIND=8) zsolsw0(kdlon), zsollw0(kdlon)
356  REAL(KIND=8) zznormcp
357  REAL(KIND=8) tauaero(kdlon,kflev,naero_grp,2)                     ! aer opt properties
358  REAL(KIND=8) pizaero(kdlon,kflev,naero_grp,2)
359  REAL(KIND=8) cgaero(kdlon,kflev,naero_grp,2)
360  REAL(KIND=8) PTAUA(kdlon,2,kflev)                         ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
361  REAL(KIND=8) POMEGAA(kdlon,2,kflev)                       ! dito for single scatt albedo
362  REAL(KIND=8) ztopswadaero(kdlon), zsolswadaero(kdlon)     ! Aerosol direct forcing at TOAand surface
363  REAL(KIND=8) ztopswad0aero(kdlon), zsolswad0aero(kdlon)   ! Aerosol direct forcing at TOAand surface
364  REAL(KIND=8) ztopswaiaero(kdlon), zsolswaiaero(kdlon)     ! dito, indirect
365!--NL
366  REAL(KIND=8) zswadaero(kdlon,kflev+1)                       ! SW Aerosol direct forcing
367  REAL(KIND=8) zlwadaero(kdlon,kflev+1)                       ! LW Aerosol direct forcing
368!-LW by CK
369  REAL(KIND=8) ztoplwadaero(kdlon), zsollwadaero(kdlon)     ! LW Aerosol direct forcing at TOAand surface
370  REAL(KIND=8) ztoplwad0aero(kdlon), zsollwad0aero(kdlon)   ! LW Aerosol direct forcing at TOAand surface
371  REAL(KIND=8) ztoplwaiaero(kdlon), zsollwaiaero(kdlon)     ! dito, indirect
372!-end
373  REAL(KIND=8) ztopsw_aero(kdlon,9), ztopsw0_aero(kdlon,9)
374  REAL(KIND=8) zsolsw_aero(kdlon,9), zsolsw0_aero(kdlon,9)
375  REAL(KIND=8) ztopswcf_aero(kdlon,3), zsolswcf_aero(kdlon,3)     
376! real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2 deje declare dans physiq.F MPL 20130618
377!MPL input supplementaires pour RECMWFL
378! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
379  REAL(KIND=8) GEMU(klon)
380!MPL input RECMWFL:
381! Tableaux aux niveaux inverses pour respecter convention Arpege
382  REAL(KIND=8) ref_liq_i(klon,klev) ! cloud droplet radius present-day from newmicro (inverted)
383  REAL(KIND=8) ref_ice_i(klon,klev) ! ice crystal radius present-day from newmicro (inverted)
384!--OB
385  REAL(KIND=8) ref_liq_pi_i(klon,klev) ! cloud droplet radius pre-industrial from newmicro (inverted)
386  REAL(KIND=8) ref_ice_pi_i(klon,klev) ! ice crystal radius pre-industrial from newmicro (inverted)
387!--end OB
388  REAL(KIND=8) paprs_i(klon,klev+1)
389  REAL(KIND=8) pplay_i(klon,klev)
390  REAL(KIND=8) cldfra_i(klon,klev)
391  REAL(KIND=8) POZON_i(kdlon,kflev, size(wo, 3)) ! mass fraction of ozone
392  ! "POZON(:, :, 1)" is for the average day-night field,
393  ! "POZON(:, :, 2)" is for daylight time.
394!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6     
395  REAL(KIND=8) PAER_i(kdlon,kflev,6)
396  REAL(KIND=8) PDP_i(klon,klev)
397  REAL(KIND=8) t_i(klon,klev),q_i(klon,klev),qsat_i(klon,klev)
398  REAL(KIND=8) flwc_i(klon,klev),fiwc_i(klon,klev)
399!MPL output RECMWFL:
400  REAL(KIND=8) ZEMTD (klon,klev+1),ZEMTD_i (klon,klev+1)       
401  REAL(KIND=8) ZEMTU (klon,klev+1),ZEMTU_i (klon,klev+1)     
402  REAL(KIND=8) ZTRSO (klon,klev+1),ZTRSO_i (klon,klev+1)   
403  REAL(KIND=8) ZTH   (klon,klev+1),ZTH_i   (klon,klev+1)   
404  REAL(KIND=8) ZCTRSO(klon,2)       
405  REAL(KIND=8) ZCEMTR(klon,2)     
406  REAL(KIND=8) ZTRSOD(klon)       
407  REAL(KIND=8) ZLWFC (klon,2)     
408  REAL(KIND=8) ZLWFT (klon,klev+1),ZLWFT_i (klon,klev+1)   
409  REAL(KIND=8) ZSWFC (klon,2)     
410  REAL(KIND=8) ZSWFT (klon,klev+1),ZSWFT_i (klon,klev+1)
411  REAL(KIND=8) ZFLUCDWN_i(klon,klev+1),ZFLUCUP_i(klon,klev+1)
412  REAL(KIND=8) PPIZA_TOT(klon,klev,NSW)
413  REAL(KIND=8) PCGA_TOT(klon,klev,NSW)
414  REAL(KIND=8) PTAU_TOT(klon,klev,NSW)
415  REAL(KIND=8) PPIZA_NAT(klon,klev,NSW)
416  REAL(KIND=8) PCGA_NAT(klon,klev,NSW)
417  REAL(KIND=8) PTAU_NAT(klon,klev,NSW)
418#ifdef CPP_RRTM
419  REAL(KIND=8) PTAU_LW_TOT(klon,klev,NLW)
420  REAL(KIND=8) PTAU_LW_NAT(klon,klev,NLW)
421#endif
422  REAL(KIND=8) PSFSWDIR(klon,NSW)
423  REAL(KIND=8) PSFSWDIF(klon,NSW)
424  REAL(KIND=8) PFSDNN(klon)
425  REAL(KIND=8) PFSDNV(klon)
426!MPL On ne redefinit pas les tableaux ZFLUX,ZFLUC,
427!MPL ZFSDWN,ZFCDWN,ZFSUP,ZFCUP car ils existent deja
428!MPL sous les noms de ZFLDN,ZFLDN0,ZFLUP,ZFLUP0,
429!MPL ZFSDN,ZFSDN0,ZFSUP,ZFSUP0
430  REAL(KIND=8) ZFLUX_i (klon,2,klev+1)
431  REAL(KIND=8) ZFLUC_i (klon,2,klev+1)
432  REAL(KIND=8) ZFSDWN_i (klon,klev+1)
433  REAL(KIND=8) ZFCDWN_i (klon,klev+1)
434  REAL(KIND=8) ZFCCDWN_i (klon,klev+1)
435  REAL(KIND=8) ZFSUP_i (klon,klev+1)
436  REAL(KIND=8) ZFCUP_i (klon,klev+1)
437  REAL(KIND=8) ZFCCUP_i (klon,klev+1)
438  REAL(KIND=8) ZFLCCDWN_i (klon,klev+1)
439  REAL(KIND=8) ZFLCCUP_i (klon,klev+1)
440! 3 lignes suivantes a activer pour CCMVAL (MPL 20100412)
441!      REAL(KIND=8) RSUN(3,2)
442!      REAL(KIND=8) SUN(3)
443!      REAL(KIND=8) SUN_FRACT(2)
444  REAL, PARAMETER:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
445  CHARACTER (LEN=80) :: abort_message
446  CHARACTER (LEN=80) :: modname='radlwsw_m'
447
448  REAL zdir, zdif
449
450! =========  INITIALISATIONS ==============================================
451 IF (lldebug) THEN
452  print*,'Entree dans radlwsw '
453  print*,'************* INITIALISATIONS *****************************'
454  print*,'klon, kdlon, klev, kflev =',klon, kdlon, klev, kflev
455 ENDIF
456
457  CALL assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
458 
459  ist=1
460  iend=klon
461  ktdia=1
462  kmode=ist
463! Aeros
464  tauaero(:,:,:,:)=0.
465  pizaero(:,:,:,:)=0.
466  cgaero(:,:,:,:)=0.
467  lldebug=.FALSE.
468
469  ztopsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
470  ztopsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
471  zsolsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
472  zsolsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
473
474  ZTOPSWADAERO(:)  = 0. !ym missing init
475  ZSOLSWADAERO(:)  = 0. !ym missing init
476  ZTOPSWAD0AERO(:) = 0. !ym missing init
477  ZSOLSWAD0AERO(:) = 0. !ym missing init
478  ZTOPSWAIAERO(:)  = 0. !ym missing init
479  ZSOLSWAIAERO(:)  = 0. !ym missing init 
480  ZTOPSWCF_AERO(:,:)= 0.!ym missing init 
481  ZSOLSWCF_AERO(:,:) =0. !ym missing init 
482
483  !
484  !-------------------------------------------
485  nb_gr = KLON / kdlon
486  IF (nb_gr*kdlon .NE. KLON) THEN
487      PRINT*, "kdlon mauvais:", KLON, kdlon, nb_gr
488      call abort_physic("radlwsw", "", 1)
489  ENDIF
490  IF (kflev .NE. KLEV) THEN
491      PRINT*, "kflev differe de KLEV, kflev, KLEV"
492      call abort_physic("radlwsw", "", 1)
493  ENDIF
494  !-------------------------------------------
495  DO k = 1, KLEV
496    DO i = 1, KLON
497      heat(i,k)=0.
498      cool(i,k)=0.
499      heat_volc(i,k)=0. !NL
500      cool_volc(i,k)=0. !NL
501      heat0(i,k)=0.
502      cool0(i,k)=0.
503    ENDDO
504  ENDDO
505  !
506  zdist = dist
507  !
508  PSCT = solaire/zdist/zdist
509
510  IF (type_trac == 'repr') THEN
511#ifdef REPROBUS
512    IF (iflag_rrtm==0) THEN
513      IF (ok_SUNTIME) PSCT = solaireTIME/zdist/zdist
514      print*,'Constante solaire: ',PSCT*zdist*zdist
515    ENDIF
516#endif
517  ENDIF
518
519 IF (lldebug) THEN
520  print*,'************** Debut boucle de 1 a ', nb_gr
521 ENDIF
522
523  DO j = 1, nb_gr
524    iof = kdlon*(j-1)
525    DO i = 1, kdlon
526      zfract(i) = fract(iof+i)
527      zrmu0(i) = rmu0(iof+i)
528
529
530      IF (iflag_rrtm==0) THEN
531!     Albedo
532        PALBD(i,1)=alb_dif(iof+i,1)
533        PALBD(i,2)=alb_dif(iof+i,2)
534        PALBP(i,1)=alb_dir(iof+i,1)
535        PALBP(i,2)=alb_dir(iof+i,2)
536! AI 02.2021 cas iflag_rrtm=1 et 2
537       ELSEIF (iflag_rrtm==1.OR.iflag_rrtm==2) THEN
538        DO kk=1,NSW
539          PALBD_NEW(i,kk)=alb_dif(iof+i,kk)
540          PALBP_NEW(i,kk)=alb_dir(iof+i,kk)
541        ENDDO
542!
543      ENDIF
544!albedo SB <<<
545
546      PEMIS(i) = 1.0    !!!!! A REVOIR (MPL)
547      PVIEW(i) = 1.66
548      PPSOL(i) = paprs(iof+i,1)
549      zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2))
550      zx_alpha2 = 1.0 - zx_alpha1
551      PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
552      PTL(i,KLEV+1) = t(iof+i,KLEV)
553      PDT0(i) = tsol(iof+i) - PTL(i,1)
554    ENDDO
555    DO k = 2, kflev
556      DO i = 1, kdlon
557        PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
558      ENDDO
559    ENDDO
560    DO k = 1, kflev
561      DO i = 1, kdlon
562        PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
563        PTAVE(i,k) = t(iof+i,k)
564        PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
565        PQS(i,k) = PWV(i,k)
566!       Confert from  column density of ozone in a cell, in kDU, to a mass fraction
567        POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
568             / (paprs(iof+i, k) - paprs(iof+i, k+1))
569        ZO3_DP(i,k) = wo(iof+i, k, 1) * RG * dobson_u * 1e3
570!       A activer pour CCMVAL on prend l'ozone impose (MPL 07042010)
571!       POZON(i,k,:) = wo(i,k,:) 
572!       print *,'RADLWSW: POZON',k, POZON(i,k,1)
573        PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
574        PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
575        PCLDSW(i,k) = cldfra(iof+i,k)
576        PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
577        PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
578        POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
579        POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
580        PCG(i,1,k) = 0.865
581        PCG(i,2,k) = 0.910
582        !-
583        ! Introduced for aerosol indirect forcings.
584        ! The following values use the cloud optical thickness calculated from
585        ! present-day aerosol concentrations whereas the quantities without the
586        ! "A" at the end are for pre-industial (natural-only) aerosol concentrations
587        !
588        PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
589        PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
590        POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
591        POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
592      ENDDO
593    ENDDO
594
595    IF (type_trac == 'repr') THEN
596#ifdef REPROBUS
597       ndimozon = size(wo, 3)
598       CALL RAD_INTERACTIF(POZON,iof)
599#endif
600    ENDIF
601    !
602    DO k = 1, kflev+1
603      DO i = 1, kdlon
604        PPMB(i,k) = paprs(iof+i,k)/100.0
605      ENDDO
606    ENDDO
607    !
608!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
609    DO kk = 1, 6
610      DO k = 1, kflev
611        DO i = 1, kdlon
612          PAER(i,k,kk) = 1.0E-15   !!!!! A REVOIR (MPL)
613        ENDDO
614      ENDDO
615    ENDDO
616    DO k = 1, kflev
617      DO i = 1, kdlon
618        tauaero(i,k,:,1)=tau_aero(iof+i,k,:,1)
619        pizaero(i,k,:,1)=piz_aero(iof+i,k,:,1)
620        cgaero(i,k,:,1) =cg_aero(iof+i,k,:,1)
621        tauaero(i,k,:,2)=tau_aero(iof+i,k,:,2)
622        pizaero(i,k,:,2)=piz_aero(iof+i,k,:,2)
623        cgaero(i,k,:,2) =cg_aero(iof+i,k,:,2)
624      ENDDO
625    ENDDO
626!
627! AI 02.2021
628#ifdef CPP_ECRAD
629  ZEMIS = 1.0
630  ZEMISW = 1.0
631  ZGELAM = longitude
632  ZGEMU = sin(latitude)
633  ZCO2 = RCO2
634  ZCH4 = RCH4
635  ZN2O = RN2O
636  ZNO2 = 0.0
637  ZCFC11 = RCFC11
638  ZCFC12 = RCFC12
639  ZHCFC22 = 0.0
640  ZCCL4 = 0.0
641  ZQ_RAIN = 0.0
642  ZQ_SNOW = 0.0
643  ZAEROSOL_OLD = 0.0
644  ZAEROSOL = 0.0
645#endif
646!
647!===== iflag_rrtm ================================================
648!     
649    IF (iflag_rrtm == 0) THEN       !!!! remettre 0 juste pour tester l'ancien rayt via rrtm
650!
651!--- Mise a zero des tableaux output du rayonnement LW-AR4 ----------             
652      DO k = 1, kflev+1
653      DO i = 1, kdlon
654!     print *,'RADLWSW: boucle mise a zero i k',i,k
655      ZFLUP(i,k)=0.
656      ZFLDN(i,k)=0.
657      ZFLUP0(i,k)=0.
658      ZFLDN0(i,k)=0.
659      ZLWFT0_i(i,k)=0.
660      ZFLUCUP_i(i,k)=0.
661      ZFLUCDWN_i(i,k)=0.
662      ENDDO
663      ENDDO
664      DO k = 1, kflev
665         DO i = 1, kdlon
666            zcool(i,k)=0.
667            zcool_volc(i,k)=0. !NL
668            zcool0(i,k)=0.
669         ENDDO
670      ENDDO
671      DO i = 1, kdlon
672      ztoplw(i)=0.
673      zsollw(i)=0.
674      ztoplw0(i)=0.
675      zsollw0(i)=0.
676      zsollwdown(i)=0.
677      ENDDO
678       ! Old radiation scheme, used for AR4 runs
679       ! average day-night ozone for longwave
680       CALL LW_LMDAR4(&
681            PPMB, PDP,&
682            PPSOL,PDT0,PEMIS,&
683            PTL, PTAVE, PWV, POZON(:, :, 1), PAER,&
684            PCLDLD,PCLDLU,&
685            PVIEW,&
686            zcool, zcool0,&
687            ztoplw,zsollw,ztoplw0,zsollw0,&
688            zsollwdown,&
689            ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
690!----- Mise a zero des tableaux output du rayonnement SW-AR4
691      DO k = 1, kflev+1
692         DO i = 1, kdlon
693            ZFSUP(i,k)=0.
694            ZFSDN(i,k)=0.
695            ZFSUP0(i,k)=0.
696            ZFSDN0(i,k)=0.
697            ZFSUPC0(i,k)=0.
698            ZFSDNC0(i,k)=0.
699            ZFLUPC0(i,k)=0.
700            ZFLDNC0(i,k)=0.
701            ZSWFT0_i(i,k)=0.
702            ZFCUP_i(i,k)=0.
703            ZFCDWN_i(i,k)=0.
704            ZFCCUP_i(i,k)=0.
705            ZFCCDWN_i(i,k)=0.
706            ZFLCCUP_i(i,k)=0.
707            ZFLCCDWN_i(i,k)=0.
708            zswadaero(i,k)=0. !--NL
709         ENDDO
710      ENDDO
711      DO k = 1, kflev
712         DO i = 1, kdlon
713            zheat(i,k)=0.
714            zheat_volc(i,k)=0.
715            zheat0(i,k)=0.
716         ENDDO
717      ENDDO
718      DO i = 1, kdlon
719      zalbpla(i)=0.
720      ztopsw(i)=0.
721      zsolsw(i)=0.
722      ztopsw0(i)=0.
723      zsolsw0(i)=0.
724      ztopswadaero(i)=0.
725      zsolswadaero(i)=0.
726      ztopswaiaero(i)=0.
727      zsolswaiaero(i)=0.
728      ENDDO
729
730      !--fraction of diffuse radiation in surface SW downward radiation
731      !--not computed with old radiation scheme
732      zsolswfdiff(:) = -999.999
733
734!     print *,'Avant SW_LMDAR4: PSCT zrmu0 zfract',PSCT, zrmu0, zfract
735       ! daylight ozone, if we have it, for short wave
736      CALL SW_AEROAR4(PSCT, zrmu0, zfract,&
737               PPMB, PDP,&
738               PPSOL, PALBD, PALBP,&
739               PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
740               PCLDSW, PTAU, POMEGA, PCG,&
741               zheat, zheat0,&
742               zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
743               ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
744               tauaero, pizaero, cgaero, &
745               PTAUA, POMEGAA,&
746               ztopswadaero,zsolswadaero,&
747               ztopswad0aero,zsolswad0aero,&
748               ztopswaiaero,zsolswaiaero, &
749               ztopsw_aero,ztopsw0_aero,&
750               zsolsw_aero,zsolsw0_aero,&
751               ztopswcf_aero,zsolswcf_aero, &
752               ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat)
753
754       ZSWFT0_i(:,:) = ZFSDN0(:,:)-ZFSUP0(:,:)
755       ZLWFT0_i(:,:) =-ZFLDN0(:,:)-ZFLUP0(:,:)
756
757       DO i=1,kdlon
758       DO k=1,kflev+1
759         lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
760         lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
761         lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
762         lwup  ( iof+i,k)   = ZFLUP  ( i,k)
763         swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
764         swdn  ( iof+i,k)   = ZFSDN  ( i,k)
765         swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
766         swup  ( iof+i,k)   = ZFSUP  ( i,k)
767       ENDDO 
768       ENDDO 
769!
770    ELSE IF (iflag_rrtm == 1) then 
771#ifdef CPP_RRTM
772!      if (prt_level.gt.10)write(lunout,*)'CPP_RRTM=.T.'
773!===== iflag_rrtm=1, on passe dans SW via RECMWFL ===============
774
775      DO k = 1, kflev+1
776      DO i = 1, kdlon
777        ZEMTD_i(i,k)=0.
778        ZEMTU_i(i,k)=0.
779        ZTRSO_i(i,k)=0.
780        ZTH_i(i,k)=0.
781        ZLWFT_i(i,k)=0.
782        ZSWFT_i(i,k)=0.
783        ZFLUX_i(i,1,k)=0.
784        ZFLUX_i(i,2,k)=0.
785        ZFLUC_i(i,1,k)=0.
786        ZFLUC_i(i,2,k)=0.
787        ZFSDWN_i(i,k)=0.
788        ZFCDWN_i(i,k)=0.
789        ZFCCDWN_i(i,k)=0.
790        ZFSUP_i(i,k)=0.
791        ZFCUP_i(i,k)=0.
792        ZFCCUP_i(i,k)=0.
793        ZFLCCDWN_i(i,k)=0.
794        ZFLCCUP_i(i,k)=0.
795      ENDDO
796      ENDDO
797!
798!--OB
799!--aerosol TOT  - anthropogenic+natural - index 2
800!--aerosol NAT  - natural only          - index 1
801!
802      DO i = 1, kdlon
803      DO k = 1, kflev
804      DO kk=1, NSW
805!
806      PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
807      PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
808      PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
809!
810      PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
811      PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
812      PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
813!
814      ENDDO
815      ENDDO
816      ENDDO
817!-end OB
818!
819!--C. Kleinschmitt
820!--aerosol TOT  - anthropogenic+natural - index 2
821!--aerosol NAT  - natural only          - index 1
822!
823      DO i = 1, kdlon
824      DO k = 1, kflev
825      DO kk=1, NLW
826!
827      PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
828      PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
829!
830      ENDDO
831      ENDDO
832      ENDDO
833!-end C. Kleinschmitt
834!     
835      DO i = 1, kdlon
836      ZCTRSO(i,1)=0.
837      ZCTRSO(i,2)=0.
838      ZCEMTR(i,1)=0.
839      ZCEMTR(i,2)=0.
840      ZTRSOD(i)=0.
841      ZLWFC(i,1)=0.
842      ZLWFC(i,2)=0.
843      ZSWFC(i,1)=0.
844      ZSWFC(i,2)=0.
845      PFSDNN(i)=0.
846      PFSDNV(i)=0.
847      DO kk = 1, NSW
848        PSFSWDIR(i,kk)=0.
849        PSFSWDIF(i,kk)=0.
850      ENDDO
851      ENDDO
852!----- Fin des mises a zero des tableaux output de RECMWF -------------------             
853!        GEMU(1:klon)=sin(rlatd(1:klon))
854! On met les donnees dans l'ordre des niveaux arpege
855         paprs_i(:,1)=paprs(:,klev+1)
856         DO k=1,klev
857            paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
858            pplay_i(1:klon,k)   =pplay(1:klon,klev+1-k)
859            cldfra_i(1:klon,k)  =cldfra(1:klon,klev+1-k)
860            PDP_i(1:klon,k)     =PDP(1:klon,klev+1-k)
861            t_i(1:klon,k)       =t(1:klon,klev+1-k)
862            q_i(1:klon,k)       =q(1:klon,klev+1-k)
863            qsat_i(1:klon,k)    =qsat(1:klon,klev+1-k)
864            flwc_i(1:klon,k)    =flwc(1:klon,klev+1-k)
865            fiwc_i(1:klon,k)    =fiwc(1:klon,klev+1-k)
866            ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)
867            ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)
868!-OB
869            ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
870            ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
871         ENDDO
872         DO k=1,kflev
873           POZON_i(1:klon,k,:)=POZON(1:klon,kflev+1-k,:)
874!!!            POZON_i(1:klon,k)=POZON(1:klon,k)            !!! on laisse 1=sol et klev=top
875!          print *,'Juste avant RECMWFL: k tsol temp',k,tsol,t(1,k)
876!!!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6     
877            DO i=1,6
878            PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
879            ENDDO
880         ENDDO
881
882!       print *,'RADLWSW: avant RECMWFL, RI0,rmu0=',solaire,rmu0
883
884!  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
885! La version ARPEGE1D utilise differentes valeurs de la constante
886! solaire suivant le rayonnement utilise.
887! A controler ...
888! SOLAR FLUX AT THE TOP (/YOMPHY3/)
889! introduce season correction
890!--------------------------------------
891! RII0 = RIP0
892! IF(LRAYFM)
893! RII0 = RIP0M   ! =rip0m if Morcrette non-each time step call.
894! IF(LRAYFM15)
895! RII0 = RIP0M15 ! =rip0m if Morcrette non-each time step call.
896         RII0=solaire/zdist/zdist
897!print*,'+++ radlwsw: solaire ,RII0',solaire,RII0
898!  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
899! Ancien appel a RECMWF (celui du cy25)
900!        CALL RECMWF (ist , iend, klon , ktdia , klev   , kmode ,
901!    s   PALBD    , PALBP   , paprs_i , pplay_i , RCO2   , cldfra_i,
902!    s   POZON_i  , PAER_i  , PDP_i   , PEMIS   , GEMU   , rmu0,
903!    s    q_i     , qsat_i  , fiwc_i  , flwc_i  , zmasq  , t_i  ,tsol,
904!    s   ZEMTD_i  , ZEMTU_i , ZTRSO_i ,
905!    s   ZTH_i    , ZCTRSO  , ZCEMTR  , ZTRSOD  ,
906!    s   ZLWFC    , ZLWFT_i , ZSWFC   , ZSWFT_i ,
907!    s   ZFLUX_i  , ZFLUC_i , ZFSDWN_i, ZFSUP_i , ZFCDWN_i,ZFCUP_i)
908!    s   'RECMWF ')
909!
910      IF (lldebug) THEN
911        CALL writefield_phy('paprs_i',paprs_i,klev+1)
912        CALL writefield_phy('pplay_i',pplay_i,klev)
913        CALL writefield_phy('cldfra_i',cldfra_i,klev)
914        CALL writefield_phy('pozon_i',POZON_i,klev)
915        CALL writefield_phy('paer_i',PAER_i,klev)
916        CALL writefield_phy('pdp_i',PDP_i,klev)
917        CALL writefield_phy('q_i',q_i,klev)
918        CALL writefield_phy('qsat_i',qsat_i,klev)
919        CALL writefield_phy('fiwc_i',fiwc_i,klev)
920        CALL writefield_phy('flwc_i',flwc_i,klev)
921        CALL writefield_phy('t_i',t_i,klev)
922        CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
923        CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
924      ENDIF
925
926! Nouvel appel a RECMWF (celui du cy32t0)
927         CALL RECMWF_AERO (ist , iend, klon , ktdia  , klev   , kmode ,&
928!                          KST,  KEND, KPROMA, KTDIA , KLEV,  KMODE
929         PALBD_NEW,PALBP_NEW, paprs_i , pplay_i , RCO2   , cldfra_i,&
930!        PALBD ,   PALBP ,    PAPRS ,   PAPRSF ,  PCCO2 , PCLFR
931         POZON_i  , PAER_i  , PDP_i   , PEMIS   , rmu0   ,&
932!        PQO3  ,    PAER  ,   PDP   ,   PEMIS  , PMU0
933          q_i     , qsat_i  , fiwc_i  , flwc_i  , zmasq  , t_i  ,tsol,&
934!         PQ       , PQS   , PQIWP ,     PQLWP , PSLM   , PT    , PTS,
935         ref_liq_i, ref_ice_i, &
936!        PREF_LIQ, PREF_ICE
937         ref_liq_pi_i, ref_ice_pi_i, &   ! rajoute par OB pour diagnostiquer effet indirect
938!        PREF_LIQ_PI, PREF_ICE_PI
939         ZEMTD_i  , ZEMTU_i , ZTRSO_i ,&
940!        PEMTD , PEMTU , PTRSO
941         ZTH_i    , ZCTRSO  , ZCEMTR  , ZTRSOD  ,&
942!        PTH   ,     PCTRSO,   PCEMTR, PTRSOD
943         ZLWFC    , ZLWFT_i , ZSWFC   , ZSWFT_i ,&
944!        PLWFC,     PLWFT,    PSWFC,    PSWFT,
945         PSFSWDIR , PSFSWDIF, PFSDNN  , PFSDNV  ,&
946!        PSFSWDIR,  PSFSWDIF, PFSDNN,    PFSDNV
947         PPIZA_TOT, PCGA_TOT,PTAU_TOT,&
948!        PPIZA_TOT, PCGA_TOT,PTAU_TOT
949         PPIZA_NAT, PCGA_NAT,PTAU_NAT,           &  ! rajoute par OB pour diagnostiquer effet direct
950!        PPIZA_NAT,PCGA_NAT,PTAU_NAT
951         PTAU_LW_TOT, PTAU_LW_NAT,               &  ! rajoute par C. Kleinschmitt
952!        PTAU_LW_TOT, PTAU_LW_NAT,
953         ZFLUX_i  , ZFLUC_i ,&
954!        PFLUX,      PFLUC
955         ZFSDWN_i , ZFSUP_i , ZFCDWN_i, ZFCUP_i, ZFCCDWN_i, ZFCCUP_i, ZFLCCDWN_i, ZFLCCUP_i, &
956!        PFSDN ,    PFSUP ,   PFSCDN , PFSCUP,   PFSCCDN,   PFSCCUP,  PFLCCDN,    PFLCCUP
957         ZTOPSWADAERO,ZSOLSWADAERO,&  ! rajoute par OB pour diagnostics
958!        PTOPSWADAERO,PSOLSWADAERO
959         ZTOPSWAD0AERO,ZSOLSWAD0AERO,&
960         ZTOPSWAIAERO,ZSOLSWAIAERO, &
961         ZTOPSWCF_AERO,ZSOLSWCF_AERO, &
962         ZSWADAERO, & !--NL
963         ZTOPLWADAERO,ZSOLLWADAERO,&  ! rajoute par C. Kleinscmitt pour LW diagnostics
964         ZTOPLWAD0AERO,ZSOLLWAD0AERO,&
965         ZTOPLWAIAERO,ZSOLLWAIAERO, &
966         ZLWADAERO, & !--NL
967         ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback) ! flags aerosols
968
969!--OB diagnostics
970! & PTOPSWAIAERO,PSOLSWAIAERO,&
971! & PTOPSWCFAERO,PSOLSWCFAERO,&
972! & PSWADAERO,& !--NL
973!!--LW diagnostics CK
974! & PTOPLWADAERO,PSOLLWADAERO,&
975! & PTOPLWAD0AERO,PSOLLWAD0AERO,&
976! & PTOPLWAIAERO,PSOLLWAIAERO,&
977! & PLWADAERO,& !--NL
978!!..end
979! & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat,&
980! & flag_aer_feedback)
981
982           
983!        print *,'RADLWSW: apres RECMWF'
984      IF (lldebug) THEN
985        CALL writefield_phy('zemtd_i',ZEMTD_i,klev+1)
986        CALL writefield_phy('zemtu_i',ZEMTU_i,klev+1)
987        CALL writefield_phy('ztrso_i',ZTRSO_i,klev+1)
988        CALL writefield_phy('zth_i',ZTH_i,klev+1)
989        CALL writefield_phy('zctrso',ZCTRSO,2)
990        CALL writefield_phy('zcemtr',ZCEMTR,2)
991        CALL writefield_phy('ztrsod',ZTRSOD,1)
992        CALL writefield_phy('zlwfc',ZLWFC,2)
993        CALL writefield_phy('zlwft_i',ZLWFT_i,klev+1)
994        CALL writefield_phy('zswfc',ZSWFC,2)
995        CALL writefield_phy('zswft_i',ZSWFT_i,klev+1)
996        CALL writefield_phy('psfswdir',PSFSWDIR,6)
997        CALL writefield_phy('psfswdif',PSFSWDIF,6)
998        CALL writefield_phy('pfsdnn',PFSDNN,1)
999        CALL writefield_phy('pfsdnv',PFSDNV,1)
1000        CALL writefield_phy('ppiza_dst',PPIZA_TOT,klev)
1001        CALL writefield_phy('pcga_dst',PCGA_TOT,klev)
1002        CALL writefield_phy('ptaurel_dst',PTAU_TOT,klev)
1003        CALL writefield_phy('zflux_i',ZFLUX_i,klev+1)
1004        CALL writefield_phy('zfluc_i',ZFLUC_i,klev+1)
1005        CALL writefield_phy('zfsdwn_i',ZFSDWN_i,klev+1)
1006        CALL writefield_phy('zfsup_i',ZFSUP_i,klev+1)
1007        CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
1008        CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
1009      ENDIF
1010! --------- output RECMWFL
1011!  ZEMTD        (KPROMA,KLEV+1)  ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
1012!  ZEMTU        (KPROMA,KLEV+1)  ; TOTAL UPWARD   LONGWAVE EMISSIVITY
1013!  ZTRSO        (KPROMA,KLEV+1)  ; TOTAL SHORTWAVE TRANSMISSIVITY
1014!  ZTH          (KPROMA,KLEV+1)  ; HALF LEVEL TEMPERATURE
1015!  ZCTRSO       (KPROMA,2)       ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
1016!  ZCEMTR       (KPROMA,2)       ; CLEAR-SKY NET LONGWAVE EMISSIVITY
1017!  ZTRSOD       (KPROMA)         ; TOTAL-SKY SURFACE SW TRANSMISSITY
1018!  ZLWFC        (KPROMA,2)       ; CLEAR-SKY LONGWAVE FLUXES
1019!  ZLWFT        (KPROMA,KLEV+1)  ; TOTAL-SKY LONGWAVE FLUXES
1020!  ZSWFC        (KPROMA,2)       ; CLEAR-SKY SHORTWAVE FLUXES
1021!  ZSWFT        (KPROMA,KLEV+1)  ; TOTAL-SKY SHORTWAVE FLUXES
1022!  PPIZA_TOT    (KPROMA,KLEV,NSW); Single scattering albedo of total aerosols
1023!  PCGA_TOT     (KPROMA,KLEV,NSW); Assymetry factor for total aerosols
1024!  PTAU_TOT     (KPROMA,KLEV,NSW); Optical depth of total aerosols
1025!  PPIZA_NAT    (KPROMA,KLEV,NSW); Single scattering albedo of natural aerosols
1026!  PCGA_NAT     (KPROMA,KLEV,NSW); Assymetry factor for natural aerosols
1027!  PTAU_NAT     (KPROMA,KLEV,NSW); Optical depth of natiral aerosols
1028!  PTAU_LW_TOT  (KPROMA,KLEV,NLW); LW Optical depth of total aerosols 
1029!  PTAU_LW_NAT  (KPROMA,KLEV,NLW); LW Optical depth of natural aerosols 
1030!  PSFSWDIR     (KPROMA,NSW)     ;
1031!  PSFSWDIF     (KPROMA,NSW)     ;
1032!  PFSDNN       (KPROMA)         ;
1033!  PFSDNV       (KPROMA)         ;
1034! ---------
1035! ---------
1036! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
1037! D autre part, on multiplie les resultats SW par fract pour etre coherent
1038! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
1039! rayonnement SW. (MPL 260609)
1040      DO k=0,klev
1041         DO i=1,klon
1042         ZEMTD(i,k+1)  = ZEMTD_i(i,k+1)
1043         ZEMTU(i,k+1)  = ZEMTU_i(i,k+1)
1044         ZTRSO(i,k+1)  = ZTRSO_i(i,k+1)
1045         ZTH(i,k+1)    = ZTH_i(i,k+1)
1046!        ZLWFT(i,k+1)  = ZLWFT_i(i,klev+1-k)
1047!        ZSWFT(i,k+1)  = ZSWFT_i(i,klev+1-k)
1048         ZFLUP(i,k+1)  = ZFLUX_i(i,1,k+1)
1049         ZFLDN(i,k+1)  = ZFLUX_i(i,2,k+1)
1050         ZFLUP0(i,k+1) = ZFLUC_i(i,1,k+1)
1051         ZFLDN0(i,k+1) = ZFLUC_i(i,2,k+1)
1052         ZFSDN(i,k+1)  = ZFSDWN_i(i,k+1)*fract(i)
1053         ZFSDN0(i,k+1) = ZFCDWN_i(i,k+1)*fract(i)
1054         ZFSDNC0(i,k+1)= ZFCCDWN_i(i,k+1)*fract(i)
1055         ZFSUP (i,k+1) = ZFSUP_i(i,k+1)*fract(i)
1056         ZFSUP0(i,k+1) = ZFCUP_i(i,k+1)*fract(i)
1057         ZFSUPC0(i,k+1)= ZFCCUP_i(i,k+1)*fract(i)
1058         ZFLDNC0(i,k+1)= ZFLCCDWN_i(i,k+1)
1059         ZFLUPC0(i,k+1)= ZFLCCUP_i(i,k+1)
1060         IF (ok_volcan) THEN
1061            ZSWADAERO(i,k+1)=ZSWADAERO(i,k+1)*fract(i) !--NL
1062         ENDIF
1063         
1064!   Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
1065!   en sortie de radlsw.F90 - MPL 7.01.09
1066         ZSWFT(i,k+1)  = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
1067         ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
1068!        WRITE(*,'("FSDN FSUP FCDN FCUP: ",4E12.5)') ZFSDWN_i(i,k+1),&
1069!        ZFSUP_i(i,k+1),ZFCDWN_i(i,k+1),ZFCUP_i(i,k+1)
1070         ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
1071         ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
1072!        print *,'FLUX2 FLUX1 FLUC2 FLUC1',ZFLUX_i(i,2,k+1),&
1073!    & ZFLUX_i(i,1,k+1),ZFLUC_i(i,2,k+1),ZFLUC_i(i,1,k+1)
1074         ENDDO
1075      ENDDO
1076
1077!--ajout OB
1078      ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
1079      ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
1080      ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
1081      ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
1082      ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
1083      ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
1084      ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
1085      ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
1086      ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
1087      ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
1088      ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
1089      ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
1090
1091!     print*,'SW_RRTM ZFSDN0 1 , klev:',ZFSDN0(1:klon,1),ZFSDN0(1:klon,klev)
1092!     print*,'SW_RRTM ZFSUP0 1 , klev:',ZFSUP0(1:klon,1),ZFSUP0(1:klon,klev)
1093!     print*,'SW_RRTM ZFSDN  1 , klev:',ZFSDN(1:klon,1),ZFSDN(1:klon,klev)
1094!     print*,'SW_RRTM ZFSUP  1 , klev:',ZFSUP(1:klon,1),ZFSUP(1:klon,klev)     
1095!     print*,'OK1'
1096! ---------
1097! ---------
1098! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
1099! LW_LMDAR4 et SW_LMDAR4
1100
1101      !--fraction of diffuse radiation in surface SW downward radiation
1102      DO i = 1, kdlon
1103       IF (fract(i).GT.0.0) THEN
1104         zdir=SUM(PSFSWDIR(i,:))
1105         zdif=SUM(PSFSWDIF(i,:))
1106         zsolswfdiff(i) = zdif/(zdir+zdif)
1107       ELSE  !--night
1108         zsolswfdiff(i) = 1.0
1109       ENDIF
1110      ENDDO
1111!
1112      DO i = 1, kdlon
1113         zsolsw(i)    = ZSWFT(i,1)
1114         zsolsw0(i)   = ZSWFT0_i(i,1)
1115!        zsolsw0(i)   = ZFSDN0(i,1)     -ZFSUP0(i,1)
1116         ztopsw(i)    = ZSWFT(i,klev+1)
1117         ztopsw0(i)   = ZSWFT0_i(i,klev+1)
1118!        ztopsw0(i)   = ZFSDN0(i,klev+1)-ZFSUP0(i,klev+1)
1119!         
1120!        zsollw(i)    = ZFLDN(i,1)      -ZFLUP(i,1)
1121!        zsollw0(i)   = ZFLDN0(i,1)     -ZFLUP0(i,1)
1122!        ztoplw(i)    = ZFLDN(i,klev+1) -ZFLUP(i,klev+1)
1123!        ztoplw0(i)   = ZFLDN0(i,klev+1)-ZFLUP0(i,klev+1)
1124         zsollw(i)    = ZLWFT(i,1)
1125         zsollw0(i)   = ZLWFT0_i(i,1)
1126         ztoplw(i)    = ZLWFT(i,klev+1)*(-1)
1127         ztoplw0(i)   = ZLWFT0_i(i,klev+1)*(-1)
1128!         
1129         IF (fract(i) == 0.) THEN
1130!!!!! A REVOIR MPL (20090630) ca n a pas de sens quand fract=0
1131! pas plus que dans le sw_AR4
1132          zalbpla(i)   = 1.0e+39
1133         ELSE
1134          zalbpla(i)   = ZFSUP(i,klev+1)/ZFSDN(i,klev+1)
1135         ENDIF
1136!!! 5 juin 2015
1137!!! Correction MP bug RRTM
1138         zsollwdown(i)= -1.*ZFLDN(i,1)
1139      ENDDO
1140!     print*,'OK2'
1141
1142! extrait de SW_AR4
1143!     DO k = 1, KFLEV
1144!        kpl1 = k+1
1145!        DO i = 1, KDLON
1146!           PHEAT(i,k) = -(ZFSUP(i,kpl1)-ZFSUP(i,k)) -(ZFSDN(i,k)-ZFSDN(i,kpl1))
1147!           PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k)
1148! ZLWFT(klon,k),ZSWFT
1149
1150      DO k=1,kflev
1151         DO i=1,kdlon
1152           zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))*RDAY*RG/RCPD/PDP(i,k)
1153           zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))*RDAY*RG/RCPD/PDP(i,k)
1154           zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
1155           zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
1156           IF (ok_volcan) THEN
1157              zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
1158              zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
1159           ENDIF
1160!          print *,'heat cool heat0 cool0 ',zheat(i,k),zcool(i,k),zheat0(i,k),zcool0(i,k)
1161!          ZFLUCUP_i(i,k)=ZFLUC_i(i,1,k)
1162!          ZFLUCDWN_i(i,k)=ZFLUC_i(i,2,k)         
1163         ENDDO
1164      ENDDO
1165#else
1166    abort_message="You should compile with -rrtm if running with iflag_rrtm=1"
1167    call abort_physic(modname, abort_message, 1)
1168#endif
1169!======================================================================
1170! AI fev 2021
1171    ELSE IF(iflag_rrtm == 2) THEN
1172    print*,'Traitement cas iflag_rrtm = ',iflag_rrtm
1173!    print*,'Mise a zero des flux '
1174#ifdef CPP_ECRAD
1175      DO k = 1, kflev+1
1176      DO i = 1, kdlon
1177        ZEMTD_i(i,k)=0.
1178        ZEMTU_i(i,k)=0.
1179        ZTRSO_i(i,k)=0.
1180        ZTH_i(i,k)=0.
1181        ZLWFT_i(i,k)=0.
1182        ZSWFT_i(i,k)=0.
1183        ZFLUX_i(i,1,k)=0.
1184        ZFLUX_i(i,2,k)=0.
1185        ZFLUC_i(i,1,k)=0.
1186        ZFLUC_i(i,2,k)=0.
1187        ZFSDWN_i(i,k)=0.
1188        ZFCDWN_i(i,k)=0.
1189        ZFCCDWN_i(i,k)=0.
1190        ZFSUP_i(i,k)=0.
1191        ZFCUP_i(i,k)=0.
1192        ZFCCUP_i(i,k)=0.
1193        ZFLCCDWN_i(i,k)=0.
1194        ZFLCCUP_i(i,k)=0.
1195      ENDDO
1196      ENDDO
1197!
1198! Aerosols A REVOIR
1199!      DO i = 1, kdlon
1200!      DO k = 1, kflev
1201!      DO kk=1, NSW
1202!
1203!      PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
1204!      PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
1205!      PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
1206!
1207!      PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
1208!      PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
1209!      PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
1210!
1211!      ENDDO
1212!      ENDDO
1213!      ENDDO
1214!-end OB
1215!
1216!      DO i = 1, kdlon
1217!      DO k = 1, kflev
1218!      DO kk=1, NLW
1219!
1220!      PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
1221!      PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
1222!
1223!      ENDDO
1224!      ENDDO
1225!      ENDDO
1226!-end C. Kleinschmitt
1227!     
1228      DO i = 1, kdlon
1229      ZCTRSO(i,1)=0.
1230      ZCTRSO(i,2)=0.
1231      ZCEMTR(i,1)=0.
1232      ZCEMTR(i,2)=0.
1233      ZTRSOD(i)=0.
1234      ZLWFC(i,1)=0.
1235      ZLWFC(i,2)=0.
1236      ZSWFC(i,1)=0.
1237      ZSWFC(i,2)=0.
1238      PFSDNN(i)=0.
1239      PFSDNV(i)=0.
1240      DO kk = 1, NSW
1241        PSFSWDIR(i,kk)=0.
1242        PSFSWDIF(i,kk)=0.
1243      ENDDO
1244      ENDDO
1245!----- Fin des mises a zero des tableaux output -------------------             
1246
1247! On met les donnees dans l'ordre des niveaux ecrad
1248!         print*,'On inverse sur la verticale '
1249         paprs_i(:,1)=paprs(:,klev+1)
1250         DO k=1,klev
1251            paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
1252            pplay_i(1:klon,k)   =pplay(1:klon,klev+1-k)
1253            cldfra_i(1:klon,k)  =cldfra(1:klon,klev+1-k)
1254            PDP_i(1:klon,k)     =PDP(1:klon,klev+1-k)
1255            t_i(1:klon,k)       =t(1:klon,klev+1-k)
1256            q_i(1:klon,k)       =q(1:klon,klev+1-k)
1257            qsat_i(1:klon,k)    =qsat(1:klon,klev+1-k)
1258            flwc_i(1:klon,k)    =flwc(1:klon,klev+1-k)
1259            fiwc_i(1:klon,k)    =fiwc(1:klon,klev+1-k)
1260            ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)
1261            ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)
1262!-OB
1263            ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
1264            ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
1265         ENDDO
1266         DO k=1,kflev
1267            ZO3_DP_i(1:klon,k)=ZO3_DP(1:klon,kflev+1-k)
1268            DO i=1,6
1269            PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
1270            ENDDO
1271         ENDDO
1272! AI 02.2021
1273! Calcul of ZTH_i
1274         DO K=2,KLEV
1275            ZTH_i(:,K)=&
1276              & (t_i(:,K-1)*pplay_i(:,K-1)*(pplay_i(:,K)-paprs_i(:,K))&
1277              & +t_i(:,K)*pplay_i(:,K)*(paprs_i(:,K)-pplay_i(:,K-1)))&
1278              & *(1.0/(paprs_i(:,K)*(pplay_i(:,K)-pplay_i(:,K-1))))
1279         ENDDO
1280            ZTH_i(:,KLEV+1)=tsol(:)
1281            ZTH_i(:,1)=t_i(:,1)-pplay_i(:,1)*(t_i(:,1)-ZTH_i(:,2))&
1282                      & /(pplay_i(:,1)-paprs_i(:,2))
1283
1284      print *,'RADLWSW: avant RADIATION_SCHEME '
1285      IF (lldebug) THEN
1286        CALL writefield_phy('rmu0',rmu0,1)
1287        CALL writefield_phy('tsol',tsol,1)
1288        CALL writefield_phy('emissiv_out',ZEMIS,1)
1289        CALL writefield_phy('emissiv_in',ZEMISW,1)
1290        CALL writefield_phy('pctsrf_ter',pctsrf(:,is_ter),1)
1291        CALL writefield_phy('pctsrf_oce',pctsrf(:,is_oce),1)
1292        CALL writefield_phy('ZGELAM',ZGELAM,1)
1293        CALL writefield_phy('ZGEMU',ZGEMU,1)
1294        CALL writefield_phy('zmasq',zmasq,1)
1295        CALL writefield_phy('paprs_i',paprs_i,klev+1)
1296        CALL writefield_phy('pplay_i',pplay_i,klev)
1297        CALL writefield_phy('t_i',t_i,klev)
1298        CALL writefield_phy('ZTH_i',ZTH_i,klev+1)
1299        CALL writefield_phy('cldfra_i',cldfra_i,klev)
1300        CALL writefield_phy('paer_i',PAER_i,klev)
1301        CALL writefield_phy('q_i',q_i,klev)
1302        CALL writefield_phy('fiwc_i',fiwc_i,klev)
1303        CALL writefield_phy('flwc_i',flwc_i,klev)
1304        CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
1305        CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
1306        CALL writefield_phy('ZO3_DP',ZO3_DP,klev)
1307      ENDIF
1308
1309      CALL RADIATION_SCHEME &
1310      & (ist, iend, klon, klev, naero_tot, NSW, &
1311! ??? naero_tot
1312      & day_cur, current_time, &
1313      & solaire, &
1314      & rmu0, tsol, PALBD_NEW,PALBP_NEW, &   
1315!       PEMIS_WINDOW (???), &
1316      &  ZEMIS, ZEMISW, &
1317!       PCCN_LAND, PCCN_SEA, & ???
1318      & pctsrf(:,is_ter), pctsrf(:,is_oce), &
1319!       longitude(rad), sin(latitude), PMASQ_ ???
1320      & ZGELAM, ZGEMU, zmasq, &
1321      & pplay_i, t_i, &
1322!       PTEMPERATURE_H ?, &
1323      & paprs_i, ZTH_i, q_i, qsat_i, &
1324      & ZCO2, ZCH4, ZN2O, ZNO2, ZCFC11, ZCFC12, ZHCFC22, ZCCL4, ZO3_DP_i, &
1325      & cldfra_i, flwc_i, fiwc_i, ZQ_RAIN, ZQ_SNOW, & 
1326      & ref_liq_i, ref_ice_i, &
1327!       aerosols
1328      & ZAEROSOL_OLD, ZAEROSOL, &
1329! Outputs
1330!       Net flux
1331      & ZSWFT_i, ZLWFT_i, ZSWFT0_i, ZLWFT0_i, &
1332!       DN flux
1333      & ZFSDWN_i, ZFLUX_i(:,2,:), ZFCDWN_i, ZFLUC_i(:,2,:), &
1334!       UP flux
1335      & ZFSUP_i, ZFLUX_i(:,1,:), ZFCUP_i, ZFLUC_i(:,1,:), &
1336!       Surf Direct flux
1337      & ZFLUX_DIR, ZFLUX_DIR_CLEAR, ZFLUX_DIR_INTO_SUN, &
1338!       UV and para flux
1339      & ZFLUX_UV, ZFLUX_PAR, ZFLUX_PAR_CLEAR, &
1340!      & ZFLUX_SW_DN_TOA,
1341      & ZEMIS_OUT, ZLWDERIVATIVE, &
1342      & PSFSWDIF, PSFSWDIR)
1343
1344      print *,'========= RADLWSW: apres RADIATION_SCHEME ==================== '
1345  IF (lldebug) THEN
1346    if (0.eq.0) then
1347      print *,' Net Flux '
1348      print *,'ZSWFT_i =', ZSWFT_i
1349      print *,'ZLWFT_i =', ZLWFT_i
1350      print *,'ZSWFT0_i =', ZSWFT0_i
1351      print*,'ZLWFT0_i =', ZLWFT0_i
1352
1353      print*,'DN Flux '
1354      print*,'ZFSDWN_i =', ZFSDWN_i
1355      print*,'ZFLUX_i(:,2,:)', ZFLUX_i(:,2,:)
1356      print*,'ZFCDWN_i =', ZFCDWN_i
1357      print*,'ZFLUC_i(:,2,:) =', ZFLUC_i(:,2,:)
1358
1359      print*,'UP Flux '
1360      print*,'ZFSUP_i =', ZFSUP_i
1361      print*,'ZFLUX_i(:,1,:) =', ZFLUX_i(:,1,:)
1362      print*,'ZFCUP_i =', ZFCUP_i
1363      print*,'ZFLUC_i(:,1,:) =', ZFLUC_i(:,1,:)
1364
1365      print*,'UV and para flux '
1366      print*,'ZFLUX_DIR =', ZFLUX_DIR
1367      print*,'ZFLUX_DIR_CLEAR', ZFLUX_DIR_CLEAR
1368      print*,'ZFLUX_DIR_INTO_SUN', ZFLUX_DIR_INTO_SUN
1369    endif
1370   ENDIF
1371
1372      IF (lldebug) THEN
1373        CALL writefield_phy('zlwft_i',ZLWFT_i,klev+1)
1374        CALL writefield_phy('zlwft0_i',ZLWFT0_i,klev+1)
1375        CALL writefield_phy('zswft_i',ZSWFT_i,klev+1)
1376        CALL writefield_phy('zswft0_i',ZSWFT0_i,klev+1)
1377        CALL writefield_phy('psfswdir',PSFSWDIR,6)
1378        CALL writefield_phy('psfswdif',PSFSWDIF,6)
1379        CALL writefield_phy('zflux_i',ZFLUX_i,klev+1)
1380        CALL writefield_phy('zfluc_i',ZFLUC_i,klev+1)
1381        CALL writefield_phy('zfsdwn_i',ZFSDWN_i,klev+1)
1382        CALL writefield_phy('zfsup_i',ZFSUP_i,klev+1)
1383        CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
1384        CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
1385      ENDIF
1386! ---------
1387! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
1388! D autre part, on multiplie les resultats SW par fract pour etre coherent
1389! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
1390! rayonnement SW. (MPL 260609)
1391      print*,'On retablit ordre des niveaux verticaux'
1392      print*,'On multiplie les flux SW par fract ?'
1393      DO k=0,klev
1394         DO i=1,klon
1395         ZEMTD(i,k+1)  = ZEMTD_i(i,k+1)
1396         ZEMTU(i,k+1)  = ZEMTU_i(i,k+1)
1397         ZTRSO(i,k+1)  = ZTRSO_i(i,k+1)
1398         ZTH(i,k+1)    = ZTH_i(i,k+1)
1399         ZFLUP(i,k+1)  = ZFLUX_i(i,1,k+1)
1400         ZFLDN(i,k+1)  = ZFLUX_i(i,2,k+1)
1401         ZFLUP0(i,k+1) = ZFLUC_i(i,1,k+1)
1402         ZFLDN0(i,k+1) = ZFLUC_i(i,2,k+1)
1403         ZFSDN(i,k+1)  = ZFSDWN_i(i,k+1)*fract(i)
1404         ZFSDN0(i,k+1) = ZFCDWN_i(i,k+1)*fract(i)
1405         ZFSDNC0(i,k+1)= ZFCCDWN_i(i,k+1)*fract(i)
1406         ZFSUP (i,k+1) = ZFSUP_i(i,k+1)*fract(i)
1407         ZFSUP0(i,k+1) = ZFCUP_i(i,k+1)*fract(i)
1408         ZFSUPC0(i,k+1)= ZFCCUP_i(i,k+1)*fract(i)
1409         ZFLDNC0(i,k+1)= ZFLCCDWN_i(i,k+1)
1410         ZFLUPC0(i,k+1)= ZFLCCUP_i(i,k+1)
1411         IF (ok_volcan) THEN
1412            ZSWADAERO(i,k+1)=ZSWADAERO(i,k+1)*fract(i) !--NL
1413         ENDIF
1414         
1415!   Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
1416!   en sortie de radlsw.F90 - MPL 7.01.09
1417         ZSWFT(i,k+1)  = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
1418         ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
1419         ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
1420         ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
1421         ENDDO
1422      ENDDO
1423
1424!--ajout OB
1425      ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
1426      ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
1427      ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
1428      ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
1429      ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
1430      ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
1431      ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
1432      ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
1433      ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
1434      ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
1435      ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
1436      ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
1437
1438! ---------
1439! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
1440! LW_LMDAR4 et SW_LMDAR4
1441
1442      !--fraction of diffuse radiation in surface SW downward radiation
1443      DO i = 1, kdlon
1444       IF (fract(i).GT.0.0) THEN
1445         zdir=SUM(PSFSWDIR(i,:))
1446         zdif=SUM(PSFSWDIF(i,:))
1447         zsolswfdiff(i) = zdif/(zdir+zdif)
1448       ELSE  !--night
1449         zsolswfdiff(i) = 1.0
1450       ENDIF
1451      ENDDO
1452!
1453      DO i = 1, kdlon
1454         zsolsw(i)    = ZSWFT(i,1)
1455         zsolsw0(i)   = ZSWFT0_i(i,1)
1456         ztopsw(i)    = ZSWFT(i,klev+1)
1457         ztopsw0(i)   = ZSWFT0_i(i,klev+1)
1458         zsollw(i)    = ZLWFT(i,1)
1459         zsollw0(i)   = ZLWFT0_i(i,1)
1460         ztoplw(i)    = ZLWFT(i,klev+1)*(-1)
1461         ztoplw0(i)   = ZLWFT0_i(i,klev+1)*(-1)
1462!         
1463         zsollwdown(i)= -1.*ZFLDN(i,1)
1464      ENDDO
1465
1466      DO k=1,kflev
1467         DO i=1,kdlon
1468           zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))*RDAY*RG/RCPD/PDP(i,k)
1469           zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))*RDAY*RG/RCPD/PDP(i,k)
1470           zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
1471           zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
1472           IF (ok_volcan) THEN
1473              zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
1474              zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
1475           ENDIF
1476         ENDDO
1477      ENDDO
1478#endif 
1479  print*,'Fin traitement ECRAD'
1480! Fin ECRAD
1481  ENDIF        ! iflag_rrtm
1482! ecrad
1483!======================================================================
1484
1485    DO i = 1, kdlon
1486      topsw(iof+i) = ztopsw(i)
1487      toplw(iof+i) = ztoplw(i)
1488      solsw(iof+i) = zsolsw(i)
1489      solswfdiff(iof+i) = zsolswfdiff(i)
1490      sollw(iof+i) = zsollw(i)
1491      sollwdown(iof+i) = zsollwdown(i)
1492      DO k = 1, kflev+1
1493        lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
1494        lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
1495        lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
1496        lwup  ( iof+i,k)   = ZFLUP  ( i,k)
1497      ENDDO
1498      topsw0(iof+i) = ztopsw0(i)
1499      toplw0(iof+i) = ztoplw0(i)
1500      solsw0(iof+i) = zsolsw0(i)
1501      sollw0(iof+i) = zsollw0(i)
1502      albpla(iof+i) = zalbpla(i)
1503
1504      DO k = 1, kflev+1
1505        swdnc0( iof+i,k)   = ZFSDNC0( i,k)
1506        swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
1507        swdn  ( iof+i,k)   = ZFSDN  ( i,k)
1508        swupc0( iof+i,k)   = ZFSUPC0( i,k)
1509        swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
1510        swup  ( iof+i,k)   = ZFSUP  ( i,k)
1511        lwdnc0( iof+i,k)   = ZFLDNC0( i,k)
1512        lwupc0( iof+i,k)   = ZFLUPC0( i,k)
1513      ENDDO
1514    ENDDO
1515    !-transform the aerosol forcings, if they have
1516    ! to be calculated
1517    IF (ok_ade) THEN
1518        DO i = 1, kdlon
1519          topswad_aero(iof+i) = ztopswadaero(i)
1520          topswad0_aero(iof+i) = ztopswad0aero(i)
1521          solswad_aero(iof+i) = zsolswadaero(i)
1522          solswad0_aero(iof+i) = zsolswad0aero(i)
1523          topsw_aero(iof+i,:) = ztopsw_aero(i,:)
1524          topsw0_aero(iof+i,:) = ztopsw0_aero(i,:)
1525          solsw_aero(iof+i,:) = zsolsw_aero(i,:)
1526          solsw0_aero(iof+i,:) = zsolsw0_aero(i,:)
1527          topswcf_aero(iof+i,:) = ztopswcf_aero(i,:)
1528          solswcf_aero(iof+i,:) = zsolswcf_aero(i,:)   
1529          !-LW
1530          toplwad_aero(iof+i) = ztoplwadaero(i)
1531          toplwad0_aero(iof+i) = ztoplwad0aero(i)
1532          sollwad_aero(iof+i) = zsollwadaero(i)
1533          sollwad0_aero(iof+i) = zsollwad0aero(i)   
1534        ENDDO
1535    ELSE
1536        DO i = 1, kdlon
1537          topswad_aero(iof+i) = 0.0
1538          solswad_aero(iof+i) = 0.0
1539          topswad0_aero(iof+i) = 0.0
1540          solswad0_aero(iof+i) = 0.0
1541          topsw_aero(iof+i,:) = 0.
1542          topsw0_aero(iof+i,:) =0.
1543          solsw_aero(iof+i,:) = 0.
1544          solsw0_aero(iof+i,:) = 0.
1545          !-LW
1546          toplwad_aero(iof+i) = 0.0
1547          sollwad_aero(iof+i) = 0.0
1548          toplwad0_aero(iof+i) = 0.0
1549          sollwad0_aero(iof+i) = 0.0
1550        ENDDO
1551    ENDIF
1552    IF (ok_aie) THEN
1553        DO i = 1, kdlon
1554          topswai_aero(iof+i) = ztopswaiaero(i)
1555          solswai_aero(iof+i) = zsolswaiaero(i)
1556          !-LW
1557          toplwai_aero(iof+i) = ztoplwaiaero(i)
1558          sollwai_aero(iof+i) = zsollwaiaero(i)
1559        ENDDO
1560    ELSE
1561        DO i = 1, kdlon
1562          topswai_aero(iof+i) = 0.0
1563          solswai_aero(iof+i) = 0.0
1564          !-LW
1565          toplwai_aero(iof+i) = 0.0
1566          sollwai_aero(iof+i) = 0.0
1567        ENDDO
1568    ENDIF
1569    DO k = 1, kflev
1570      DO i = 1, kdlon
1571        !        scale factor to take into account the difference between
1572        !        dry air and watter vapour scpecifi! heat capacity
1573        zznormcp=1.0+RVTMP2*PWV(i,k)
1574        heat(iof+i,k) = zheat(i,k)/zznormcp
1575        cool(iof+i,k) = zcool(i,k)/zznormcp
1576        heat0(iof+i,k) = zheat0(i,k)/zznormcp
1577        cool0(iof+i,k) = zcool0(i,k)/zznormcp
1578        IF(ok_volcan) THEN !NL
1579           heat_volc(iof+i,k) = zheat_volc(i,k)/zznormcp
1580           cool_volc(iof+i,k) = zcool_volc(i,k)/zznormcp
1581        ENDIF
1582      ENDDO
1583    ENDDO
1584
1585 ENDDO ! j = 1, nb_gr
1586
1587END SUBROUTINE radlwsw
1588
1589end module radlwsw_m
Note: See TracBrowser for help on using the repository browser.