source: LMDZ6/branches/Optimisation_LMDZ/libf/phylmd/radlwsw_m.F90 @ 3758

Last change on this file since 3758 was 3630, checked in by Laurent Fairhead, 5 years ago

Parameter new_aod is not needed anymore as it is assumed to be true
all the time. This means that we cannot replay AR4 simulations with new
LMDZ sources (we probably couldn't anyway)
LF, OB

  • 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: 47.1 KB
Line 
1!
2! $Id: radlwsw_m.F90 3630 2020-02-10 10:04:40Z adurocher $
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 pour 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,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
48
49  USE DIMPHY
50  USE assert_m, ONLY : assert
51  USE infotrac_phy, ONLY : type_trac
52  USE write_field_phy
53#ifdef REPROBUS
54  USE CHEM_REP, ONLY : solaireTIME, ok_SUNTIME, ndimozon
55#endif
56#ifdef CPP_RRTM
57!    modules necessaires au rayonnement
58!    -----------------------------------------
59!     USE YOMCST   , ONLY : RG       ,RD       ,RTT      ,RPI
60!     USE YOERAD   , ONLY : NSW      ,LRRTM    ,LINHOM   , LCCNL,LCCNO,
61!     USE YOERAD   , ONLY : NSW      ,LRRTM    ,LCCNL    ,LCCNO ,&
62! NSW mis dans .def MPL 20140211
63! NLW ajoute par OB
64      USE YOERAD   , ONLY : NLW, LRRTM    ,LCCNL    ,LCCNO ,&
65          NRADIP   , NRADLP , NICEOPT, NLIQOPT ,RCCNLND  , RCCNSEA
66      USE YOELW    , ONLY : NSIL     ,NTRA     ,NUA      ,TSTAND   ,XP
67      USE YOESW    , ONLY : RYFWCA   ,RYFWCB   ,RYFWCC   ,RYFWCD,&   
68          RYFWCE   ,RYFWCF   ,REBCUA   ,REBCUB   ,REBCUC,&   
69          REBCUD   ,REBCUE   ,REBCUF   ,REBCUI   ,REBCUJ,& 
70          REBCUG   ,REBCUH   ,RHSAVI   ,RFULIO   ,RFLAA0,& 
71          RFLAA1   ,RFLBB0   ,RFLBB1   ,RFLBB2   ,RFLBB3,& 
72          RFLCC0   ,RFLCC1   ,RFLCC2   ,RFLCC3   ,RFLDD0,& 
73          RFLDD1   ,RFLDD2   ,RFLDD3   ,RFUETA   ,RASWCA,&
74          RASWCB   ,RASWCC   ,RASWCD   ,RASWCE   ,RASWCF
75!    &    RASWCB   ,RASWCC   ,RASWCD   ,RASWCE   ,RASWCF, RLINLI
76      USE YOERDU   , ONLY : NUAER  ,NTRAER ,REPLOG ,REPSC  ,REPSCW ,DIFF
77!      USE YOETHF   , ONLY : RTICE
78      USE YOERRTWN , ONLY : DELWAVE   ,TOTPLNK     
79      USE YOMPHY3  , ONLY : RII0
80#endif
81      USE aero_mod
82
83  !======================================================================
84  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
85  ! Objet: interface entre le modele et les rayonnements
86  ! Arguments:
87  ! dist-----input-R- distance astronomique terre-soleil
88  ! rmu0-----input-R- cosinus de l'angle zenithal
89  ! fract----input-R- duree d'ensoleillement normalisee
90  ! co2_ppm--input-R- concentration du gaz carbonique (en ppm)
91  ! paprs----input-R- pression a inter-couche (Pa)
92  ! pplay----input-R- pression au milieu de couche (Pa)
93  ! tsol-----input-R- temperature du sol (en K)
94  ! alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible
95  ! alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge   
96  ! t--------input-R- temperature (K)
97  ! q--------input-R- vapeur d'eau (en kg/kg)
98  ! cldfra---input-R- fraction nuageuse (entre 0 et 1)
99  ! cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
100  ! cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
101  ! ok_ade---input-L- apply the Aerosol Direct Effect or not?
102  ! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
103  ! ok_volcan-input-L- activate volcanic diags (SW heat & LW cool rate, SW & LW flux)
104  ! flag_aerosol-input-I- aerosol flag from 0 to 6
105  ! flag_aerosol_strat-input-I- use stratospheric aerosols flag (0, 1, 2)
106  ! flag_aer_feedback-input-I- activate aerosol radiative feedback (T, F)
107  ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
108  ! cldtaupi-input-R- epaisseur optique des nuages dans le visible
109  !                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
110  !                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
111  !                   it is needed for the diagnostics of the aerosol indirect radiative forcing     
112  !
113  ! heat-----output-R- echauffement atmospherique (visible) (K/jour)
114  ! cool-----output-R- refroidissement dans l'IR (K/jour)
115  ! albpla---output-R- albedo planetaire (entre 0 et 1)
116  ! topsw----output-R- flux solaire net au sommet de l'atm.
117  ! toplw----output-R- ray. IR montant au sommet de l'atmosphere
118  ! solsw----output-R- flux solaire net a la surface
119  ! sollw----output-R- ray. IR montant a la surface
120  ! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
121  ! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
122  ! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
123  ! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
124  !
125  ! heat_volc-----output-R- echauffement atmospherique  du au forcage volcanique (visible) (K/s)
126  ! cool_volc-----output-R- refroidissement dans l'IR du au forcage volcanique (K/s)
127  !
128  ! ATTENTION: swai and swad have to be interpreted in the following manner:
129  ! ---------
130  ! ok_ade=F & ok_aie=F -both are zero
131  ! ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
132  !                        indirect is zero
133  ! ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
134  !                        direct is zero
135  ! ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
136  !                        aerosol direct forcing is F_{AD} = topswai-topswad
137  !
138  ! --------- RRTM: output RECMWFL
139  ! ZEMTD (KPROMA,KLEV+1)         ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
140  ! ZEMTU (KPROMA,KLEV+1)         ; TOTAL UPWARD   LONGWAVE EMISSIVITY
141  ! ZTRSO (KPROMA,KLEV+1)         ; TOTAL SHORTWAVE TRANSMISSIVITY
142  ! ZTH   (KPROMA,KLEV+1)         ; HALF LEVEL TEMPERATURE
143  ! ZCTRSO(KPROMA,2)              ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
144  ! ZCEMTR(KPROMA,2)              ; CLEAR-SKY NET LONGWAVE EMISSIVITY
145  ! ZTRSOD(KPROMA)                ; TOTAL-SKY SURFACE SW TRANSMISSITY
146  ! ZLWFC (KPROMA,2)              ; CLEAR-SKY LONGWAVE FLUXES
147  ! ZLWFT (KPROMA,KLEV+1)         ; TOTAL-SKY LONGWAVE FLUXES
148  ! ZLWFT0(KPROMA,KLEV+1)         ; CLEAR-SKY LONGWAVE FLUXES      ! added by MPL 090109
149  ! ZSWFC (KPROMA,2)              ; CLEAR-SKY SHORTWAVE FLUXES
150  ! ZSWFT (KPROMA,KLEV+1)         ; TOTAL-SKY SHORTWAVE FLUXES
151  ! ZSWFT0(KPROMA,KLEV+1)         ; CLEAR-SKY SHORTWAVE FLUXES     ! added by MPL 090109
152  ! ZFLUX (KLON,2,KLEV+1)         ; TOTAL LW FLUXES  1=up, 2=DWN   ! added by MPL 080411
153  ! ZFLUC (KLON,2,KLEV+1)         ; CLEAR SKY LW FLUXES            ! added by MPL 080411
154  ! ZFSDWN(klon,KLEV+1)           ; TOTAL SW  DWN FLUXES           ! added by MPL 080411
155  ! ZFCDWN(klon,KLEV+1)           ; CLEAR SKY SW  DWN FLUXES       ! added by MPL 080411
156  ! ZFCCDWN(klon,KLEV+1)          ; CLEAR SKY CLEAN (NO AEROSOL) SW  DWN FLUXES      ! added by OB 211117
157  ! ZFSUP (klon,KLEV+1)           ; TOTAL SW  UP  FLUXES           ! added by MPL 080411
158  ! ZFCUP (klon,KLEV+1)           ; CLEAR SKY SW  UP  FLUXES       ! added by MPL 080411
159  ! ZFCCUP (klon,KLEV+1)          ; CLEAR SKY CLEAN (NO AEROSOL) SW  UP  FLUXES      ! added by OB 211117
160  ! ZFLCCDWN(klon,KLEV+1)         ; CLEAR SKY CLEAN (NO AEROSOL) LW  DWN FLUXES      ! added by OB 211117
161  ! ZFLCCUP (klon,KLEV+1)         ; CLEAR SKY CLEAN (NO AEROSOL) LW  UP  FLUXES      ! added by OB 211117
162 
163  !======================================================================
164 
165  ! ====================================================================
166  ! Adapte au modele de chimie INCA par Celine Deandreis & Anne Cozic -- 2009
167  ! 1 = ZERO   
168  ! 2 = AER total   
169  ! 3 = NAT   
170  ! 4 = BC   
171  ! 5 = SO4   
172  ! 6 = POM   
173  ! 7 = DUST   
174  ! 8 = SS   
175  ! 9 = NO3   
176  !
177  ! ====================================================================
178  include "YOETHF.h"
179  include "YOMCST.h"
180  include "clesphys.h"
181
182! Input arguments
183  REAL,    INTENT(in)  :: dist
184  REAL,    INTENT(in)  :: rmu0(KLON), fract(KLON)
185  REAL,    INTENT(in)  :: paprs(KLON,KLEV+1), pplay(KLON,KLEV)
186!albedo SB >>>
187! REAL,    INTENT(in)  :: alb1(KLON), alb2(KLON), tsol(KLON)
188  REAL,    INTENT(in)  :: tsol(KLON)
189  REAL,    INTENT(in) :: alb_dir(KLON,NSW),alb_dif(KLON,NSW)
190  real, intent(in) :: SFRWL(6)
191!albedo SB <<<
192  REAL,    INTENT(in)  :: t(KLON,KLEV), q(KLON,KLEV)
193
194  REAL, INTENT(in):: wo(:, :, :) ! dimension(KLON,KLEV, 1 or 2)
195  ! column-density of ozone in a layer, in kilo-Dobsons
196  ! "wo(:, :, 1)" is for the average day-night field,
197  ! "wo(:, :, 2)" is for daylight time.
198
199  LOGICAL, INTENT(in)  :: ok_ade, ok_aie                                 ! switches whether to use aerosol direct (indirect) effects or not
200  LOGICAL, INTENT(in)  :: ok_volcan                                      ! produce volcanic diags (SW/LW heat flux and rate)
201  LOGICAL              :: lldebug
202  INTEGER, INTENT(in)  :: flag_aerosol                                   ! takes value 0 (no aerosol) or 1 to 6 (aerosols)
203  INTEGER, INTENT(in)  :: flag_aerosol_strat                             ! use stratospheric aerosols
204  LOGICAL, INTENT(in)  :: flag_aer_feedback                              ! activate aerosol radiative feedback
205  REAL,    INTENT(in)  :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV)
206  REAL,    INTENT(in)  :: tau_aero(KLON,KLEV,naero_grp,2)                        ! aerosol optical properties (see aeropt.F)
207  REAL,    INTENT(in)  :: piz_aero(KLON,KLEV,naero_grp,2)                        ! aerosol optical properties (see aeropt.F)
208  REAL,    INTENT(in)  :: cg_aero(KLON,KLEV,naero_grp,2)                         ! aerosol optical properties (see aeropt.F)
209!--OB
210  REAL,    INTENT(in)  :: tau_aero_sw_rrtm(KLON,KLEV,2,NSW)                 ! aerosol optical properties RRTM
211  REAL,    INTENT(in)  :: piz_aero_sw_rrtm(KLON,KLEV,2,NSW)                 ! aerosol optical properties RRTM
212  REAL,    INTENT(in)  :: cg_aero_sw_rrtm(KLON,KLEV,2,NSW)                  ! aerosol optical properties RRTM
213!--OB fin
214
215!--C. Kleinschmitt
216#ifdef CPP_RRTM
217  REAL,    INTENT(in)  :: tau_aero_lw_rrtm(KLON,KLEV,2,NLW)                 ! LW aerosol optical properties RRTM
218#else
219  REAL,    INTENT(in)  :: tau_aero_lw_rrtm(KLON,KLEV,2,nbands_lw_rrtm)
220#endif
221!--C. Kleinschmitt end
222
223  REAL,    INTENT(in)  :: cldtaupi(KLON,KLEV)                            ! cloud optical thickness for pre-industrial aerosol concentrations
224  REAL,    INTENT(in)  :: qsat(klon,klev) ! Variable pour iflag_rrtm=1
225  REAL,    INTENT(in)  :: flwc(klon,klev) ! Variable pour iflag_rrtm=1
226  REAL,    INTENT(in)  :: fiwc(klon,klev) ! Variable pour iflag_rrtm=1
227  REAL,    INTENT(in)  :: ref_liq(klon,klev) ! cloud droplet radius present-day from newmicro
228  REAL,    INTENT(in)  :: ref_ice(klon,klev) ! ice crystal radius   present-day from newmicro
229  REAL,    INTENT(in)  :: ref_liq_pi(klon,klev) ! cloud droplet radius pre-industrial from newmicro
230  REAL,    INTENT(in)  :: ref_ice_pi(klon,klev) ! ice crystal radius   pre-industrial from newmicro
231
232! Output arguments
233  REAL,    INTENT(out) :: heat(KLON,KLEV), cool(KLON,KLEV)
234  REAL,    INTENT(out) :: heat0(KLON,KLEV), cool0(KLON,KLEV)
235  REAL,    INTENT(out) :: heat_volc(KLON,KLEV), cool_volc(KLON,KLEV) !NL
236  REAL,    INTENT(out) :: topsw(KLON), toplw(KLON)
237  REAL,    INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON)
238  REAL,    INTENT(out) :: topsw0(KLON), toplw0(KLON), solsw0(KLON), sollw0(KLON)
239  REAL,    INTENT(out) :: sollwdown(KLON)
240  REAL,    INTENT(out) :: swdn(KLON,kflev+1),swdn0(KLON,kflev+1), swdnc0(KLON,kflev+1)
241  REAL,    INTENT(out) :: swup(KLON,kflev+1),swup0(KLON,kflev+1), swupc0(KLON,kflev+1)
242  REAL,    INTENT(out) :: lwdn(KLON,kflev+1),lwdn0(KLON,kflev+1), lwdnc0(KLON,kflev+1)
243  REAL,    INTENT(out) :: lwup(KLON,kflev+1),lwup0(KLON,kflev+1), lwupc0(KLON,kflev+1)
244  REAL,    INTENT(out) :: topswad_aero(KLON), solswad_aero(KLON)         ! output: aerosol direct forcing at TOA and surface
245  REAL,    INTENT(out) :: topswai_aero(KLON), solswai_aero(KLON)         ! output: aerosol indirect forcing atTOA and surface
246  REAL,    INTENT(out) :: toplwad_aero(KLON), sollwad_aero(KLON)         ! output: LW aerosol direct forcing at TOA and surface
247  REAL,    INTENT(out) :: toplwai_aero(KLON), sollwai_aero(KLON)         ! output: LW aerosol indirect forcing atTOA and surface
248  REAL, DIMENSION(klon), INTENT(out)    :: topswad0_aero
249  REAL, DIMENSION(klon), INTENT(out)    :: solswad0_aero
250  REAL, DIMENSION(klon), INTENT(out)    :: toplwad0_aero
251  REAL, DIMENSION(klon), INTENT(out)    :: sollwad0_aero
252  REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw_aero
253  REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw0_aero
254  REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw_aero
255  REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw0_aero
256  REAL, DIMENSION(kdlon,3), INTENT(out) :: topswcf_aero
257  REAL, DIMENSION(kdlon,3), INTENT(out) :: solswcf_aero
258  REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZSWFT0_i
259  REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZLWFT0_i
260
261! Local variables
262  REAL(KIND=8) ZFSUP(KDLON,KFLEV+1)
263  REAL(KIND=8) ZFSDN(KDLON,KFLEV+1)
264  REAL(KIND=8) ZFSUP0(KDLON,KFLEV+1)
265  REAL(KIND=8) ZFSDN0(KDLON,KFLEV+1)
266  REAL(KIND=8) ZFSUPC0(KDLON,KFLEV+1)
267  REAL(KIND=8) ZFSDNC0(KDLON,KFLEV+1)
268  REAL(KIND=8) ZFLUP(KDLON,KFLEV+1)
269  REAL(KIND=8) ZFLDN(KDLON,KFLEV+1)
270  REAL(KIND=8) ZFLUP0(KDLON,KFLEV+1)
271  REAL(KIND=8) ZFLDN0(KDLON,KFLEV+1)
272  REAL(KIND=8) ZFLUPC0(KDLON,KFLEV+1)
273  REAL(KIND=8) ZFLDNC0(KDLON,KFLEV+1)
274  REAL(KIND=8) zx_alpha1, zx_alpha2
275  INTEGER k, kk, i, j, iof, nb_gr
276  INTEGER ist,iend,ktdia,kmode
277  REAL(KIND=8) PSCT
278  REAL(KIND=8) PALBD(kdlon,2), PALBP(kdlon,2)
279!  MPL 06.01.09: pour RRTM, creation de PALBD_NEW et PALBP_NEW
280! avec NSW en deuxieme dimension       
281  REAL(KIND=8) PALBD_NEW(kdlon,NSW), PALBP_NEW(kdlon,NSW)
282  REAL(KIND=8) PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
283  REAL(KIND=8) PPSOL(kdlon), PDP(kdlon,KLEV)
284  REAL(KIND=8) PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
285  REAL(KIND=8) PTAVE(kdlon,kflev)
286  REAL(KIND=8) PWV(kdlon,kflev), PQS(kdlon,kflev)
287
288  real(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
289  ! "POZON(:, :, 1)" is for the average day-night field,
290  ! "POZON(:, :, 2)" is for daylight time.
291!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6 
292  REAL(KIND=8) PAER(kdlon,kflev,6)
293  REAL(KIND=8) PCLDLD(kdlon,kflev)
294  REAL(KIND=8) PCLDLU(kdlon,kflev)
295  REAL(KIND=8) PCLDSW(kdlon,kflev)
296  REAL(KIND=8) PTAU(kdlon,2,kflev)
297  REAL(KIND=8) POMEGA(kdlon,2,kflev)
298  REAL(KIND=8) PCG(kdlon,2,kflev)
299  REAL(KIND=8) zfract(kdlon), zrmu0(kdlon), zdist
300  REAL(KIND=8) zheat(kdlon,kflev), zcool(kdlon,kflev)
301  REAL(KIND=8) zheat0(kdlon,kflev), zcool0(kdlon,kflev)
302  REAL(KIND=8) zheat_volc(kdlon,kflev), zcool_volc(kdlon,kflev) !NL
303  REAL(KIND=8) ztopsw(kdlon), ztoplw(kdlon)
304  REAL(KIND=8) zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
305  REAL(KIND=8) zsollwdown(kdlon)
306  REAL(KIND=8) ztopsw0(kdlon), ztoplw0(kdlon)
307  REAL(KIND=8) zsolsw0(kdlon), zsollw0(kdlon)
308  REAL(KIND=8) zznormcp
309  REAL(KIND=8) tauaero(kdlon,kflev,naero_grp,2)                     ! aer opt properties
310  REAL(KIND=8) pizaero(kdlon,kflev,naero_grp,2)
311  REAL(KIND=8) cgaero(kdlon,kflev,naero_grp,2)
312  REAL(KIND=8) PTAUA(kdlon,2,kflev)                         ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
313  REAL(KIND=8) POMEGAA(kdlon,2,kflev)                       ! dito for single scatt albedo
314  REAL(KIND=8) ztopswadaero(kdlon), zsolswadaero(kdlon)     ! Aerosol direct forcing at TOAand surface
315  REAL(KIND=8) ztopswad0aero(kdlon), zsolswad0aero(kdlon)   ! Aerosol direct forcing at TOAand surface
316  REAL(KIND=8) ztopswaiaero(kdlon), zsolswaiaero(kdlon)     ! dito, indirect
317!--NL
318  REAL(KIND=8) zswadaero(kdlon,kflev+1)                       ! SW Aerosol direct forcing
319  REAL(KIND=8) zlwadaero(kdlon,kflev+1)                       ! LW Aerosol direct forcing
320!-LW by CK
321  REAL(KIND=8) ztoplwadaero(kdlon), zsollwadaero(kdlon)     ! LW Aerosol direct forcing at TOAand surface
322  REAL(KIND=8) ztoplwad0aero(kdlon), zsollwad0aero(kdlon)   ! LW Aerosol direct forcing at TOAand surface
323  REAL(KIND=8) ztoplwaiaero(kdlon), zsollwaiaero(kdlon)     ! dito, indirect
324!-end
325  REAL(KIND=8) ztopsw_aero(kdlon,9), ztopsw0_aero(kdlon,9)
326  REAL(KIND=8) zsolsw_aero(kdlon,9), zsolsw0_aero(kdlon,9)
327  REAL(KIND=8) ztopswcf_aero(kdlon,3), zsolswcf_aero(kdlon,3)     
328! real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2 deje declare dans physiq.F MPL 20130618
329!MPL input supplementaires pour RECMWFL
330! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
331      REAL(KIND=8) GEMU(klon)
332!MPL input RECMWFL:
333! Tableaux aux niveaux inverses pour respecter convention Arpege
334      REAL(KIND=8) ref_liq_i(klon,klev) ! cloud droplet radius present-day from newmicro (inverted)
335      REAL(KIND=8) ref_ice_i(klon,klev) ! ice crystal radius present-day from newmicro (inverted)
336!--OB
337      REAL(KIND=8) ref_liq_pi_i(klon,klev) ! cloud droplet radius pre-industrial from newmicro (inverted)
338      REAL(KIND=8) ref_ice_pi_i(klon,klev) ! ice crystal radius pre-industrial from newmicro (inverted)
339!--end OB
340      REAL(KIND=8) paprs_i(klon,klev+1)
341      REAL(KIND=8) pplay_i(klon,klev)
342      REAL(KIND=8) cldfra_i(klon,klev)
343      REAL(KIND=8) POZON_i(kdlon,kflev, size(wo, 3)) ! mass fraction of ozone
344  ! "POZON(:, :, 1)" is for the average day-night field,
345  ! "POZON(:, :, 2)" is for daylight time.
346!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6     
347      REAL(KIND=8) PAER_i(kdlon,kflev,6)
348      REAL(KIND=8) PDP_i(klon,klev)
349      REAL(KIND=8) t_i(klon,klev),q_i(klon,klev),qsat_i(klon,klev)
350      REAL(KIND=8) flwc_i(klon,klev),fiwc_i(klon,klev)
351!MPL output RECMWFL:
352      REAL(KIND=8) ZEMTD (klon,klev+1),ZEMTD_i (klon,klev+1)       
353      REAL(KIND=8) ZEMTU (klon,klev+1),ZEMTU_i (klon,klev+1)     
354      REAL(KIND=8) ZTRSO (klon,klev+1),ZTRSO_i (klon,klev+1)   
355      REAL(KIND=8) ZTH   (klon,klev+1),ZTH_i   (klon,klev+1)   
356      REAL(KIND=8) ZCTRSO(klon,2)       
357      REAL(KIND=8) ZCEMTR(klon,2)     
358      REAL(KIND=8) ZTRSOD(klon)       
359      REAL(KIND=8) ZLWFC (klon,2)     
360      REAL(KIND=8) ZLWFT (klon,klev+1),ZLWFT_i (klon,klev+1)   
361      REAL(KIND=8) ZSWFC (klon,2)     
362      REAL(KIND=8) ZSWFT (klon,klev+1),ZSWFT_i (klon,klev+1)
363      REAL(KIND=8) ZFLUCDWN_i(klon,klev+1),ZFLUCUP_i(klon,klev+1)
364      REAL(KIND=8) PPIZA_TOT(klon,klev,NSW)
365      REAL(KIND=8) PCGA_TOT(klon,klev,NSW)
366      REAL(KIND=8) PTAU_TOT(klon,klev,NSW)
367      REAL(KIND=8) PPIZA_NAT(klon,klev,NSW)
368      REAL(KIND=8) PCGA_NAT(klon,klev,NSW)
369      REAL(KIND=8) PTAU_NAT(klon,klev,NSW)
370#ifdef CPP_RRTM
371      REAL(KIND=8) PTAU_LW_TOT(klon,klev,NLW)
372      REAL(KIND=8) PTAU_LW_NAT(klon,klev,NLW)
373#endif
374      REAL(KIND=8) PSFSWDIR(klon,NSW)
375      REAL(KIND=8) PSFSWDIF(klon,NSW)
376      REAL(KIND=8) PFSDNN(klon)
377      REAL(KIND=8) PFSDNV(klon)
378!MPL On ne redefinit pas les tableaux ZFLUX,ZFLUC,
379!MPL ZFSDWN,ZFCDWN,ZFSUP,ZFCUP car ils existent deja
380!MPL sous les noms de ZFLDN,ZFLDN0,ZFLUP,ZFLUP0,
381!MPL ZFSDN,ZFSDN0,ZFSUP,ZFSUP0
382      REAL(KIND=8) ZFLUX_i (klon,2,klev+1)
383      REAL(KIND=8) ZFLUC_i (klon,2,klev+1)
384      REAL(KIND=8) ZFSDWN_i (klon,klev+1)
385      REAL(KIND=8) ZFCDWN_i (klon,klev+1)
386      REAL(KIND=8) ZFCCDWN_i (klon,klev+1)
387      REAL(KIND=8) ZFSUP_i (klon,klev+1)
388      REAL(KIND=8) ZFCUP_i (klon,klev+1)
389      REAL(KIND=8) ZFCCUP_i (klon,klev+1)
390      REAL(KIND=8) ZFLCCDWN_i (klon,klev+1)
391      REAL(KIND=8) ZFLCCUP_i (klon,klev+1)
392! 3 lignes suivantes a activer pour CCMVAL (MPL 20100412)
393!      REAL(KIND=8) RSUN(3,2)
394!      REAL(KIND=8) SUN(3)
395!      REAL(KIND=8) SUN_FRACT(2)
396  real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
397  CHARACTER (LEN=80) :: abort_message
398  CHARACTER (LEN=80) :: modname='radlwsw_m'
399
400  call assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
401  ! initialisation
402  ist=1
403  iend=klon
404  ktdia=1
405  kmode=ist
406  tauaero(:,:,:,:)=0.
407  pizaero(:,:,:,:)=0.
408  cgaero(:,:,:,:)=0.
409  lldebug=.FALSE.
410
411  ztopsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
412  ztopsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
413  zsolsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
414  zsolsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
415
416
417   ZTOPSWADAERO(:)  = 0. !ym missing init
418   ZSOLSWADAERO(:)  = 0. !ym missing init
419   ZTOPSWAD0AERO(:) = 0. !ym missing init
420   ZSOLSWAD0AERO(:) = 0. !ym missing init
421   ZTOPSWAIAERO(:)  = 0. !ym missing init
422   ZSOLSWAIAERO(:)  = 0. !ym missing init 
423   ZTOPSWCF_AERO(:,:)= 0.!ym missing init 
424   ZSOLSWCF_AERO(:,:) =0. !ym missing init 
425
426  !
427  !-------------------------------------------
428  nb_gr = KLON / kdlon
429  IF (nb_gr*kdlon .NE. KLON) THEN
430      PRINT*, "kdlon mauvais:", KLON, kdlon, nb_gr
431      call abort_physic("radlwsw", "", 1)
432  ENDIF
433  IF (kflev .NE. KLEV) THEN
434      PRINT*, "kflev differe de KLEV, kflev, KLEV"
435      call abort_physic("radlwsw", "", 1)
436  ENDIF
437  !-------------------------------------------
438  DO k = 1, KLEV
439    DO i = 1, KLON
440      heat(i,k)=0.
441      cool(i,k)=0.
442      heat_volc(i,k)=0. !NL
443      cool_volc(i,k)=0. !NL
444      heat0(i,k)=0.
445      cool0(i,k)=0.
446    ENDDO
447  ENDDO
448  !
449  zdist = dist
450  !
451  PSCT = solaire/zdist/zdist
452
453  IF (type_trac == 'repr') THEN
454#ifdef REPROBUS
455     if(ok_SUNTIME) PSCT = solaireTIME/zdist/zdist
456     print*,'Constante solaire: ',PSCT*zdist*zdist
457#endif
458  END IF
459
460  DO j = 1, nb_gr
461    iof = kdlon*(j-1)
462    DO i = 1, kdlon
463      zfract(i) = fract(iof+i)
464!     zfract(i) = 1.     !!!!!!  essai MPL 19052010
465      zrmu0(i) = rmu0(iof+i)
466
467
468!albedo SB >>>
469!
470      IF (iflag_rrtm==0) THEN
471!
472        PALBD(i,1)=alb_dif(iof+i,1)
473        PALBD(i,2)=alb_dif(iof+i,2)
474        PALBP(i,1)=alb_dir(iof+i,1)
475        PALBP(i,2)=alb_dir(iof+i,2)
476!
477      ELSEIF (iflag_rrtm==1) THEn
478!
479        DO kk=1,NSW
480          PALBD_NEW(i,kk)=alb_dif(iof+i,kk)
481          PALBP_NEW(i,kk)=alb_dir(iof+i,kk)
482        ENDDO
483!
484      ENDIF
485!albedo SB <<<
486
487
488      PEMIS(i) = 1.0    !!!!! A REVOIR (MPL)
489      PVIEW(i) = 1.66
490      PPSOL(i) = paprs(iof+i,1)
491      zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2))
492      zx_alpha2 = 1.0 - zx_alpha1
493      PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
494      PTL(i,KLEV+1) = t(iof+i,KLEV)
495      PDT0(i) = tsol(iof+i) - PTL(i,1)
496    ENDDO
497    DO k = 2, kflev
498      DO i = 1, kdlon
499        PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
500      ENDDO
501    ENDDO
502    DO k = 1, kflev
503      DO i = 1, kdlon
504        PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
505        PTAVE(i,k) = t(iof+i,k)
506        PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
507        PQS(i,k) = PWV(i,k)
508!       Confert from  column density of ozone in a cell, in kDU, to a mass fraction
509        POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
510             / (paprs(iof+i, k) - paprs(iof+i, k+1))
511!       A activer pour CCMVAL on prend l'ozone impose (MPL 07042010)
512!       POZON(i,k,:) = wo(i,k,:) 
513!       print *,'RADLWSW: POZON',k, POZON(i,k,1)
514        PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
515        PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
516        PCLDSW(i,k) = cldfra(iof+i,k)
517        PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
518        PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
519        POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
520        POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
521        PCG(i,1,k) = 0.865
522        PCG(i,2,k) = 0.910
523        !-
524        ! Introduced for aerosol indirect forcings.
525        ! The following values use the cloud optical thickness calculated from
526        ! present-day aerosol concentrations whereas the quantities without the
527        ! "A" at the end are for pre-industial (natural-only) aerosol concentrations
528        !
529        PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
530        PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
531        POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
532        POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
533      ENDDO
534    ENDDO
535
536    IF (type_trac == 'repr') THEN
537#ifdef REPROBUS
538       ndimozon = size(wo, 3)
539       CALL RAD_INTERACTIF(POZON,iof)
540#endif
541    END IF
542
543    !
544    DO k = 1, kflev+1
545      DO i = 1, kdlon
546        PPMB(i,k) = paprs(iof+i,k)/100.0
547      ENDDO
548    ENDDO
549    !
550!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6
551    DO kk = 1, 6
552      DO k = 1, kflev
553        DO i = 1, kdlon
554          PAER(i,k,kk) = 1.0E-15   !!!!! A REVOIR (MPL)
555        ENDDO
556      ENDDO
557    ENDDO
558    DO k = 1, kflev
559      DO i = 1, kdlon
560        tauaero(i,k,:,1)=tau_aero(iof+i,k,:,1)
561        pizaero(i,k,:,1)=piz_aero(iof+i,k,:,1)
562        cgaero(i,k,:,1) =cg_aero(iof+i,k,:,1)
563        tauaero(i,k,:,2)=tau_aero(iof+i,k,:,2)
564        pizaero(i,k,:,2)=piz_aero(iof+i,k,:,2)
565        cgaero(i,k,:,2) =cg_aero(iof+i,k,:,2)
566      ENDDO
567    ENDDO
568
569!
570!===== iflag_rrtm ================================================
571!     
572    IF (iflag_rrtm == 0) THEN       !!!! remettre 0 juste pour tester l'ancien rayt via rrtm
573!--- Mise a zero des tableaux output du rayonnement LW-AR4 ----------             
574      DO k = 1, kflev+1
575      DO i = 1, kdlon
576!     print *,'RADLWSW: boucle mise a zero i k',i,k
577      ZFLUP(i,k)=0.
578      ZFLDN(i,k)=0.
579      ZFLUP0(i,k)=0.
580      ZFLDN0(i,k)=0.
581      ZLWFT0_i(i,k)=0.
582      ZFLUCUP_i(i,k)=0.
583      ZFLUCDWN_i(i,k)=0.
584      ENDDO
585      ENDDO
586      DO k = 1, kflev
587         DO i = 1, kdlon
588            zcool(i,k)=0.
589            zcool_volc(i,k)=0. !NL
590            zcool0(i,k)=0.
591         ENDDO
592      ENDDO
593      DO i = 1, kdlon
594      ztoplw(i)=0.
595      zsollw(i)=0.
596      ztoplw0(i)=0.
597      zsollw0(i)=0.
598      zsollwdown(i)=0.
599      ENDDO
600       ! Old radiation scheme, used for AR4 runs
601       ! average day-night ozone for longwave
602       CALL LW_LMDAR4(&
603            PPMB, PDP,&
604            PPSOL,PDT0,PEMIS,&
605            PTL, PTAVE, PWV, POZON(:, :, 1), PAER,&
606            PCLDLD,PCLDLU,&
607            PVIEW,&
608            zcool, zcool0,&
609            ztoplw,zsollw,ztoplw0,zsollw0,&
610            zsollwdown,&
611            ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
612!----- Mise a zero des tableaux output du rayonnement SW-AR4
613      DO k = 1, kflev+1
614         DO i = 1, kdlon
615            ZFSUP(i,k)=0.
616            ZFSDN(i,k)=0.
617            ZFSUP0(i,k)=0.
618            ZFSDN0(i,k)=0.
619            ZFSUPC0(i,k)=0.
620            ZFSDNC0(i,k)=0.
621            ZFLUPC0(i,k)=0.
622            ZFLDNC0(i,k)=0.
623            ZSWFT0_i(i,k)=0.
624            ZFCUP_i(i,k)=0.
625            ZFCDWN_i(i,k)=0.
626            ZFCCUP_i(i,k)=0.
627            ZFCCDWN_i(i,k)=0.
628            ZFLCCUP_i(i,k)=0.
629            ZFLCCDWN_i(i,k)=0.
630            zswadaero(i,k)=0. !--NL
631         ENDDO
632      ENDDO
633      DO k = 1, kflev
634         DO i = 1, kdlon
635            zheat(i,k)=0.
636            zheat_volc(i,k)=0.
637            zheat0(i,k)=0.
638         ENDDO
639      ENDDO
640      DO i = 1, kdlon
641      zalbpla(i)=0.
642      ztopsw(i)=0.
643      zsolsw(i)=0.
644      ztopsw0(i)=0.
645      zsolsw0(i)=0.
646      ztopswadaero(i)=0.
647      zsolswadaero(i)=0.
648      ztopswaiaero(i)=0.
649      zsolswaiaero(i)=0.
650      ENDDO
651!     print *,'Avant SW_LMDAR4: PSCT zrmu0 zfract',PSCT, zrmu0, zfract
652       ! daylight ozone, if we have it, for short wave
653      CALL SW_AEROAR4(PSCT, zrmu0, zfract,&
654               PPMB, PDP,&
655               PPSOL, PALBD, PALBP,&
656               PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
657               PCLDSW, PTAU, POMEGA, PCG,&
658               zheat, zheat0,&
659               zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
660               ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
661               tauaero, pizaero, cgaero, &
662               PTAUA, POMEGAA,&
663               ztopswadaero,zsolswadaero,&
664               ztopswad0aero,zsolswad0aero,&
665               ztopswaiaero,zsolswaiaero, &
666               ztopsw_aero,ztopsw0_aero,&
667               zsolsw_aero,zsolsw0_aero,&
668               ztopswcf_aero,zsolswcf_aero, &
669               ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat)
670
671       ZSWFT0_i(:,:) = ZFSDN0(:,:)-ZFSUP0(:,:)
672       ZLWFT0_i(:,:) =-ZFLDN0(:,:)-ZFLUP0(:,:)
673
674       DO i=1,kdlon
675       DO k=1,kflev+1
676!        print *,'iof i k klon klev=',iof,i,k,klon,klev
677         lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
678         lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
679         lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
680         lwup  ( iof+i,k)   = ZFLUP  ( i,k)
681         swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
682         swdn  ( iof+i,k)   = ZFSDN  ( i,k)
683         swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
684         swup  ( iof+i,k)   = ZFSUP  ( i,k)
685       ENDDO 
686       ENDDO 
687!          print*,'SW_AR4 ZFSDN0 1 , klev:',ZFSDN0(1:klon,1),ZFSDN0(1:klon,klev)
688!          print*,'SW_AR4 swdn0  1 , klev:',swdn0(1:klon,1),swdn0(1:klon,klev)
689!          print*,'SW_AR4 ZFSUP0 1 , klev:',ZFSUP0(1:klon,1),ZFSUP0(1:klon,klev)
690!          print*,'SW_AR4 swup0  1 , klev:',swup0(1:klon,1),swup0(1:klon,klev)
691!          print*,'SW_AR4 ZFSDN  1 , klev:',ZFSDN(1:klon,1) ,ZFSDN(1:klon,klev)
692!          print*,'SW_AR4 ZFSUP  1 , klev:',ZFSUP(1:klon,1) ,ZFSUP(1:klon,klev)
693    ELSE 
694#ifdef CPP_RRTM
695!      if (prt_level.gt.10)write(lunout,*)'CPP_RRTM=.T.'
696!===== iflag_rrtm=1, on passe dans SW via RECMWFL ===============
697
698      DO k = 1, kflev+1
699      DO i = 1, kdlon
700      ZEMTD_i(i,k)=0.
701      ZEMTU_i(i,k)=0.
702      ZTRSO_i(i,k)=0.
703      ZTH_i(i,k)=0.
704      ZLWFT_i(i,k)=0.
705      ZSWFT_i(i,k)=0.
706      ZFLUX_i(i,1,k)=0.
707      ZFLUX_i(i,2,k)=0.
708      ZFLUC_i(i,1,k)=0.
709      ZFLUC_i(i,2,k)=0.
710      ZFSDWN_i(i,k)=0.
711      ZFCDWN_i(i,k)=0.
712      ZFCCDWN_i(i,k)=0.
713      ZFSUP_i(i,k)=0.
714      ZFCUP_i(i,k)=0.
715      ZFCCUP_i(i,k)=0.
716      ZFLCCDWN_i(i,k)=0.
717      ZFLCCUP_i(i,k)=0.
718      ENDDO
719      ENDDO
720!
721!--OB
722!--aerosol TOT  - anthropogenic+natural - index 2
723!--aerosol NAT  - natural only          - index 1
724!
725      DO i = 1, kdlon
726      DO k = 1, kflev
727      DO kk=1, NSW
728!
729      PTAU_TOT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,2,kk)
730      PPIZA_TOT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,2,kk)
731      PCGA_TOT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,2,kk)
732!
733      PTAU_NAT(i,kflev+1-k,kk)=tau_aero_sw_rrtm(i,k,1,kk)
734      PPIZA_NAT(i,kflev+1-k,kk)=piz_aero_sw_rrtm(i,k,1,kk)
735      PCGA_NAT(i,kflev+1-k,kk)=cg_aero_sw_rrtm(i,k,1,kk)
736!
737      ENDDO
738      ENDDO
739      ENDDO
740!-end OB
741!
742!--C. Kleinschmitt
743!--aerosol TOT  - anthropogenic+natural - index 2
744!--aerosol NAT  - natural only          - index 1
745!
746      DO i = 1, kdlon
747      DO k = 1, kflev
748      DO kk=1, NLW
749!
750      PTAU_LW_TOT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,2,kk)
751      PTAU_LW_NAT(i,kflev+1-k,kk)=tau_aero_lw_rrtm(i,k,1,kk)
752!
753      ENDDO
754      ENDDO
755      ENDDO
756!-end C. Kleinschmitt
757!     
758      DO i = 1, kdlon
759      ZCTRSO(i,1)=0.
760      ZCTRSO(i,2)=0.
761      ZCEMTR(i,1)=0.
762      ZCEMTR(i,2)=0.
763      ZTRSOD(i)=0.
764      ZLWFC(i,1)=0.
765      ZLWFC(i,2)=0.
766      ZSWFC(i,1)=0.
767      ZSWFC(i,2)=0.
768      PFSDNN(i)=0.
769      PFSDNV(i)=0.
770      DO kk = 1, NSW
771      PSFSWDIR(i,kk)=0.
772      PSFSWDIF(i,kk)=0.
773      ENDDO
774      ENDDO
775!----- Fin des mises a zero des tableaux output de RECMWF -------------------             
776!        GEMU(1:klon)=sin(rlatd(1:klon))
777! On met les donnees dans l'ordre des niveaux arpege
778         paprs_i(:,1)=paprs(:,klev+1)
779         do k=1,klev
780            paprs_i(1:klon,k+1) =paprs(1:klon,klev+1-k)
781            pplay_i(1:klon,k)   =pplay(1:klon,klev+1-k)
782            cldfra_i(1:klon,k)  =cldfra(1:klon,klev+1-k)
783            PDP_i(1:klon,k)     =PDP(1:klon,klev+1-k)
784            t_i(1:klon,k)       =t(1:klon,klev+1-k)
785            q_i(1:klon,k)       =q(1:klon,klev+1-k)
786            qsat_i(1:klon,k)    =qsat(1:klon,klev+1-k)
787            flwc_i(1:klon,k)    =flwc(1:klon,klev+1-k)
788            fiwc_i(1:klon,k)    =fiwc(1:klon,klev+1-k)
789            ref_liq_i(1:klon,k) =ref_liq(1:klon,klev+1-k)
790            ref_ice_i(1:klon,k) =ref_ice(1:klon,klev+1-k)
791!-OB
792            ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
793            ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
794         enddo
795         do k=1,kflev
796           POZON_i(1:klon,k,:)=POZON(1:klon,kflev+1-k,:)
797!!!            POZON_i(1:klon,k)=POZON(1:klon,k)            !!! on laisse 1=sol et klev=top
798!          print *,'Juste avant RECMWFL: k tsol temp',k,tsol,t(1,k)
799!!!!!!! Modif MPL 6.01.09 avec RRTM, on passe de 5 a 6     
800            do i=1,6
801            PAER_i(1:klon,k,i)=PAER(1:klon,kflev+1-k,i)
802            enddo
803         enddo
804!       print *,'RADLWSW: avant RECMWFL, RI0,rmu0=',solaire,rmu0
805
806!  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
807! La version ARPEGE1D utilise differentes valeurs de la constante
808! solaire suivant le rayonnement utilise.
809! A controler ...
810! SOLAR FLUX AT THE TOP (/YOMPHY3/)
811! introduce season correction
812!--------------------------------------
813! RII0 = RIP0
814! IF(LRAYFM)
815! RII0 = RIP0M   ! =rip0m if Morcrette non-each time step call.
816! IF(LRAYFM15)
817! RII0 = RIP0M15 ! =rip0m if Morcrette non-each time step call.
818         RII0=solaire/zdist/zdist
819!print*,'+++ radlwsw: solaire ,RII0',solaire,RII0
820!  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
821! Ancien appel a RECMWF (celui du cy25)
822!        CALL RECMWF (ist , iend, klon , ktdia , klev   , kmode ,
823!    s   PALBD    , PALBP   , paprs_i , pplay_i , RCO2   , cldfra_i,
824!    s   POZON_i  , PAER_i  , PDP_i   , PEMIS   , GEMU   , rmu0,
825!    s    q_i     , qsat_i  , fiwc_i  , flwc_i  , zmasq  , t_i  ,tsol,
826!    s   ZEMTD_i  , ZEMTU_i , ZTRSO_i ,
827!    s   ZTH_i    , ZCTRSO  , ZCEMTR  , ZTRSOD  ,
828!    s   ZLWFC    , ZLWFT_i , ZSWFC   , ZSWFT_i ,
829!    s   ZFLUX_i  , ZFLUC_i , ZFSDWN_i, ZFSUP_i , ZFCDWN_i,ZFCUP_i)
830!    s   'RECMWF ')
831!
832      if(lldebug) then
833        CALL writefield_phy('paprs_i',paprs_i,klev+1)
834        CALL writefield_phy('pplay_i',pplay_i,klev)
835        CALL writefield_phy('cldfra_i',cldfra_i,klev)
836        CALL writefield_phy('pozon_i',POZON_i,klev)
837        CALL writefield_phy('paer_i',PAER_i,klev)
838        CALL writefield_phy('pdp_i',PDP_i,klev)
839        CALL writefield_phy('q_i',q_i,klev)
840        CALL writefield_phy('qsat_i',qsat_i,klev)
841        CALL writefield_phy('fiwc_i',fiwc_i,klev)
842        CALL writefield_phy('flwc_i',flwc_i,klev)
843        CALL writefield_phy('t_i',t_i,klev)
844        CALL writefield_phy('palbd_new',PALBD_NEW,NSW)
845        CALL writefield_phy('palbp_new',PALBP_NEW,NSW)
846      endif
847
848! Nouvel appel a RECMWF (celui du cy32t0)
849         CALL RECMWF_AERO (ist , iend, klon , ktdia  , klev   , kmode ,&
850         PALBD_NEW,PALBP_NEW, paprs_i , pplay_i , RCO2   , cldfra_i,&
851         POZON_i  , PAER_i  , PDP_i   , PEMIS   , rmu0   ,&
852          q_i     , qsat_i  , fiwc_i  , flwc_i  , zmasq  , t_i  ,tsol,&
853         ref_liq_i, ref_ice_i, &
854         ref_liq_pi_i, ref_ice_pi_i, &   ! rajoute par OB pour diagnostiquer effet indirect
855         ZEMTD_i  , ZEMTU_i , ZTRSO_i ,&
856         ZTH_i    , ZCTRSO  , ZCEMTR  , ZTRSOD  ,&
857         ZLWFC    , ZLWFT_i , ZSWFC   , ZSWFT_i ,&
858         PSFSWDIR , PSFSWDIF, PFSDNN  , PFSDNV  ,&
859         PPIZA_TOT, PCGA_TOT,PTAU_TOT,&
860         PPIZA_NAT, PCGA_NAT,PTAU_NAT,           &  ! rajoute par OB pour diagnostiquer effet direct
861         PTAU_LW_TOT, PTAU_LW_NAT,               &  ! rajoute par C. Kleinschmitt
862         ZFLUX_i  , ZFLUC_i ,&
863         ZFSDWN_i , ZFSUP_i , ZFCDWN_i, ZFCUP_i, ZFCCDWN_i, ZFCCUP_i, ZFLCCDWN_i, ZFLCCUP_i, &
864         ZTOPSWADAERO,ZSOLSWADAERO,&  ! rajoute par OB pour diagnostics
865         ZTOPSWAD0AERO,ZSOLSWAD0AERO,&
866         ZTOPSWAIAERO,ZSOLSWAIAERO, &
867         ZTOPSWCF_AERO,ZSOLSWCF_AERO, &
868         ZSWADAERO, & !--NL
869         ZTOPLWADAERO,ZSOLLWADAERO,&  ! rajoute par C. Kleinscmitt pour LW diagnostics
870         ZTOPLWAD0AERO,ZSOLLWAD0AERO,&
871         ZTOPLWAIAERO,ZSOLLWAIAERO, &
872         ZLWADAERO, & !--NL
873         ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback) ! flags aerosols
874           
875!        print *,'RADLWSW: apres RECMWF'
876      if(lldebug) then
877        CALL writefield_phy('zemtd_i',ZEMTD_i,klev+1)
878        CALL writefield_phy('zemtu_i',ZEMTU_i,klev+1)
879        CALL writefield_phy('ztrso_i',ZTRSO_i,klev+1)
880        CALL writefield_phy('zth_i',ZTH_i,klev+1)
881        CALL writefield_phy('zctrso',ZCTRSO,2)
882        CALL writefield_phy('zcemtr',ZCEMTR,2)
883        CALL writefield_phy('ztrsod',ZTRSOD,1)
884        CALL writefield_phy('zlwfc',ZLWFC,2)
885        CALL writefield_phy('zlwft_i',ZLWFT_i,klev+1)
886        CALL writefield_phy('zswfc',ZSWFC,2)
887        CALL writefield_phy('zswft_i',ZSWFT_i,klev+1)
888        CALL writefield_phy('psfswdir',PSFSWDIR,6)
889        CALL writefield_phy('psfswdif',PSFSWDIF,6)
890        CALL writefield_phy('pfsdnn',PFSDNN,1)
891        CALL writefield_phy('pfsdnv',PFSDNV,1)
892        CALL writefield_phy('ppiza_dst',PPIZA_TOT,klev)
893        CALL writefield_phy('pcga_dst',PCGA_TOT,klev)
894        CALL writefield_phy('ptaurel_dst',PTAU_TOT,klev)
895        CALL writefield_phy('zflux_i',ZFLUX_i,klev+1)
896        CALL writefield_phy('zfluc_i',ZFLUC_i,klev+1)
897        CALL writefield_phy('zfsdwn_i',ZFSDWN_i,klev+1)
898        CALL writefield_phy('zfsup_i',ZFSUP_i,klev+1)
899        CALL writefield_phy('zfcdwn_i',ZFCDWN_i,klev+1)
900        CALL writefield_phy('zfcup_i',ZFCUP_i,klev+1)
901      endif
902! --------- output RECMWFL
903!  ZEMTD        (KPROMA,KLEV+1)  ; TOTAL DOWNWARD LONGWAVE EMISSIVITY
904!  ZEMTU        (KPROMA,KLEV+1)  ; TOTAL UPWARD   LONGWAVE EMISSIVITY
905!  ZTRSO        (KPROMA,KLEV+1)  ; TOTAL SHORTWAVE TRANSMISSIVITY
906!  ZTH          (KPROMA,KLEV+1)  ; HALF LEVEL TEMPERATURE
907!  ZCTRSO       (KPROMA,2)       ; CLEAR-SKY SHORTWAVE TRANSMISSIVITY
908!  ZCEMTR       (KPROMA,2)       ; CLEAR-SKY NET LONGWAVE EMISSIVITY
909!  ZTRSOD       (KPROMA)         ; TOTAL-SKY SURFACE SW TRANSMISSITY
910!  ZLWFC        (KPROMA,2)       ; CLEAR-SKY LONGWAVE FLUXES
911!  ZLWFT        (KPROMA,KLEV+1)  ; TOTAL-SKY LONGWAVE FLUXES
912!  ZSWFC        (KPROMA,2)       ; CLEAR-SKY SHORTWAVE FLUXES
913!  ZSWFT        (KPROMA,KLEV+1)  ; TOTAL-SKY SHORTWAVE FLUXES
914!  PPIZA_TOT    (KPROMA,KLEV,NSW); Single scattering albedo of total aerosols
915!  PCGA_TOT     (KPROMA,KLEV,NSW); Assymetry factor for total aerosols
916!  PTAU_TOT     (KPROMA,KLEV,NSW); Optical depth of total aerosols
917!  PPIZA_NAT    (KPROMA,KLEV,NSW); Single scattering albedo of natural aerosols
918!  PCGA_NAT     (KPROMA,KLEV,NSW); Assymetry factor for natural aerosols
919!  PTAU_NAT     (KPROMA,KLEV,NSW); Optical depth of natiral aerosols
920!  PTAU_LW_TOT  (KPROMA,KLEV,NLW); LW Optical depth of total aerosols 
921!  PTAU_LW_NAT  (KPROMA,KLEV,NLW); LW Optical depth of natural aerosols 
922!  PSFSWDIR     (KPROMA,NSW)     ;
923!  PSFSWDIF     (KPROMA,NSW)     ;
924!  PFSDNN       (KPROMA)         ;
925!  PFSDNV       (KPROMA)         ;
926! ---------
927! ---------
928! On retablit l'ordre des niveaux lmd pour les tableaux de sortie
929! D autre part, on multiplie les resultats SW par fract pour etre coherent
930! avec l ancien rayonnement AR4. Si nuit, fract=0 donc pas de
931! rayonnement SW. (MPL 260609)
932      DO k=0,klev
933         DO i=1,klon
934         ZEMTD(i,k+1)  = ZEMTD_i(i,k+1)
935         ZEMTU(i,k+1)  = ZEMTU_i(i,k+1)
936         ZTRSO(i,k+1)  = ZTRSO_i(i,k+1)
937         ZTH(i,k+1)    = ZTH_i(i,k+1)
938!        ZLWFT(i,k+1)  = ZLWFT_i(i,klev+1-k)
939!        ZSWFT(i,k+1)  = ZSWFT_i(i,klev+1-k)
940         ZFLUP(i,k+1)  = ZFLUX_i(i,1,k+1)
941         ZFLDN(i,k+1)  = ZFLUX_i(i,2,k+1)
942         ZFLUP0(i,k+1) = ZFLUC_i(i,1,k+1)
943         ZFLDN0(i,k+1) = ZFLUC_i(i,2,k+1)
944         ZFSDN(i,k+1)  = ZFSDWN_i(i,k+1)*fract(i)
945         ZFSDN0(i,k+1) = ZFCDWN_i(i,k+1)*fract(i)
946         ZFSDNC0(i,k+1)= ZFCCDWN_i(i,k+1)*fract(i)
947         ZFSUP (i,k+1) = ZFSUP_i(i,k+1)*fract(i)
948         ZFSUP0(i,k+1) = ZFCUP_i(i,k+1)*fract(i)
949         ZFSUPC0(i,k+1)= ZFCCUP_i(i,k+1)*fract(i)
950         ZFLDNC0(i,k+1)= ZFLCCDWN_i(i,k+1)
951         ZFLUPC0(i,k+1)= ZFLCCUP_i(i,k+1)
952         IF(ok_volcan) THEN
953            ZSWADAERO(i,k+1)=ZSWADAERO(i,k+1)*fract(i) !--NL
954         ENDIF
955         
956!   Nouveau calcul car visiblement ZSWFT et ZSWFC sont nuls dans RRTM cy32
957!   en sortie de radlsw.F90 - MPL 7.01.09
958         ZSWFT(i,k+1)  = (ZFSDWN_i(i,k+1)-ZFSUP_i(i,k+1))*fract(i)
959         ZSWFT0_i(i,k+1) = (ZFCDWN_i(i,k+1)-ZFCUP_i(i,k+1))*fract(i)
960!        WRITE(*,'("FSDN FSUP FCDN FCUP: ",4E12.5)') ZFSDWN_i(i,k+1),&
961!        ZFSUP_i(i,k+1),ZFCDWN_i(i,k+1),ZFCUP_i(i,k+1)
962         ZLWFT(i,k+1) =-ZFLUX_i(i,2,k+1)-ZFLUX_i(i,1,k+1)
963         ZLWFT0_i(i,k+1)=-ZFLUC_i(i,2,k+1)-ZFLUC_i(i,1,k+1)
964!        print *,'FLUX2 FLUX1 FLUC2 FLUC1',ZFLUX_i(i,2,k+1),&
965!    & ZFLUX_i(i,1,k+1),ZFLUC_i(i,2,k+1),ZFLUC_i(i,1,k+1)
966         ENDDO
967      ENDDO
968
969!--ajout OB
970      ZTOPSWADAERO(:) =ZTOPSWADAERO(:) *fract(:)
971      ZSOLSWADAERO(:) =ZSOLSWADAERO(:) *fract(:)
972      ZTOPSWAD0AERO(:)=ZTOPSWAD0AERO(:)*fract(:)
973      ZSOLSWAD0AERO(:)=ZSOLSWAD0AERO(:)*fract(:)
974      ZTOPSWAIAERO(:) =ZTOPSWAIAERO(:) *fract(:)
975      ZSOLSWAIAERO(:) =ZSOLSWAIAERO(:) *fract(:)
976      ZTOPSWCF_AERO(:,1)=ZTOPSWCF_AERO(:,1)*fract(:)
977      ZTOPSWCF_AERO(:,2)=ZTOPSWCF_AERO(:,2)*fract(:)
978      ZTOPSWCF_AERO(:,3)=ZTOPSWCF_AERO(:,3)*fract(:)
979      ZSOLSWCF_AERO(:,1)=ZSOLSWCF_AERO(:,1)*fract(:)
980      ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
981      ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
982
983!     print*,'SW_RRTM ZFSDN0 1 , klev:',ZFSDN0(1:klon,1),ZFSDN0(1:klon,klev)
984!     print*,'SW_RRTM ZFSUP0 1 , klev:',ZFSUP0(1:klon,1),ZFSUP0(1:klon,klev)
985!     print*,'SW_RRTM ZFSDN  1 , klev:',ZFSDN(1:klon,1),ZFSDN(1:klon,klev)
986!     print*,'SW_RRTM ZFSUP  1 , klev:',ZFSUP(1:klon,1),ZFSUP(1:klon,klev)     
987!     print*,'OK1'
988! ---------
989! ---------
990! On renseigne les champs LMDz, pour avoir la meme chose qu'en sortie de
991! LW_LMDAR4 et SW_LMDAR4
992      DO i = 1, kdlon
993         zsolsw(i)    = ZSWFT(i,1)
994         zsolsw0(i)   = ZSWFT0_i(i,1)
995!        zsolsw0(i)   = ZFSDN0(i,1)     -ZFSUP0(i,1)
996         ztopsw(i)    = ZSWFT(i,klev+1)
997         ztopsw0(i)   = ZSWFT0_i(i,klev+1)
998!        ztopsw0(i)   = ZFSDN0(i,klev+1)-ZFSUP0(i,klev+1)
999!         
1000!        zsollw(i)    = ZFLDN(i,1)      -ZFLUP(i,1)
1001!        zsollw0(i)   = ZFLDN0(i,1)     -ZFLUP0(i,1)
1002!        ztoplw(i)    = ZFLDN(i,klev+1) -ZFLUP(i,klev+1)
1003!        ztoplw0(i)   = ZFLDN0(i,klev+1)-ZFLUP0(i,klev+1)
1004         zsollw(i)    = ZLWFT(i,1)
1005         zsollw0(i)   = ZLWFT0_i(i,1)
1006         ztoplw(i)    = ZLWFT(i,klev+1)*(-1)
1007         ztoplw0(i)   = ZLWFT0_i(i,klev+1)*(-1)
1008!         
1009           IF (fract(i) == 0.) THEN
1010!!!!! A REVOIR MPL (20090630) ca n a pas de sens quand fract=0
1011! pas plus que dans le sw_AR4
1012          zalbpla(i)   = 1.0e+39
1013         ELSE
1014          zalbpla(i)   = ZFSUP(i,klev+1)/ZFSDN(i,klev+1)
1015         ENDIF
1016!!! 5 juin 2015
1017!!! Correction MP bug RRTM
1018         zsollwdown(i)= -1.*ZFLDN(i,1)
1019      ENDDO
1020!     print*,'OK2'
1021
1022! extrait de SW_AR4
1023!     DO k = 1, KFLEV
1024!        kpl1 = k+1
1025!        DO i = 1, KDLON
1026!           PHEAT(i,k) = -(ZFSUP(i,kpl1)-ZFSUP(i,k)) -(ZFSDN(i,k)-ZFSDN(i,kpl1))
1027!           PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k)
1028! ZLWFT(klon,k),ZSWFT
1029
1030      do k=1,kflev
1031         do i=1,kdlon
1032           zheat(i,k)=(ZSWFT(i,k+1)-ZSWFT(i,k))*RDAY*RG/RCPD/PDP(i,k)
1033           zheat0(i,k)=(ZSWFT0_i(i,k+1)-ZSWFT0_i(i,k))*RDAY*RG/RCPD/PDP(i,k)
1034           zcool(i,k)=(ZLWFT(i,k)-ZLWFT(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
1035           zcool0(i,k)=(ZLWFT0_i(i,k)-ZLWFT0_i(i,k+1))*RDAY*RG/RCPD/PDP(i,k)
1036           IF(ok_volcan) THEN
1037              zheat_volc(i,k)=(ZSWADAERO(i,k+1)-ZSWADAERO(i,k))*RG/RCPD/PDP(i,k) !NL
1038              zcool_volc(i,k)=(ZLWADAERO(i,k)-ZLWADAERO(i,k+1))*RG/RCPD/PDP(i,k) !NL
1039           ENDIF
1040!          print *,'heat cool heat0 cool0 ',zheat(i,k),zcool(i,k),zheat0(i,k),zcool0(i,k)
1041!          ZFLUCUP_i(i,k)=ZFLUC_i(i,1,k)
1042!          ZFLUCDWN_i(i,k)=ZFLUC_i(i,2,k)         
1043         enddo
1044      enddo
1045#else
1046    abort_message="You should compile with -rrtm if running with iflag_rrtm=1"
1047    call abort_physic(modname, abort_message, 1)
1048#endif
1049    ENDIF ! iflag_rrtm
1050!======================================================================
1051
1052    DO i = 1, kdlon
1053      topsw(iof+i) = ztopsw(i)
1054      toplw(iof+i) = ztoplw(i)
1055      solsw(iof+i) = zsolsw(i)
1056      sollw(iof+i) = zsollw(i)
1057      sollwdown(iof+i) = zsollwdown(i)
1058      DO k = 1, kflev+1
1059        lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
1060        lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
1061        lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
1062        lwup  ( iof+i,k)   = ZFLUP  ( i,k)
1063      ENDDO
1064      topsw0(iof+i) = ztopsw0(i)
1065      toplw0(iof+i) = ztoplw0(i)
1066      solsw0(iof+i) = zsolsw0(i)
1067      sollw0(iof+i) = zsollw0(i)
1068      albpla(iof+i) = zalbpla(i)
1069
1070      DO k = 1, kflev+1
1071        swdnc0( iof+i,k)   = ZFSDNC0( i,k)
1072        swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
1073        swdn  ( iof+i,k)   = ZFSDN  ( i,k)
1074        swupc0( iof+i,k)   = ZFSUPC0( i,k)
1075        swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
1076        swup  ( iof+i,k)   = ZFSUP  ( i,k)
1077        lwdnc0( iof+i,k)   = ZFLDNC0( i,k)
1078        lwupc0( iof+i,k)   = ZFLUPC0( i,k)
1079      ENDDO
1080    ENDDO
1081    !-transform the aerosol forcings, if they have
1082    ! to be calculated
1083    IF (ok_ade) THEN
1084        DO i = 1, kdlon
1085          topswad_aero(iof+i) = ztopswadaero(i)
1086          topswad0_aero(iof+i) = ztopswad0aero(i)
1087          solswad_aero(iof+i) = zsolswadaero(i)
1088          solswad0_aero(iof+i) = zsolswad0aero(i)
1089! MS the following lines seem to be wrong, why is iof on right hand side???
1090!          topsw_aero(iof+i,:) = ztopsw_aero(iof+i,:)
1091!          topsw0_aero(iof+i,:) = ztopsw0_aero(iof+i,:)
1092!          solsw_aero(iof+i,:) = zsolsw_aero(iof+i,:)
1093!          solsw0_aero(iof+i,:) = zsolsw0_aero(iof+i,:)
1094          topsw_aero(iof+i,:) = ztopsw_aero(i,:)
1095          topsw0_aero(iof+i,:) = ztopsw0_aero(i,:)
1096          solsw_aero(iof+i,:) = zsolsw_aero(i,:)
1097          solsw0_aero(iof+i,:) = zsolsw0_aero(i,:)
1098          topswcf_aero(iof+i,:) = ztopswcf_aero(i,:)
1099          solswcf_aero(iof+i,:) = zsolswcf_aero(i,:)   
1100          !-LW
1101          toplwad_aero(iof+i) = ztoplwadaero(i)
1102          toplwad0_aero(iof+i) = ztoplwad0aero(i)
1103          sollwad_aero(iof+i) = zsollwadaero(i)
1104          sollwad0_aero(iof+i) = zsollwad0aero(i)   
1105        ENDDO
1106    ELSE
1107        DO i = 1, kdlon
1108          topswad_aero(iof+i) = 0.0
1109          solswad_aero(iof+i) = 0.0
1110          topswad0_aero(iof+i) = 0.0
1111          solswad0_aero(iof+i) = 0.0
1112          topsw_aero(iof+i,:) = 0.
1113          topsw0_aero(iof+i,:) =0.
1114          solsw_aero(iof+i,:) = 0.
1115          solsw0_aero(iof+i,:) = 0.
1116          !-LW
1117          toplwad_aero(iof+i) = 0.0
1118          sollwad_aero(iof+i) = 0.0
1119          toplwad0_aero(iof+i) = 0.0
1120          sollwad0_aero(iof+i) = 0.0
1121        ENDDO
1122    ENDIF
1123    IF (ok_aie) THEN
1124        DO i = 1, kdlon
1125          topswai_aero(iof+i) = ztopswaiaero(i)
1126          solswai_aero(iof+i) = zsolswaiaero(i)
1127          !-LW
1128          toplwai_aero(iof+i) = ztoplwaiaero(i)
1129          sollwai_aero(iof+i) = zsollwaiaero(i)
1130        ENDDO
1131    ELSE
1132        DO i = 1, kdlon
1133          topswai_aero(iof+i) = 0.0
1134          solswai_aero(iof+i) = 0.0
1135          !-LW
1136          toplwai_aero(iof+i) = 0.0
1137          sollwai_aero(iof+i) = 0.0
1138        ENDDO
1139    ENDIF
1140    DO k = 1, kflev
1141      DO i = 1, kdlon
1142        !        scale factor to take into account the difference between
1143        !        dry air and watter vapour scpecifi! heat capacity
1144        zznormcp=1.0+RVTMP2*PWV(i,k)
1145        heat(iof+i,k) = zheat(i,k)/zznormcp
1146        cool(iof+i,k) = zcool(i,k)/zznormcp
1147        heat0(iof+i,k) = zheat0(i,k)/zznormcp
1148        cool0(iof+i,k) = zcool0(i,k)/zznormcp
1149        IF(ok_volcan) THEN !NL
1150           heat_volc(iof+i,k) = zheat_volc(i,k)/zznormcp
1151           cool_volc(iof+i,k) = zcool_volc(i,k)/zznormcp
1152        ENDIF
1153      ENDDO
1154    ENDDO
1155
1156 ENDDO ! j = 1, nb_gr
1157
1158END SUBROUTINE radlwsw
1159
1160end module radlwsw_m
Note: See TracBrowser for help on using the repository browser.