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

Contents of /trunk/phylmd/Radlwsw/radlwsw.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 53 - (show annotations)
Fri Oct 7 13:11:58 2011 UTC (12 years, 7 months ago) by guez
Original Path: trunk/libf/phylmd/Radlwsw/radlwsw.f90
File size: 14116 byte(s)


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

  ViewVC Help
Powered by ViewVC 1.1.21