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

  ViewVC Help
Powered by ViewVC 1.1.21