My Project
 All Classes Files Functions Variables Macros
sw_aeroAR4.F90
Go to the documentation of this file.
1 !
2 ! $Id$
3 !
4 SUBROUTINE sw_aeroar4(PSCT, PRMU0, PFRAC, &
5  ppmb, pdp, &
6  ppsol, palbd, palbp,&
7  ptave, pwv, pqs, pozon, paer,&
8  pcldsw, ptau, pomega, pcg,&
9  pheat, pheat0,&
10  palbpla,ptopsw,psolsw,ptopsw0,psolsw0,&
11  zfsup,zfsdn,zfsup0,zfsdn0,&
12  tauaero, pizaero, cgaero,&
13  ptaua, pomegaa,&
14  ptopswadaero,psolswadaero,&
15  ptopswad0aero,psolswad0aero,&
16  ptopswaiaero,psolswaiaero,&
17  ptopswaero,ptopsw0aero,&
18  psolswaero,psolsw0aero,&
19  ptopswcfaero,psolswcfaero,&
20  ok_ade, ok_aie, flag_aerosol, flag_aerosol_strat )
21 
22  USE dimphy
23  USE phys_output_mod, ONLY : swaero_diag
24  IMPLICIT NONE
25 
26 #include "YOMCST.h"
27 #include "clesphys.h"
28 #include "iniprint.h"
29  !
30  ! ------------------------------------------------------------------
31  !
32  ! PURPOSE.
33  ! --------
34  !
35  ! THIS ROUTINE COMPUTES THE SHORTWAVE RADIATION FLUXES IN TWO
36  ! SPECTRAL INTERVALS FOLLOWING FOUQUART AND BONNEL (1980).
37  !
38  ! METHOD.
39  ! -------
40  !
41  ! 1. COMPUTES ABSORBER AMOUNTS (SWU)
42  ! 2. COMPUTES FLUXES IN 1ST SPECTRAL INTERVAL (SW1S)
43  ! 3. COMPUTES FLUXES IN 2ND SPECTRAL INTERVAL (SW2S)
44  !
45  ! REFERENCE.
46  ! ----------
47  !
48  ! SEE RADIATION'S PART OF THE ECMWF RESEARCH DEPARTMENT
49  ! DOCUMENTATION, AND FOUQUART AND BONNEL (1980)
50  !
51  ! AUTHOR.
52  ! -------
53  ! JEAN-JACQUES MORCRETTE *ECMWF*
54  !
55  ! MODIFICATIONS.
56  ! --------------
57  ! ORIGINAL : 89-07-14
58  ! 1995-01-01 J.-J. MORCRETTE Direct/Diffuse Albedo
59  ! 2003-11-27 J. QUAAS Introduce aerosol forcings (based on BOUCHER)
60  ! 2009-04 A. COZIC - C.DEANDREIS Indroduce NAT/BC/POM/DUST/SS aerosol forcing
61  ! 2012-09 O. BOUCHER - reorganise aerosol cases with ok_ade, ok_aie, flag_aerosol
62  ! ------------------------------------------------------------------
63  !
64  !* ARGUMENTS:
65  !
66  REAL(KIND=8) psct ! constante solaire (valeur conseillee: 1370)
67 
68  REAL(KIND=8) ppsol(kdlon) ! SURFACE PRESSURE (PA)
69  REAL(KIND=8) pdp(kdlon,kflev) ! LAYER THICKNESS (PA)
70  REAL(KIND=8) ppmb(kdlon,kflev+1) ! HALF-LEVEL PRESSURE (MB)
71 
72  REAL(KIND=8) prmu0(kdlon) ! COSINE OF ZENITHAL ANGLE
73  REAL(KIND=8) pfrac(kdlon) ! fraction de la journee
74 
75  REAL(KIND=8) ptave(kdlon,kflev) ! LAYER TEMPERATURE (K)
76  REAL(KIND=8) pwv(kdlon,kflev) ! SPECIFI! HUMIDITY (KG/KG)
77  REAL(KIND=8) pqs(kdlon,kflev) ! SATURATED WATER VAPOUR (KG/KG)
78  REAL(KIND=8) pozon(kdlon,kflev) ! OZONE CONCENTRATION (KG/KG)
79  REAL(KIND=8) paer(kdlon,kflev,5) ! AEROSOLS' OPTICAL THICKNESS
80 
81  REAL(KIND=8) palbd(kdlon,2) ! albedo du sol (lumiere diffuse)
82  REAL(KIND=8) palbp(kdlon,2) ! albedo du sol (lumiere parallele)
83 
84  REAL(KIND=8) pcldsw(kdlon,kflev) ! CLOUD FRACTION
85  REAL(KIND=8) ptau(kdlon,2,kflev) ! CLOUD OPTICAL THICKNESS (pre-industrial value)
86  REAL(KIND=8) pcg(kdlon,2,kflev) ! ASYMETRY FACTOR
87  REAL(KIND=8) pomega(kdlon,2,kflev) ! SINGLE SCATTERING ALBEDO
88 
89  REAL(KIND=8) pheat(kdlon,kflev) ! SHORTWAVE HEATING (K/DAY)
90  REAL(KIND=8) pheat0(kdlon,kflev)! SHORTWAVE HEATING (K/DAY) clear-sky
91  REAL(KIND=8) palbpla(kdlon) ! PLANETARY ALBEDO
92  REAL(KIND=8) ptopsw(kdlon) ! SHORTWAVE FLUX AT T.O.A.
93  REAL(KIND=8) psolsw(kdlon) ! SHORTWAVE FLUX AT SURFACE
94  REAL(KIND=8) ptopsw0(kdlon) ! SHORTWAVE FLUX AT T.O.A. (CLEAR-SKY)
95  REAL(KIND=8) psolsw0(kdlon) ! SHORTWAVE FLUX AT SURFACE (CLEAR-SKY)
96  !
97  !* LOCAL VARIABLES:
98  !
99  real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
100 
101  REAL(KIND=8) zoz(kdlon,kflev)
102  ! column-density of ozone in layer, in kilo-Dobsons
103 
104  REAL(KIND=8) zaki(kdlon,2)
105  REAL(KIND=8) zcld(kdlon,kflev)
106  REAL(KIND=8) zclear(kdlon)
107  REAL(KIND=8) zdsig(kdlon,kflev)
108  REAL(KIND=8) zfact(kdlon)
109  REAL(KIND=8) zfd(kdlon,kflev+1)
110  REAL(KIND=8) zfdown(kdlon,kflev+1)
111  REAL(KIND=8) zfu(kdlon,kflev+1)
112  REAL(KIND=8) zfup(kdlon,kflev+1)
113  REAL(KIND=8) zrmu(kdlon)
114  REAL(KIND=8) zsec(kdlon)
115  REAL(KIND=8) zud(kdlon,5,kflev+1)
116  REAL(KIND=8) zcldsw0(kdlon,kflev)
117 
118  REAL(KIND=8) zfsup(kdlon,kflev+1)
119  REAL(KIND=8) zfsdn(kdlon,kflev+1)
120  REAL(KIND=8) zfsup0(kdlon,kflev+1)
121  REAL(KIND=8) zfsdn0(kdlon,kflev+1)
122 
123  INTEGER inu, jl, jk, i, k, kpl1
124 
125  INTEGER swpas ! Every swpas steps, sw is calculated
126  parameter(swpas=1)
127 
128  INTEGER, SAVE :: itapsw = 0
129  !$OMP THREADPRIVATE(itapsw)
130  LOGICAL, SAVE :: appel1er = .true.
131  !$OMP THREADPRIVATE(appel1er)
132  LOGICAL, SAVE :: initialized = .false.
133  !$OMP THREADPRIVATE(initialized)
134 
135  !jq-local flag introduced for aerosol forcings
136  REAL(KIND=8), SAVE :: flag_aer
137  !$OMP THREADPRIVATE(flag_aer)
138 
139  LOGICAL ok_ade, ok_aie ! use aerosol forcings or not?
140  LOGICAL flag_aerosol_strat ! use stratospehric aerosols
141  INTEGER flag_aerosol ! global flag for aerosol 0 (no aerosol) or 1-5 (aerosols)
142  REAL(KIND=8) tauaero(kdlon,kflev,9,2) ! aerosol optical properties
143  REAL(KIND=8) pizaero(kdlon,kflev,9,2) ! (see aeropt.F)
144  REAL(KIND=8) cgaero(kdlon,kflev,9,2) ! -"-
145  REAL(KIND=8) ptaua(kdlon,2,kflev) ! CLOUD OPTICAL THICKNESS (present-day value)
146  REAL(KIND=8) pomegaa(kdlon,2,kflev) ! SINGLE SCATTERING ALBEDO
147  REAL(KIND=8) ptopswadaero(kdlon) ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
148  REAL(KIND=8) psolswadaero(kdlon) ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
149  REAL(KIND=8) ptopswad0aero(kdlon) ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
150  REAL(KIND=8) psolswad0aero(kdlon) ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
151  REAL(KIND=8) ptopswaiaero(kdlon) ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL IND)
152  REAL(KIND=8) psolswaiaero(kdlon) ! SHORTWAVE FLUX AT SURFACE(+AEROSOL IND)
153  REAL(KIND=8) ptopswaero(kdlon,9) ! SW TOA AS DRF nat & ant
154  REAL(KIND=8) ptopsw0aero(kdlon,9) ! SW SRF AS DRF nat & ant
155  REAL(KIND=8) psolswaero(kdlon,9) ! SW TOA CS DRF nat & ant
156  REAL(KIND=8) psolsw0aero(kdlon,9) ! SW SRF CS DRF nat & ant
157  REAL(KIND=8) ptopswcfaero(kdlon,3) ! SW TOA AS cloudRF nat & ant
158  REAL(KIND=8) psolswcfaero(kdlon,3) ! SW SRF AS cloudRF nat & ant
159 
160  !jq - Fluxes including aerosol effects
161  REAL(KIND=8),ALLOCATABLE,SAVE :: zfsupad_aero(:,:)
162  !$OMP THREADPRIVATE(ZFSUPAD_AERO)
163  REAL(KIND=8),ALLOCATABLE,SAVE :: zfsdnad_aero(:,:)
164  !$OMP THREADPRIVATE(ZFSDNAD_AERO)
165  !jq - Fluxes including aerosol effects
166  REAL(KIND=8),ALLOCATABLE,SAVE :: zfsupad0_aero(:,:)
167  !$OMP THREADPRIVATE(ZFSUPAD0_AERO)
168  REAL(KIND=8),ALLOCATABLE,SAVE :: zfsdnad0_aero(:,:)
169  !$OMP THREADPRIVATE(ZFSDNAD0_AERO)
170  REAL(KIND=8),ALLOCATABLE,SAVE :: zfsupai_aero(:,:)
171  !$OMP THREADPRIVATE(ZFSUPAI_AERO)
172  REAL(KIND=8),ALLOCATABLE,SAVE :: zfsdnai_aero(:,:)
173  !$OMP THREADPRIVATE(ZFSDNAI_AERO)
174  REAL(KIND=8),ALLOCATABLE,SAVE :: zfsup_aero(:,:,:)
175  !$OMP THREADPRIVATE(ZFSUP_AERO)
176  REAL(KIND=8),ALLOCATABLE,SAVE :: zfsdn_aero(:,:,:)
177  !$OMP THREADPRIVATE(ZFSDN_AERO)
178  REAL(KIND=8),ALLOCATABLE,SAVE :: zfsup0_aero(:,:,:)
179  !$OMP THREADPRIVATE(ZFSUP0_AERO)
180  REAL(KIND=8),ALLOCATABLE,SAVE :: zfsdn0_aero(:,:,:)
181  !$OMP THREADPRIVATE(ZFSDN0_AERO)
182 
183 ! Key to define the aerosol effect acting on climate
184 ! OB: AEROSOLFEEDBACK_ACTIVE is now a LOGICAL
185 ! TRUE: fluxes use natural and/or anthropogenic aerosols according to ok_ade and ok_aie, DEFAULT
186 ! FALSE: fluxes use no aerosols (case 1)
187 
188  LOGICAL,SAVE :: aerosolfeedback_active = .true.
189 !$OMP THREADPRIVATE(AEROSOLFEEDBACK_ACTIVE)
190 
191  CHARACTER (LEN=20) :: modname='sw_aeroAR4'
192  CHARACTER (LEN=80) :: abort_message
193 
194  IF(.NOT.initialized) THEN
195  flag_aer=0.
196  initialized=.true.
197  ALLOCATE(zfsupad_aero(kdlon,kflev+1))
198  ALLOCATE(zfsdnad_aero(kdlon,kflev+1))
199  ALLOCATE(zfsupad0_aero(kdlon,kflev+1))
200  ALLOCATE(zfsdnad0_aero(kdlon,kflev+1))
201  ALLOCATE(zfsupai_aero(kdlon,kflev+1))
202  ALLOCATE(zfsdnai_aero(kdlon,kflev+1))
203 !-OB decrease size of these arrays to what is needed
204 ! | direct effect
205 !ind effect | no aerosol natural total
206 !natural (PTAU) | 1 3 2 --ZFSUP/ZFSDN
207 !total (PTAUA) | 5 4 --ZFSUP/ZFSDN
208 !no cloud | 1 3 2 --ZFSUP0/ZFSDN0
209 ! so we need which case when ?
210 ! ok_ade and ok_aie = 4-5, 4-2 and 2
211 ! ok_ade and not ok_aie = 2-3 and 2
212 ! not ok_ade and ok_aie = 5-3 and 5
213 ! not ok_ade and not ok_aie = 3
214 ! therefore the cases have the folliwng switches
215 ! 3 = not ok_ade or not ok_aie
216 ! 4 = ok_ade and ok_aie
217 ! 2 = ok_ade
218 ! 5 = ok_aie
219  ALLOCATE(zfsup_aero(kdlon,kflev+1,5))
220  ALLOCATE(zfsdn_aero(kdlon,kflev+1,5))
221  ALLOCATE(zfsup0_aero(kdlon,kflev+1,3))
222  ALLOCATE(zfsdn0_aero(kdlon,kflev+1,3))
223 ! end OB modif
224  zfsupad_aero(:,:)=0.
225  zfsdnad_aero(:,:)=0.
226  zfsupad0_aero(:,:)=0.
227  zfsdnad0_aero(:,:)=0.
228  zfsupai_aero(:,:)=0.
229  zfsdnai_aero(:,:)=0.
230  zfsup_aero(:,:,:)=0.
231  zfsdn_aero(:,:,:)=0.
232  zfsup0_aero(:,:,:)=0.
233  zfsdn0_aero(:,:,:)=0.
234  ENDIF
235 
236  IF (appel1er) THEN
237  WRITE(lunout,*)'SW calling frequency : ', swpas
238  WRITE(lunout,*) " In general, it should be 1"
239  appel1er = .false.
240  ENDIF
241  ! ------------------------------------------------------------------
242  IF (mod(itapsw,swpas).EQ.0) THEN
243 
244  DO jk = 1 , kflev
245  DO jl = 1, kdlon
246  zcldsw0(jl,jk) = 0.0
247  zoz(jl,jk) = pozon(jl,jk)*46.6968/rg &
248  *pdp(jl,jk)*(101325.0/ppsol(jl))
249  ENDDO
250  ENDDO
251 
252 ! clear sky with no aerosols at all is computed IF ACTIVEFEEDBACK_ACTIVE is false or for extended diag
253  IF ( swaero_diag .or. .not. aerosolfeedback_active .OR. flag_aerosol .EQ. 0 ) THEN
254 
255  ! clear-sky: zero aerosol effect
256  flag_aer=0.0
257  CALL swu_lmdar4(psct,zcldsw0,ppmb,ppsol,&
258  prmu0,pfrac,ptave,pwv,&
259  zaki,zcld,zclear,zdsig,zfact,zrmu,zsec,zud)
260  inu = 1
261  CALL sw1s_lmdar4(inu,paer, flag_aer, &
262  tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
263  palbd, palbp, pcg, zcld, zclear, zcldsw0,&
264  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
265  zfd, zfu)
266  inu = 2
267  CALL sw2s_lmdar4(inu, paer, flag_aer, &
268  tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
269  zaki, palbd, palbp, pcg, zcld, zclear, zcldsw0,&
270  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
271  pwv, pqs,&
272  zfdown, zfup)
273  DO jk = 1 , kflev+1
274  DO jl = 1, kdlon
275  zfsup0_aero(jl,jk,1) = (zfup(jl,jk) + zfu(jl,jk)) * zfact(jl)
276  zfsdn0_aero(jl,jk,1) = (zfdown(jl,jk) + zfd(jl,jk)) * zfact(jl)
277  ENDDO
278  ENDDO
279  ENDIF ! swaero_diag .or. .not. AEROSOLFEEDBACK_ACTIVE
280 
281 ! cloudy sky with no aerosols at all is either computed IF no indirect effect is asked for, or for extended diag
282  IF ( swaero_diag .or. .not. aerosolfeedback_active .OR. flag_aerosol .EQ. 0 ) THEN
283  ! cloudy-sky: zero aerosol effect
284  flag_aer=0.0
285  CALL swu_lmdar4(psct,pcldsw,ppmb,ppsol,&
286  prmu0,pfrac,ptave,pwv,&
287  zaki,zcld,zclear,zdsig,zfact,zrmu,zsec,zud)
288  inu = 1
289  CALL sw1s_lmdar4(inu, paer, flag_aer, &
290  tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
291  palbd, palbp, pcg, zcld, zclear, pcldsw,&
292  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
293  zfd, zfu)
294  inu = 2
295  CALL sw2s_lmdar4(inu, paer, flag_aer, &
296  tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
297  zaki, palbd, palbp, pcg, zcld, zclear, pcldsw,&
298  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
299  pwv, pqs,&
300  zfdown, zfup)
301 
302  DO jk = 1 , kflev+1
303  DO jl = 1, kdlon
304  zfsup_aero(jl,jk,1) = (zfup(jl,jk) + zfu(jl,jk)) * zfact(jl)
305  zfsdn_aero(jl,jk,1) = (zfdown(jl,jk) + zfd(jl,jk)) * zfact(jl)
306  ENDDO
307  ENDDO
308  ENDIF ! swaero_diag .or. .not. AEROSOLFEEDBACK_ACTIVE
309 
310  IF (flag_aerosol .GT. 0 .OR. flag_aerosol_strat) THEN
311 
312  IF (ok_ade.and.swaero_diag .or. .not. ok_ade) THEN
313 
314  ! clear sky direct effect natural aerosol
315  ! CAS AER (3)
316  flag_aer=1.0
317  CALL swu_lmdar4(psct,zcldsw0,ppmb,ppsol,&
318  prmu0,pfrac,ptave,pwv,&
319  zaki,zcld,zclear,zdsig,zfact,zrmu,zsec,zud)
320  inu = 1
321  CALL sw1s_lmdar4(inu, paer, flag_aer,&
322  tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
323  palbd, palbp, pcg, zcld, zclear, pcldsw,&
324  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
325  zfd, zfu)
326  inu = 2
327  CALL sw2s_lmdar4(inu, paer, flag_aer,&
328  tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
329  zaki, palbd, palbp, pcg, zcld, zclear, pcldsw,&
330  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
331  pwv, pqs,&
332  zfdown, zfup)
333 
334  DO jk = 1 , kflev+1
335  DO jl = 1, kdlon
336  zfsup0_aero(jl,jk,3) = (zfup(jl,jk) + zfu(jl,jk)) * zfact(jl)
337  zfsdn0_aero(jl,jk,3) = (zfdown(jl,jk) + zfd(jl,jk)) * zfact(jl)
338  ENDDO
339  ENDDO
340  ENDIF !--end not swaero_diag or not ok_ade
341 
342  IF (ok_ade) THEN
343 
344  ! clear sky direct effect of total aerosol
345  ! CAS AER (2)
346  flag_aer=1.0
347  CALL swu_lmdar4(psct,zcldsw0,ppmb,ppsol,&
348  prmu0,pfrac,ptave,pwv,&
349  zaki,zcld,zclear,zdsig,zfact,zrmu,zsec,zud)
350  inu = 1
351  CALL sw1s_lmdar4(inu, paer, flag_aer,&
352  tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
353  palbd, palbp, pcg, zcld, zclear, pcldsw,&
354  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
355  zfd, zfu)
356  inu = 2
357  CALL sw2s_lmdar4(inu, paer, flag_aer,&
358  tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
359  zaki, palbd, palbp, pcg, zcld, zclear, pcldsw,&
360  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
361  pwv, pqs,&
362  zfdown, zfup)
363 
364  DO jk = 1 , kflev+1
365  DO jl = 1, kdlon
366  zfsup0_aero(jl,jk,2) = (zfup(jl,jk) + zfu(jl,jk)) * zfact(jl)
367  zfsdn0_aero(jl,jk,2) = (zfdown(jl,jk) + zfd(jl,jk)) * zfact(jl)
368  ENDDO
369  ENDDO
370 
371  ! cloudy-sky with natural aerosols for indirect effect
372  ! but total aerosols for direct effect
373  ! PTAU
374  ! CAS AER (2)
375  flag_aer=1.0
376  CALL swu_lmdar4(psct,pcldsw,ppmb,ppsol,&
377  prmu0,pfrac,ptave,pwv,&
378  zaki,zcld,zclear,zdsig,zfact,zrmu,zsec,zud)
379  inu = 1
380  CALL sw1s_lmdar4(inu, paer, flag_aer,&
381  tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
382  palbd, palbp, pcg, zcld, zclear, pcldsw,&
383  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
384  zfd, zfu)
385  inu = 2
386  CALL sw2s_lmdar4(inu, paer, flag_aer,&
387  tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
388  zaki, palbd, palbp, pcg, zcld, zclear, pcldsw,&
389  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
390  pwv, pqs,&
391  zfdown, zfup)
392 
393  DO jk = 1 , kflev+1
394  DO jl = 1, kdlon
395  zfsup_aero(jl,jk,2) = (zfup(jl,jk) + zfu(jl,jk)) * zfact(jl)
396  zfsdn_aero(jl,jk,2) = (zfdown(jl,jk) + zfd(jl,jk)) * zfact(jl)
397  ENDDO
398  ENDDO
399 
400  ENDIF !-end ok_ade
401 
402  IF ( .not. ok_ade .or. .not. ok_aie ) THEN
403 
404  ! cloudy-sky with natural aerosols for indirect effect
405  ! and natural aerosols for direct effect
406  ! PTAU
407  ! CAS AER (3)
408  ! cloudy-sky direct effect natural aerosol
409  flag_aer=1.0
410  CALL swu_lmdar4(psct,pcldsw,ppmb,ppsol,&
411  prmu0,pfrac,ptave,pwv,&
412  zaki,zcld,zclear,zdsig,zfact,zrmu,zsec,zud)
413  inu = 1
414  CALL sw1s_lmdar4(inu, paer, flag_aer,&
415  tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
416  palbd, palbp, pcg, zcld, zclear, pcldsw,&
417  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
418  zfd, zfu)
419  inu = 2
420  CALL sw2s_lmdar4(inu, paer, flag_aer,&
421  tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
422  zaki, palbd, palbp, pcg, zcld, zclear, pcldsw,&
423  zdsig, pomega, zoz, zrmu, zsec, ptau, zud,&
424  pwv, pqs,&
425  zfdown, zfup)
426 
427  DO jk = 1 , kflev+1
428  DO jl = 1, kdlon
429  zfsup_aero(jl,jk,3) = (zfup(jl,jk) + zfu(jl,jk)) * zfact(jl)
430  zfsdn_aero(jl,jk,3) = (zfdown(jl,jk) + zfd(jl,jk)) * zfact(jl)
431  ENDDO
432  ENDDO
433 
434  ENDIF !--true/false or false/true
435 
436  IF (ok_ade .and. ok_aie) THEN
437 
438  ! cloudy-sky with total aerosols for indirect effect
439  ! and total aerosols for direct effect
440  ! PTAUA
441  ! CAS AER (2)
442  flag_aer=1.0
443  CALL swu_lmdar4(psct,pcldsw,ppmb,ppsol,&
444  prmu0,pfrac,ptave,pwv,&
445  zaki,zcld,zclear,zdsig,zfact,zrmu,zsec,zud)
446  inu = 1
447  CALL sw1s_lmdar4(inu, paer, flag_aer,&
448  tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
449  palbd, palbp, pcg, zcld, zclear, pcldsw,&
450  zdsig, pomegaa, zoz, zrmu, zsec, ptaua, zud,&
451  zfd, zfu)
452  inu = 2
453  CALL sw2s_lmdar4(inu, paer, flag_aer,&
454  tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
455  zaki, palbd, palbp, pcg, zcld, zclear, pcldsw,&
456  zdsig, pomegaa, zoz, zrmu, zsec, ptaua, zud,&
457  pwv, pqs,&
458  zfdown, zfup)
459 
460  DO jk = 1 , kflev+1
461  DO jl = 1, kdlon
462  zfsup_aero(jl,jk,4) = (zfup(jl,jk) + zfu(jl,jk)) * zfact(jl)
463  zfsdn_aero(jl,jk,4) = (zfdown(jl,jk) + zfd(jl,jk)) * zfact(jl)
464  ENDDO
465  ENDDO
466 
467  ENDIF ! ok_ade .and. ok_aie
468 
469  IF (ok_aie) THEN
470  ! cloudy-sky with total aerosols for indirect effect
471  ! and natural aerosols for direct effect
472  ! PTAUA
473  ! CAS AER (3)
474  flag_aer=1.0
475  CALL swu_lmdar4(psct,pcldsw,ppmb,ppsol,&
476  prmu0,pfrac,ptave,pwv,&
477  zaki,zcld,zclear,zdsig,zfact,zrmu,zsec,zud)
478  inu = 1
479  CALL sw1s_lmdar4(inu, paer, flag_aer,&
480  tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
481  palbd, palbp, pcg, zcld, zclear, pcldsw,&
482  zdsig, pomegaa, zoz, zrmu, zsec, ptaua, zud,&
483  zfd, zfu)
484  inu = 2
485  CALL sw2s_lmdar4(inu, paer, flag_aer,&
486  tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
487  zaki, palbd, palbp, pcg, zcld, zclear, pcldsw,&
488  zdsig, pomegaa, zoz, zrmu, zsec, ptaua, zud,&
489  pwv, pqs,&
490  zfdown, zfup)
491 
492  DO jk = 1 , kflev+1
493  DO jl = 1, kdlon
494  zfsup_aero(jl,jk,5) = (zfup(jl,jk) + zfu(jl,jk)) * zfact(jl)
495  zfsdn_aero(jl,jk,5) = (zfdown(jl,jk) + zfd(jl,jk)) * zfact(jl)
496  ENDDO
497  ENDDO
498 
499  ENDIF ! ok_aie
500 
501  ENDIF !--if flag_aerosol GT 0 OR flag_aerosol_strat
502 
503  itapsw = 0
504  ENDIF
505  itapsw = itapsw + 1
506 
507  IF ( aerosolfeedback_active .AND. (flag_aerosol .GT. 0 .OR. flag_aerosol_strat) ) THEN
508  IF ( ok_ade .and. ok_aie ) THEN
509  zfsup(:,:) = zfsup_aero(:,:,4)
510  zfsdn(:,:) = zfsdn_aero(:,:,4)
511  zfsup0(:,:) = zfsup0_aero(:,:,2)
512  zfsdn0(:,:) = zfsdn0_aero(:,:,2)
513  ENDIF
514 
515  IF ( ok_ade .and. (.not. ok_aie) ) THEN
516  zfsup(:,:) = zfsup_aero(:,:,2)
517  zfsdn(:,:) = zfsdn_aero(:,:,2)
518  zfsup0(:,:) = zfsup0_aero(:,:,2)
519  zfsdn0(:,:) = zfsdn0_aero(:,:,2)
520  ENDIF
521 
522  IF ( (.not. ok_ade) .and. ok_aie ) THEN
523  zfsup(:,:) = zfsup_aero(:,:,5)
524  zfsdn(:,:) = zfsdn_aero(:,:,5)
525  zfsup0(:,:) = zfsup0_aero(:,:,3)
526  zfsdn0(:,:) = zfsdn0_aero(:,:,3)
527  ENDIF
528 
529  IF ((.not. ok_ade) .and. (.not. ok_aie)) THEN
530  zfsup(:,:) = zfsup_aero(:,:,3)
531  zfsdn(:,:) = zfsdn_aero(:,:,3)
532  zfsup0(:,:) = zfsup0_aero(:,:,3)
533  zfsdn0(:,:) = zfsdn0_aero(:,:,3)
534  ENDIF
535 
536 ! MS the following allows to compute the forcing diagostics without
537 ! letting the aerosol forcing act on the meteorology
538 ! SEE logic above
539  ELSE
540  zfsup(:,:) = zfsup_aero(:,:,1)
541  zfsdn(:,:) = zfsdn_aero(:,:,1)
542  zfsup0(:,:) = zfsup0_aero(:,:,1)
543  zfsdn0(:,:) = zfsdn0_aero(:,:,1)
544  ENDIF
545 
546 ! Now computes heating rates
547  DO k = 1, kflev
548  kpl1 = k+1
549  DO i = 1, kdlon
550  pheat(i,k) = -(zfsup(i,kpl1)-zfsup(i,k))-(zfsdn(i,k)-zfsdn(i,kpl1))
551  pheat(i,k) = pheat(i,k) * rday*rg/rcpd / pdp(i,k)
552  pheat0(i,k) = -(zfsup0(i,kpl1)-zfsup0(i,k))-(zfsdn0(i,k)-zfsdn0(i,kpl1))
553  pheat0(i,k) = pheat0(i,k) * rday*rg/rcpd / pdp(i,k)
554  ENDDO
555  ENDDO
556 
557  DO i = 1, kdlon
558 ! effective SW surface albedo calculation
559  palbpla(i) = zfsup(i,kflev+1)/(zfsdn(i,kflev+1)+1.0e-20)
560 
561 ! clear sky net fluxes at TOA and SRF
562  psolsw0(i) = zfsdn0(i,1) - zfsup0(i,1)
563  ptopsw0(i) = zfsdn0(i,kflev+1) - zfsup0(i,kflev+1)
564 
565 ! cloudy sky net fluxes at TOA and SRF
566  psolsw(i) = zfsdn(i,1) - zfsup(i,1)
567  ptopsw(i) = zfsdn(i,kflev+1) - zfsup(i,kflev+1)
568 
569 ! net anthropogenic forcing direct and 1st indirect effect diagnostics
570 ! requires a natural aerosol field read and used
571 ! Difference of net fluxes from double call to radiation
572 
573 IF (ok_ade) THEN
574 
575 ! indices 1: natural; 2 anthropogenic
576 
577 ! TOA/SRF all sky natural forcing
578  psolswaero(i,1) = (zfsdn_aero(i,1,3) - zfsup_aero(i,1,3))-(zfsdn_aero(i,1,1) - zfsup_aero(i,1,1))
579  ptopswaero(i,1) = (zfsdn_aero(i,kflev+1,3) - zfsup_aero(i,kflev+1,3))- (zfsdn_aero(i,kflev+1,1) - zfsup_aero(i,kflev+1,1))
580 
581 ! TOA/SRF clear sky natural forcing
582  psolsw0aero(i,1) = (zfsdn0_aero(i,1,3) - zfsup0_aero(i,1,3))-(zfsdn0_aero(i,1,1) - zfsup0_aero(i,1,1))
583  ptopsw0aero(i,1) = (zfsdn0_aero(i,kflev+1,3) - zfsup0_aero(i,kflev+1,3))-(zfsdn0_aero(i,kflev+1,1) - zfsup0_aero(i,kflev+1,1))
584 
585  IF (ok_aie) THEN
586 
587 ! TOA/SRF all sky anthropogenic forcing
588  psolswaero(i,2) = (zfsdn_aero(i,1,4) - zfsup_aero(i,1,4))-(zfsdn_aero(i,1,5) - zfsup_aero(i,1,5))
589  ptopswaero(i,2) = (zfsdn_aero(i,kflev+1,4) - zfsup_aero(i,kflev+1,4))- (zfsdn_aero(i,kflev+1,5) - zfsup_aero(i,kflev+1,5))
590 
591  ELSE
592 
593 ! TOA/SRF all sky anthropogenic forcing
594  psolswaero(i,2) = (zfsdn_aero(i,1,2) - zfsup_aero(i,1,2))-(zfsdn_aero(i,1,3) - zfsup_aero(i,1,3))
595  ptopswaero(i,2) = (zfsdn_aero(i,kflev+1,2) - zfsup_aero(i,kflev+1,2))- (zfsdn_aero(i,kflev+1,3) - zfsup_aero(i,kflev+1,3))
596 
597  ENDIF
598 
599 ! TOA/SRF clear sky anthropogenic forcing
600  psolsw0aero(i,2) = (zfsdn0_aero(i,1,2) - zfsup0_aero(i,1,2))-(zfsdn0_aero(i,1,3) - zfsup0_aero(i,1,3))
601  ptopsw0aero(i,2) = (zfsdn0_aero(i,kflev+1,2) - zfsup0_aero(i,kflev+1,2))-(zfsdn0_aero(i,kflev+1,3) - zfsup0_aero(i,kflev+1,3))
602 
603 ! direct anthropogenic forcing , as in old LMDzT, however differences of net fluxes
604  psolswadaero(i) = psolswaero(i,2)
605  ptopswadaero(i) = ptopswaero(i,2)
606  psolswad0aero(i) = psolsw0aero(i,2)
607  ptopswad0aero(i) = ptopsw0aero(i,2)
608 
609 ! OB: these diagnostics may not always work but who need them
610 ! Cloud forcing indices 1: natural; 2 anthropogenic; 3: zero aerosol direct effect
611 ! Instantaneously computed cloudy sky direct aerosol effect, cloud forcing due to aerosols above clouds
612 ! natural
613  psolswcfaero(i,1) = psolswaero(i,1) - psolsw0aero(i,1)
614  ptopswcfaero(i,1) = ptopswaero(i,1) - ptopsw0aero(i,1)
615 
616 ! Instantaneously computed cloudy SKY DIRECT aerosol effect, cloud forcing due to aerosols above clouds
617 ! anthropogenic
618  psolswcfaero(i,2) = psolswaero(i,2) - psolsw0aero(i,2)
619  ptopswcfaero(i,2) = ptopswaero(i,2) - ptopsw0aero(i,2)
620 
621 ! Cloudforcing without aerosol
622 ! zero
623  psolswcfaero(i,3) = (zfsdn_aero(i,1,1) - zfsup_aero(i,1,1))-(zfsdn0_aero(i,1,1) - zfsup0_aero(i,1,1))
624  ptopswcfaero(i,3) = (zfsdn_aero(i,kflev+1,1) - zfsup_aero(i,kflev+1,1))- (zfsdn0_aero(i,kflev+1,1) - zfsup0_aero(i,kflev+1,1))
625 
626 ENDIF
627 
628 IF (ok_aie) THEN
629  IF (ok_ade) THEN
630  psolswaiaero(i) = (zfsdn_aero(i,1,4) - zfsup_aero(i,1,4))-(zfsdn_aero(i,1,2) - zfsup_aero(i,1,2))
631  ptopswaiaero(i) = (zfsdn_aero(i,kflev+1,4) - zfsup_aero(i,kflev+1,4))-(zfsdn_aero(i,kflev+1,2) - zfsup_aero(i,kflev+1,2))
632  ELSE
633  psolswaiaero(i) = (zfsdn_aero(i,1,5) - zfsup_aero(i,1,5))-(zfsdn_aero(i,1,3) - zfsup_aero(i,1,3))
634  ptopswaiaero(i) = (zfsdn_aero(i,kflev+1,5) - zfsup_aero(i,kflev+1,5))-(zfsdn_aero(i,kflev+1,3) - zfsup_aero(i,kflev+1,3))
635  ENDIF
636 ENDIF
637 
638 ENDDO
639 
640 END SUBROUTINE sw_aeroar4