source: LMDZ4/trunk/libf/phylmd/radlwsw.F @ 1085

Last change on this file since 1085 was 998, checked in by Laurent Fairhead, 16 years ago

Modifications necessaires a la preparation au passage au nouveau rayonnement
RRTM MPL
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.7 KB
RevLine 
[524]1!
2! $Header$
3!
4      SUBROUTINE radlwsw(dist, rmu0, fract,
[888]5     .                  paprs, pplay,tsol,alb1, alb2, t,q,wo,
[524]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,
[998]16     .                  cldtaupi, topswai, solswai,qsat,flwc,fiwc)
[524]17c     
[766]18      USE dimphy
[524]19      IMPLICIT none
20c======================================================================
21c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
22c Objet: interface entre le modele et les rayonnements
23c Arguments:
24c dist-----input-R- distance astronomique terre-soleil
25c rmu0-----input-R- cosinus de l'angle zenithal
26c fract----input-R- duree d'ensoleillement normalisee
27c co2_ppm--input-R- concentration du gaz carbonique (en ppm)
28c solaire--input-R- constante solaire (W/m**2)
29c paprs----input-R- pression a inter-couche (Pa)
30c pplay----input-R- pression au milieu de couche (Pa)
31c tsol-----input-R- temperature du sol (en K)
[888]32c alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible
33c alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge
[524]34c t--------input-R- temperature (K)
35c q--------input-R- vapeur d'eau (en kg/kg)
[652]36c wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505
[524]37c cldfra---input-R- fraction nuageuse (entre 0 et 1)
38c cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value)
39c cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1)
40c ok_ade---input-L- apply the Aerosol Direct Effect or not?
41c ok_aie---input-L- apply the Aerosol Indirect Effect or not?
42c tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F)
43c cldtaupi-input-R- epaisseur optique des nuages dans le visible
44c                   calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller
45c                   droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd
46c                   it is needed for the diagnostics of the aerosol indirect radiative forcing     
47c
48c heat-----output-R- echauffement atmospherique (visible) (K/jour)
49c cool-----output-R- refroidissement dans l'IR (K/jour)
50c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
51c albpla---output-R- albedo planetaire (entre 0 et 1)
52c topsw----output-R- flux solaire net au sommet de l'atm.
53c toplw----output-R- ray. IR montant au sommet de l'atmosphere
54c solsw----output-R- flux solaire net a la surface
55c sollw----output-R- ray. IR montant a la surface
56c solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir)
57c topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir)
58c solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind)
59c topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind)
60c
61c ATTENTION: swai and swad have to be interpreted in the following manner:
62c ---------
63c ok_ade=F & ok_aie=F -both are zero
64c ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad
65c                        indirect is zero
66c ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
67c                        direct is zero
68c ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai
69c                        aerosol direct forcing is F_{AD} = topswai-topswad
70c
71     
72c======================================================================
[766]73cym#include "dimensions.h"
74cym#include "dimphy.h"
75cym#include "raddim.h"
[524]76#include "YOETHF.h"
77c
78      real rmu0(klon), fract(klon), dist
79cIM   real co2_ppm
80cIM   real solaire
81#include "clesphys.h"
82c
83      real paprs(klon,klev+1), pplay(klon,klev)
[888]84      real alb1(klon), alb2(klon), tsol(klon)
[524]85      real t(klon,klev), q(klon,klev), 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)
[644]93cIM output 3D
[524]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)
[644]98c
[524]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)
103c
104      REAL*8 zx_alpha1, zx_alpha2
105c
106#include "YOMCST.h"
107c
108      INTEGER k, kk, i, j, iof, nb_gr
[998]109      EXTERNAL LW_LMDAR4,SW_LMDAR4
[524]110c
111cIM ctes ds clesphys.h  REAL*8 RCO2, RCH4, RN2O, RCFC11, RCFC12
112      REAL*8 PSCT
113c
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)
127c
128      REAL*8 zfract(kdlon), zrmu0(kdlon), zdist
129c
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)
[644]134cIM
[524]135      REAL*8 zsollwdown(kdlon)
[644]136c
[524]137      REAL*8 ztopsw0(kdlon), ztoplw0(kdlon)
138      REAL*8 zsolsw0(kdlon), zsollw0(kdlon)
139      REAL*8 zznormcp
[644]140cIM output 3D : SWup, SWdn, LWup, LWdn
[524]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)
[998]145      REAL qsat(klon,klev),flwc(klon,klev),fiwc(klon,klev)
[524]146c-OB
147cjq the following quantities are needed for the aerosol radiative forcings
148
149      real topswad(klon), solswad(klon) ! output: aerosol direct forcing at TOA and surface
150      real topswai(klon), solswai(klon) ! output: aerosol indirect forcing atTOA and surface
151      real tau_ae(klon,klev,2), piz_ae(klon,klev,2), cg_ae(klon,klev,2) ! aerosol optical properties (see aeropt.F)
152      real cldtaupi(klon,klev)  ! cloud optical thickness for pre-industrial aerosol concentrations
153                                ! (i.e., with a smaller droplet concentrationand thus larger droplet radii)
154      logical ok_ade, ok_aie    ! switches whether to use aerosol direct (indirect) effects or not
155      real*8 tauae(kdlon,kflev,2) ! aer opt properties
156      real*8 pizae(kdlon,kflev,2)
157      real*8 cgae(kdlon,kflev,2)
158      REAL*8 PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use
159      REAL*8 POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo
160      REAL*8 ztopswad(kdlon), zsolswad(kdlon) ! Aerosol direct forcing at TOAand surface
161      REAL*8 ztopswai(kdlon), zsolswai(kdlon) ! dito, indirect
162cjq-end
[557]163!rv
164      tauae(:,:,:)=0.
165      pizae(:,:,:)=0.
166      cgae(:,:,:)=0.
167!rv
[524]168     
169c
170c-------------------------------------------
171      nb_gr = klon / kdlon
172      IF (nb_gr*kdlon .NE. klon) THEN
173         PRINT*, "kdlon mauvais:", klon, kdlon, nb_gr
174         CALL abort
175      ENDIF
176      IF (kflev .NE. klev) THEN
177          PRINT*, "kflev differe de klev, kflev, klev"
178          CALL abort
179      ENDIF
180c-------------------------------------------
181      DO k = 1, klev
182      DO i = 1, klon
183         heat(i,k)=0.
184         cool(i,k)=0.
185         heat0(i,k)=0.
186         cool0(i,k)=0.
187      ENDDO
188      ENDDO
189c
190      zdist = dist
191c
192cIM anciennes valeurs
193c     RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97
194c
195cIM : on met RCO2, RCH4, RN2O, RCFC11 et RCFC12 dans clesphys.h /lecture ds conf_phys.F90
196c     RCH4 = 1.65E-06* 16.043/28.97
197c     RN2O = 306.E-09* 44.013/28.97
198c     RCFC11 = 280.E-12* 137.3686/28.97
199c     RCFC12 = 484.E-12* 120.9140/28.97
200cIM anciennes valeurs
201c     RCH4 = 1.72E-06* 16.043/28.97
202c     RN2O = 310.E-09* 44.013/28.97
203c
204c     PRINT*,'IMradlwsw : solaire, co2= ', solaire, co2_ppm
205      PSCT = solaire/zdist/zdist
206c
207      DO 99999 j = 1, nb_gr
208      iof = kdlon*(j-1)
209c
210      DO i = 1, kdlon
211         zfract(i) = fract(iof+i)
212         zrmu0(i) = rmu0(iof+i)
[888]213         PALBD(i,1) = alb1(iof+i)
214!         PALBD(i,2) = alb1(iof+i)
215         PALBD(i,2) = alb2(iof+i)
216         PALBP(i,1) = alb1(iof+i)
217!         PALBP(i,2) = alb1(iof+i)
218         PALBP(i,2) = alb2(iof+i)
[524]219cIM cf. JLD pour etre en accord avec ORCHIDEE il faut mettre PEMIS(i) = 0.96
220         PEMIS(i) = 1.0
221         PVIEW(i) = 1.66
222         PPSOL(i) = paprs(iof+i,1)
223         zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))
224     .             / (pplay(iof+i,1)-pplay(iof+i,2))
225         zx_alpha2 = 1.0 - zx_alpha1
226         PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2
227         PTL(i,klev+1) = t(iof+i,klev)
228         PDT0(i) = tsol(iof+i) - PTL(i,1)
229      ENDDO
230      DO k = 2, kflev
231      DO i = 1, kdlon
232         PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5
233      ENDDO
234      ENDDO
235      DO k = 1, kflev
236      DO i = 1, kdlon
237         PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1)
238         PTAVE(i,k) = t(iof+i,k)
239         PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
240         PQS(i,k) = PWV(i,k)
241c wo:    cm.atm (epaisseur en cm dans la situation standard)
242c POZON: kg/kg
[699]243         POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968
[524]244     .               /(paprs(iof+i,k)-paprs(iof+i,k+1))
245     .               *(paprs(iof+i,1)/101325.0)
246         PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
247         PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
248         PCLDSW(i,k) = cldfra(iof+i,k)
249         PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable
250         PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines
251         POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k))
252         POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k))
253         PCG(i,1,k) = 0.865
254         PCG(i,2,k) = 0.910
255c-OB
256cjq Introduced for aerosol indirect forcings.
257cjq The following values use the cloud optical thickness calculated from
258cjq present-day aerosol concentrations whereas the quantities without the
259cjq "A" at the end are for pre-industial (natural-only) aerosol concentrations
260cjq
261         PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable
262         PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines
263         POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k))
264         POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k))
265cjq-end
266      ENDDO
267      ENDDO
268c
269      DO k = 1, kflev+1
270      DO i = 1, kdlon
271         PPMB(i,k) = paprs(iof+i,k)/100.0
272      ENDDO
273      ENDDO
274c
275      DO kk = 1, 5
276      DO k = 1, kflev
277      DO i = 1, kdlon
278         PAER(i,k,kk) = 1.0E-15
279      ENDDO
280      ENDDO
281      ENDDO
282c-OB
283      DO k = 1, kflev
284      DO i = 1, kdlon
285        tauae(i,k,1)=tau_ae(iof+i,k,1)
286        pizae(i,k,1)=piz_ae(iof+i,k,1)
287        cgae(i,k,1) =cg_ae(iof+i,k,1)
288        tauae(i,k,2)=tau_ae(iof+i,k,2)
289        pizae(i,k,2)=piz_ae(iof+i,k,2)
290        cgae(i,k,2) =cg_ae(iof+i,k,2)
291      ENDDO
292      ENDDO
293c
[998]294c===== si iflag_rrtm=0 ================================================
[524]295cIM ctes ds clesphys.h   CALL LW(RCO2,RCH4,RN2O,RCFC11,RCFC12,
[998]296cIM ctes ds clesphys.h   CALL SW(PSCT, RCO2, zrmu0, zfract,
297c     
298      if (iflag_rrtm.eq.0) then
299         CALL LW_LMDAR4(
[524]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)
[998]309         CALL SW_LMDAR4(PSCT, zrmu0, zfract,
[524]310     S        PPMB, PDP,
311     S        PPSOL, PALBD, PALBP,
312     S        PTAVE, PWV, PQS, POZON, PAER,
313     S        PCLDSW, PTAU, POMEGA, PCG,
314     S        zheat, zheat0,
315     S        zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,
316     S        ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,
317     S        tauae, pizae, cgae, ! aerosol optical properties
318     s        PTAUA, POMEGAA,
319     s        ztopswad,zsolswad,ztopswai,zsolswai, ! diagnosed aerosol forcing
320     J        ok_ade, ok_aie) ! apply aerosol effects or not?
[998]321      else
322c===== si iflag_rrtm=1, on passe dans SW via RECMWFL ===============
323          PRINT*, "Cette option ne fonctionne pas encore !!!"
324         CALL abort
325         endif   ! if(iflag_rrtm=0)
[524]326
327c======================================================================
328      DO i = 1, kdlon
329         radsol(iof+i) = zsolsw(i) + zsollw(i)
330         topsw(iof+i) = ztopsw(i)
331         toplw(iof+i) = ztoplw(i)
332         solsw(iof+i) = zsolsw(i)
333         sollw(iof+i) = zsollw(i)
334         sollwdown(iof+i) = zsollwdown(i)
335cIM
336         DO k = 1, kflev+1
337         lwdn0 ( iof+i,k)   = ZFLDN0 ( i,k)
338         lwdn  ( iof+i,k)   = ZFLDN  ( i,k)
339         lwup0 ( iof+i,k)   = ZFLUP0 ( i,k)
340         lwup  ( iof+i,k)   = ZFLUP  ( i,k)
341         ENDDO
[644]342c
[524]343         topsw0(iof+i) = ztopsw0(i)
344         toplw0(iof+i) = ztoplw0(i)
345         solsw0(iof+i) = zsolsw0(i)
346         sollw0(iof+i) = zsollw0(i)
347         albpla(iof+i) = zalbpla(i)
[644]348cIM
[524]349         DO k = 1, kflev+1
350         swdn0 ( iof+i,k)   = ZFSDN0 ( i,k)
351         swdn  ( iof+i,k)   = ZFSDN  ( i,k)
352         swup0 ( iof+i,k)   = ZFSUP0 ( i,k)
353         swup  ( iof+i,k)   = ZFSUP  ( i,k)
354         ENDDO !k=1, kflev+1
355      ENDDO
356cjq-transform the aerosol forcings, if they have
357cjq to be calculated
358      IF (ok_ade) THEN
359      DO i = 1, kdlon
360         topswad(iof+i) = ztopswad(i)
361         solswad(iof+i) = zsolswad(i)
362      ENDDO
363      ELSE
364      DO i = 1, kdlon
365         topswad(iof+i) = 0.0
366         solswad(iof+i) = 0.0
367      ENDDO
368      ENDIF
369      IF (ok_aie) THEN
370      DO i = 1, kdlon
371         topswai(iof+i) = ztopswai(i)
372         solswai(iof+i) = zsolswai(i)
373      ENDDO
374      ELSE
375      DO i = 1, kdlon
376         topswai(iof+i) = 0.0
377         solswai(iof+i) = 0.0
378      ENDDO
379      ENDIF
380cjq-end
381      DO k = 1, kflev
382c      DO i = 1, kdlon
383c         heat(iof+i,k) = zheat(i,k)
384c         cool(iof+i,k) = zcool(i,k)
385c         heat0(iof+i,k) = zheat0(i,k)
386c         cool0(iof+i,k) = zcool0(i,k)
387c      ENDDO
388      DO i = 1, kdlon
389C        scale factor to take into account the difference between
390C        dry air and watter vapour scpecific heat capacity
391         zznormcp=1.0+RVTMP2*PWV(i,k)
392         heat(iof+i,k) = zheat(i,k)/zznormcp
393         cool(iof+i,k) = zcool(i,k)/zznormcp
394         heat0(iof+i,k) = zheat0(i,k)/zznormcp
395         cool0(iof+i,k) = zcool0(i,k)/zznormcp
396      ENDDO
397      ENDDO
398c
39999999 CONTINUE
400      RETURN
401      END
Note: See TracBrowser for help on using the repository browser.