source: LMDZ5/branches/LMDZ6_rc0/libf/phylmd/sw_aeroAR4.F90

Last change on this file was 2575, checked in by acozic, 8 years ago

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