/[lmdze]/trunk/Sources/phylmd/Radlwsw/radlwsw.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/Radlwsw/radlwsw.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/Radlwsw/radlwsw.f revision 52 by guez, Fri Sep 23 12:28:01 2011 UTC trunk/libf/phylmd/Radlwsw/radlwsw.f90 revision 62 by guez, Thu Jul 26 14:37:37 2012 UTC
# Line 1  Line 1 
1  !  module radlwsw_m
2  ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/radlwsw.F,v 1.4 2005/06/06 13:16:33 fairhead Exp $  
3  !    IMPLICIT none
4        SUBROUTINE radlwsw(dist, rmu0, fract,  
5       .                  paprs, pplay,tsol,albedo, alblw, t,q,wo,  contains
6       .                  cldfra, cldemi, cldtaupd,  
7       .                  heat,heat0,cool,cool0,radsol,albpla,    SUBROUTINE radlwsw(dist, rmu0, fract, paprs, pplay, tsol, albedo, alblw, &
8       .                  topsw,toplw,solsw,sollw,         t, q, wo, cldfra, cldemi, cldtaupd, heat, heat0, cool, cool0, radsol, &
9       .                  sollwdown,         albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, toplw0, solsw0, &
10       .                  topsw0,toplw0,solsw0,sollw0,         sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, swup0, swup, ok_ade, &
11       .                  lwdn0, lwdn, lwup0, lwup,         ok_aie, tau_ae, piz_ae, cg_ae, topswad, solswad, cldtaupi, topswai, &
12       .                  swdn0, swdn, swup0, swup,         solswai)
13       .                  ok_ade, ok_aie,  
14       .                  tau_ae, piz_ae, cg_ae,      ! From LMDZ4/libf/phylmd/radlwsw.F, version 1.4 2005/06/06 13:16:33
15       .                  topswad, solswad,      ! Author: Z. X. Li (LMD/CNRS) date: 1996/07/19
16       .                  cldtaupi, topswai, solswai)      ! Objet : interface entre le modèle et les rayonnements
17  c            ! Rayonnements solaire et infrarouge
18        use dimphy  
19        use clesphys      USE dimphy, ONLY: klev, klon
20        use SUPHEC_M      USE clesphys, ONLY: bug_ozone, solaire
21        use raddim, only: kflev, kdlon      USE suphec_m, ONLY: rg
22        use yoethf_m      USE raddim, ONLY: kdlon
23        IMPLICIT none      USE yoethf_m, ONLY: rvtmp2
24  c======================================================================      use sw_m, only: sw
25  c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719          
26  c Objet: interface entre le modele et les rayonnements      ! Arguments:
27  c Arguments:      ! dist-----input-R- distance astronomique terre-soleil
28  c dist-----input-R- distance astronomique terre-soleil      ! rmu0-----input-R- cosinus de l'angle zenithal
29  c rmu0-----input-R- cosinus de l'angle zenithal      ! fract----input-R- duree d'ensoleillement normalisee
30  c fract----input-R- duree d'ensoleillement normalisee      ! co2_ppm--input-R- concentration du gaz carbonique (en ppm)
31  c co2_ppm--input-R- concentration du gaz carbonique (en ppm)      ! solaire--input-R- constante solaire (W/m**2)
32  c solaire--input-R- constante solaire (W/m**2)      ! paprs----input-R- pression a inter-couche (Pa)
33  c paprs----input-R- pression a inter-couche (Pa)      ! pplay----input-R- pression au milieu de couche (Pa)
34  c pplay----input-R- pression au milieu de couche (Pa)      ! tsol-----input-R- temperature du sol (en K)
35  c tsol-----input-R- temperature du sol (en K)      ! albedo---input-R- albedo du sol (entre 0 et 1)
36  c albedo---input-R- albedo du sol (entre 0 et 1)      ! t--------input-R- temperature (K)
37  c t--------input-R- temperature (K)      ! q--------input-R- vapeur d'eau (en kg/kg)
38  c q--------input-R- vapeur d'eau (en kg/kg)      ! wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505
39  c wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505      ! cldfra---input-R- fraction nuageuse (entre 0 et 1)
40  c cldfra---input-R- fraction nuageuse (entre 0 et 1)      ! cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
41  c cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)      ! cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
42  c cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)      ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
43  c ok_ade---input-L- apply the Aerosol Direct Effect or not?      ! cldtaupi-input-R- epaisseur optique des nuages dans le visible
44  c ok_aie---input-L- apply the Aerosol Indirect Effect or not?      !                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
45  c tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)      !                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
46  c cldtaupi-input-R- epaisseur optique des nuages dans le visible      !                   it is needed for the diagnostics of the aerosol indirect radiative forcing      
47  c                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller  
48  c                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd      ! cool-----output-R- refroidissement dans l'IR (K/jour)
49  c                   it is needed for the diagnostics of the aerosol indirect radiative forcing            ! radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
50  c      ! albpla---output-R- albedo planetaire (entre 0 et 1)
51  c heat-----output-R- echauffement atmospherique (visible) (K/jour)      ! topsw----output-R- flux solaire net au sommet de l'atm.
52  c cool-----output-R- refroidissement dans l'IR (K/jour)      ! solsw----output-R- flux solaire net a la surface
53  c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)      ! sollw----output-R- ray. IR montant a la surface
54  c albpla---output-R- albedo planetaire (entre 0 et 1)      ! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
55  c topsw----output-R- flux solaire net au sommet de l'atm.      ! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
56  c toplw----output-R- ray. IR montant au sommet de l'atmosphere      ! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
57  c solsw----output-R- flux solaire net a la surface      ! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
58  c sollw----output-R- ray. IR montant a la surface  
59  c solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)      ! ATTENTION: swai and swad have to be interpreted in the following manner:
60  c topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)      ! ok_ade = F & ok_aie = F -both are zero
61  c solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)      ! ok_ade = T & ok_aie = F -aerosol direct forcing is F_{AD} = topsw-topswad
62  c topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)      !                        indirect is zero
63  c      ! ok_ade = F & ok_aie = T -aerosol indirect forcing is F_{AI} = topsw-topswai
64  c ATTENTION: swai and swad have to be interpreted in the following manner:      !                        direct is zero
65  c ---------      ! ok_ade = T & ok_aie = T -aerosol indirect forcing is F_{AI} = topsw-topswai
66  c ok_ade=F & ok_aie=F -both are zero      !                        aerosol direct forcing is F_{AD} = topswai-topswad
67  c ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad  
68  c                        indirect is zero      real rmu0(klon), fract(klon), dist
69  c ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai  
70  c                        direct is zero      real, intent(in):: paprs(klon, klev+1)
71  c ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai      real, intent(in):: pplay(klon, klev)
72  c                        aerosol direct forcing is F_{AD} = topswai-topswad      real albedo(klon), alblw(klon), tsol(klon)
73  c      real, intent(in):: t(klon, klev)
74              real q(klon, klev)
75  c======================================================================      real, intent(in):: wo(klon, klev)
76  c      real cldfra(klon, klev), cldemi(klon, klev), cldtaupd(klon, klev)
77        real rmu0(klon), fract(klon), dist  
78  cIM   real co2_ppm      real, intent(out):: heat(klon, klev)
79  cIM   real solaire      ! échauffement atmosphérique (visible) (K/jour)
80  c  
81        real, intent(in):: paprs(klon,klev+1)      real cool(klon, klev)
82        real, intent(in):: pplay(klon,klev)      real heat0(klon, klev), cool0(klon, klev)
83        real albedo(klon), alblw(klon), tsol(klon)      real radsol(klon), topsw(klon)
84        real, intent(in):: t(klon,klev)  
85        real q(klon,klev)      real, intent(out):: toplw(klon)
86        real, intent(in):: wo(klon,klev)      ! rayonnement infrarouge montant au sommet de l'atmosphère
87        real cldfra(klon,klev), cldemi(klon,klev), cldtaupd(klon,klev)  
88        real heat(klon,klev), cool(klon,klev)      real solsw(klon), sollw(klon), albpla(klon)
89        real heat0(klon,klev), cool0(klon,klev)      real topsw0(klon), solsw0(klon), sollw0(klon)
90        real radsol(klon), topsw(klon), toplw(klon)      real, intent(out):: toplw0(klon)
91        real solsw(klon), sollw(klon), albpla(klon)      real sollwdown(klon)
92        real topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)      !IM output 3D
93        real sollwdown(klon)      DOUBLE PRECISION ZFSUP(KDLON, KLEV+1)
94  cIM output 3D      DOUBLE PRECISION ZFSDN(KDLON, KLEV+1)
95        REAL*8 ZFSUP(KDLON,KFLEV+1)      DOUBLE PRECISION ZFSUP0(KDLON, KLEV+1)
96        REAL*8 ZFSDN(KDLON,KFLEV+1)      DOUBLE PRECISION ZFSDN0(KDLON, KLEV+1)
97        REAL*8 ZFSUP0(KDLON,KFLEV+1)  
98        REAL*8 ZFSDN0(KDLON,KFLEV+1)      DOUBLE PRECISION ZFLUP(KDLON, KLEV+1)
99  c      DOUBLE PRECISION ZFLDN(KDLON, KLEV+1)
100        REAL*8 ZFLUP(KDLON,KFLEV+1)      DOUBLE PRECISION ZFLUP0(KDLON, KLEV+1)
101        REAL*8 ZFLDN(KDLON,KFLEV+1)      DOUBLE PRECISION ZFLDN0(KDLON, KLEV+1)
102        REAL*8 ZFLUP0(KDLON,KFLEV+1)  
103        REAL*8 ZFLDN0(KDLON,KFLEV+1)      DOUBLE PRECISION zx_alpha1, zx_alpha2
104  c  
105        REAL*8 zx_alpha1, zx_alpha2      INTEGER k, kk, i, iof, nb_gr
106  c      EXTERNAL lw
107  c  
108        INTEGER k, kk, i, j, iof, nb_gr      DOUBLE PRECISION PSCT
109        EXTERNAL lw, sw  
110  c      DOUBLE PRECISION PALBD(kdlon, 2), PALBP(kdlon, 2)
111  cIM ctes ds clesphys.h  REAL*8 RCO2, RCH4, RN2O, RCFC11, RCFC12      DOUBLE PRECISION PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)
112        REAL*8 PSCT      DOUBLE PRECISION PPSOL(kdlon), PDP(kdlon, klev)
113  c      DOUBLE PRECISION PTL(kdlon, klev+1), PPMB(kdlon, klev+1)
114        REAL*8 PALBD(kdlon,2), PALBP(kdlon,2)      DOUBLE PRECISION PTAVE(kdlon, klev)
115        REAL*8 PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon)      DOUBLE PRECISION PWV(kdlon, klev), PQS(kdlon, klev), POZON(kdlon, klev)
116        REAL*8 PPSOL(kdlon), PDP(kdlon,klev)      DOUBLE PRECISION PAER(kdlon, klev, 5)
117        REAL*8 PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1)      DOUBLE PRECISION PCLDLD(kdlon, klev)
118        REAL*8 PTAVE(kdlon,kflev)      DOUBLE PRECISION PCLDLU(kdlon, klev)
119        REAL*8 PWV(kdlon,kflev), PQS(kdlon,kflev), POZON(kdlon,kflev)      DOUBLE PRECISION PCLDSW(kdlon, klev)
120        REAL*8 PAER(kdlon,kflev,5)      DOUBLE PRECISION PTAU(kdlon, 2, klev)
121        REAL*8 PCLDLD(kdlon,kflev)      DOUBLE PRECISION POMEGA(kdlon, 2, klev)
122        REAL*8 PCLDLU(kdlon,kflev)      DOUBLE PRECISION PCG(kdlon, 2, klev)
123        REAL*8 PCLDSW(kdlon,kflev)  
124        REAL*8 PTAU(kdlon,2,kflev)      DOUBLE PRECISION zfract(kdlon), zrmu0(kdlon), zdist
125        REAL*8 POMEGA(kdlon,2,kflev)  
126        REAL*8 PCG(kdlon,2,kflev)      DOUBLE PRECISION zheat(kdlon, klev), zcool(kdlon, klev)
127  c      DOUBLE PRECISION zheat0(kdlon, klev), zcool0(kdlon, klev)
128        REAL*8 zfract(kdlon), zrmu0(kdlon), zdist      DOUBLE PRECISION ztopsw(kdlon), ztoplw(kdlon)
129  c      DOUBLE PRECISION zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)
130        REAL*8 zheat(kdlon,kflev), zcool(kdlon,kflev)      DOUBLE PRECISION zsollwdown(kdlon)
131        REAL*8 zheat0(kdlon,kflev), zcool0(kdlon,kflev)  
132        REAL*8 ztopsw(kdlon), ztoplw(kdlon)      DOUBLE PRECISION ztopsw0(kdlon), ztoplw0(kdlon)
133        REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon)      DOUBLE PRECISION zsolsw0(kdlon), zsollw0(kdlon)
134  cIM      DOUBLE PRECISION zznormcp
135        REAL*8 zsollwdown(kdlon)      !IM output 3D: SWup, SWdn, LWup, LWdn
136  c      REAL swdn(klon, klev+1), swdn0(klon, klev+1)
137        REAL*8 ztopsw0(kdlon), ztoplw0(kdlon)      REAL swup(klon, klev+1), swup0(klon, klev+1)
138        REAL*8 zsolsw0(kdlon), zsollw0(kdlon)      REAL lwdn(klon, klev+1), lwdn0(klon, klev+1)
139        REAL*8 zznormcp      REAL lwup(klon, klev+1), lwup0(klon, klev+1)
140  cIM output 3D : SWup, SWdn, LWup, LWdn  
141        REAL swdn(klon,kflev+1),swdn0(klon,kflev+1)      !jq the following quantities are needed for the aerosol radiative forcings
142        REAL swup(klon,kflev+1),swup0(klon,kflev+1)  
143        REAL lwdn(klon,kflev+1),lwdn0(klon,kflev+1)      real topswad(klon), solswad(klon)
144        REAL lwup(klon,kflev+1),lwup0(klon,kflev+1)      ! output: aerosol direct forcing at TOA and surface
145  c-OB  
146  cjq the following quantities are needed for the aerosol radiative forcings      real topswai(klon), solswai(klon)
147        ! output: aerosol indirect forcing atTOA and surface
148        real topswad(klon), solswad(klon) ! output: aerosol direct forcing at TOA and surface  
149        real topswai(klon), solswai(klon) ! output: aerosol indirect forcing atTOA and surface      real tau_ae(klon, klev, 2), piz_ae(klon, klev, 2), cg_ae(klon, klev, 2)
150        real tau_ae(klon,klev,2), piz_ae(klon,klev,2), cg_ae(klon,klev,2) ! aerosol optical properties (see aeropt.F)      ! aerosol optical properties (see aeropt.F)
151        real cldtaupi(klon,klev)  ! cloud optical thickness for pre-industrial aerosol concentrations  
152                                  ! (i.e., with a smaller droplet concentrationand thus larger droplet radii)      real cldtaupi(klon, klev)
153        logical ok_ade, ok_aie    ! switches whether to use aerosol direct (indirect) effects or not      ! cloud optical thickness for pre-industrial aerosol concentrations
154        real*8 tauae(kdlon,kflev,2) ! aer opt properties      ! (i.e., with a smaller droplet concentrationand thus larger droplet radii)
155        real*8 pizae(kdlon,kflev,2)  
156        real*8 cgae(kdlon,kflev,2)      logical ok_ade, ok_aie
157        REAL*8 PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use      ! switches whether to use aerosol direct (indirect) effects or not
158        REAL*8 POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo      ! ok_ade---input-L- apply the Aerosol Direct Effect or not?
159        REAL*8 ztopswad(kdlon), zsolswad(kdlon) ! Aerosol direct forcing at TOAand surface      ! ok_aie---input-L- apply the Aerosol Indirect Effect or not?
160        REAL*8 ztopswai(kdlon), zsolswai(kdlon) ! dito, indirect  
161  cjq-end      double precision tauae(kdlon, klev, 2) ! aer opt properties
162  !rv      double precision pizae(kdlon, klev, 2)
163        tauae(:,:,:)=0.      double precision cgae(kdlon, klev, 2)
164        pizae(:,:,:)=0.  
165        cgae(:,:,:)=0.      DOUBLE PRECISION PTAUA(kdlon, 2, klev)
166  !rv      ! present-day value of cloud opt thickness (PTAU is pre-industrial
167              ! value), local use
168  c  
169  c-------------------------------------------      DOUBLE PRECISION POMEGAA(kdlon, 2, klev) ! dito for single scatt albedo
170        nb_gr = klon / kdlon  
171        IF (nb_gr*kdlon .NE. klon) THEN      DOUBLE PRECISION ztopswad(kdlon), zsolswad(kdlon)
172           PRINT*, "kdlon mauvais:", klon, kdlon, nb_gr      ! Aerosol direct forcing at TOAand surface
173           stop 1  
174        ENDIF      DOUBLE PRECISION ztopswai(kdlon), zsolswai(kdlon) ! dito, indirect
175        IF (kflev .NE. klev) THEN  
176            PRINT*, "kflev differe de klev, kflev, klev"      !----------------------------------------------------------------------
177            stop 1  
178        ENDIF      tauae = 0.
179  c-------------------------------------------      pizae = 0.
180        DO k = 1, klev      cgae = 0.
181        DO i = 1, klon  
182           heat(i,k)=0.      nb_gr = klon / kdlon
183           cool(i,k)=0.      IF (nb_gr * kdlon /= klon) THEN
184           heat0(i,k)=0.         PRINT *, "kdlon mauvais :", klon, kdlon, nb_gr
185           cool0(i,k)=0.         stop 1
186        ENDDO      ENDIF
187        ENDDO      
188  c      heat = 0.
189        zdist = dist      cool = 0.
190  c      heat0 = 0.
191  cIM anciennes valeurs      cool0 = 0.
192  c     RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97      zdist = dist
193  c      PSCT = solaire / zdist / zdist
194  cIM : on met RCO2, RCH4, RN2O, RCFC11 et RCFC12 dans clesphys.h /lecture ds conf_phys.F90  
195  c     RCH4 = 1.65E-06* 16.043/28.97      loop_iof: DO iof = 0, klon - kdlon, kdlon
196  c     RN2O = 306.E-09* 44.013/28.97         DO i = 1, kdlon
197  c     RCFC11 = 280.E-12* 137.3686/28.97            zfract(i) = fract(iof+i)
198  c     RCFC12 = 484.E-12* 120.9140/28.97            zrmu0(i) = rmu0(iof+i)
199  cIM anciennes valeurs            PALBD(i, 1) = albedo(iof+i)
200  c     RCH4 = 1.72E-06* 16.043/28.97            PALBD(i, 2) = alblw(iof+i)
201  c     RN2O = 310.E-09* 44.013/28.97            PALBP(i, 1) = albedo(iof+i)
202  c            PALBP(i, 2) = alblw(iof+i)
203  c     PRINT*,'IMradlwsw : solaire, co2= ', solaire, co2_ppm            ! cf. JLD pour etre en accord avec ORCHIDEE il faut mettre
204        PSCT = solaire/zdist/zdist            ! PEMIS(i) = 0.96
205  c            PEMIS(i) = 1.0
206        DO 99999 j = 1, nb_gr            PVIEW(i) = 1.66
207        iof = kdlon*(j-1)            PPSOL(i) = paprs(iof+i, 1)
208  c            zx_alpha1 = (paprs(iof+i, 1)-pplay(iof+i, 2))  &
209        DO i = 1, kdlon                 / (pplay(iof+i, 1)-pplay(iof+i, 2))
210           zfract(i) = fract(iof+i)            zx_alpha2 = 1.0 - zx_alpha1
211           zrmu0(i) = rmu0(iof+i)            PTL(i, 1) = t(iof+i, 1) * zx_alpha1 + t(iof+i, 2) * zx_alpha2
212           PALBD(i,1) = albedo(iof+i)            PTL(i, klev+1) = t(iof+i, klev)
213  !         PALBD(i,2) = albedo(iof+i)            PDT0(i) = tsol(iof+i) - PTL(i, 1)
214           PALBD(i,2) = alblw(iof+i)         ENDDO
215           PALBP(i,1) = albedo(iof+i)         DO k = 2, klev
216  !         PALBP(i,2) = albedo(iof+i)            DO i = 1, kdlon
217           PALBP(i,2) = alblw(iof+i)               PTL(i, k) = (t(iof+i, k)+t(iof+i, k-1))*0.5
218  cIM cf. JLD pour etre en accord avec ORCHIDEE il faut mettre PEMIS(i) = 0.96            ENDDO
219           PEMIS(i) = 1.0         ENDDO
220           PVIEW(i) = 1.66         DO k = 1, klev
221           PPSOL(i) = paprs(iof+i,1)            DO i = 1, kdlon
222           zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))               PDP(i, k) = paprs(iof+i, k)-paprs(iof+i, k+1)
223       .             / (pplay(iof+i,1)-pplay(iof+i,2))               PTAVE(i, k) = t(iof+i, k)
224           zx_alpha2 = 1.0 - zx_alpha1               PWV(i, k) = MAX (q(iof+i, k), 1.0e-12)
225           PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2               PQS(i, k) = PWV(i, k)
226           PTL(i,klev+1) = t(iof+i,klev)               ! wo:    cm.atm (epaisseur en cm dans la situation standard)
227           PDT0(i) = tsol(iof+i) - PTL(i,1)               ! POZON: kg/kg
228        ENDDO               IF (bug_ozone) then
229        DO k = 2, kflev                  POZON(i, k) = MAX(wo(iof+i, k), 1.0e-12)*RG/46.6968 &
230        DO i = 1, kdlon                       /(paprs(iof+i, k)-paprs(iof+i, k+1)) &
231           PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5                       *(paprs(iof+i, 1)/101325.0)
232        ENDDO               ELSE
233        ENDDO                  ! le calcul qui suit est maintenant fait dans ozonecm (MPL)
234        DO k = 1, kflev                  POZON(i, k) = wo(i, k)
235        DO i = 1, kdlon               ENDIF
236           PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)               PCLDLD(i, k) = cldfra(iof+i, k)*cldemi(iof+i, k)
237           PTAVE(i,k) = t(iof+i,k)               PCLDLU(i, k) = cldfra(iof+i, k)*cldemi(iof+i, k)
238           PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)               PCLDSW(i, k) = cldfra(iof+i, k)
239           PQS(i,k) = PWV(i,k)               PTAU(i, 1, k) = MAX(cldtaupi(iof+i, k), 1.0e-05)
240  c wo:    cm.atm (epaisseur en cm dans la situation standard)               ! (1e-12 serait instable)
241  c POZON: kg/kg               PTAU(i, 2, k) = MAX(cldtaupi(iof+i, k), 1.0e-05)
242           IF (bug_ozone) then               ! (pour 32-bit machines)
243             POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968               POMEGA(i, 1, k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i, 1, k))
244       .               /(paprs(iof+i,k)-paprs(iof+i,k+1))               POMEGA(i, 2, k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i, 2, k))
245       .               *(paprs(iof+i,1)/101325.0)               PCG(i, 1, k) = 0.865
246           ELSE               PCG(i, 2, k) = 0.910
247  c le calcul qui suit est maintenant fait dans ozonecm (MPL)  
248             POZON(i,k) = wo(i,k)               ! Introduced for aerosol indirect forcings.  The
249           ENDIF               ! following values use the cloud optical thickness
250           PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)               ! calculated from present-day aerosol concentrations
251           PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)               ! whereas the quantities without the "A" at the end are
252           PCLDSW(i,k) = cldfra(iof+i,k)               ! for pre-industial (natural-only) aerosol concentrations
253           PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable               PTAUA(i, 1, k) = MAX(cldtaupd(iof+i, k), 1.0e-05)
254           PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines               ! (1e-12 serait instable)
255           POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))               PTAUA(i, 2, k) = MAX(cldtaupd(iof+i, k), 1.0e-05)
256           POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))               ! (pour 32-bit machines)
257           PCG(i,1,k) = 0.865               POMEGAA(i, 1, k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i, 1, k))
258           PCG(i,2,k) = 0.910               POMEGAA(i, 2, k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i, 2, k))
259  c-OB               !jq-end
260  cjq Introduced for aerosol indirect forcings.            ENDDO
261  cjq The following values use the cloud optical thickness calculated from         ENDDO
262  cjq present-day aerosol concentrations whereas the quantities without the  
263  cjq "A" at the end are for pre-industial (natural-only) aerosol concentrations         DO k = 1, klev+1
264  cjq            DO i = 1, kdlon
265           PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable               PPMB(i, k) = paprs(iof+i, k)/100.0
266           PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines            ENDDO
267           POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))         ENDDO
268           POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))  
269  cjq-end         DO kk = 1, 5
270        ENDDO            DO k = 1, klev
271        ENDDO               DO i = 1, kdlon
272  c                  PAER(i, k, kk) = 1.0E-15
273        DO k = 1, kflev+1               ENDDO
274        DO i = 1, kdlon            ENDDO
275           PPMB(i,k) = paprs(iof+i,k)/100.0         ENDDO
276        ENDDO  
277        ENDDO         DO k = 1, klev
278  c            DO i = 1, kdlon
279        DO kk = 1, 5               tauae(i, k, 1) = tau_ae(iof+i, k, 1)
280        DO k = 1, kflev               pizae(i, k, 1) = piz_ae(iof+i, k, 1)
281        DO i = 1, kdlon               cgae(i, k, 1) =cg_ae(iof+i, k, 1)
282           PAER(i,k,kk) = 1.0E-15               tauae(i, k, 2) = tau_ae(iof+i, k, 2)
283        ENDDO               pizae(i, k, 2) = piz_ae(iof+i, k, 2)
284        ENDDO               cgae(i, k, 2) =cg_ae(iof+i, k, 2)
285        ENDDO            ENDDO
286  c-OB         ENDDO
287        DO k = 1, kflev  
288        DO i = 1, kdlon         CALL LW(PPMB, PDP, PPSOL, PDT0, PEMIS, PTL, PTAVE, PWV, POZON, PAER, &
289          tauae(i,k,1)=tau_ae(iof+i,k,1)              PCLDLD, PCLDLU, PVIEW, zcool, zcool0, ztoplw, zsollw, ztoplw0, &
290          pizae(i,k,1)=piz_ae(iof+i,k,1)              zsollw0, zsollwdown, ZFLUP, ZFLDN, ZFLUP0, ZFLDN0)
291          cgae(i,k,1) =cg_ae(iof+i,k,1)         CALL SW(PSCT, zrmu0, zfract, PPMB, PDP, PPSOL, PALBD, PALBP, PTAVE, &
292          tauae(i,k,2)=tau_ae(iof+i,k,2)              PWV, PQS, POZON, PAER, PCLDSW, PTAU, POMEGA, PCG, zheat, zheat0, &
293          pizae(i,k,2)=piz_ae(iof+i,k,2)              zalbpla, ztopsw, zsolsw, ztopsw0, zsolsw0, ZFSUP, ZFSDN, ZFSUP0, &
294          cgae(i,k,2) =cg_ae(iof+i,k,2)              ZFSDN0, tauae, pizae, cgae, PTAUA, POMEGAA, ztopswad, zsolswad, &
295        ENDDO              ztopswai, zsolswai, ok_ade, ok_aie)
296        ENDDO  
297  c         DO i = 1, kdlon
298  c======================================================================            radsol(iof+i) = zsolsw(i) + zsollw(i)
299  cIM ctes ds clesphys.h   CALL LW(RCO2,RCH4,RN2O,RCFC11,RCFC12,            topsw(iof+i) = ztopsw(i)
300        CALL LW(            toplw(iof+i) = ztoplw(i)
301       .        PPMB, PDP,            solsw(iof+i) = zsolsw(i)
302       .        PPSOL,PDT0,PEMIS,            sollw(iof+i) = zsollw(i)
303       .        PTL, PTAVE, PWV, POZON, PAER,            sollwdown(iof+i) = zsollwdown(i)
304       .        PCLDLD,PCLDLU,  
305       .        PVIEW,            DO k = 1, klev+1
306       .        zcool, zcool0,               lwdn0 ( iof+i, k)   = ZFLDN0 ( i, k)
307       .        ztoplw,zsollw,ztoplw0,zsollw0,               lwdn  ( iof+i, k)   = ZFLDN  ( i, k)
308       .        zsollwdown,               lwup0 ( iof+i, k)   = ZFLUP0 ( i, k)
309       .        ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)               lwup  ( iof+i, k)   = ZFLUP  ( i, k)
310  cIM ctes ds clesphys.h   CALL SW(PSCT, RCO2, zrmu0, zfract,            ENDDO
311        CALL SW(PSCT, zrmu0, zfract,  
312       S        PPMB, PDP,            topsw0(iof+i) = ztopsw0(i)
313       S        PPSOL, PALBD, PALBP,            toplw0(iof+i) = ztoplw0(i)
314       S        PTAVE, PWV, PQS, POZON, PAER,            solsw0(iof+i) = zsolsw0(i)
315       S        PCLDSW, PTAU, POMEGA, PCG,            sollw0(iof+i) = zsollw0(i)
316       S        zheat, zheat0,            albpla(iof+i) = zalbpla(i)
317       S        zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,  
318       S        ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,            DO k = 1, klev+1
319       S        tauae, pizae, cgae, ! aerosol optical properties               swdn0 ( iof+i, k)   = ZFSDN0 ( i, k)
320       s        PTAUA, POMEGAA,               swdn  ( iof+i, k)   = ZFSDN  ( i, k)
321       s        ztopswad,zsolswad,ztopswai,zsolswai, ! diagnosed aerosol forcing               swup0 ( iof+i, k)   = ZFSUP0 ( i, k)
322       J        ok_ade, ok_aie) ! apply aerosol effects or not?               swup  ( iof+i, k)   = ZFSUP  ( i, k)
323              ENDDO
324  c======================================================================         ENDDO
325        DO i = 1, kdlon         ! transform the aerosol forcings, if they have to be calculated
326           radsol(iof+i) = zsolsw(i) + zsollw(i)         IF (ok_ade) THEN
327           topsw(iof+i) = ztopsw(i)            DO i = 1, kdlon
328           toplw(iof+i) = ztoplw(i)               topswad(iof+i) = ztopswad(i)
329           solsw(iof+i) = zsolsw(i)               solswad(iof+i) = zsolswad(i)
330           sollw(iof+i) = zsollw(i)            ENDDO
331           sollwdown(iof+i) = zsollwdown(i)         ELSE
332  cIM            DO i = 1, kdlon
333           DO k = 1, kflev+1               topswad(iof+i) = 0.0
334           lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)               solswad(iof+i) = 0.0
335           lwdn  ( iof+i,k)   = ZFLDN  ( i,k)            ENDDO
336           lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)         ENDIF
337           lwup  ( iof+i,k)   = ZFLUP  ( i,k)         IF (ok_aie) THEN
338           ENDDO            DO i = 1, kdlon
339  c               topswai(iof+i) = ztopswai(i)
340           topsw0(iof+i) = ztopsw0(i)               solswai(iof+i) = zsolswai(i)
341           toplw0(iof+i) = ztoplw0(i)            ENDDO
342           solsw0(iof+i) = zsolsw0(i)         ELSE
343           sollw0(iof+i) = zsollw0(i)            DO i = 1, kdlon
344           albpla(iof+i) = zalbpla(i)               topswai(iof+i) = 0.0
345  cIM               solswai(iof+i) = 0.0
346           DO k = 1, kflev+1            ENDDO
347           swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)         ENDIF
348           swdn  ( iof+i,k)   = ZFSDN  ( i,k)  
349           swup0 ( iof+i,k)   = ZFSUP0 ( i,k)         DO k = 1, klev
350           swup  ( iof+i,k)   = ZFSUP  ( i,k)            DO i = 1, kdlon
351           ENDDO !k=1, kflev+1               ! scale factor to take into account the difference
352        ENDDO               ! between dry air and water vapour specific heat capacity
353  cjq-transform the aerosol forcings, if they have               zznormcp = 1. + RVTMP2 * PWV(i, k)
354  cjq to be calculated               heat(iof+i, k) = zheat(i, k) / zznormcp
355        IF (ok_ade) THEN               cool(iof+i, k) = zcool(i, k)/zznormcp
356        DO i = 1, kdlon               heat0(iof+i, k) = zheat0(i, k)/zznormcp
357           topswad(iof+i) = ztopswad(i)               cool0(iof+i, k) = zcool0(i, k)/zznormcp
358           solswad(iof+i) = zsolswad(i)            ENDDO
359        ENDDO         ENDDO
360        ELSE      end DO loop_iof
361        DO i = 1, kdlon  
362           topswad(iof+i) = 0.0    END SUBROUTINE radlwsw
363           solswad(iof+i) = 0.0  
364        ENDDO  end module radlwsw_m
       ENDIF  
       IF (ok_aie) THEN  
       DO i = 1, kdlon  
          topswai(iof+i) = ztopswai(i)  
          solswai(iof+i) = zsolswai(i)  
       ENDDO  
       ELSE  
       DO i = 1, kdlon  
          topswai(iof+i) = 0.0  
          solswai(iof+i) = 0.0  
       ENDDO  
       ENDIF  
 cjq-end  
       DO k = 1, kflev  
 c      DO i = 1, kdlon  
 c         heat(iof+i,k) = zheat(i,k)  
 c         cool(iof+i,k) = zcool(i,k)  
 c         heat0(iof+i,k) = zheat0(i,k)  
 c         cool0(iof+i,k) = zcool0(i,k)  
 c      ENDDO  
       DO i = 1, kdlon  
 C        scale factor to take into account the difference between  
 C        dry air and watter vapour scpecific heat capacity  
          zznormcp=1.0+RVTMP2*PWV(i,k)  
          heat(iof+i,k) = zheat(i,k)/zznormcp  
          cool(iof+i,k) = zcool(i,k)/zznormcp  
          heat0(iof+i,k) = zheat0(i,k)/zznormcp  
          cool0(iof+i,k) = zcool0(i,k)/zznormcp  
       ENDDO  
       ENDDO  
 c  
 99999 CONTINUE  
       RETURN  
       END  

Legend:
Removed from v.52  
changed lines
  Added in v.62

  ViewVC Help
Powered by ViewVC 1.1.21