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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 71 - (hide annotations)
Mon Jul 8 18:12:18 2013 UTC (10 years, 10 months ago) by guez
Original Path: trunk/libf/phylmd/Radlwsw/radlwsw.f90
File size: 13846 byte(s)
No reason to call inidissip in ce0l.

In inidissip, set random seed to 1 beacuse PGI compiler does not
accept all zeros.

dq was computed needlessly in caladvtrac. Arguments masse and dq of
calfis not used.

Replaced real*8 by double precision.

Pass arrays with inverted order of vertical levels to conflx instead
of creating local variables for this inside conflx.

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

  ViewVC Help
Powered by ViewVC 1.1.21