GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/rrtm/srtm_reftra.F90 Lines: 0 62 0.0 %
Date: 2023-06-30 12:51:15 Branches: 0 14 0.0 %

Line Branch Exec Source
1
#ifdef RS6K
2
@PROCESS HOT NOSTRICT
3
#endif
4
SUBROUTINE 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
69
USE PARKIND1  ,ONLY : JPIM     ,JPRB
70
USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
71
72
USE PARSRTM , ONLY : JPLAY
73
!USE YOESW   , ONLY : NDBUG
74
USE YOERDU  , ONLY : REPLOG
75
76
IMPLICIT NONE
77
78
INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV
79
INTEGER(KIND=JPIM),INTENT(OUT)   :: KMODTS
80
LOGICAL           ,INTENT(IN)    :: LDRTCHK(JPLAY)
81
REAL(KIND=JPRB)   ,INTENT(IN)    :: PGG(JPLAY)
82
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRMUZ
83
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTAU(JPLAY)
84
REAL(KIND=JPRB)   ,INTENT(IN)    :: PW(JPLAY)
85
REAL(KIND=JPRB)   ,INTENT(INOUT)   :: PREF(JPLAY)
86
REAL(KIND=JPRB)   ,INTENT(INOUT)   :: PREFD(JPLAY)
87
REAL(KIND=JPRB)   ,INTENT(INOUT)   :: PTRA(JPLAY)
88
REAL(KIND=JPRB)   ,INTENT(INOUT)   :: PTRAD(JPLAY)
89
!     ------------------------------------------------------------------
90
91
INTEGER(KIND=JPIM) :: JK
92
93
REAL(KIND=JPRB) :: ZA, ZA1, ZA2
94
REAL(KIND=JPRB) :: ZBETA, ZDEND, ZDENR, ZDENT
95
REAL(KIND=JPRB) :: ZE1, ZE2, ZEM1, ZEM2, ZEMM, ZEP1, ZEP2
96
REAL(KIND=JPRB) :: ZG, ZG3, ZGAMMA1, ZGAMMA2, ZGAMMA3, ZGAMMA4, ZGT
97
REAL(KIND=JPRB) :: ZR1, ZR2, ZR3, ZR4, ZR5, ZRK, ZRK2, ZRKG, ZRM1, ZRP, ZRP1, ZRPP
98
REAL(KIND=JPRB) :: ZSR3, ZT1, ZT2, ZT3, ZT4, ZT5, ZTO1
99
REAL(KIND=JPRB) :: ZW, ZWCRIT, ZWO
100
REAL(KIND=JPRB) :: ZHOOK_HANDLE,EXP500,ZTEMP
101
102
!     ------------------------------------------------------------------
103
104
IF (LHOOK) CALL DR_HOOK('SRTM_REFTRA',0,ZHOOK_HANDLE)
105
EXP500=EXP(500.0_JPRB)
106
ZSR3=SQRT(3._JPRB)
107
ZWCRIT=0.9995_JPRB
108
KMODTS=2
109
110
!NDBUG=3
111
112
DO 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
260
ENDDO
261
262
!     ------------------------------------------------------------------
263
IF (LHOOK) CALL DR_HOOK('SRTM_REFTRA',1,ZHOOK_HANDLE)
264
END SUBROUTINE SRTM_REFTRA