source: trunk/LMDZ.MARS/libf/phymars/swr_fouquart.F @ 3727

Last change on this file since 3727 was 3727, checked in by emillour, 3 months ago

Mars PCM:
Code tidying: put routines in modules, remove useless "return" statements and
remove obsolete and unused scopyi.F
EM
File size: 16.5 KB
Line 
1      MODULE swr_fouquart_mod
2     
3      IMPLICIT NONE
4     
5      CONTAINS
6     
7      SUBROUTINE SWR_FOUQUART ( KDLON, KFLEV, KNU
8     S     ,  aerosol,QVISsQREF3d,omegaVIS3d,gVIS3d
9     &     ,  albedo,PDSIG,PPSOL,PRMU,PSEC
10     S     ,  PFD,PFU )
11
12      use dimradmars_mod, only: sunfr, ndlo2, nsun,
13     &                          ndlon, nflev, naerkind
14      use yomlw_h, only: nlaylte
15      use callkeys_mod, only: rayleigh
16      IMPLICIT NONE
17C     
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_mod , 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
80      REAL QVISsQREF3d(NDLO2,KFLEV,nsun,naerkind)
81      REAL omegaVIS3d(NDLO2,KFLEV,nsun,naerkind)
82      REAL gVIS3d(NDLO2,KFLEV,nsun,naerkind)
83
84      REAL PPSOL(NDLO2)
85      REAL PFD(NDLO2,KFLEV+1),PFU(NDLO2,KFLEV+1)
86      REAL PRMU(NDLO2)
87
88C     LOCAL ARRAYS
89C     ------------
90 
91      INTEGER jk,ja,jl,jae, jkl,jklp1,jkm1,jaj
92      REAL ZTRAY, ZRATIO,ZGAR, ZFF
93      real zfacoa,zcorae
94      real ZMUE, zgap,zbmu0, zww,zto,zden,zmu1,zbmu1,zden1,zre11
95
96      REAL ZC1I(NDLON,NFLEV+1), ZGG(NDLON), ZREF(NDLON)
97     S ,  ZRE1(NDLON), ZRE2(NDLON)
98     S ,  ZRMUZ(NDLON), ZRNEB(NDLON), ZR21(NDLON)
99     S ,  ZR23(NDLON),  ZSS1(NDLON), ZTO1(NDLON), ZTR(NDLON,2,NFLEV+1)
100     S ,  ZTR1(NDLON), ZTR2(NDLON), ZW(NDLON)
101
102      REAL ZRAY1(NDLO2,NFLEV+1), ZRAY2(NDLO2,NFLEV+1)
103     s   ,  ZREFZ(NDLO2,2,NFLEV+1)
104     S   ,  ZRMUE(NDLO2,NFLEV+1)
105     S   ,  ZCGAZ(NDLO2,NFLEV),ZPIZAZ(NDLO2,NFLEV),ZTAUAZ(NDLO2,NFLEV)
106
107      REAL  ZRAYL(NDLON)
108     S     ,   ZRJ(NDLON,6,NFLEV+1)
109     S     ,  ZRK(NDLON,6,NFLEV+1)
110     S     ,  ZTRA1(NDLON,NFLEV+1), ZTRA2(NDLON,NFLEV+1)
111
112c     Function
113c     --------
114      real CVMGT
115
116C    --------------------------------
117C     OPTICAL PARAMETERS FOR AEROSOLS
118C     -------------------------------
119C     
120      DO  JK = 1 , nlaylte+1
121         DO  JA = 1 , 6
122            DO JL = 1 , KDLON
123               ZRJ(JL,JA,JK) = 0.
124               ZRK(JL,JA,JK) = 0.
125            END DO
126         END DO
127      END DO
128
129c Computing TOTAL single scattering parameters by adding
130c  properties of all the NAERKIND kind of aerosols
131
132      DO JK = 1 , nlaylte
133         DO  JL = 1 , KDLON
134            ZCGAZ(JL,JK) = 0.
135            ZPIZAZ(JL,JK) =  0.
136            ZTAUAZ(JL,JK) = 0.
137         END DO
138         DO 106 JAE=1,naerkind
139            DO 105 JL = 1 , KDLON
140c              Mean Extinction optical depth in the spectral band
141c              ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
142               ZTAUAZ(JL,JK)=ZTAUAZ(JL,JK)
143     S              +aerosol(JL,JK,JAE)*QVISsQREF3d(JL,JK,KNU,JAE)
144c              Single scattering albedo
145c              ~~~~~~~~~~~~~~~~~~~~~~~~
146               ZPIZAZ(JL,JK)=ZPIZAZ(JL,JK)+aerosol(JL,JK,JAE)*
147     S           QVISsQREF3d(JL,JK,KNU,JAE)*
148     &           omegaVIS3d(JL,JK,KNU,JAE)
149c              Assymetry factor
150c              ~~~~~~~~~~~~~~~~
151               ZCGAZ(JL,JK) =  ZCGAZ(JL,JK) +aerosol(JL,JK,JAE)*
152     S           QVISsQREF3d(JL,JK,KNU,JAE)*
153     &           omegaVIS3d(JL,JK,KNU,JAE)*gVIS3d(JL,JK,KNU,JAE)
154 105        CONTINUE
155 106     CONTINUE
156      END DO
157C     
158      DO JK = 1 , nlaylte
159         DO JL = 1 , KDLON
160            ZCGAZ(JL,JK) = CVMGT( 0., ZCGAZ(JL,JK) / ZPIZAZ(JL,JK),
161     S            (ZPIZAZ(JL,JK).EQ.0) )
162            ZPIZAZ(JL,JK) = CVMGT( 1., ZPIZAZ(JL,JK) / ZTAUAZ(JL,JK),
163     S           (ZTAUAZ(JL,JK).EQ.0) )
164         END DO
165      END DO
166
167C     --------------------------------
168C     INCLUDING RAYLEIGH SCATERRING
169C     -------------------------------
170      if (rayleigh) then
171
172        call swrayleigh(kdlon,knu,ppsol,prmu,ZRAYL)
173
174c       Modifying mean aerosol parameters to account rayleigh scat by gas:
175
176        DO JK = 1 , nlaylte
177           DO JL = 1 , KDLON
178c             Rayleigh opacity in each layer :
179              ZTRAY = ZRAYL(JL) * PDSIG(JL,JK)
180c             ratio Tau(rayleigh) / Tau (total)
181              ZRATIO = ZTRAY / (ZTRAY + ZTAUAZ(JL,JK))
182              ZGAR = ZCGAZ(JL,JK)
183              ZFF = ZGAR * ZGAR
184                ZTAUAZ(JL,JK)=ZTRAY+ZTAUAZ(JL,JK)*(1.-ZPIZAZ(JL,JK)*ZFF)
185              ZCGAZ(JL,JK) = ZGAR * (1. - ZRATIO) / (1. + ZGAR)
186              ZPIZAZ(JL,JK) =ZRATIO+(1.-ZRATIO)*ZPIZAZ(JL,JK)*(1.-ZFF)
187     S           / (1. -ZPIZAZ(JL,JK) * ZFF)
188           END DO
189        END DO
190      end if
191
192     
193C     ----------------------------------------------
194C     TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL
195C     ----------------------------------------------
196C     
197 200  CONTINUE
198     
199      DO JL = 1 , KDLON
200         ZR23(JL) = 0.
201         ZC1I(JL,nlaylte+1) = 0.
202      END DO
203     
204      DO JK = 1 , nlaylte
205         JKL = nlaylte+1 - JK
206         JKLP1 = JKL + 1
207         DO JL = 1 , KDLON
208            ZFACOA = 1.-ZPIZAZ(JL,JKL)*ZCGAZ(JL,JKL)*ZCGAZ(JL,JKL)
209            ZCORAE = ZFACOA * ZTAUAZ(JL,JKL) * PSEC(JL)
210            ZR21(JL) = EXP(-ZCORAE   )
211            ZSS1(JL) =  1.0-ZR21(JL)
212            ZC1I(JL,JKL) = 1.0-(1.0-ZSS1(JL))*(1.0-ZC1I(JL,JKLP1))
213         END DO
214      END DO
215
216C     -----------------------------------------------
217C     REFLECTIVITY/TRANSMISSIVITY FOR PURE SCATTERING
218C     -----------------------------------------------
219C     
220      DO JL = 1 , KDLON
221         ZRAY1(JL,nlaylte+1) = 0.
222         ZRAY2(JL,nlaylte+1) = 0.
223         ZREFZ(JL,2,1) = albedo(JL,KNU)
224         ZREFZ(JL,1,1) = albedo(JL,KNU)
225         ZTRA1(JL,nlaylte+1) = 1.
226         ZTRA2(JL,nlaylte+1) = 1.
227      END DO
228
229      DO JK = 2 , nlaylte+1
230         JKM1 = JK-1
231         DO 342 JL = 1 , KDLON
232            ZRNEB(JL)= 1.e-5   ! used to be "cloudiness" (PCLDSW in Morcrette)
233
234            ZRE1(JL)=0.
235            ZTR1(JL)=0.
236            ZRE2(JL)=0.
237            ZTR2(JL)=0.
238     
239C           EQUIVALENT ZENITH ANGLE
240c           ~~~~~~~~~~~~~~~~~~~~~~~
241            ZMUE = (1.-ZC1I(JL,JK)) * PSEC(JL)
242     S           + ZC1I(JL,JK) * 1.66
243            ZRMUE(JL,JK) = 1./ZMUE
244
245C     ------------------------------------------------------------------
246C          REFLECT./TRANSMISSIVITY DUE TO AEROSOLS (and rayleigh ?)
247C     ------------------------------------------------------------------
248
249            ZGAP = ZCGAZ(JL,JKM1)
250            ZBMU0 = 0.5 - 0.75 * ZGAP / ZMUE
251            ZWW =ZPIZAZ(JL,JKM1)
252            ZTO = ZTAUAZ(JL,JKM1)
253            ZDEN = 1. + (1. - ZWW + ZBMU0 * ZWW) * ZTO * ZMUE
254     S           + (1-ZWW) * (1. - ZWW +2.*ZBMU0*ZWW)*ZTO*ZTO*ZMUE*ZMUE
255            ZRAY1(JL,JKM1) = ZBMU0 * ZWW * ZTO * ZMUE / ZDEN
256            ZTRA1(JL,JKM1) = 1. / ZDEN
257C     
258            ZMU1 = 0.5
259            ZBMU1 = 0.5 - 0.75 * ZGAP * ZMU1
260            ZDEN1= 1. + (1. - ZWW + ZBMU1 * ZWW) * ZTO / ZMU1
261     S         + (1-ZWW) * (1. - ZWW +2.*ZBMU1*ZWW)*ZTO*ZTO/ZMU1/ZMU1
262            ZRAY2(JL,JKM1) = ZBMU1 * ZWW * ZTO / ZMU1 / ZDEN1
263            ZTRA2(JL,JKM1) = 1. / ZDEN1
264
265            ZGG(JL) =  ZCGAZ(JL,JKM1)
266            ZW(JL) =ZPIZAZ(JL,JKM1)
267            ZREF(JL) = ZREFZ(JL,1,JKM1)
268            ZRMUZ(JL) = ZRMUE(JL,JK)
269            ZTO1(JL) =  ZTAUAZ(JL,JKM1)/ZPIZAZ(JL,JKM1)
270
271 342     CONTINUE
272
273C     
274         CALL DEDD ( KDLON
275     S        , ZGG,ZREF,ZRMUZ,ZTO1,ZW
276     S        , ZRE1,ZRE2,ZTR1,ZTR2     )
277C     
278         DO JL = 1 , KDLON
279C     
280            ZREFZ(JL,1,JK) = (1.-ZRNEB(JL)) * (ZRAY1(JL,JKM1)
281     S           + ZREFZ(JL,1,JKM1) * ZTRA1(JL,JKM1)
282     S           * ZTRA2(JL,JKM1)
283     S           / (1.-ZRAY2(JL,JKM1)*ZREFZ(JL,1,JKM1)))
284     S           + ZRNEB(JL) * ZRE2(JL)
285C     
286            ZTR(JL,1,JKM1) = ZRNEB(JL) * ZTR2(JL) + (ZTRA1(JL,JKM1)
287     S           / (1.-ZRAY2(JL,JKM1)*ZREFZ(JL,1,JKM1)))
288     S           * (1.-ZRNEB(JL))
289C     
290            ZREFZ(JL,2,JK) = (1.-ZRNEB(JL)) * (ZRAY1(JL,JKM1)
291     S           + ZREFZ(JL,2,JKM1) * ZTRA1(JL,JKM1)
292     S           * ZTRA2(JL,JKM1) )
293     S           + ZRNEB(JL) * ZRE1(JL)
294C     
295            ZTR(JL,2,JKM1) = ZRNEB(JL) * ZTR1(JL)
296     S           + ZTRA1(JL,JKM1) * (1.-ZRNEB(JL))
297C     
298         END DO
299      END DO
300C     
301C     
302C     ------------------------------------------------------------------
303C     
304C     *         3.5    REFLECT./TRANSMISSIVITY BETWEEN SURFACE AND LEVEL
305C     -------------------------------------------------
306C     
307 350  CONTINUE
308C     
309      IF (KNU.EQ.1) THEN
310         JAJ = 2
311         DO 351 JL = 1 , KDLON
312            ZRJ(JL,JAJ,nlaylte+1) = 1.
313            ZRK(JL,JAJ,nlaylte+1) = ZREFZ(JL, 1,nlaylte+1)
314 351     CONTINUE
315C     
316         DO 353 JK = 1 , nlaylte
317            JKL = nlaylte+1 - JK
318            JKLP1 = JKL + 1
319            DO 352 JL = 1 , KDLON
320               ZRE11= ZRJ(JL,JAJ,JKLP1) * ZTR(JL,  1,JKL)
321               ZRJ(JL,JAJ,JKL) = ZRE11
322               ZRK(JL,JAJ,JKL) = ZRE11 * ZREFZ(JL,  1,JKL)
323 352        CONTINUE
324 353     CONTINUE
325 354     CONTINUE
326C     
327      ELSE
328C     
329         DO 358 JAJ = 1 , 2
330            DO 355 JL = 1 , KDLON
331               ZRJ(JL,JAJ,nlaylte+1) = 1.
332               ZRK(JL,JAJ,nlaylte+1) = ZREFZ(JL,JAJ,nlaylte+1)
333 355        CONTINUE
334C     
335            DO 357 JK = 1 , nlaylte
336               JKL = nlaylte+1 - JK
337               JKLP1 = JKL + 1
338               DO 356 JL = 1 , KDLON
339                  ZRE11= ZRJ(JL,JAJ,JKLP1) * ZTR(JL,JAJ,JKL)
340                  ZRJ(JL,JAJ,JKL) = ZRE11
341                  ZRK(JL,JAJ,JKL) = ZRE11 * ZREFZ(JL,JAJ,JKL)
342 356           CONTINUE
343 357        CONTINUE
344 358     CONTINUE
345      END IF
346
347C     
348C     
349C     
350C     ------------------------------------------------------------------
351C     ---------------
352C     DOWNWARD FLUXES
353C     ---------------
354C   
355      JAJ = 2
356
357      do JK = 1 , nlaylte+1
358        JKL = nlaylte+1 - JK + 1
359        DO  JL = 1 , KDLON
360            PFD(JL,JKL) =   ZRJ(JL,JAJ,JKL) * sunfr(KNU)
361        end do
362      end do
363C   
364C  -------------
365C  UPWARD FLUXES
366C  -------------
367      DO JK = 1 , nlaylte+1
368         DO  JL = 1 , KDLON
369c           ZRK = upward flux / incident top flux
370            PFU(JL,JK) =    ZRK(JL,JAJ,JK) * sunfr(KNU)
371         END DO
372      END DO
373
374C     
375
376      END SUBROUTINE SWR_FOUQUART
377
378CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
379
380      SUBROUTINE DEDD (KDLON,PGG,PREF,PRMUZ,PTO1,PW
381     S                ,      PRE1,PRE2,PTR1,PTR2         )
382      use dimradmars_mod, only: ndlo2
383      implicit none
384C
385C
386C**** *DEDD* - DELTA-EDDINGTON IN A CLOUDY LAYER
387C
388C     PURPOSE.
389C     --------
390C           COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY OF A CLOUDY
391C     LAYER USING THE DELTA-EDDINGTON'S APPROXIMATION.
392C
393C**   INTERFACE.
394C     ----------
395C          *DEDD* IS CALLED BY *SW*.
396C
397C     SUBROUTINE DEDD (KDLON,PGG,PREF,PRMUZ,PTO1,PW
398C    S                ,      PRE1,PRE2,PTR1,PTR2         )
399C
400C        EXPLICIT ARGUMENTS :
401C        --------------------
402C PGG    : (NDLON)             ; ASSYMETRY FACTOR
403C PREF   : (NDLON)             ; REFLECTIVITY OF THE UNDERLYING LAYER
404C PRMUZ  : (NDLON)             ; COSINE OF SOLAR ZENITH ANGLE
405C PTO1   : (NDLON)             ; OPTICAL THICKNESS
406C PW     : (NDLON)             ; SINGLE SCATTERING ALBEDO
407C     ==== OUTPUTS ===
408C PRE1   : (NDLON)             ; LAYER REFLECTIVITY ASSUMING NO
409C                              ; REFLECTION FROM UNDERLYING LAYER
410C PTR1   : (NDLON)             ; LAYER TRANSMISSIVITY ASSUMING NO
411C                              ; REFLECTION FROM UNDERLYING LAYER
412C PRE2   : (NDLON)             ; LAYER REFLECTIVITY ASSUMING
413C                              ; REFLECTION FROM UNDERLYING LAYER
414C PTR2   : (NDLON)             ; LAYER TRANSMISSIVITY ASSUMING
415C                              ; REFLECTION FROM UNDERLYING LAYER
416C
417C        IMPLICIT ARGUMENTS :   NONE
418C        --------------------
419C
420C     METHOD.
421C     -------
422C
423C          STANDARD DELTA-EDDINGTON LAYER CALCULATIONS.
424C
425C     EXTERNALS.
426C     ----------
427C
428C          NONE
429C
430C     REFERENCE.
431C     ----------
432C
433C        SEE RADIATION'S PART OF THE MODEL'S DOCUMENTATION AND
434C        ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "IN CORE MODEL"
435C
436C     AUTHOR.
437C     -------
438C        JEAN-JACQUES MORCRETTE  *ECMWF*
439C
440C     MODIFICATIONS.
441C     --------------
442C        ORIGINAL : 88-12-15
443C     ------------------------------------------------------------------
444C
445C*       0.1   ARGUMENTS
446C              ---------
447      INTEGER KDLON
448C
449      REAL PGG(NDLO2),PREF(NDLO2),PRMUZ(NDLO2),PTO1(NDLO2),PW(NDLO2)
450      REAL PRE1(NDLO2),PRE2(NDLO2),PTR1(NDLO2),PTR2(NDLO2)
451
452c   local
453      integer jl
454      real*8 ZFF,ZGP,ZTOP,ZWCP,ZDT,ZX1,ZWM,ZRM2,ZRK,ZX2,ZRP,ZALPHA
455      real*8 ZBETA,ZEXMU0,ZEXKP,ZEXKM,ZXP2P,ZXM2P,ZAP2B,ZAM2B
456      real*8 ZA11,ZA12,ZA13,ZA22,ZA21,ZA23,ZDENA,ZC1A,ZC2A
457      real*8 ZRI0A,ZRI1A,ZRI0B,ZRI1B
458      real*8 ZB21,ZB22,ZB23,ZDENB,ZC1B,ZC2B
459      real*8 ZRI0C,ZRI1C,ZRI0D,ZRI1D
460C
461C     ------------------------------------------------------------------
462C
463C*         1.      DELTA-EDDINGTON CALCULATIONS
464C
465 100  CONTINUE
466C
467      DO 131 JL   =   1 , KDLON
468C
469C*         1.1     SET UP THE DELTA-MODIFIED PARAMETERS
470C
471 110  CONTINUE
472C
473      ZFF = PGG(JL)*PGG(JL)
474      ZGP = PGG(JL)/(1.+PGG(JL))
475      ZTOP = (1.- PW(JL) * ZFF) * PTO1(JL)
476      ZWCP = (1-ZFF)* PW(JL) /(1.- PW(JL) * ZFF)
477      ZDT = 2./3.
478      ZX1 = 1.-ZWCP*ZGP
479      ZWM = 1.-ZWCP
480      ZRM2 =  PRMUZ(JL) * PRMUZ(JL)
481      ZRK = SQRT(3.*ZWM*ZX1)
482      ZX2 = 4.*(1.-ZRK*ZRK*ZRM2)
483      ZRP = SQRT(3.*ZWM/ZX1)
484      ZALPHA = 3.*ZWCP*ZRM2*(1.+ZGP*ZWM)/ZX2
485      ZBETA = 3.*ZWCP* PRMUZ(JL) *(1.+3.*ZGP*ZRM2*ZWM)/ZX2
486      ZEXMU0 = EXP(-ZTOP/ PRMUZ(JL) )
487      ZEXKP = EXP(ZRK*ZTOP)
488      ZEXKM = 1./ZEXKP
489      ZXP2P = 1.+ZDT*ZRP
490      ZXM2P = 1.-ZDT*ZRP
491      ZAP2B = ZALPHA+ZDT*ZBETA
492      ZAM2B = ZALPHA-ZDT*ZBETA
493C
494C*         1.2     WITHOUT REFLECTION FROM THE UNDERLYING LAYER
495C
496 120  CONTINUE
497C
498      ZA11 = ZXP2P
499      ZA12 = ZXM2P
500      ZA13 = ZAP2B
501      ZA22 = ZXP2P*ZEXKP
502      ZA21 = ZXM2P*ZEXKM
503      ZA23 = ZAM2B*ZEXMU0
504      ZDENA = ZA11 * ZA22 - ZA21 * ZA12
505      ZC1A = (ZA22*ZA13-ZA12*ZA23)/ZDENA
506      ZC2A = (ZA11*ZA23-ZA21*ZA13)/ZDENA
507      ZRI0A = ZC1A+ZC2A-ZALPHA
508      ZRI1A = ZRP*(ZC1A-ZC2A)-ZBETA
509      PRE1(JL) = (ZRI0A-ZDT*ZRI1A)/ PRMUZ(JL)
510      ZRI0B = ZC1A*ZEXKM+ZC2A*ZEXKP-ZALPHA*ZEXMU0
511      ZRI1B = ZRP*(ZC1A*ZEXKM-ZC2A*ZEXKP)-ZBETA*ZEXMU0
512      PTR1(JL) = ZEXMU0+(ZRI0B+ZDT*ZRI1B)/ PRMUZ(JL)
513C
514C*         1.3     WITH REFLECTION FROM THE UNDERLYING LAYER
515C
516 130  CONTINUE
517C
518      ZB21 = ZA21- PREF(JL) *ZXP2P*ZEXKM
519      ZB22 = ZA22- PREF(JL) *ZXM2P*ZEXKP
520      ZB23 = ZA23- PREF(JL) *ZEXMU0*(ZAP2B - PRMUZ(JL) )
521      ZDENB = ZA11 * ZB22 - ZB21 * ZA12
522      ZC1B = (ZB22*ZA13-ZA12*ZB23)/ZDENB
523      ZC2B = (ZA11*ZB23-ZB21*ZA13)/ZDENB
524      ZRI0C = ZC1B+ZC2B-ZALPHA
525      ZRI1C = ZRP*(ZC1B-ZC2B)-ZBETA
526      PRE2(JL) = (ZRI0C-ZDT*ZRI1C) / PRMUZ(JL)
527      ZRI0D = ZC1B*ZEXKM + ZC2B*ZEXKP - ZALPHA*ZEXMU0
528      ZRI1D = ZRP * (ZC1B*ZEXKM - ZC2B*ZEXKP) - ZBETA*ZEXMU0
529      PTR2(JL) = ZEXMU0 + (ZRI0D + ZDT*ZRI1D) / PRMUZ(JL)
530C
531 131  CONTINUE
532
533      END SUBROUTINE DEDD
534
535      END MODULE swr_fouquart_mod
Note: See TracBrowser for help on using the repository browser.