source: LMDZ5/trunk/libf/phylmd/sw_aeroAR4.F90 @ 2251

Last change on this file since 2251 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.
  • 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
File size: 23.8 KB
Line 
1!
2! $Id$
3!
4SUBROUTINE 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
573IF (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
626ENDIF
627
628IF (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
636ENDIF
637
638ENDDO
639
640END SUBROUTINE SW_AEROAR4
Note: See TracBrowser for help on using the repository browser.