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

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

  ViewVC Help
Powered by ViewVC 1.1.21