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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 38 - (hide 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 guez 3 !
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 guez 38 use SUPHEC_M
21 guez 3 use raddim, only: kflev, kdlon
22 guez 38 use yoethf_m
23 guez 3 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 guez 10 real, intent(in):: pplay(klon,klev)
83 guez 3 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