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

Last change on this file since 1584 was 1575, checked in by jghattas, 13 years ago
  • Added suffix _mpi_rank to lmdz.out text file. Each processus now write into seperate file but only if lunout/=6 (lunout is set in run.def).
  • Change some print* into write(lunout,*)
  • phytrac.F90 : always include ini_histrac and write_histrac. The file histrac.nc is written if ecrit_tra> 0 (set in physiq.def). Change default value of ecrit_tra into 0.
File size: 21.4 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 )
21
22  USE dimphy
23
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  !        95-01-01   J.-J. MORCRETTE  Direct/Diffuse Albedo
59  !        03-11-27   J. QUAAS Introduce aerosol forcings (based on BOUCHER)
60  !        09-04      A. COZIC - C.DEANDREIS Indroduce NAT/BC/POM/DUST/SS aerosol forcing
61  !     ------------------------------------------------------------------
62  !
63  !* ARGUMENTS:
64  !
65  REAL(KIND=8) PSCT  ! constante solaire (valeur conseillee: 1370)
66
67  REAL(KIND=8) PPSOL(KDLON)        ! SURFACE PRESSURE (PA)
68  REAL(KIND=8) PDP(KDLON,KFLEV)    ! LAYER THICKNESS (PA)
69  REAL(KIND=8) PPMB(KDLON,KFLEV+1) ! HALF-LEVEL PRESSURE (MB)
70
71  REAL(KIND=8) PRMU0(KDLON)  ! COSINE OF ZENITHAL ANGLE
72  REAL(KIND=8) PFRAC(KDLON)  ! fraction de la journee
73
74  REAL(KIND=8) PTAVE(KDLON,KFLEV)  ! LAYER TEMPERATURE (K)
75  REAL(KIND=8) PWV(KDLON,KFLEV)    ! SPECIFI! HUMIDITY (KG/KG)
76  REAL(KIND=8) PQS(KDLON,KFLEV)    ! SATURATED WATER VAPOUR (KG/KG)
77  REAL(KIND=8) POZON(KDLON,KFLEV)  ! OZONE CONCENTRATION (KG/KG)
78  REAL(KIND=8) PAER(KDLON,KFLEV,5) ! AEROSOLS' OPTICAL THICKNESS
79
80  REAL(KIND=8) PALBD(KDLON,2)  ! albedo du sol (lumiere diffuse)
81  REAL(KIND=8) PALBP(KDLON,2)  ! albedo du sol (lumiere parallele)
82
83  REAL(KIND=8) PCLDSW(KDLON,KFLEV)    ! CLOUD FRACTION
84  REAL(KIND=8) PTAU(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS
85  REAL(KIND=8) PCG(KDLON,2,KFLEV)     ! ASYMETRY FACTOR
86  REAL(KIND=8) POMEGA(KDLON,2,KFLEV)  ! SINGLE SCATTERING ALBEDO
87
88  REAL(KIND=8) PHEAT(KDLON,KFLEV) ! SHORTWAVE HEATING (K/DAY)
89  REAL(KIND=8) PHEAT0(KDLON,KFLEV)! SHORTWAVE HEATING (K/DAY) clear-sky
90  REAL(KIND=8) PALBPLA(KDLON)     ! PLANETARY ALBEDO
91  REAL(KIND=8) PTOPSW(KDLON)      ! SHORTWAVE FLUX AT T.O.A.
92  REAL(KIND=8) PSOLSW(KDLON)      ! SHORTWAVE FLUX AT SURFACE
93  REAL(KIND=8) PTOPSW0(KDLON)     ! SHORTWAVE FLUX AT T.O.A. (CLEAR-SKY)
94  REAL(KIND=8) PSOLSW0(KDLON)     ! SHORTWAVE FLUX AT SURFACE (CLEAR-SKY)
95  !
96  !* LOCAL VARIABLES:
97  !
98  real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
99
100  REAL(KIND=8) ZOZ(KDLON,KFLEV)
101  ! column-density of ozone in layer, in kilo-Dobsons
102
103  REAL(KIND=8) ZAKI(KDLON,2)     
104  REAL(KIND=8) ZCLD(KDLON,KFLEV)
105  REAL(KIND=8) ZCLEAR(KDLON)
106  REAL(KIND=8) ZDSIG(KDLON,KFLEV)
107  REAL(KIND=8) ZFACT(KDLON)
108  REAL(KIND=8) ZFD(KDLON,KFLEV+1)
109  REAL(KIND=8) ZFDOWN(KDLON,KFLEV+1)
110  REAL(KIND=8) ZFU(KDLON,KFLEV+1)
111  REAL(KIND=8) ZFUP(KDLON,KFLEV+1)
112  REAL(KIND=8) ZRMU(KDLON)
113  REAL(KIND=8) ZSEC(KDLON)
114  REAL(KIND=8) ZUD(KDLON,5,KFLEV+1)
115  REAL(KIND=8) ZCLDSW0(KDLON,KFLEV)
116
117  REAL(KIND=8) ZFSUP(KDLON,KFLEV+1)
118  REAL(KIND=8) ZFSDN(KDLON,KFLEV+1)
119  REAL(KIND=8) ZFSUP0(KDLON,KFLEV+1)
120  REAL(KIND=8) ZFSDN0(KDLON,KFLEV+1)
121
122  INTEGER inu, jl, jk, i, k, kpl1
123
124  INTEGER swpas  ! Every swpas steps, sw is calculated
125  PARAMETER(swpas=1)
126
127  INTEGER, SAVE :: itapsw = 0
128  !$OMP THREADPRIVATE(itapsw)
129  LOGICAL, SAVE :: appel1er = .TRUE.
130  !$OMP THREADPRIVATE(appel1er)
131  LOGICAL, SAVE :: initialized = .FALSE.
132  !$OMP THREADPRIVATE(initialized)
133
134  !jq-Introduced for aerosol forcings
135  REAL(KIND=8), SAVE :: flag_aer
136  !$OMP THREADPRIVATE(flag_aer)
137
138  LOGICAL ok_ade, ok_aie    ! use aerosol forcings or not?
139  REAL(KIND=8) tauaero(kdlon,kflev,9,2)  ! aerosol optical properties
140  REAL(KIND=8) pizaero(kdlon,kflev,9,2)  ! (see aeropt.F)
141  REAL(KIND=8) cgaero(kdlon,kflev,9,2)   ! -"-
142  REAL(KIND=8) PTAUA(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS (pre-industrial value)
143  REAL(KIND=8) POMEGAA(KDLON,2,KFLEV)  ! SINGLE SCATTERING ALBEDO
144  REAL(KIND=8) PTOPSWADAERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
145  REAL(KIND=8) PSOLSWADAERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
146  REAL(KIND=8) PTOPSWAD0AERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
147  REAL(KIND=8) PSOLSWAD0AERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
148  REAL(KIND=8) PTOPSWAIAERO(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL IND)
149  REAL(KIND=8) PSOLSWAIAERO(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL IND)
150  REAL(KIND=8) PTOPSWAERO(KDLON,9)       ! SW TOA AS DRF nat & ant
151  REAL(KIND=8) PTOPSW0AERO(KDLON,9)      ! SW SRF AS DRF nat & ant
152  REAL(KIND=8) PSOLSWAERO(KDLON,9)       ! SW TOA CS DRF nat & ant
153  REAL(KIND=8) PSOLSW0AERO(KDLON,9)      ! SW SRF CS DRF nat & ant
154  REAL(KIND=8) PTOPSWCFAERO(KDLON,3)   !  SW TOA AS cloudRF nat & ant
155  REAL(KIND=8) PSOLSWCFAERO(KDLON,3)   !  SW SRF AS cloudRF nat & ant
156
157  !jq - Fluxes including aerosol effects
158  REAL(KIND=8),ALLOCATABLE,SAVE :: ZFSUPAD_AERO(:,:)
159  !$OMP THREADPRIVATE(ZFSUPAD_AERO)
160  REAL(KIND=8),ALLOCATABLE,SAVE :: ZFSDNAD_AERO(:,:)
161  !$OMP THREADPRIVATE(ZFSDNAD_AERO)
162  !jq - Fluxes including aerosol effects
163  REAL(KIND=8),ALLOCATABLE,SAVE :: ZFSUPAD0_AERO(:,:)
164  !$OMP THREADPRIVATE(ZFSUPAD0_AERO)
165  REAL(KIND=8),ALLOCATABLE,SAVE :: ZFSDNAD0_AERO(:,:)
166  !$OMP THREADPRIVATE(ZFSDNAD0_AERO)
167  REAL(KIND=8),ALLOCATABLE,SAVE :: ZFSUPAI_AERO(:,:)
168  !$OMP THREADPRIVATE(ZFSUPAI_AERO)
169  REAL(KIND=8),ALLOCATABLE,SAVE :: ZFSDNAI_AERO(:,:)
170  !$OMP THREADPRIVATE(ZFSDNAI_AERO)
171  REAL(KIND=8),ALLOCATABLE,SAVE ::  ZFSUP_AERO(:,:,:)
172  !$OMP THREADPRIVATE(ZFSUP_AERO)
173  REAL(KIND=8),ALLOCATABLE,SAVE ::  ZFSDN_AERO(:,:,:)
174  !$OMP THREADPRIVATE(ZFSDN_AERO)
175  REAL(KIND=8),ALLOCATABLE,SAVE ::  ZFSUP0_AERO(:,:,:)
176  !$OMP THREADPRIVATE(ZFSUP0_AERO)
177  REAL(KIND=8),ALLOCATABLE,SAVE ::  ZFSDN0_AERO(:,:,:)
178  !$OMP THREADPRIVATE(ZFSDN0_AERO)
179
180! Key to define the aerosol effect acting on climate
181! 0: aerosol feedback active according to ok_ade, ok_aie  DEFAULT
182! 1: no feedback , zero aerosol fluxes are used for climate, diagnostics according to ok_ade_ok_aie
183! 2: feedback according to total aerosol direct effect used for climate, diagnostics according to ok_ade, ok_aie
184! 3: feedback according to natural aerosol direct effect used for climate, diagnostics according to ok_ade_ok_aie
185
186  INTEGER,SAVE :: AEROSOLFEEDBACK_ACTIVE = 0
187!$OMP THREADPRIVATE(AEROSOLFEEDBACK_ACTIVE) 
188
189      CHARACTER (LEN=20) :: modname='sw_aeroAR4'
190      CHARACTER (LEN=80) :: abort_message
191
192  IF ((.not. ok_ade) .and. (AEROSOLFEEDBACK_ACTIVE .ge. 2)) THEN
193     abort_message ='Error: direct effect is not activated but assumed to be active - see sw_aeroAR4.F90'
194     CALL abort_gcm (modname,abort_message,1)
195  ENDIF
196  AEROSOLFEEDBACK_ACTIVE=MIN(MAX(AEROSOLFEEDBACK_ACTIVE,0),3)
197  IF  (AEROSOLFEEDBACK_ACTIVE .gt. 3) THEN
198     abort_message ='Error: AEROSOLFEEDBACK_ACTIVE options go only until 3'
199     CALL abort_gcm (modname,abort_message,1)
200  ENDIF
201
202  IF(.NOT.initialized) THEN
203     flag_aer=0.
204     initialized=.TRUE.
205     ALLOCATE(ZFSUPAD_AERO(KDLON,KFLEV+1))
206     ALLOCATE(ZFSDNAD_AERO(KDLON,KFLEV+1))
207     ALLOCATE(ZFSUPAD0_AERO(KDLON,KFLEV+1))
208     ALLOCATE(ZFSDNAD0_AERO(KDLON,KFLEV+1))
209     ALLOCATE(ZFSUPAI_AERO(KDLON,KFLEV+1))
210     ALLOCATE(ZFSDNAI_AERO(KDLON,KFLEV+1))
211     ALLOCATE(ZFSUP_AERO (KDLON,KFLEV+1,9))
212     ALLOCATE(ZFSDN_AERO (KDLON,KFLEV+1,9))
213     ALLOCATE(ZFSUP0_AERO(KDLON,KFLEV+1,9))
214     ALLOCATE(ZFSDN0_AERO(KDLON,KFLEV+1,9))
215     ZFSUPAD_AERO(:,:)=0.
216     ZFSDNAD_AERO(:,:)=0.
217     ZFSUPAD0_AERO(:,:)=0.
218     ZFSDNAD0_AERO(:,:)=0.
219     ZFSUPAI_AERO(:,:)=0.
220     ZFSDNAI_AERO(:,:)=0.
221     ZFSUP_AERO (:,:,:)=0.
222     ZFSDN_AERO (:,:,:)=0.
223     ZFSUP0_AERO(:,:,:)=0.
224     ZFSDN0_AERO(:,:,:)=0.
225  ENDIF
226
227  IF (appel1er) THEN
228     WRITE(lunout,*) 'SW calling frequency : ', swpas
229     WRITE(lunout,*) "   In general, it should be 1"
230     appel1er = .FALSE.
231  ENDIF
232  !     ------------------------------------------------------------------
233  IF (MOD(itapsw,swpas).EQ.0) THEN
234
235     DO JK = 1 , KFLEV
236        DO JL = 1, KDLON
237           ZCLDSW0(JL,JK) = 0.0
238           ZOZ(JL,JK) = POZON(JL,JK)*46.6968/RG &
239                *PDP(JL,JK)*(101325.0/PPSOL(JL))
240        ENDDO
241     ENDDO
242
243! clear sky is either computed IF no direct effect is asked for, or for extended diag)
244     IF (( lev_histmth .ge. 4 ) .or. ( .not. ok_ade )) THEN   
245
246     ! clear-sky: zero aerosol effect
247     flag_aer=0.0
248     CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
249          PRMU0,PFRAC,PTAVE,PWV,&
250          ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
251     INU = 1
252     CALL SW1S_LMDAR4(INU,PAER, flag_aer, &
253          tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
254          PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,&
255          ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
256          ZFD, ZFU)
257     INU = 2
258     CALL SW2S_LMDAR4(INU, PAER, flag_aer, &
259          tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
260          ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,&
261          ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
262          PWV, PQS,&
263          ZFDOWN, ZFUP)
264     DO JK = 1 , KFLEV+1
265        DO JL = 1, KDLON
266           ZFSUP0_AERO(JL,JK,1) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
267           ZFSDN0_AERO(JL,JK,1) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
268        ENDDO
269     ENDDO
270     ENDIF
271
272! cloudy sky is either computed IF no indirect effect is asked for, or for extended diag)
273     IF (( lev_histmth .ge. 4 ) .or. ( .not. ok_aie )) THEN   
274     ! cloudy-sky: zero aerosol effect
275     flag_aer=0.0
276     CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
277          PRMU0,PFRAC,PTAVE,PWV,&
278          ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
279     INU = 1
280     CALL SW1S_LMDAR4(INU, PAER, flag_aer, &
281          tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
282          PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
283          ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
284          ZFD, ZFU)
285     INU = 2
286     CALL SW2S_LMDAR4(INU, PAER, flag_aer, &
287          tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),&
288          ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
289          ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
290          PWV, PQS,&
291          ZFDOWN, ZFUP)
292
293     DO JK = 1 , KFLEV+1
294        DO JL = 1, KDLON
295           ZFSUP_AERO(JL,JK,1) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
296           ZFSDN_AERO(JL,JK,1) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
297        ENDDO
298     ENDDO
299     ENDIF
300
301
302     IF (ok_ade) THEN
303
304        ! clear sky (Anne Cozic 03/07/2007) direct effect of total aerosol
305        ! CAS AER (2)
306        flag_aer=1.0
307        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
308             PRMU0,PFRAC,PTAVE,PWV,&
309             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
310        INU = 1
311        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
312             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
313             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
314             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
315             ZFD, ZFU)
316        INU = 2
317        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
318             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
319             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
320             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
321             PWV, PQS,&
322             ZFDOWN, ZFUP)
323
324        DO JK = 1 , KFLEV+1
325           DO JL = 1, KDLON
326              ZFSUP0_AERO(JL,JK,2) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
327              ZFSDN0_AERO(JL,JK,2) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
328           ENDDO
329        ENDDO
330
331! cloudy sky is either computed IF no indirect effect is asked for, or for extended diag)
332        IF (( lev_histmth .ge. 2 ) .or. (.not. ok_aie)) THEN 
333        ! cloudy-sky aerosol direct effect of total aerosol
334        flag_aer=1.0
335        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
336             PRMU0,PFRAC,PTAVE,PWV,&
337             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
338        INU = 1
339        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
340             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
341             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
342             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
343             ZFD, ZFU)
344        INU = 2
345        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
346             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
347             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
348             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
349             PWV, PQS,&
350             ZFDOWN, ZFUP)
351
352        DO JK = 1 , KFLEV+1
353           DO JL = 1, KDLON
354              ZFSUP_AERO(JL,JK,2) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
355              ZFSDN_AERO(JL,JK,2) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
356           ENDDO
357        ENDDO
358        ENDIF
359
360! natural aeroosl clear sky is  computed  for extended diag)
361        IF ( lev_histmth .ge. 4 ) THEN           
362        ! clear sky direct effect natural aerosol
363        flag_aer=1.0
364        CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
365             PRMU0,PFRAC,PTAVE,PWV,&
366             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
367        INU = 1
368        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
369             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
370             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
371             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
372             ZFD, ZFU)
373        INU = 2
374        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
375             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
376             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
377             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
378             PWV, PQS,&
379             ZFDOWN, ZFUP)
380
381        DO JK = 1 , KFLEV+1
382           DO JL = 1, KDLON
383              ZFSUP0_AERO(JL,JK,3) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
384              ZFSDN0_AERO(JL,JK,3) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
385           ENDDO
386        ENDDO
387        ENDIF
388
389! cloud sky natural is for extended diagnostics
390        IF ( lev_histmth .ge. 2 ) THEN
391        ! cloudy-sky direct effect natural aerosol
392        flag_aer=1.0
393        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
394             PRMU0,PFRAC,PTAVE,PWV,&
395             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
396        INU = 1
397        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
398             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
399             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
400             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
401             ZFD, ZFU)
402        INU = 2
403        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
404             tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),&
405             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
406             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
407             PWV, PQS,&
408             ZFDOWN, ZFUP)
409
410        DO JK = 1 , KFLEV+1
411           DO JL = 1, KDLON
412              ZFSUP_AERO(JL,JK,3) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
413              ZFSDN_AERO(JL,JK,3) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
414           ENDDO
415        ENDDO
416        ENDIF
417
418     ENDIF ! ok_ade
419
420! cloudy sky needs to be computed in all cases IF ok_aie is activated
421     IF (ok_aie) THEN
422        !jq   cloudy-sky + aerosol direct + aerosol indirect of total aerosol
423        flag_aer=1.0
424        CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
425             PRMU0,PFRAC,PTAVE,PWV,&
426             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
427        INU = 1
428        CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
429             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
430             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
431             ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,&
432             ZFD, ZFU)
433        INU = 2
434        CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
435             tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),&
436             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
437             ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,&
438             PWV, PQS,&
439             ZFDOWN, ZFUP)
440        DO JK = 1 , KFLEV+1
441           DO JL = 1, KDLON
442              ZFSUP_AERO(JL,JK,4) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
443              ZFSDN_AERO(JL,JK,4) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
444           ENDDO
445        ENDDO
446     ENDIF ! ok_aie     
447
448     itapsw = 0
449  ENDIF
450  itapsw = itapsw + 1
451
452  IF  ( AEROSOLFEEDBACK_ACTIVE .eq. 0) THEN
453  IF ( ok_ade .and. ok_aie  ) THEN
454    ZFSUP(:,:) =    ZFSUP_AERO(:,:,4)
455    ZFSDN(:,:) =    ZFSDN_AERO(:,:,4)
456    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,2)
457    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,2)
458  ENDIF
459  IF ( ok_ade .and. (.not. ok_aie) )  THEN
460    ZFSUP(:,:) =    ZFSUP_AERO(:,:,2)
461    ZFSDN(:,:) =    ZFSDN_AERO(:,:,2)
462    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,2)
463    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,2)
464  ENDIF
465
466  IF ( (.not. ok_ade) .and. ok_aie  )  THEN
467    print*,'Warning: indirect effect in cloudy regions includes direct aerosol effect'
468    ZFSUP(:,:) =    ZFSUP_AERO(:,:,4)
469    ZFSDN(:,:) =    ZFSDN_AERO(:,:,4)
470    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,1)
471    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,1)
472  ENDIF
473  IF ((.not. ok_ade) .and. (.not. ok_aie)) THEN
474    ZFSUP(:,:) =    ZFSUP_AERO(:,:,1)
475    ZFSDN(:,:) =    ZFSDN_AERO(:,:,1)
476    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,1)
477    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,1)
478  ENDIF
479
480! MS the following allows to compute the forcing diagostics without
481! letting the aerosol forcing act on the meteorology
482! SEE logic above
483  ELSEIF  ( AEROSOLFEEDBACK_ACTIVE .gt. 0) THEN
484    ZFSUP(:,:) =    ZFSUP_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
485    ZFSDN(:,:) =    ZFSDN_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
486    ZFSUP0(:,:) =   ZFSUP0_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
487    ZFSDN0(:,:) =   ZFSDN0_AERO(:,:,AEROSOLFEEDBACK_ACTIVE)
488  ENDIF
489 
490
491  DO k = 1, KFLEV
492     kpl1 = k+1
493     DO i = 1, KDLON
494        PHEAT(i,k) = -(ZFSUP(i,kpl1)-ZFSUP(i,k))-(ZFSDN(i,k)-ZFSDN(i,kpl1))
495        PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k)
496        PHEAT0(i,k) = -(ZFSUP0(i,kpl1)-ZFSUP0(i,k))-(ZFSDN0(i,k)-ZFSDN0(i,kpl1))
497        PHEAT0(i,k) = PHEAT0(i,k) * RDAY*RG/RCPD / PDP(i,k)
498     ENDDO
499  ENDDO
500
501  DO i = 1, KDLON
502! effective SW surface albedo calculation
503     PALBPLA(i) = ZFSUP(i,KFLEV+1)/(ZFSDN(i,KFLEV+1)+1.0e-20)
504     
505! clear sky net fluxes at TOA and SRF
506     PSOLSW0(i) = ZFSDN0(i,1) - ZFSUP0(i,1)
507     PTOPSW0(i) = ZFSDN0(i,KFLEV+1) - ZFSUP0(i,KFLEV+1)
508
509! cloudy sky net fluxes at TOA and SRF
510     PSOLSW(i) = ZFSDN(i,1) - ZFSUP(i,1)
511     PTOPSW(i) = ZFSDN(i,KFLEV+1) - ZFSUP(i,KFLEV+1)
512
513
514! net anthropogenic forcing direct and 1st indirect effect diagnostics
515! requires a natural aerosol field read and used
516! Difference of net fluxes from double call to radiation
517
518
519IF (ok_ade) THEN
520
521! indices 1: natural; 2 anthropogenic
522! TOA/SRF all sky natural forcing
523     PSOLSWAERO(i,1) = (ZFSDN_AERO(i,1,3) - ZFSUP_AERO(i,1,3))-(ZFSDN_AERO(i,1,1) - ZFSUP_AERO(i,1,1))
524     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))
525
526! TOA/SRF all sky anthropogenic forcing
527     PSOLSWAERO(i,2) = (ZFSDN_AERO(i,1,2) - ZFSUP_AERO(i,1,2))-(ZFSDN_AERO(i,1,3) - ZFSUP_AERO(i,1,3))
528     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))
529
530! TOA/SRF clear sky natural forcing
531     PSOLSW0AERO(i,1) = (ZFSDN0_AERO(i,1,3) - ZFSUP0_AERO(i,1,3))-(ZFSDN0_AERO(i,1,1) - ZFSUP0_AERO(i,1,1))
532     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))
533
534! TOA/SRF clear sky anthropogenic forcing
535     PSOLSW0AERO(i,2) = (ZFSDN0_AERO(i,1,2) - ZFSUP0_AERO(i,1,2))-(ZFSDN0_AERO(i,1,3) - ZFSUP0_AERO(i,1,3))
536     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))
537
538! Cloud forcing indices 1: natural; 2 anthropogenic; 3: zero aerosol direct effect
539! Instantaneously computed cloudy sky direct aerosol effect, cloud forcing due to aerosols above clouds
540! natural
541     PSOLSWCFAERO(i,1) = PSOLSWAERO(i,1) - PSOLSW0AERO(i,1)
542     PTOPSWCFAERO(i,1) = PTOPSWAERO(i,1) - PTOPSW0AERO(i,1)
543
544! Instantaneously computed cloudy SKY DIRECT aerosol effect, cloud forcing due to aerosols above clouds
545! anthropogenic
546     PSOLSWCFAERO(i,2) = PSOLSWAERO(i,2) - PSOLSW0AERO(i,2)
547     PTOPSWCFAERO(i,2) = PTOPSWAERO(i,2) - PTOPSW0AERO(i,2)
548
549! Cloudforcing without aerosol
550! zero
551     PSOLSWCFAERO(i,3) = (ZFSDN_AERO(i,1,1) - ZFSUP_AERO(i,1,1))-(ZFSDN0_AERO(i,1,1) - ZFSUP0_AERO(i,1,1))
552     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))
553
554! direct anthropogenic forcing , as in old LMDzT, however differences of net fluxes
555     PSOLSWADAERO(i) = PSOLSWAERO(i,2)
556     PTOPSWADAERO(i) = PTOPSWAERO(i,2)
557     PSOLSWAD0AERO(i) = PSOLSW0AERO(i,2)
558     PTOPSWAD0AERO(i) = PTOPSW0AERO(i,2)
559
560ENDIF
561
562
563IF (ok_aie) THEN
564     PSOLSWAIAERO(i) = (ZFSDN_AERO(i,1,4) - ZFSUP_AERO(i,1,4))-(ZFSDN_AERO(i,1,2) - ZFSUP_AERO(i,1,2))
565     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))
566ENDIF
567
568  ENDDO
569END SUBROUTINE SW_AEROAR4
Note: See TracBrowser for help on using the repository browser.