source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/radlwsw_m.F90 @ 5441

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