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

Last change on this file since 3807 was 3759, checked in by jbclement, 5 weeks ago

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