LMDZ
srtm_reftra.F90
Go to the documentation of this file.
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
srtm_reftra
subroutine srtm_reftra(KLEV, KMODTS, LDRTCHK, PGG, PRMUZ, PTAU, PW, PREF, PREFD, PTRA, PTRAD)
Definition: srtm_reftra.F90:10
parsrtm::jplay
integer(kind=jpim), parameter jplay
Definition: parsrtm.F90:19
dimphy::klev
integer, save klev
Definition: dimphy.F90:7
parkind1::jprb
integer, parameter jprb
Definition: parkind1.F90:31
x
!$Header!c c INCLUDE fxyprim h c c c Fonctions in line c c REAL fyprim REAL rj c c il faut la calculer avant d appeler ces fonctions c c c Fonctions a changer selon x(x) et y(y) choisis.c-----------------------------------------------------------------c c.....ici
yoerdu::replog
real(kind=jprb) replog
Definition: yoerdu.F90:19
yomhook::lhook
logical lhook
Definition: yomhook.F90:12
yomhook::dr_hook
subroutine dr_hook(CDNAME, KSWITCH, PKEY)
Definition: yomhook.F90:17
parkind1::jpim
integer, parameter jpim
Definition: parkind1.F90:13
yoerdu
Definition: yoerdu.F90:1