source: LMDZ5/branches/testing/libf/phylmd/rrtm/srtm_reftra.F90

Last change on this file was 1999, checked in by Laurent Fairhead, 11 years ago

Merged trunk changes r1920:1997 into testing branch
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 7.5 KB
Line 
1#ifdef RS6K
2@PROCESS HOT NOSTRICT
3#endif
4SUBROUTINE SRTM_REFTRA &
5 & ( KLEV  , KMODTS, &
6 &   LDRTCHK, &
7 &   PGG   , PRMUZ, PTAU , PW, &
8 &   PREF  , PREFD, PTRA , PTRAD &
9 & ) 
10 
11!**** *SRTM_REFTRA* - REFLECTIVITY AND TRANSMISSIVITY
12
13!     PURPOSE.
14!     --------
15!           COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY OF A CLEAR OR
16!     CLOUDY LAYER USING A CHOICE OF VARIOUS APPROXIMATIONS.
17
18!**   INTERFACE.
19!     ----------
20!          *SRTM_REFTRA* IS CALLED BY *SRTM_SPCVRT*
21
22!        EXPLICIT ARGUMENTS :
23!        --------------------
24! INPUTS
25! ------
26!      KMODTS  = 1 EDDINGTON (JOSEPH ET AL., 1976)
27!              = 2 PIFM (ZDUNKOWSKI ET AL., 1980)
28!              = 3 DISCRETE ORDINATES (LIOU, 1973)
29!      LDRTCHK = .T. IF CLOUDY
30!              = .F. IF CLEAR-SKY
31!      PGG     = ASSYMETRY FACTOR
32!      PRMUZ   = COSINE SOLAR ZENITH ANGLE
33!      PTAU    = OPTICAL THICKNESS
34!      PW      = SINGLE SCATTERING ALBEDO
35
36! OUTPUTS
37! -------
38!      PREF    : COLLIMATED BEAM REFLECTIVITY
39!      PREFD   : DIFFUSE BEAM REFLECTIVITY
40!      PTRA    : COLLIMATED BEAM TRANSMISSIVITY
41!      PTRAD   : DIFFUSE BEAM TRANSMISSIVITY
42
43!     METHOD.
44!     -------
45!          STANDARD DELTA-EDDINGTON, P.I.F.M., OR D.O.M. LAYER CALCULATIONS.
46
47!     EXTERNALS.
48!     ----------
49!          NONE
50
51!     REFERENCE.
52!     ----------
53
54!     AUTHOR.
55!     -------
56!        JEAN-JACQUES MORCRETTE  *ECMWF*
57
58!     MODIFICATIONS.
59!     --------------
60!        ORIGINAL : 03-02-27
61!        M.Hamrud   01-Oct-2003      CY28 Cleaning
62!        Mike Iacono, AER, Mar 2004: bug fix
63
64!     ------------------------------------------------------------------
65
66!*       0.1   ARGUMENTS
67!              ---------
68
69USE PARKIND1  ,ONLY : JPIM     ,JPRB
70USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
71
72USE PARSRTM , ONLY : JPLAY
73!USE YOESW   , ONLY : NDBUG
74USE YOERDU  , ONLY : REPLOG
75
76IMPLICIT NONE
77
78INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV
79INTEGER(KIND=JPIM),INTENT(OUT)   :: KMODTS
80LOGICAL           ,INTENT(IN)    :: LDRTCHK(JPLAY)
81REAL(KIND=JPRB)   ,INTENT(IN)    :: PGG(JPLAY)
82REAL(KIND=JPRB)   ,INTENT(IN)    :: PRMUZ
83REAL(KIND=JPRB)   ,INTENT(IN)    :: PTAU(JPLAY)
84REAL(KIND=JPRB)   ,INTENT(IN)    :: PW(JPLAY)
85REAL(KIND=JPRB)   ,INTENT(INOUT)   :: PREF(JPLAY)
86REAL(KIND=JPRB)   ,INTENT(INOUT)   :: PREFD(JPLAY)
87REAL(KIND=JPRB)   ,INTENT(INOUT)   :: PTRA(JPLAY)
88REAL(KIND=JPRB)   ,INTENT(INOUT)   :: PTRAD(JPLAY)
89!     ------------------------------------------------------------------
90
91INTEGER(KIND=JPIM) :: JK
92
93REAL(KIND=JPRB) :: ZA, ZA1, ZA2
94REAL(KIND=JPRB) :: ZBETA, ZDEND, ZDENR, ZDENT
95REAL(KIND=JPRB) :: ZE1, ZE2, ZEM1, ZEM2, ZEMM, ZEP1, ZEP2
96REAL(KIND=JPRB) :: ZG, ZG3, ZGAMMA1, ZGAMMA2, ZGAMMA3, ZGAMMA4, ZGT
97REAL(KIND=JPRB) :: ZR1, ZR2, ZR3, ZR4, ZR5, ZRK, ZRK2, ZRKG, ZRM1, ZRP, ZRP1, ZRPP
98REAL(KIND=JPRB) :: ZSR3, ZT1, ZT2, ZT3, ZT4, ZT5, ZTO1
99REAL(KIND=JPRB) :: ZW, ZWCRIT, ZWO
100REAL(KIND=JPRB) :: ZHOOK_HANDLE,EXP500,ZTEMP
101
102!     ------------------------------------------------------------------
103
104IF (LHOOK) CALL DR_HOOK('SRTM_REFTRA',0,ZHOOK_HANDLE)
105EXP500=EXP(500.0_JPRB)
106ZSR3=SQRT(3._JPRB)
107ZWCRIT=0.9995_JPRB
108KMODTS=2
109
110!NDBUG=3
111
112DO JK=1,KLEV
113!  if (NDBUG < 2) then
114!    print 9000,JK,LRTCHK(JK),PTAU(JK),PW(JK),PGG(JK),PRMUZ
115  9000 format(1x,'SRTM_REFTRA:inputs:',I3,L8,4E13.6)
116!  end if
117  IF (.NOT.LDRTCHK(JK)) THEN
118    PREF(JK) =0.0_JPRB
119    PTRA(JK) =1.0_JPRB
120    PREFD(JK)=0.0_JPRB
121    PTRAD(JK)=1.0_JPRB
122!    if (NDBUG < 2) then
123!      print 9001,JL,JK,PREF(JK),PTRA(JK),PREFD(JK),PTRAD(JK)
124    9001  format(1x,'SRTM_REFTRA:not.LRTCKH:',2I3,4F10.6)
125!    end if
126  ELSE
127    ZTO1=PTAU(JK)
128    ZW  =PW(JK)
129    ZG  =PGG(JK) 
130
131!-- GENERAL TWO-STREAM EXPRESSIONS
132
133    ZG3= 3._JPRB * ZG
134    IF (KMODTS == 1) THEN
135      ZGAMMA1= (7._JPRB - ZW * (4._JPRB + ZG3)) * 0.25_JPRB
136      ZGAMMA2=-(1._JPRB - ZW * (4._JPRB - ZG3)) * 0.25_JPRB
137      ZGAMMA3= (2._JPRB - ZG3 * PRMUZ ) * 0.25_JPRB
138    ELSEIF (KMODTS == 2) THEN 
139      ZGAMMA1= (8._JPRB - ZW * (5._JPRB + ZG3)) * 0.25_JPRB
140      ZGAMMA2=  3._JPRB *(ZW * (1._JPRB - ZG )) * 0.25_JPRB
141      ZGAMMA3= (2._JPRB - ZG3 * PRMUZ ) * 0.25_JPRB
142    ELSEIF (KMODTS == 3) THEN 
143      ZGAMMA1= ZSR3 * (2._JPRB - ZW * (1._JPRB + ZG)) * 0.5_JPRB
144      ZGAMMA2= ZSR3 * ZW * (1._JPRB - ZG ) * 0.5_JPRB
145      ZGAMMA3= (1._JPRB - ZSR3 * ZG * PRMUZ ) * 0.5_JPRB
146    ENDIF
147    ZGAMMA4= 1._JPRB - ZGAMMA3
148   
149!-- RECOMPUTE ORIGINAL S.S.A. TO TEST FOR CONSERVATIVE SOLUTION
150    ZWO= ZW / (1._JPRB - (1._JPRB - ZW) * (ZG / (1._JPRB - ZG))**2)
151!   ZTEMP=(1._JPRB - ZG)**2
152!   ZWO= ZW*ZTEMP/ (ZTEMP - (1._JPRB - ZW)*(ZG **2))
153   
154    IF (ZWO >= ZWCRIT) THEN
155!!-- conservative scattering
156
157      ZA  = ZGAMMA1 * PRMUZ
158      ZA1 = ZA - ZGAMMA3
159      ZGT = ZGAMMA1 * ZTO1
160       
161!-- Homogeneous reflectance and transmittance
162
163! collimated beam
164
165      ZE1 = MIN ( ZTO1 / PRMUZ , 500._JPRB)
166      ZE2 = EXP ( - ZE1 )
167      ZTEMP=1.0_JPRB/(1._JPRB + ZGT)
168      PREF(JK) = (ZGT - ZA1 * (1._JPRB - ZE2)) *ZTEMP
169      PTRA(JK) = 1._JPRB - PREF(JK)
170
171! isotropic incidence
172
173      PREFD(JK) = ZGT *ZTEMP
174      PTRAD(JK) = 1._JPRB - PREFD(JK)       
175
176!      if (NDBUG < 2) then
177!        print 9002,JL,JK,PREF(JK),PTRA(JK),PREFD(JK),PTRAD(JK)
178      9002  format(1x,'SRTM_REFTRA: consrv: LDRTCHK:',2I3,4F10.6)
179!      end if
180
181    ELSE
182
183!-- non-conservative scattering
184
185      ZA1 = ZGAMMA1 * ZGAMMA4 + ZGAMMA2 * ZGAMMA3
186      ZA2 = ZGAMMA1 * ZGAMMA3 + ZGAMMA2 * ZGAMMA4
187!      ZRK = SQRT ( ZGAMMA1**2 - ZGAMMA2**2)
188      ZRK = SQRT ( MAX ( REPLOG, ZGAMMA1**2 - ZGAMMA2**2) )
189      ZRP = ZRK * PRMUZ               
190      ZRP1 = 1._JPRB + ZRP
191      ZRM1 = 1._JPRB - ZRP
192      ZRK2 = 2._JPRB * ZRK
193      ZRPP = 1._JPRB - ZRP*ZRP
194      ZRKG = ZRK + ZGAMMA1
195      ZR1  = ZRM1 * (ZA2 + ZRK * ZGAMMA3)
196      ZR2  = ZRP1 * (ZA2 - ZRK * ZGAMMA3)
197      ZR3  = ZRK2 * (ZGAMMA3 - ZA2 * PRMUZ )
198      ZR4  = ZRPP * ZRKG
199      ZR5  = ZRPP * (ZRK - ZGAMMA1)
200      ZT1  = ZRP1 * (ZA1 + ZRK * ZGAMMA4)
201      ZT2  = ZRM1 * (ZA1 - ZRK * ZGAMMA4)
202      ZT3  = ZRK2 * (ZGAMMA4 + ZA1 * PRMUZ )
203      ZT4  = ZR4
204      ZT5  = ZR5
205      ZBETA = - ZR5 / ZR4
206       
207!-- Homogeneous reflectance and transmittance
208
209      IF(ZRK * ZTO1 > 500._JPRB)THEN
210        ZEP1=EXP500
211      ELSE
212        ZEP1=EXP(ZRK * ZTO1)
213      ENDIF
214      ZEM1=1.0_JPRB/ZEP1
215      IF(ZTO1 > 500._JPRB*PRMUZ)THEN
216        ZEP2=EXP500
217      ELSE
218        ZEP2=EXP(ZTO1 / PRMUZ)
219      ENDIF
220      ZEM2=1.0_JPRB/ZEP2
221
222!     ZE1 = MIN ( ZRK * ZTO1, 500._JPRB)
223!     ZE2 = MIN ( ZTO1 / PRMUZ , 500._JPRB)
224
225!     ZEP1 = EXP( ZE1 )
226!      ZEM1 = EXP(-ZE1 )
227!     ZEM1=1.0_JPRB/ZEP1
228
229!     ZEP2 = EXP( ZE2 )
230!      ZEM2 = EXP(-ZE2 )
231!     ZEM2=1.0_JPRB/ZEP2
232
233! collimated beam
234
235      ZDENR = ZR4*ZEP1 + ZR5*ZEM1
236!- bug noticed by Mike Iacono
237!      PREF(JK) = ZWO * (ZR1*ZEP1 - ZR2*ZEM1 - ZR3*ZEM2) / ZDENR
238      PREF(JK) = ZW  * (ZR1*ZEP1 - ZR2*ZEM1 - ZR3*ZEM2) / ZDENR
239
240      ZDENT = ZT4*ZEP1 + ZT5*ZEM1
241!- bug noticed by Mike Iacono
242!      PTRA(JK) = ZEM2 * (1._JPRB - ZWO * (ZT1*ZEP1 - ZT2*ZEM1 - ZT3*ZEP2) / ZDENT)
243      PTRA(JK) = ZEM2 * (1._JPRB - ZW  * (ZT1*ZEP1 - ZT2*ZEM1 - ZT3*ZEP2) / ZDENT)
244
245! diffuse beam
246
247      ZEMM = ZEM1*ZEM1
248      ZDEND = 1._JPRB / ( (1._JPRB - ZBETA*ZEMM ) * ZRKG)
249      PREFD(JK) =  ZGAMMA2 * (1._JPRB - ZEMM) * ZDEND
250      PTRAD(JK) =  ZRK2*ZEM1*ZDEND
251
252!      if (NDBUG < 2) then       
253!        print 9003,JL,JK,PREF(JK),PTRA(JK),PREFD(JK),PTRAD(JK)
254      9003  format(1x,'SRTM_REFTRA: OMG<1:  LDRTCHK:',2I3,4F10.6)
255!      end if
256    ENDIF
257
258  ENDIF         
259
260ENDDO   
261
262!     ------------------------------------------------------------------
263IF (LHOOK) CALL DR_HOOK('SRTM_REFTRA',1,ZHOOK_HANDLE)
264END SUBROUTINE SRTM_REFTRA     
Note: See TracBrowser for help on using the repository browser.