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

Annotation of /trunk/Sources/phylmd/Radlwsw/radlwsw.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 53 - (hide 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 guez 53 module radlwsw_m
2 guez 3
3 guez 53 IMPLICIT none
4 guez 3
5 guez 53 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