source: CONFIG/publications/ICOLMDZORINCA_CO2_Transport_GMD_2023/INCA/src/INCA_SRC/radlwsw_inca.F90 @ 6610

Last change on this file since 6610 was 6610, checked in by acosce, 10 months ago

INCA used for ICOLMDZORINCA_CO2_Transport_GMD_2023
File size: 40.8 KB
Line 
1! $Id: radlwsw_inca.F90 163 2010-02-22 15:41:45Z acosce $
2!! =========================================================================
3!! INCA - INteraction with Chemistry and Aerosols
4!!
5!! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE)
6!!           Unite mixte CEA-CNRS-UVSQ
7!!
8!! Contributors to this INCA subroutine:
9!!
10!! Celine Deandreis, LSCE
11!!
12!! Anne Cozic, LSCE, anne.cozic@cea.fr
13!!
14!! This software is a computer program whose purpose is to simulate the
15!! atmospheric gas phase and aerosol composition. The model is designed to be
16!! used within a transport model or a general circulation model. This version
17!! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts
18!! for emissions, transport (resolved and sub-grid scale), photochemical
19!! transformations, and scavenging (dry deposition and washout) of chemical
20!! species and aerosols interactively in the GCM. Several versions of the INCA
21!! model are currently used depending on the envisaged applications with the
22!! chemistry-climate model.
23!!
24!! This software is governed by the CeCILL  license under French law and
25!! abiding by the rules of distribution of free software.  You can  use,
26!! modify and/ or redistribute the software under the terms of the CeCILL
27!! license as circulated by CEA, CNRS and INRIA at the following URL
28!! "http://www.cecill.info".
29!!
30!! As a counterpart to the access to the source code and  rights to copy,
31!! modify and redistribute granted by the license, users are provided only
32!! with a limited warranty  and the software's author,  the holder of the
33!! economic rights,  and the successive licensors  have only  limited
34!! liability.
35!!
36!! In this respect, the user's attention is drawn to the risks associated
37!! with loading,  using,  modifying and/or developing or reproducing the
38!! software by the user in light of its specific status of free software,
39!! that may mean  that it is complicated to manipulate,  and  that  also
40!! therefore means  that it is reserved for developers  and  experienced
41!! professionals having in-depth computer knowledge. Users are therefore
42!! encouraged to load and test the software's suitability as regards their
43!! requirements in conditions enabling the security of their systems and/or
44!! data to be ensured and,  more generally, to use and operate it in the
45!! same conditions as regards security.
46!!
47!! The fact that you are presently reading this means that you have had
48!! knowledge of the CeCILL license and that you accept its terms.
49!! =========================================================================
50
51#include <inca_define.h>
52
53SUBROUTINE radlwsw_inca(chemistry_couple, kdlon,kflev,dist, rmu0, fract, &
54   solaire,paprs, pplay,tsol,albedo, alblw, t,q,size_wo, wo,&
55   cldfra, cldemi, cldtaupd,&
56   heat,heat0,cool,cool0,albpla,&
57   topsw,toplw,solsw,sollw,&
58   sollwdown,&
59   topsw0,toplw0,solsw0,sollw0,&
60   lwdn0, lwdn, lwup0, lwup,&
61   swdn0, swdn, swup0, swup,&
62   ok_ade, ok_aie,&
63   tau_inca, piz_inca, cg_inca,&
64   topswad_inca, solswad_inca,&
65   topswad0_inca, solswad0_inca,&
66   topsw_inca, topsw0_inca,&
67   solsw_inca, solsw0_inca,&
68   cldtaupi, topswai_inca, solswai_inca)
69
70  USE INCA_DIM
71  USE PRINT_INCA
72  USE AEROSOL_DIAG 
73  USE CHEM_MODS, ONLY : o3_inca
74
75  IMPLICIT NONE
76  !======================================================================
77  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
78  ! Objet: interface entre le modele et les rayonnements
79  ! Arguments:
80  ! dist-----input-R- distance astronomique terre-soleil
81  ! rmu0-----input-R- cosinus de l'angle zenithal
82  ! fract----input-R- duree d'ensoleillement normalisee
83  ! co2_ppm--input-R- concentration du gaz carbonique (en ppm)
84  ! solaire--input-R- constante solaire (W/m**2)
85  ! paprs----input-R- pression a inter-couche (Pa)
86  ! pplay----input-R- pression au milieu de couche (Pa)
87  ! tsol-----input-R- temperature du sol (en K)
88  ! albedo---input-R- albedo du sol (entre 0 et 1)
89  ! t--------input-R- temperature (K)
90  ! q--------input-R- vapeur d'eau (en kg/kg)
91  ! wo-------input-R- contenu en ozone (en kilo-Dobsons 21/10/2016) correction MPL 100505
92  ! cldfra---input-R- fraction nuageuse (entre 0 et 1)
93  ! cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
94  ! cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
95  ! ok_ade---input-L- apply the Aerosol Direct Effect or not?
96  ! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
97  ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
98  ! cldtaupi-input-R- epaisseur optique des nuages dans le visible
99  !                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
100  !                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
101  !                   it is needed for the diagnostics of the aerosol indirect radiative forcing     
102  !
103  ! heat-----output-R- echauffement atmospherique (visible) (K/jour)
104  ! cool-----output-R- refroidissement dans l'IR (K/jour)
105  ! albpla---output-R- albedo planetaire (entre 0 et 1)
106  ! topsw----output-R- flux solaire net au sommet de l'atm.
107  ! toplw----output-R- ray. IR montant au sommet de l'atmosphere
108  ! solsw----output-R- flux solaire net a la surface
109  ! sollw----output-R- ray. IR montant a la surface
110  ! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
111  ! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
112  ! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
113  ! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
114  !
115  ! ATTENTION: swai and swad have to be interpreted in the following manner:
116  ! ---------
117  ! ok_ade=F & ok_aie=F -both are zero
118  ! ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
119  !                        indirect is zero
120  ! ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
121  !                        direct is zero
122  ! ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
123  !                        aerosol direct forcing is F_{AD} = topswai-topswad
124  !
125 
126  !======================================================================
127 
128  ! ====================================================================
129  ! Adapte au modele de chimie INCA par Celine Deandreis -- 2007
130  ! 1 = ZERO   
131  ! 2 = AER total   
132  ! 3 = NAT   
133  ! 4 = BC   
134  ! 5 = SO4   
135  ! 6 = POM   
136  ! 7 = DUST
137  ! 8 = SS
138  ! 9 = FNO3
139  ! 10 = DNO3   
140  ! 11 = SNO3
141  !
142  ! ====================================================================
143#include "YOETHF_I.h"
144#include "YOMCST_I.h"
145  EXTERNAL lw_LMDAR4, sw_LMDAR4
146
147
148  LOGICAL, INTENT(in) :: chemistry_couple
149  INTEGER, INTENT(in) :: kdlon,kflev
150  REAL, INTENT(in) :: solaire
151  REAL, INTENT(in) :: dist
152  REAL, INTENT(in) :: rmu0(PLON), fract(PLON)
153  REAL,INTENT(in) :: paprs(PLON,PLEV+1), pplay(PLON,PLEV)
154  REAL,INTENT(in) ::albedo(PLON), alblw(PLON), tsol(PLON)
155  REAL,INTENT(in) :: t(PLON,PLEV), q(PLON,PLEV)
156  INTEGER, INTENT(in) :: size_wo
157  REAL, INTENT(in):: wo(PLON,PLEV,size_wo)  ! column-density of ozone in a layer, in kilo-Dobsons
158                                            ! "wo(:, :, 1)" is for the average day-night field,
159                                            ! "wo(:, :, 2)" is for daylight time.
160  LOGICAL, INTENT(in) :: ok_ade, ok_aie     ! switches whether to use aerosol direct (indirect) effects or not
161  REAL, INTENT(in) :: cldfra(PLON,PLEV), cldemi(PLON,PLEV), cldtaupd(PLON,PLEV)
162  REAL, INTENT(in) :: tau_inca(PLON,PLEV,naero_grp,2) ! aerosol optical properties (see aeropt.F)
163  REAL, INTENT(in) :: piz_inca(PLON,PLEV,naero_grp,2) ! aerosol optical properties (see aeropt.F)
164  REAL, INTENT(in) :: cg_inca(PLON,PLEV,naero_grp,2)        ! aerosol optical properties (see aeropt.F)
165  REAL, INTENT(in) :: cldtaupi(PLON,PLEV)  ! cloud optical thickness for pre-industrial aerosol concentrations
166
167
168  REAL, INTENT(out):: heat(PLON,PLEV), cool(PLON,PLEV)
169  REAL, INTENT(out):: heat0(PLON,PLEV), cool0(PLON,PLEV)
170  REAL, INTENT(out) :: topsw(PLON), toplw(PLON)
171  REAL, INTENT(out) :: solsw(PLON), sollw(PLON), albpla(PLON)
172  REAL, INTENT(out) :: topsw0(PLON), toplw0(PLON), solsw0(PLON), sollw0(PLON)
173  REAL, INTENT(out) :: sollwdown(PLON)
174  REAL, INTENT(out) :: swdn(PLON,kflev+1),swdn0(PLON,kflev+1)
175  REAL, INTENT(out) :: swup(PLON,kflev+1),swup0(PLON,kflev+1)
176  REAL, INTENT(out) :: lwdn(PLON,kflev+1),lwdn0(PLON,kflev+1)
177  REAL, INTENT(out) :: lwup(PLON,kflev+1),lwup0(PLON,kflev+1)
178  REAL, INTENT(out) :: topswad_inca(PLON), solswad_inca(PLON) ! output: aerosol direct forcing at TOA and surface
179  REAL, INTENT(out) :: topswad0_inca(PLON), solswad0_inca(PLON) ! output: aerosol direct forcing at TOA and surface
180  REAL, INTENT(out) :: topswai_inca(PLON), solswai_inca(PLON) ! output: aerosol indirect forcing atTOA and surface
181  REAL*8, INTENT(out) :: topsw_inca(kdlon,naero_grp), topsw0_inca(kdlon,naero_grp)
182  REAL*8, INTENT(out) :: solsw_inca(kdlon,naero_grp), solsw0_inca(kdlon,naero_grp)
183
184#ifdef AER
185
186
187  REAL*8 ZFSUP(KDLON,KFLEV+1)
188  REAL*8 ZFSDN(KDLON,KFLEV+1)
189  REAL*8 ZFSUP0(KDLON,KFLEV+1)
190  REAL*8 ZFSDN0(KDLON,KFLEV+1)
191  REAL*8 ZFLUP(KDLON,KFLEV+1)
192  REAL*8 ZFLDN(KDLON,KFLEV+1)
193  REAL*8 ZFLUP0(KDLON,KFLEV+1)
194  REAL*8 ZFLDN0(KDLON,KFLEV+1)
195  REAL*8 zx_alpha1, zx_alpha2
196  INTEGER k, kk, i, j, iof, nb_gr
197  REAL*8 PSCT
198  REAL*8 PALBD(kdlon,2), PALBP(kdlon,2)
199  REAL*8 PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
200  REAL*8 PPSOL(kdlon), PDP(kdlon,PLEV)
201  REAL*8 PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)
202  REAL*8 PTAVE(kdlon,kflev)
203  REAL*8 PWV(kdlon,kflev), PQS(kdlon,kflev)
204  REAL*8 POZON(kdlon,kflev, size(wo,3))! mass fraction of ozone
205  ! "POZON(:, :, 1)" is for the average day-night field,
206  ! "POZON(:, :, 2)" is for daylight time.
207
208  REAL*8 PAER(kdlon,kflev,5)
209  REAL*8 PCLDLD(kdlon,kflev)
210  REAL*8 PCLDLU(kdlon,kflev)
211  REAL*8 PCLDSW(kdlon,kflev)
212  REAL*8 PTAU(kdlon,2,kflev)
213  REAL*8 POMEGA(kdlon,2,kflev)
214  REAL*8 PCG(kdlon,2,kflev)
215  REAL*8 zfract(kdlon), zrmu0(kdlon), zdist
216  REAL*8 zheat(kdlon,kflev), zcool(kdlon,kflev)
217  REAL*8 zheat0(kdlon,kflev), zcool0(kdlon,kflev)
218  REAL*8 ztopsw(kdlon), ztoplw(kdlon)
219  REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
220  REAL*8 zsollwdown(kdlon)
221  REAL*8 ztopsw0(kdlon), ztoplw0(kdlon)
222  REAL*8 zsolsw0(kdlon), zsollw0(kdlon)
223  REAL*8 zznormcp
224  REAL*8 tauinca(kdlon,kflev,naero_grp,2) ! aer opt properties
225  REAL*8 pizinca(kdlon,kflev,naero_grp,2)
226  REAL*8 cginca(kdlon,kflev,naero_grp,2)
227  REAL*8 PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
228  REAL*8 POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo
229  REAL*8 ztopswadinca(kdlon), zsolswadinca(kdlon) ! Aerosol direct forcing at TOAand surface
230  REAL*8 ztopswad0inca(kdlon), zsolswad0inca(kdlon) ! Aerosol direct forcing at TOAand surface
231  REAL*8 ztopswaiinca(kdlon), zsolswaiinca(kdlon) ! dito, indirect
232  REAL*8 ztopsw_inca(kdlon,naero_grp), ztopsw0_inca(kdlon,naero_grp)
233  REAL*8 zsolsw_inca(kdlon,naero_grp), zsolsw0_inca(kdlon,naero_grp)
234  REAL*8 ZCLEAR(KDLON) 
235  real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
236
237  ! initialisation
238  tauinca(:,:,:,:)=0.
239  pizinca(:,:,:,:)=0.
240  cginca(:,:,:,:)=0.
241 
242  !
243  !-------------------------------------------
244  nb_gr = PLON / kdlon
245  IF (nb_gr*kdlon .NE. PLON) THEN
246      WRITE(lunout,*)  "kdlon mauvais:", PLON, kdlon, nb_gr
247      CALL abort
248  ENDIF
249  IF (kflev .NE. PLEV) THEN
250      WRITE(lunout,*) "kflev differe de PLEV, kflev, PLEV"
251      CALL abort
252  ENDIF
253  !-------------------------------------------
254  DO k = 1, PLEV
255    DO i = 1, PLON
256      heat(i,k)=0.
257      cool(i,k)=0.
258      heat0(i,k)=0.
259      cool0(i,k)=0.
260    ENDDO
261  ENDDO
262  !
263  zdist = dist
264  !
265  PSCT = solaire/zdist/zdist
266  DO j = 1, nb_gr
267    iof = kdlon*(j-1)
268    DO i = 1, kdlon
269      zfract(i) = fract(iof+i)
270      zrmu0(i) = rmu0(iof+i)
271      PALBD(i,1) = albedo(iof+i)
272      PALBD(i,2) = alblw(iof+i)
273      PALBP(i,1) = albedo(iof+i)
274      PALBP(i,2) = alblw(iof+i)
275      PEMIS(i) = 1.0 
276      PVIEW(i) = 1.66
277      PPSOL(i) = paprs(iof+i,1)
278      zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2))
279      zx_alpha2 = 1.0 - zx_alpha1
280      PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
281      PTL(i,PLEV+1) = t(iof+i,PLEV)
282      PDT0(i) = tsol(iof+i) - PTL(i,1)
283    ENDDO
284    DO k = 2, kflev
285      DO i = 1, kdlon
286        PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
287      ENDDO
288    ENDDO
289    DO k = 1, kflev
290      DO i = 1, kdlon
291        PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
292        PTAVE(i,k) = t(iof+i,k)
293        PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
294        PQS(i,k) = PWV(i,k)
295!       Confert from  column density of ozone in a cell, in kDU, to a mass fraction
296        POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
297             / (paprs(iof+i, k) - paprs(iof+i, k+1))
298
299!Old version : 21/10/2016 - correction d'une modification d'unite datant de 06/2009 dans lmdz
300        ! wo:    cm.atm (epaisseur en cm dans la situation standard)
301        ! POZON: kg/kg
302!        POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968 &
303!           /(paprs(iof+i,k)-paprs(iof+i,k+1))&
304!           *(paprs(iof+i,1)/101325.0)
305        PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
306        PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
307        PCLDSW(i,k) = cldfra(iof+i,k)
308        PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
309        PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
310        POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
311        POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
312        PCG(i,1,k) = 0.865
313        PCG(i,2,k) = 0.910
314        !-
315        ! Introduced for aerosol indirect forcings.
316        ! The following values use the cloud optical thickness calculated from
317        ! present-day aerosol concentrations whereas the quantities without the
318        ! "A" at the end are for pre-industial (natural-only) aerosol concentrations
319        !
320        PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
321        PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
322        POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
323        POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
324      ENDDO
325    ENDDO
326
327    if (chemistry_couple) then
328       POZON(:,:,1) = o3_inca(:,:) 
329    endif
330    !
331    DO k = 1, kflev+1
332      DO i = 1, kdlon
333        PPMB(i,k) = paprs(iof+i,k)/100.0
334      ENDDO
335    ENDDO
336    !
337    DO kk = 1, 5
338      DO k = 1, kflev
339        DO i = 1, kdlon
340          PAER(i,k,kk) = 1.0E-15
341        ENDDO
342      ENDDO
343    ENDDO
344    DO k = 1, kflev
345      DO i = 1, kdlon
346        tauinca(i,k,:,1)=tau_inca(iof+i,k,:,1)
347        pizinca(i,k,:,1)=piz_inca(iof+i,k,:,1)
348        cginca(i,k,:,1) =cg_inca(iof+i,k,:,1)
349        tauinca(i,k,:,2)=tau_inca(iof+i,k,:,2)
350        pizinca(i,k,:,2)=piz_inca(iof+i,k,:,2)
351        cginca(i,k,:,2) =cg_inca(iof+i,k,:,2)
352      ENDDO
353    ENDDO
354    !
355    !======================================================================
356    CALL LW_LMDAR4(&
357       PPMB, PDP,&
358       PPSOL,PDT0,PEMIS,&
359       PTL, PTAVE, PWV, POZON(:, :, 1), PAER,&
360       PCLDLD,PCLDLU,&
361       PVIEW,&
362       zcool, zcool0,&
363       ztoplw,zsollw,ztoplw0,zsollw0,&
364       zsollwdown,&
365       ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
366
367    CALL SW_INCA(kdlon,kflev,PSCT, zrmu0, zfract,&
368       PPMB, PDP,&
369       PPSOL, PALBD, PALBP,&
370       PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
371       PCLDSW, PTAU, POMEGA, PCG,&
372       zheat, zheat0,&
373       zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,&
374       ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
375       tauinca, pizinca, cginca, &! aerosol optical properties
376       PTAUA, POMEGAA,&
377       ztopswadinca,zsolswadinca,&
378       ztopswad0inca,zsolswad0inca,&
379       ztopswaiinca,zsolswaiinca, & ! diagnosed aerosol forcing
380       ztopsw_inca,ztopsw0_inca,&
381       zsolsw_inca,zsolsw0_inca,&
382       ok_ade, ok_aie,ZCLEAR) ! apply aerosol effects or not?
383   
384       DO i=1, kdlon
385        cldfract(i)=1-ZCLEAR(i) 
386       END DO
387
388    !======================================================================
389    DO i = 1, kdlon
390      topsw(iof+i) = ztopsw(i)
391      toplw(iof+i) = ztoplw(i)
392      solsw(iof+i) = zsolsw(i)
393      sollw(iof+i) = zsollw(i)
394      sollwdown(iof+i) = zsollwdown(i)
395      DO k = 1, kflev+1
396        lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
397        lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
398        lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
399        lwup  ( iof+i,k)   = ZFLUP  ( i,k)
400      ENDDO
401      topsw0(iof+i) = ztopsw0(i)
402      toplw0(iof+i) = ztoplw0(i)
403      solsw0(iof+i) = zsolsw0(i)
404      sollw0(iof+i) = zsollw0(i)
405      albpla(iof+i) = zalbpla(i)
406      DO k = 1, kflev+1
407        swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
408        swdn  ( iof+i,k)   = ZFSDN  ( i,k)
409        swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
410        swup  ( iof+i,k)   = ZFSUP  ( i,k)
411      ENDDO
412    ENDDO
413    !-transform the aerosol forcings, if they have
414    ! to be calculated
415    IF (ok_ade) THEN
416      DO i = 1, kdlon
417!         topsw_inca(iof+i,:) = ztopsw_inca(iof+i,:)
418!         topsw0_inca(iof+i,:) = ztopsw0_inca(iof+i,:)
419!         solsw_inca(iof+i,:) = zsolsw_inca(iof+i,:)
420!         solsw0_inca(iof+i,:) = zsolsw0_inca(iof+i,:)
421         ! le calcul de topsw_inca et cie a ete modifie en 2013 suite
422         ! a des modifications du code cote lmdz. Ce diagnostique
423         ! ne comporte plus que les flux net pour ant et nat (defini
424         ! sur 1 et 2) - verif en 2015 Yves et Anne
425
426          ! nat
427          topsw_inca(iof+i,1) = ztopsw_inca(i,3)-ztopsw_inca(i,1)
428          ! ant
429          topsw_inca(iof+i,2) = ztopsw_inca(i,2)-ztopsw_inca(i,3)
430          ! nat
431          topsw0_inca(iof+i,1) = ztopsw0_inca(i,3)-ztopsw0_inca(i,1)
432          ! ant
433          topsw0_inca(iof+i,2) = ztopsw0_inca(i,2)-ztopsw0_inca(i,3)
434          ! nat
435          solsw_inca(iof+i,1) = zsolsw_inca(i,3)-zsolsw_inca(i,1)
436          ! ant
437          solsw_inca(iof+i,2) = zsolsw_inca(i,2)-zsolsw_inca(i,3)
438          ! nat
439          solsw0_inca(iof+i,1) = zsolsw0_inca(i,3)-zsolsw0_inca(i,1)
440          ! ant
441          solsw0_inca(iof+i,2) = zsolsw0_inca(i,2)-zsolsw0_inca(i,3)
442
443
444
445!RAFF AD=direct antro
446          topswad_inca(iof+i) = ztopsw_inca(i,2)-ztopsw_inca(i,3)
447          solswad_inca(iof+i) = zsolsw_inca(i,2)-zsolsw_inca(i,3)
448          topswad0_inca(iof+i) = ztopsw0_inca(i,2)-ztopsw0_inca(i,3)
449          solswad0_inca(iof+i) = zsolsw0_inca(i,2)-zsolsw0_inca(i,3)
450
451        ENDDO
452    ELSE
453        DO i = 1, kdlon
454          topswad_inca(iof+i) = 0.0
455          solswad_inca(iof+i) = 0.0
456          topswad0_inca(iof+i) = 0.0
457          solswad0_inca(iof+i) = 0.0
458          topsw_inca(iof+i,:) = 0.
459          topsw0_inca(iof+i,:) =0.
460          solsw_inca(iof+i,:) = 0.
461          solsw0_inca(iof+i,:) = 0.
462        ENDDO
463    ENDIF
464    IF (ok_aie) THEN
465        DO i = 1, kdlon
466!          topswai_inca(iof+i) = ztopswaiinca(i)
467!          solswai_inca(iof+i) = zsolswaiinca(i)
468          topswai_inca(iof+i) = ztopsw_inca(i,2)-ztopswaiinca(i)
469          solswai_inca(iof+i) = zsolsw_inca(i,2)- zsolswaiinca(i)
470        ENDDO
471    ELSE
472        DO i = 1, kdlon
473          topswai_inca(iof+i) = 0.0
474          solswai_inca(iof+i) = 0.0
475        ENDDO
476    ENDIF
477    DO k = 1, kflev
478      DO i = 1, kdlon
479        !        scale factor to take into account the difference between
480        !        dry air and watter vapour scpecifi! heat capacity
481        zznormcp=1.0+RVTMP2*PWV(i,k)
482        heat(iof+i,k) = zheat(i,k)/zznormcp
483        cool(iof+i,k) = zcool(i,k)/zznormcp
484        heat0(iof+i,k) = zheat0(i,k)/zznormcp
485        cool0(iof+i,k) = zcool0(i,k)/zznormcp
486      ENDDO
487    ENDDO
488    !
489  ENDDO
490
491  swtoaas_ad(:) =  topswad_inca(:) 
492  swtoacs_ad(:) =  topswad0_inca(:)
493  if (ok_aie) then
494     DO i = 1, kdlon
495        swtoaas_ai(iof+i) = ztopswaiinca(i)
496        swsrfas_ai(iof+i) = zsolswaiinca(i)
497     enddo
498  else
499     swtoaas_ai(:) = 0.0
500     swsrfas_ai(:) = 0.0
501  endif
502  if (ok_ade) then
503     swtoaas(:,:)  =  ztopsw_inca(:,:) 
504     swtoacs(:,:)  =  ztopsw0_inca(:,:)
505     swsrfas(:,:)  =  zsolsw_inca(:,:) 
506     swsrfcs(:,:)  =  zsolsw0_inca(:,:)
507  else
508     swtoaas(:,:)  = 0.0
509     swtoacs(:,:)  = 0.0
510     swsrfas(:,:)  = 0.0
511     swsrfcs(:,:)  = 0.0
512  endif
513  swsrfas_ad(:) =  solswad_inca(:) 
514  swsrfcs_ad(:) =  solswad0_inca(:)
515  cld_tau(:,:)  =  cldtaupd(:,:)
516  cld_taupi(:,:)=  cldtaupi(:,:)
517  cld_emi(:,:)  =  cldemi(:,:)
518  tops(:)       = topsw(:)
519  tops0(:)      = topsw0(:)
520  topl(:)       = toplw(:)
521  topl0(:)      = toplw0(:)
522
523!RAF
524
525  dforctoaas(:,1)=swtoaas(:,2)-swtoaas(:,3)
526  dforctoacs(:,1)=swtoacs(:,2)-swtoacs(:,3)
527  dforcsrfas(:,1)=swsrfas(:,2)-swsrfas(:,3)
528  dforcsrfcs(:,1)=swsrfcs(:,2)-swsrfcs(:,3)
529  do i=2,naero_grp
530    dforctoaas(:,i)=swtoaas(:,i)-swtoaas(:,1)
531    dforctoacs(:,i)=swtoacs(:,i)-swtoacs(:,1)
532    dforcsrfas(:,i)=swsrfas(:,i)-swsrfas(:,1)
533    dforcsrfcs(:,i)=swsrfcs(:,i)-swsrfcs(:,1)
534  end do       
535  IF (ok_aie) THEN
536  iforctoaas(:)=swtoaas_ai(:)-swtoaas(:,2) 
537  iforcsrfas(:)=swsrfas_ai(:)-swsrfas(:,2) 
538  ELSE
539  iforctoaas(:)=0.
540  iforcsrfas(:)=0.
541  ENDIF
542
543!RAFF Cloud forcing+ impact of aerosols on cloud forcing
544  cforctoa_0(:) = topsw_inca(:,1)-topsw0_inca(:,1)
545  cforcsrf_0(:) = solsw_inca(:,1)-solsw0_inca(:,1)
546  dcforctoa_nat(:) = dforctoaas(:,3)-dforctoacs(:,3)
547  dcforcsrf_nat(:) =  dforcsrfas(:,3)-dforcsrfcs(:,3)
548  dcforctoa_antr(:) = dforctoaas(:,1)-dforctoacs(:,1)
549  dcforcsrf_antr(:) = dforcsrfas(:,1)-dforcsrfcs(:,1)
550
551!RAFF2 Forcing in cloudy regions
552  DO i = 1, PLON       
553    if (cldfract(i).gt.0.) then
554     cRFtoa_nat(i)=(dforctoaas(i,3)-ZCLEAR(i)*dforctoacs(i,3))/cldfract(i)
555     cRFsrf_nat(i)=(dforcsrfas(i,3)-ZCLEAR(i)*dforcsrfcs(i,3))/cldfract(i)
556     cRFtoa_antr(i)=(dforctoaas(i,1)-ZCLEAR(i)*dforctoacs(i,1))/cldfract(i)
557     cRFsrf_antr(i)=(dforcsrfas(i,1)-ZCLEAR(i)*dforcsrfcs(i,1))/cldfract(i)
558    else
559     cRFtoa_nat(i)=0.
560     cRFsrf_nat(i)=0.
561     cRFtoa_antr(i)=0.
562     cRFsrf_antr(i)=0.
563    end if
564  END DO       
565
566#endif
567ENDSUBROUTINE radlwsw_inca
568
569#ifdef AER
570
571SUBROUTINE SW_INCA(kdlon,kflev,PSCT, PRMU0, PFRAC, &
572   PPMB, PDP, &
573   PPSOL, PALBD, PALBP,&
574   PTAVE, PWV, PQS, POZON, PAER,&
575   PCLDSW, PTAU, POMEGA, PCG,&
576   PHEAT, PHEAT0,&
577   PALBPLA,PTOPSW,PSOLSW,PTOPSW0,PSOLSW0,&
578   ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,&
579   tauinca, pizinca, cginca,&
580   PTAUA, POMEGAA,&
581   PTOPSWADINCA,PSOLSWADINCA,&
582   PTOPSWAD0INCA,PSOLSWAD0INCA,&
583   PTOPSWAIINCA,PSOLSWAIINCA,&
584   PTOPSWINCA,PTOPSW0INCA,&
585   PSOLSWINCA,PSOLSW0INCA,&
586   ok_ade, ok_aie,ZCLEAR)
587
588  USE PRINT_INCA
589  USE PARAM_CHEM
590  USE AEROSOL_DIAG
591
592  IMPLICIT NONE
593 
594#include "YOMCST_I.h"
595  !
596  !     ------------------------------------------------------------------
597  !
598  !     PURPOSE.
599  !     --------
600  !
601  !          THIS ROUTINE COMPUTES THE SHORTWAVE RADIATION FLUXES IN TWO
602  !     SPECTRAL INTERVALS FOLLOWING FOUQUART AND BONNEL (1980).
603  !
604  !     METHOD.
605  !     -------
606  !
607  !          1. COMPUTES ABSORBER AMOUNTS                 (SWU)
608  !          2. COMPUTES FLUXES IN 1ST SPECTRAL INTERVAL  (SW1S)
609  !          3. COMPUTES FLUXES IN 2ND SPECTRAL INTERVAL  (SW2S)
610  !
611  !     REFERENCE.
612  !     ----------
613  !
614  !        SEE RADIATION'S PART OF THE ECMWF RESEARCH DEPARTMENT
615  !        DOCUMENTATION, AND FOUQUART AND BONNEL (1980)
616  !
617  !     AUTHOR.
618  !     -------
619  !        JEAN-JACQUES MORCRETTE  *ECMWF*
620  !
621  !     MODIFICATIONS.
622  !     --------------
623  !        ORIGINAL : 89-07-14
624  !        95-01-01   J.-J. MORCRETTE  Direct/Diffuse Albedo
625  !        03-11-27   J. QUAAS Introduce aerosol forcings (based on BOUCHER)
626  !     ------------------------------------------------------------------
627  !
628  !* ARGUMENTS:
629  !
630  REAL*8 PSCT  ! constante solaire (valeur conseillee: 1370)
631 
632  INTEGER, INTENT(in) :: kdlon,kflev
633  REAL*8 PPSOL(KDLON)        ! SURFACE PRESSURE (PA)
634  REAL*8 PDP(KDLON,KFLEV)    ! LAYER THICKNESS (PA)
635  REAL*8 PPMB(KDLON,KFLEV+1) ! HALF-LEVEL PRESSURE (MB)
636 
637  REAL*8 PRMU0(KDLON)  ! COSINE OF ZENITHAL ANGLE
638  REAL*8 PFRAC(KDLON)  ! fraction de la journee
639 
640  REAL*8 PTAVE(KDLON,KFLEV)  ! LAYER TEMPERATURE (K)
641  REAL*8 PWV(KDLON,KFLEV)    ! SPECIFI! HUMIDITY (KG/KG)
642  REAL*8 PQS(KDLON,KFLEV)    ! SATURATED WATER VAPOUR (KG/KG)
643  REAL*8 POZON(KDLON,KFLEV)  ! OZONE CONCENTRATION (KG/KG)
644  REAL*8 PAER(KDLON,KFLEV,5) ! AEROSOLS' OPTICAL THICKNESS
645 
646  REAL*8 PALBD(KDLON,2)  ! albedo du sol (lumiere diffuse)
647  REAL*8 PALBP(KDLON,2)  ! albedo du sol (lumiere parallele)
648 
649  REAL*8 PCLDSW(KDLON,KFLEV)    ! CLOUD FRACTION
650  REAL*8 PTAU(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS
651  REAL*8 PCG(KDLON,2,KFLEV)     ! ASYMETRY FACTOR
652  REAL*8 POMEGA(KDLON,2,KFLEV)  ! SINGLE SCATTERING ALBEDO
653 
654  REAL*8 PHEAT(KDLON,KFLEV) ! SHORTWAVE HEATING (K/DAY)
655  REAL*8 PHEAT0(KDLON,KFLEV)! SHORTWAVE HEATING (K/DAY) clear-sky
656  REAL*8 PALBPLA(KDLON)     ! PLANETARY ALBEDO
657  REAL*8 PTOPSW(KDLON)      ! SHORTWAVE FLUX AT T.O.A.
658  REAL*8 PSOLSW(KDLON)      ! SHORTWAVE FLUX AT SURFACE
659  REAL*8 PTOPSW0(KDLON)     ! SHORTWAVE FLUX AT T.O.A. (CLEAR-SKY)
660  REAL*8 PSOLSW0(KDLON)     ! SHORTWAVE FLUX AT SURFACE (CLEAR-SKY)
661  !
662  !* LOCAL VARIABLES:
663  !
664  REAL*8 ZOZ(KDLON,KFLEV)
665  REAL*8 ZAKI(KDLON,2)     
666  REAL*8 ZCLD(KDLON,KFLEV)
667  REAL*8 ZCLEAR(KDLON) 
668  REAL*8 ZDSIG(KDLON,KFLEV)
669  REAL*8 ZFACT(KDLON)
670  REAL*8 ZFD(KDLON,KFLEV+1)
671  REAL*8 ZFDOWN(KDLON,KFLEV+1)
672  REAL*8 ZFU(KDLON,KFLEV+1)
673  REAL*8 ZFUP(KDLON,KFLEV+1)
674  REAL*8 ZRMU(KDLON)
675  REAL*8 ZSEC(KDLON)
676  REAL*8 ZUD(KDLON,5,KFLEV+1)
677  REAL*8 ZCLDSW0(KDLON,KFLEV)
678 
679  REAL*8 ZFSUP(KDLON,KFLEV+1)
680  REAL*8 ZFSDN(KDLON,KFLEV+1)
681  REAL*8 ZFSUP0(KDLON,KFLEV+1)
682  REAL*8 ZFSDN0(KDLON,KFLEV+1)
683 
684  INTEGER inu, jl, jk, i, k, kpl1
685
686  INTEGER swpas  ! Every swpas steps, sw is calculated
687  PARAMETER(swpas=1)
688 
689  INTEGER itapsw
690  LOGICAL appel1er
691  DATA itapsw /0/
692  DATA appel1er /.TRUE./
693  SAVE itapsw,appel1er
694  !$OMP THREADPRIVATE(appel1er)
695  !$OMP THREADPRIVATE(itapsw)
696  !jq-Introduced for aerosol forcings
697  REAL*8 flag_aer
698  LOGICAL ok_ade, ok_aie    ! use aerosol forcings or not?
699 
700
701  REAL*8 tauinca(kdlon,kflev,naero_grp,2)  ! aerosol optical properties
702  REAL*8 pizinca(kdlon,kflev,naero_grp,2)  ! (see aeropt.F)
703  REAL*8 cginca(kdlon,kflev,naero_grp,2)   ! -"-
704  REAL*8 PTAUA(KDLON,2,KFLEV)    ! CLOUD OPTICAL THICKNESS (pre-industrial value)
705  REAL*8 POMEGAA(KDLON,2,KFLEV)  ! SINGLE SCATTERING ALBEDO
706  REAL*8 PTOPSWADINCA(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
707  REAL*8 PSOLSWADINCA(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
708  REAL*8 PTOPSWAD0INCA(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR)
709  REAL*8 PSOLSWAD0INCA(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR)
710  REAL*8 PTOPSWAIINCA(KDLON)     ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL IND)
711  REAL*8 PSOLSWAIINCA(KDLON)     ! SHORTWAVE FLUX AT SURFACE(+AEROSOL IND)
712  REAL*8 PTOPSWINCA(KDLON,naero_grp)
713  REAL*8 PTOPSW0INCA(KDLON,naero_grp)
714  REAL*8 PSOLSWINCA(KDLON,naero_grp)
715  REAL*8 PSOLSW0INCA(KDLON,naero_grp)
716 
717  !jq - Fluxes including aerosol effects
718  REAL*8,ALLOCATABLE,SAVE :: ZFSUPAD_INCA(:,:)
719  !$OMP THREADPRIVATE(ZFSUPAD_INCA)
720  REAL*8,ALLOCATABLE,SAVE :: ZFSDNAD_INCA(:,:)
721  !$OMP THREADPRIVATE(ZFSDNAD_INCA)
722  !jq - Fluxes including aerosol effects
723  REAL*8,ALLOCATABLE,SAVE :: ZFSUPAD0_INCA(:,:)
724  !$OMP THREADPRIVATE(ZFSUPAD0_INCA)
725  REAL*8,ALLOCATABLE,SAVE :: ZFSDNAD0_INCA(:,:)
726  !$OMP THREADPRIVATE(ZFSDNAD0_INCA)
727  REAL*8,ALLOCATABLE,SAVE :: ZFSUPAI_INCA(:,:)
728  !$OMP THREADPRIVATE(ZFSUPAI_INCA)
729  REAL*8,ALLOCATABLE,SAVE :: ZFSDNAI_INCA(:,:)
730  !$OMP THREADPRIVATE(ZFSDNAI_INCA)
731  REAL*8,ALLOCATABLE,SAVE ::  ZFSUP_INCA(:,:,:)
732  !$OMP THREADPRIVATE(ZFSUP_INCA)
733  REAL*8,ALLOCATABLE,SAVE ::  ZFSDN_INCA(:,:,:)
734  !$OMP THREADPRIVATE(ZFSDN_INCA)
735  REAL*8,ALLOCATABLE,SAVE ::  ZFSUP0_INCA(:,:,:)
736  !$OMP THREADPRIVATE(ZFSUP0_INCA)
737  REAL*8,ALLOCATABLE,SAVE ::  ZFSDN0_INCA(:,:,:)
738  !$OMP THREADPRIVATE(ZFSDN0_INCA)
739!!HEAT
740  REAL*8,ALLOCATABLE,SAVE :: ZFSUP_HEAT(:,:)
741  !$OMP THREADPRIVATE(ZFSUP_HEAT)
742  REAL*8,ALLOCATABLE,SAVE :: ZFSDN_HEAT(:,:)
743  !$OMP THREADPRIVATE(ZFSDN_HEAT)
744  REAL*8,ALLOCATABLE,SAVE :: ZFSUP0_HEAT(:,:)
745  !$OMP THREADPRIVATE(ZFSUP0_HEAT)
746  REAL*8,ALLOCATABLE,SAVE :: ZFSDN0_HEAT(:,:)
747  !$OMP THREADPRIVATE(ZFSDN0_HEAT)
748
749  INTEGER ispec 
750 
751  LOGICAL initialized
752  !rv
753  SAVE flag_aer
754  !$OMP THREADPRIVATE(flag_aer)
755  DATA initialized/.FALSE./
756  SAVE initialized
757  !$OMP THREADPRIVATE(initialized)
758 
759 
760  IF(.NOT.initialized) THEN
761      flag_aer=0.
762      initialized=.TRUE.
763      ALLOCATE(ZFSUPAD_INCA(KDLON,KFLEV+1))
764      ALLOCATE(ZFSDNAD_INCA(KDLON,KFLEV+1))
765      ALLOCATE(ZFSUPAD0_INCA(KDLON,KFLEV+1))
766      ALLOCATE(ZFSDNAD0_INCA(KDLON,KFLEV+1))
767      ALLOCATE(ZFSUPAI_INCA(KDLON,KFLEV+1))
768      ALLOCATE(ZFSDNAI_INCA(KDLON,KFLEV+1))
769      ALLOCATE(ZFSUP_INCA (KDLON,KFLEV+1,naero_grp))
770      ALLOCATE(ZFSDN_INCA (KDLON,KFLEV+1,naero_grp))
771      ALLOCATE(ZFSUP0_INCA(KDLON,KFLEV+1,naero_grp))
772      ALLOCATE(ZFSDN0_INCA(KDLON,KFLEV+1,naero_grp))
773
774      ALLOCATE(ZFSUP_HEAT(KDLON,KFLEV+1))
775      ALLOCATE(ZFSDN_HEAT(KDLON,KFLEV+1))
776      ALLOCATE(ZFSUP0_HEAT(KDLON,KFLEV+1))
777      ALLOCATE(ZFSDN0_HEAT(KDLON,KFLEV+1))
778
779      ZFSUPAD_INCA(:,:)=0.
780      ZFSDNAD_INCA(:,:)=0.
781      ZFSUPAD0_INCA(:,:)=0.
782      ZFSDNAD0_INCA(:,:)=0.
783      ZFSUPAI_INCA(:,:)=0.
784      ZFSDNAI_INCA(:,:)=0.
785      ZFSUP_INCA (:,:,:)=0.
786      ZFSDN_INCA (:,:,:)=0.
787      ZFSUP0_INCA(:,:,:)=0.
788      ZFSDN0_INCA(:,:,:)=0.
789
790      ZFSUP_HEAT(:,:)=0.
791      ZFSDN_HEAT(:,:)=0.
792      ZFSUP0_HEAT(:,:)=0.
793      ZFSDN0_HEAT(:,:)=0.
794  ENDIF
795 
796 
797  IF (appel1er) THEN
798      WRITE(lunout,*) 'SW calling frequency : ', swpas
799      WRITE(lunout,*) "   In general, it should be 1"
800      appel1er = .FALSE.
801  ENDIF
802!     ------------------------------------------------------------------
803  IF (MOD(itapsw,swpas).EQ.0) THEN
804     
805      DO JK = 1 , KFLEV
806        DO JL = 1, KDLON
807          ZCLDSW0(JL,JK) = 0.0
808          ZOZ(JL,JK) = POZON(JL,JK)*46.6968/RG &
809             *PDP(JL,JK)*(101325.0/PPSOL(JL))
810        ENDDO
811      ENDDO
812     
813     
814      ! clear-sky:
815      flag_aer=0.0
816      CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
817         PRMU0,PFRAC,PTAVE,PWV,&
818         ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
819      INU = 1
820      CALL SW1S_LMDAR4(INU,PAER, flag_aer, &
821         tauinca(:,:,1,:), pizinca(:,:,1,:), cginca(:,:,1,:),&
822         PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,&
823         ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
824         ZFD, ZFU)
825      INU = 2
826      CALL SW2S_LMDAR4(INU, PAER, flag_aer, &
827         tauinca(:,:,1,:), pizinca(:,:,1,:), cginca(:,:,1,:),&
828         ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,&
829         ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
830         PWV, PQS,&
831         ZFDOWN, ZFUP)
832      DO JK = 1 , KFLEV+1
833        DO JL = 1, KDLON
834          ZFSUP0(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
835          ZFSDN0(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
836          ZFSUP0_INCA(JL,JK,1) = ZFSUP0(JL,JK) 
837          ZFSDN0_INCA(JL,JK,1) = ZFSDN0(JL,JK)
838        ENDDO
839      ENDDO
840     
841     
842      ! cloudy-sky:
843      flag_aer=0.0
844      CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
845         PRMU0,PFRAC,PTAVE,PWV,&
846         ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
847      INU = 1
848      CALL SW1S_LMDAR4(INU, PAER, flag_aer, &
849         tauinca(:,:,1,:), pizinca(:,:,1,:), cginca(:,:,1,:),&
850         PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
851         ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
852         ZFD, ZFU)
853      INU = 2
854      CALL SW2S_LMDAR4(INU, PAER, flag_aer, &
855         tauinca(:,:,1,:), pizinca(:,:,1,:), cginca(:,:,1,:),&
856         ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
857         ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
858         PWV, PQS,&
859         ZFDOWN, ZFUP)
860     
861      DO JK = 1 , KFLEV+1
862        DO JL = 1, KDLON
863          ZFSUP(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
864          ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
865          ZFSUP_INCA(JL,JK,1) = ZFSUP(JL,JK) 
866          ZFSDN_INCA(JL,JK,1) = ZFSDN(JL,JK)
867        ENDDO
868      ENDDO
869     
870     
871      IF (ok_ade) THEN
872         
873          ! clear sky (Anne 03/07/2007)
874          ! CAS AER (2)
875          flag_aer=1.0
876          CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
877             PRMU0,PFRAC,PTAVE,PWV,&
878             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
879          INU = 1
880          CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
881             tauinca(:,:,2,:), pizinca(:,:,2,:), cginca(:,:,2,:),&
882             PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,&
883             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
884             ZFD, ZFU)
885          INU = 2
886          CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
887             tauinca(:,:,2,:), pizinca(:,:,2,:), cginca(:,:,2,:),&
888             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,&
889             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
890             PWV, PQS,&
891             ZFDOWN, ZFUP)
892         
893          DO JK = 1 , KFLEV+1
894            DO JL = 1, KDLON
895              ZFSUPAD0_INCA(JL,JK) = ZFSUP0(JL,JK) 
896              ZFSDNAD0_INCA(JL,JK) = ZFSDN0(JL,JK) 
897              ZFSUP0(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
898              ZFSDN0(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
899              ZFSUP0_INCA(JL,JK,2) = ZFSUP0(JL,JK) 
900              ZFSDN0_INCA(JL,JK,2) = ZFSDN0(JL,JK)
901            ENDDO
902          ENDDO
903         
904          ! cloudy-sky + aerosol dir OB
905          ! Anne AER
906          flag_aer=1.0
907          CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
908             PRMU0,PFRAC,PTAVE,PWV,&
909             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
910          INU = 1
911          CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
912             tauinca(:,:,2,:), pizinca(:,:,2,:), cginca(:,:,2,:),&
913             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
914             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
915             ZFD, ZFU)
916          INU = 2
917          CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
918             tauinca(:,:,2,:), pizinca(:,:,2,:), cginca(:,:,2,:),&
919             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
920             ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
921             PWV, PQS,&
922             ZFDOWN, ZFUP)
923         
924          DO JK = 1 , KFLEV+1
925            DO JL = 1, KDLON
926              ZFSUPAD_INCA(JL,JK) = ZFSUP(JL,JK) 
927              ZFSDNAD_INCA(JL,JK) = ZFSDN(JL,JK) 
928              ZFSUP(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
929              ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
930              ZFSUP_INCA(JL,JK,2) = ZFSUP(JL,JK) 
931              ZFSDN_INCA(JL,JK,2) = ZFSDN(JL,JK)
932            ENDDO
933          ENDDO
934         
935           !CAS NAT BC SO4 POM DUSS CNO3 FNO3
936          do ispec=3,naero_grp
937             ! clear sky
938             flag_aer=1.0
939             CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,&
940                PRMU0,PFRAC,PTAVE,PWV,&
941                ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
942             INU = 1
943             CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
944                tauinca(:,:,ispec,:), pizinca(:,:,ispec,:), cginca(:,:,ispec,:),&
945                PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,&
946                ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
947                ZFD, ZFU)
948             INU = 2
949             CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
950                tauinca(:,:,ispec,:), pizinca(:,:,ispec,:), cginca(:,:,ispec,:),&
951                ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,&
952                ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
953                PWV, PQS,&
954                ZFDOWN, ZFUP)
955         
956             DO JK = 1 , KFLEV+1
957               DO JL = 1, KDLON
958                 ZFSUP0_INCA(JL,JK,ispec) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
959                 ZFSDN0_INCA(JL,JK,ispec) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
960               ENDDO
961             ENDDO 
962         
963             ! cloudy-sky
964             flag_aer=1.0
965             CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
966                PRMU0,PFRAC,PTAVE,PWV,&
967                ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
968             INU = 1
969             CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
970                tauinca(:,:,ispec,:), pizinca(:,:,ispec,:), cginca(:,:,ispec,:),&
971                PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
972                ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
973                ZFD, ZFU)
974             INU = 2
975             CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
976                tauinca(:,:,ispec,:), pizinca(:,:,ispec,:), cginca(:,:,ispec,:),&
977                ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
978                ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,&
979                PWV, PQS,&
980                ZFDOWN, ZFUP)
981         
982             DO JK = 1 , KFLEV+1
983               DO JL = 1, KDLON
984                 ZFSUP_INCA(JL,JK,ispec) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
985                 ZFSDN_INCA(JL,JK,ispec) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
986               ENDDO
987             ENDDO
988         
989          end do !ispec         
990      ENDIF ! ok_ade
991     
992     
993      IF (ok_aie) THEN
994         
995          !jq   cloudy-sky + aerosol direct + aerosol indirect
996          flag_aer=1.0
997          CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,&
998             PRMU0,PFRAC,PTAVE,PWV,&
999             ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD)
1000          INU = 1
1001          CALL SW1S_LMDAR4(INU, PAER, flag_aer,&
1002             tauinca(:,:,2,:), pizinca(:,:,2,:), cginca(:,:,2,:),&
1003             PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
1004             ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,&
1005             ZFD, ZFU)
1006          INU = 2
1007          CALL SW2S_LMDAR4(INU, PAER, flag_aer,&
1008             tauinca(:,:,2,:), pizinca(:,:,2,:), cginca(:,:,2,:),&
1009             ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,&
1010             ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,&
1011             PWV, PQS,&
1012             ZFDOWN, ZFUP)
1013          DO JK = 1 , KFLEV+1
1014            DO JL = 1, KDLON
1015! attention dans zfsupai et zfdnai on stocke la valeur issue de l'effet direct
1016              ZFSUPAI_INCA(JL,JK) = ZFSUP(JL,JK) 
1017              ZFSDNAI_INCA(JL,JK) = ZFSDN(JL,JK)         
1018              ZFSUP(JL,JK) = (ZFUP(JL,JK)   + ZFU(JL,JK)) * ZFACT(JL)
1019              ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL)
1020              ZFSUP_INCA(JL,JK,2) = ZFSUP(JL,JK) 
1021              ZFSDN_INCA(JL,JK,2) = ZFSDN(JL,JK)
1022            ENDDO
1023          ENDDO
1024      ENDIF ! ok_aie     
1025     
1026     
1027      itapsw = 0
1028  ENDIF
1029  itapsw = itapsw + 1
1030 
1031
1032! HEAT COMPUTATION
1033
1034  DO i = 1, KDLON       
1035  DO k = 1, KFLEV+1
1036    if (feedb .eq. 0) then
1037       ZFSUP_HEAT(i,k)=ZFSUP_INCA(i,k,1)       
1038       ZFSDN_HEAT(i,k)=ZFSDN_INCA(i,k,1)       
1039       ZFSUP0_HEAT(i,k)=ZFSUP0_INCA(i,k,1)             
1040       ZFSDN0_HEAT(i,k)=ZFSDN0_INCA(i,k,1)
1041    elseif (feedb .eq. 1) then
1042        IF ((ok_aie) .and. (ok_ade)) THEN
1043!           ZFSUP_HEAT(i,k)=ZFSUPAI_INCA(i,k)   
1044!           ZFSDN_HEAT(i,k)=ZFSDNAI_INCA(i,k)
1045! correction par rapport au code sw_aeroAR4
1046! il faut prendre la valeur issues de ok_aie
1047! et non pas ok_ade
1048           ZFSUP_HEAT(i,k)=ZFSUP_INCA(i,k,2)   
1049           ZFSDN_HEAT(i,k)=ZFSDN_INCA(i,k,2)
1050           ZFSUP0_HEAT(i,k)=ZFSUP0_INCA(i,k,2)         
1051           ZFSDN0_HEAT(i,k)=ZFSDN0_INCA(i,k,2)
1052        ENDIF
1053        IF ((.not. ok_aie) .and. (ok_ade)) THEN
1054           ZFSUP_HEAT(i,k)=ZFSUP_INCA(i,k,2)   
1055           ZFSDN_HEAT(i,k)=ZFSDN_INCA(i,k,2)
1056           ZFSUP0_HEAT(i,k)=ZFSUP0_INCA(i,k,2)         
1057           ZFSDN0_HEAT(i,k)=ZFSDN0_INCA(i,k,2)
1058        ENDIF
1059        IF ((.not. ok_aie) .and. (.not. ok_ade)) THEN   
1060           ZFSUP_HEAT(i,k)=ZFSUP_INCA(i,k,1)   
1061           ZFSDN_HEAT(i,k)=ZFSDN_INCA(i,k,1)
1062           ZFSUP0_HEAT(i,k)=ZFSUP0_INCA(i,k,1)         
1063           ZFSDN0_HEAT(i,k)=ZFSDN0_INCA(i,k,1)
1064        ENDIF
1065        IF ((ok_aie) .and. (.not. ok_ade)) THEN 
1066!           ZFSUP_HEAT(i,k)=ZFSUPAI_INCA(i,k)   
1067!           ZFSDN_HEAT(i,k)=ZFSDNAI_INCA(i,k)
1068! correction par rapport au code sw_aeroAR4
1069! il faut prendre la valeur issues de ok_aie
1070! et non pas ok_ade
1071           ZFSUP_HEAT(i,k)=ZFSUP_INCA(i,k,2)   
1072           ZFSDN_HEAT(i,k)=ZFSDN_INCA(i,k,2)
1073           ZFSUP0_HEAT(i,k)=ZFSUP0_INCA(i,k,1)         
1074           ZFSDN0_HEAT(i,k)=ZFSDN0_INCA(i,k,1)
1075         END IF
1076    endif !feedb
1077  END DO
1078  END DO
1079
1080!raf: var _HEAT
1081  DO k = 1, KFLEV
1082    kpl1 = k+1
1083    DO i = 1, KDLON
1084    PHEAT(i,k) = -(ZFSUP_HEAT(i,kpl1)-ZFSUP_HEAT(i,k)) &
1085         -(ZFSDN_HEAT(i,k)-ZFSDN_HEAT(i,kpl1))
1086          PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k)
1087    PHEAT0(i,k) = -(ZFSUP0_HEAT(i,kpl1)-ZFSUP0_HEAT(i,k)) &
1088      -(ZFSDN0_HEAT(i,k)-ZFSDN0_HEAT(i,kpl1))
1089      PHEAT0(i,k) = PHEAT0(i,k) * RDAY*RG/RCPD / PDP(i,k)
1090    ENDDO
1091  ENDDO 
1092
1093  DO i = 1, KDLON
1094    PALBPLA(i) = ZFSUP(i,KFLEV+1)/(ZFSDN(i,KFLEV+1)+1.0e-20)
1095    ! clear sky
1096   
1097    PSOLSW0(i) = ZFSDN0_HEAT(i,1) - ZFSUP0_HEAT(i,1)
1098    PTOPSW0(i) = ZFSDN0_HEAT(i,KFLEV+1) - ZFSUP0_HEAT(i,KFLEV+1)
1099   
1100    PSOLSW(i) = ZFSDN_HEAT(i,1) - ZFSUP_HEAT(i,1)
1101    PTOPSW(i) = ZFSDN_HEAT(i,KFLEV+1) - ZFSUP_HEAT(i,KFLEV+1)
1102
1103
1104    PSOLSW0INCA(i,:) = ZFSDN0_INCA(i,1,:) - ZFSUP0_INCA(i,1,:)
1105    PTOPSW0INCA(i,:) = &
1106       ZFSDN0_INCA(i,KFLEV+1,:) - ZFSUP0_INCA(i,KFLEV+1,:)
1107   
1108    PSOLSWINCA(i,:) = ZFSDN_INCA(i,1,:) - ZFSUP_INCA(i,1,:)
1109    PTOPSWINCA(i,:) = &
1110       ZFSDN_INCA(i,KFLEV+1,:) - ZFSUP_INCA(i,KFLEV+1,:)
1111   
1112    PSOLSWADINCA(i) = ZFSDNAD_INCA(i,1) - ZFSUPAD_INCA(i,1)
1113    PTOPSWADINCA(i) = &
1114       ZFSDNAD_INCA(i,KFLEV+1) - ZFSUPAD_INCA(i,KFLEV+1)
1115   
1116    PSOLSWAD0INCA(i) = ZFSDNAD0_INCA(i,1) - ZFSUPAD0_INCA(i,1)
1117    PTOPSWAD0INCA(i) = &
1118       ZFSDNAD0_INCA(i,KFLEV+1) - ZFSUPAD0_INCA(i,KFLEV+1)
1119   
1120    PSOLSWAIINCA(i) = ZFSDNAI_INCA(i,1) - ZFSUPAI_INCA(i,1)
1121    PTOPSWAIINCA(i) = &
1122       ZFSDNAI_INCA(i,KFLEV+1) - ZFSUPAI_INCA(i,KFLEV+1)
1123   
1124  ENDDO
1125 
1126 
1127END SUBROUTINE SW_INCA
1128
1129#endif
Note: See TracBrowser for help on using the repository browser.