source: LMDZ6/branches/DYNAMICO-conv-GC/libf/phylmd/radlwsw_m.F90 @ 3603

Last change on this file since 3603 was 3406, checked in by jghattas, 6 years ago

Added all modifications in the model code that were used for the simulations with DYANMICO during the Grand Challeng 2018. Modifications done by Y. Meurdesoif, L. Fairhead and A.K. Traore

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