source: LMDZ6/branches/Amaury_dev/libf/phylmd/sw_aeroAR4.F90 @ 5218

Last change on this file since 5218 was 5160, checked in by abarral, 4 months ago

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