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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 38 - (show annotations)
Thu Jan 6 17:52:19 2011 UTC (13 years, 4 months ago) by guez
Original Path: trunk/libf/phylmd/Radlwsw/radlwsw.f
File size: 14815 byte(s)
Extracted ASCII art from "inigeom" into a separate text file in the
documentation.

"test_disvert" now creates a separate file for layer thicknesses.

Moved variables from module "yomcst" to module "suphec_m" because this
is where those variables are defined. Kept in "yomcst" only parameters
of Earth orbit. Gave the attribute "parameter" to some variables of
module "suphec_m".

Variables of module "yoethf" were defined in procedure "suphec". Moved
these definitions to a new procedure "yoethf" in module "yoethf_m".

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

  ViewVC Help
Powered by ViewVC 1.1.21