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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 52 - (show annotations)
Fri Sep 23 12:28:01 2011 UTC (12 years, 7 months ago) by guez
File size: 14839 byte(s)
Split "conflx.f" into single-procedure files in directory "Conflx".

Split "cv_routines.f" into single-procedure files in directory
"CV_routines". Made module "cvparam" from included file
"cvparam.h". No included file other than "netcdf.inc" left in LMDZE.

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

  ViewVC Help
Powered by ViewVC 1.1.21