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

  ViewVC Help
Powered by ViewVC 1.1.21