source: trunk/mesoscale/LMDZ.MARS/libf_gcm/phymars/swr.F @ 91

Last change on this file since 91 was 57, checked in by aslmd, 15 years ago

mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)
File size: 11.8 KB
RevLine 
[57]1      SUBROUTINE SWR ( KDLON, KFLEV, KNU
2     S     ,  aerosol,albedo,PDSIG,PPSOL,PRMU,PSEC
3     S     ,  PFD,PFU )
4
5      IMPLICIT NONE
6C     
7#include "dimensions.h"
8#include "dimphys.h"
9#include "dimradmars.h"
10#include "callkeys.h"
11
12#include "yomaer.h"
13#include "yomlw.h"
14
15#include "fisice.h"
16
17#include "aerice.h"
18C     
19C   SWR - Continuum scattering computations
20C     
21C     PURPOSE.
22C     --------
23C     Computes the reflectivity and transmissivity in case oF
24C     Continuum scattering
25c     F. Forget (1999)
26c
27c      BASED ON MORCRETTE EARTH MODEL
28C     (See radiation's part of the ecmwf research department
29C     documentation, and Fouquart and BonneL (1980)
30C     
31C     IMPLICIT ARGUMENTS :
32C     --------------------
33C     
34C     ==== INPUTS ===
35c
36c    KDLON :  number of horizontal grid points
37c    KFLEV :  number of vertical layers
38c    KNU   :   Solar band # (1 or 2)
39c   aerosol               aerosol extinction optical depth
40c                         at reference wavelength "longrefvis" set
41c                         in dimradmars.h , in each layer, for one of
42c                         the "naerkind" kind of aerosol optical properties.
43c    albedo   hemispheric surface albedo
44c                         albedo (i,1) : mean albedo for solar band#1
45c                                        (see below)
46c                         albedo (i,2) : mean albedo for solar band#2
47c                                        (see below)
48c    PDSIG      layer thickness in sigma coordinates
49c    PPSOL       Surface pressure (Pa)
50c    PRMU:  cos of solar zenith angle (=1 when sun at zenith)
51c           (CORRECTED for high zenith angle (atmosphere), unlike mu0)
52c    PSEC   =1./PRMU
53
54C     ==== OUTPUTS ===
55c
56c    PFD : downward flux in spectral band #INU in a given mesh
57c         (normalized to the total incident flux at the top of the atmosphere)
58c    PFU : upward flux in specatral band #INU in a given mesh
59c         (normalized to the total incident flux at the top of the atmosphere)
60C
61C     
62C     METHOD.
63C     -------
64C     
65C     Computes continuum fluxes corresponding to aerosoL
66C     Or/and rayleigh scattering (no molecular gas absorption)
67C     
68C-----------------------------------------------------------------------
69C     
70C     
71C-----------------------------------------------------------------------
72C     
73     
74C     ARGUMENTS
75C     ---------
76      INTEGER KDLON, KFLEV, KNU
77      REAL aerosol(NDLO2,KFLEV,naerkind), albedo(NDLO2,2),
78     S     PDSIG(NDLO2,KFLEV),PSEC(NDLO2)
79      REAL PPSOL(NDLO2)
80      REAL PFD(NDLO2,KFLEV+1),PFU(NDLO2,KFLEV+1)
81      REAL PRMU(NDLO2)
82
83C     LOCAL ARRAYS
84C     ------------
85 
86      INTEGER jk,ja,jl,jae, jkl,jklp1,jkm1,jaj
87      REAL ZTRAY, ZRATIO,ZGAR, ZFF
88      real zfacoa,zcorae
89      real ZMUE, zgap,zbmu0, zww,zto,zden,zmu1,zbmu1,zden1,zre11
90
91      REAL ZC1I(NDLON,NFLEV+1), ZGG(NDLON), ZREF(NDLON)
92     S ,  ZRE1(NDLON), ZRE2(NDLON)
93     S ,  ZRMUZ(NDLON), ZRNEB(NDLON), ZR21(NDLON)
94     S ,  ZR23(NDLON),  ZSS1(NDLON), ZTO1(NDLON), ZTR(NDLON,2,NFLEV+1)
95     S ,  ZTR1(NDLON), ZTR2(NDLON), ZW(NDLON)
96
97      REAL ZRAY1(NDLO2,NFLEV+1), ZRAY2(NDLO2,NFLEV+1)
98     s   ,  ZREFZ(NDLO2,2,NFLEV+1)
99     S   ,  ZRMUE(NDLO2,NFLEV+1)
100     S   ,  ZCGAZ(NDLO2,NFLEV),ZPIZAZ(NDLO2,NFLEV),ZTAUAZ(NDLO2,NFLEV)
101
102      REAL  ZRAYL(NDLON)
103     S     ,   ZRJ(NDLON,6,NFLEV+1)
104     S     ,  ZRK(NDLON,6,NFLEV+1)
105     S     ,  ZTRA1(NDLON,NFLEV+1), ZTRA2(NDLON,NFLEV+1)
106
107
108      real ray,coefsizew
109
110c     Function
111c     --------
112      real CVMGT
113
114C    --------------------------------
115C     OPTICAL PARAMETERS FOR AEROSOLS
116C     -------------------------------
117C     
118      DO  JK = 1 , nlaylte+1
119         DO  JA = 1 , 6
120            DO JL = 1 , KDLON
121               ZRJ(JL,JA,JK) = 0.
122               ZRK(JL,JA,JK) = 0.
123            END DO
124         END DO
125      END DO
126
127c Computing TOTAL single scattering parameters by adding
128c  properties of all the NAERKIND kind of aerosols
129
130      DO JK = 1 , nlaylte
131         DO  JL = 1 , KDLON
132            ZCGAZ(JL,JK) = 0.
133            ZPIZAZ(JL,JK) =  0.
134            ZTAUAZ(JL,JK) = 0.
135         END DO
136         DO 106 JAE=1,naerkind
137            DO 105 JL = 1 , KDLON
138c              Mean Extinction optical depth in the spectral band
139c              ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
140               ZTAUAZ(JL,JK)=ZTAUAZ(JL,JK)
141     S              +aerosol(JL,JK,JAE)*QVISsQREF(KNU,JAE)
142c              Single scattering albedo
143c              ~~~~~~~~~~~~~~~~~~~~~~~~
144c TEST : to account for the varying w with varying crystal size
145               if (activice.and.JAE.eq.naerkind.and.KNU.eq.2) then
146                 ray=min( max(rice(JL,JK)*1.e+6, 1.),10.)
147                 coefsizew=(0.0001417*ray**2.-0.00328*ray+0.99667)
148     &             /omegavis(KNU,JAE)
149               else
150                 coefsizew=1.
151               endif
152               ZPIZAZ(JL,JK)=ZPIZAZ(JL,JK)+aerosol(JL,JK,JAE)
153     S           * QVISsQREF(KNU,JAE)*omegavis(KNU,JAE)*coefsizew
154c              Assymetry factor
155c              ~~~~~~~~~~~~~~~~
156               ZCGAZ(JL,JK) =  ZCGAZ(JL,JK) +aerosol(JL,JK,JAE)
157     S            * QVISsQREF(KNU,JAE)*omegavis(KNU,JAE)*gvis(KNU,JAE)
158 105        CONTINUE
159 106     CONTINUE
160      END DO
161C     
162      DO JK = 1 , nlaylte
163         DO JL = 1 , KDLON
164            ZCGAZ(JL,JK) = CVMGT( 0., ZCGAZ(JL,JK) / ZPIZAZ(JL,JK),
165     S            (ZPIZAZ(JL,JK).EQ.0) )
166            ZPIZAZ(JL,JK) = CVMGT( 1., ZPIZAZ(JL,JK) / ZTAUAZ(JL,JK),
167     S           (ZTAUAZ(JL,JK).EQ.0) )
168         END DO
169      END DO
170
171C     --------------------------------
172C     INCLUDING RAYLEIGH SCATERRING
173C     -------------------------------
174      if (rayleigh) then
175
176        call swrayleigh(kdlon,knu,ppsol,prmu,ZRAYL)
177
178c       Modifying mean aerosol parameters to account rayleigh scat by gas:
179
180        DO JK = 1 , nlaylte
181           DO JL = 1 , KDLON
182c             Rayleigh opacity in each layer :
183              ZTRAY = ZRAYL(JL) * PDSIG(JL,JK)
184c             ratio Tau(rayleigh) / Tau (total)
185              ZRATIO = ZTRAY / (ZTRAY + ZTAUAZ(JL,JK))
186              ZGAR = ZCGAZ(JL,JK)
187              ZFF = ZGAR * ZGAR
188                ZTAUAZ(JL,JK)=ZTRAY+ZTAUAZ(JL,JK)*(1.-ZPIZAZ(JL,JK)*ZFF)
189              ZCGAZ(JL,JK) = ZGAR * (1. - ZRATIO) / (1. + ZGAR)
190              ZPIZAZ(JL,JK) =ZRATIO+(1.-ZRATIO)*ZPIZAZ(JL,JK)*(1.-ZFF)
191     S           / (1. -ZPIZAZ(JL,JK) * ZFF)
192           END DO
193        END DO
194      end if
195
196     
197C     ----------------------------------------------
198C     TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL
199C     ----------------------------------------------
200C     
201 200  CONTINUE
202     
203      DO JL = 1 , KDLON
204         ZR23(JL) = 0.
205         ZC1I(JL,nlaylte+1) = 0.
206      END DO
207     
208      DO JK = 1 , nlaylte
209         JKL = nlaylte+1 - JK
210         JKLP1 = JKL + 1
211         DO JL = 1 , KDLON
212            ZFACOA = 1.-ZPIZAZ(JL,JKL)*ZCGAZ(JL,JKL)*ZCGAZ(JL,JKL)
213            ZCORAE = ZFACOA * ZTAUAZ(JL,JKL) * PSEC(JL)
214            ZR21(JL) = EXP(-ZCORAE   )
215            ZSS1(JL) =  1.0-ZR21(JL)
216            ZC1I(JL,JKL) = 1.0-(1.0-ZSS1(JL))*(1.0-ZC1I(JL,JKLP1))
217         END DO
218      END DO
219
220C     -----------------------------------------------
221C     REFLECTIVITY/TRANSMISSIVITY FOR PURE SCATTERING
222C     -----------------------------------------------
223C     
224      DO JL = 1 , KDLON
225         ZRAY1(JL,nlaylte+1) = 0.
226         ZRAY2(JL,nlaylte+1) = 0.
227         ZREFZ(JL,2,1) = albedo(JL,KNU)
228         ZREFZ(JL,1,1) = albedo(JL,KNU)
229         ZTRA1(JL,nlaylte+1) = 1.
230         ZTRA2(JL,nlaylte+1) = 1.
231      END DO
232
233      DO JK = 2 , nlaylte+1
234         JKM1 = JK-1
235         DO 342 JL = 1 , KDLON
236            ZRNEB(JL)= 1.e-5   ! used to be "cloudiness" (PCLDSW in Morcrette)
237
238            ZRE1(JL)=0.
239            ZTR1(JL)=0.
240            ZRE2(JL)=0.
241            ZTR2(JL)=0.
242     
243C           EQUIVALENT ZENITH ANGLE
244c           ~~~~~~~~~~~~~~~~~~~~~~~
245            ZMUE = (1.-ZC1I(JL,JK)) * PSEC(JL)
246     S           + ZC1I(JL,JK) * 1.66
247            ZRMUE(JL,JK) = 1./ZMUE
248
249C     ------------------------------------------------------------------
250C          REFLECT./TRANSMISSIVITY DUE TO AEROSOLS (and rayleigh ?)
251C     ------------------------------------------------------------------
252
253            ZGAP = ZCGAZ(JL,JKM1)
254            ZBMU0 = 0.5 - 0.75 * ZGAP / ZMUE
255            ZWW =ZPIZAZ(JL,JKM1)
256            ZTO = ZTAUAZ(JL,JKM1)
257            ZDEN = 1. + (1. - ZWW + ZBMU0 * ZWW) * ZTO * ZMUE
258     S           + (1-ZWW) * (1. - ZWW +2.*ZBMU0*ZWW)*ZTO*ZTO*ZMUE*ZMUE
259            ZRAY1(JL,JKM1) = ZBMU0 * ZWW * ZTO * ZMUE / ZDEN
260            ZTRA1(JL,JKM1) = 1. / ZDEN
261C     
262            ZMU1 = 0.5
263            ZBMU1 = 0.5 - 0.75 * ZGAP * ZMU1
264            ZDEN1= 1. + (1. - ZWW + ZBMU1 * ZWW) * ZTO / ZMU1
265     S         + (1-ZWW) * (1. - ZWW +2.*ZBMU1*ZWW)*ZTO*ZTO/ZMU1/ZMU1
266            ZRAY2(JL,JKM1) = ZBMU1 * ZWW * ZTO / ZMU1 / ZDEN1
267            ZTRA2(JL,JKM1) = 1. / ZDEN1
268
269            ZGG(JL) =  ZCGAZ(JL,JKM1)
270            ZW(JL) =ZPIZAZ(JL,JKM1)
271            ZREF(JL) = ZREFZ(JL,1,JKM1)
272            ZRMUZ(JL) = ZRMUE(JL,JK)
273            ZTO1(JL) =  ZTAUAZ(JL,JKM1)/ZPIZAZ(JL,JKM1)
274
275 342     CONTINUE
276
277C     
278         CALL DEDD ( KDLON
279     S        , ZGG,ZREF,ZRMUZ,ZTO1,ZW
280     S        , ZRE1,ZRE2,ZTR1,ZTR2     )
281C     
282         DO JL = 1 , KDLON
283C     
284            ZREFZ(JL,1,JK) = (1.-ZRNEB(JL)) * (ZRAY1(JL,JKM1)
285     S           + ZREFZ(JL,1,JKM1) * ZTRA1(JL,JKM1)
286     S           * ZTRA2(JL,JKM1)
287     S           / (1.-ZRAY2(JL,JKM1)*ZREFZ(JL,1,JKM1)))
288     S           + ZRNEB(JL) * ZRE2(JL)
289C     
290            ZTR(JL,1,JKM1) = ZRNEB(JL) * ZTR2(JL) + (ZTRA1(JL,JKM1)
291     S           / (1.-ZRAY2(JL,JKM1)*ZREFZ(JL,1,JKM1)))
292     S           * (1.-ZRNEB(JL))
293C     
294            ZREFZ(JL,2,JK) = (1.-ZRNEB(JL)) * (ZRAY1(JL,JKM1)
295     S           + ZREFZ(JL,2,JKM1) * ZTRA1(JL,JKM1)
296     S           * ZTRA2(JL,JKM1) )
297     S           + ZRNEB(JL) * ZRE1(JL)
298C     
299            ZTR(JL,2,JKM1) = ZRNEB(JL) * ZTR1(JL)
300     S           + ZTRA1(JL,JKM1) * (1.-ZRNEB(JL))
301C     
302         END DO
303      END DO
304C     
305C     
306C     ------------------------------------------------------------------
307C     
308C     *         3.5    REFLECT./TRANSMISSIVITY BETWEEN SURFACE AND LEVEL
309C     -------------------------------------------------
310C     
311 350  CONTINUE
312C     
313      IF (KNU.EQ.1) THEN
314         JAJ = 2
315         DO 351 JL = 1 , KDLON
316            ZRJ(JL,JAJ,nlaylte+1) = 1.
317            ZRK(JL,JAJ,nlaylte+1) = ZREFZ(JL, 1,nlaylte+1)
318 351     CONTINUE
319C     
320         DO 353 JK = 1 , nlaylte
321            JKL = nlaylte+1 - JK
322            JKLP1 = JKL + 1
323            DO 352 JL = 1 , KDLON
324               ZRE11= ZRJ(JL,JAJ,JKLP1) * ZTR(JL,  1,JKL)
325               ZRJ(JL,JAJ,JKL) = ZRE11
326               ZRK(JL,JAJ,JKL) = ZRE11 * ZREFZ(JL,  1,JKL)
327 352        CONTINUE
328 353     CONTINUE
329 354     CONTINUE
330C     
331      ELSE
332C     
333         DO 358 JAJ = 1 , 2
334            DO 355 JL = 1 , KDLON
335               ZRJ(JL,JAJ,nlaylte+1) = 1.
336               ZRK(JL,JAJ,nlaylte+1) = ZREFZ(JL,JAJ,nlaylte+1)
337 355        CONTINUE
338C     
339            DO 357 JK = 1 , nlaylte
340               JKL = nlaylte+1 - JK
341               JKLP1 = JKL + 1
342               DO 356 JL = 1 , KDLON
343                  ZRE11= ZRJ(JL,JAJ,JKLP1) * ZTR(JL,JAJ,JKL)
344                  ZRJ(JL,JAJ,JKL) = ZRE11
345                  ZRK(JL,JAJ,JKL) = ZRE11 * ZREFZ(JL,JAJ,JKL)
346 356           CONTINUE
347 357        CONTINUE
348 358     CONTINUE
349      END IF
350
351C     
352C     
353C     
354C     ------------------------------------------------------------------
355C     ---------------
356C     DOWNWARD FLUXES
357C     ---------------
358C   
359      JAJ = 2
360
361      do JK = 1 , nlaylte+1
362        JKL = nlaylte+1 - JK + 1
363        DO  JL = 1 , KDLON
364            PFD(JL,JKL) =   ZRJ(JL,JAJ,JKL) * sunfr(KNU)
365        end do
366      end do
367C   
368C  -------------
369C  UPWARD FLUXES
370C  -------------
371      DO JK = 1 , nlaylte+1
372         DO  JL = 1 , KDLON
373c           ZRK = upward flux / incident top flux
374            PFU(JL,JK) =    ZRK(JL,JAJ,JK) * sunfr(KNU)
375         END DO
376      END DO
377
378C     
379      RETURN
380      END
Note: See TracBrowser for help on using the repository browser.