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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 62 - (show annotations)
Thu Jul 26 14:37:37 2012 UTC (11 years, 9 months ago) by guez
File size: 14142 byte(s)
Changed handling of compiler in compilation system.

Removed the prefix letters "y", "p", "t" or "z" in some names of variables.

Replaced calls to NetCDF by calls to NetCDF95.

Extracted "ioget_calendar" procedures from "calendar.f90" into a
separate file.

Extracted to a separate file, "mathop2.f90", procedures that were not
part of the generic interface "mathop" in "mathop.f90".

Removed computation of "dq" in "bilan_dyn", which was not used.

In "iniadvtrac", removed schemes 20 Slopes and 30 Prather. Was not
compatible with declarations of array sizes.

In "clcdrag", "ustarhb", "vdif_kcay", "yamada4" and "coefkz", changed
the size of some arrays from "klon" to "knon".

Removed possible call to "conema3" in "physiq".

Removed unused argument "cd" in "yamada".

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

  ViewVC Help
Powered by ViewVC 1.1.21