source: LMDZ5/branches/IPSLCM5A2.1_ISO/libf/phyiso/rrtm/rrtm_rtrn1a_140gp.F90

Last change on this file was 3331, checked in by acozic, 7 years ago

Add modification for isotopes

  • Property svn:executable set to *
File size: 30.8 KB
RevLine 
[3331]1SUBROUTINE RRTM_RTRN1A_140GP (KLEV,K_ISTART,K_IEND,K_ICLDLYR,P_CLDFRAC,P_TAUCLD,P_ABSS1,&
2 & P_OD,P_TAUSF1,P_CLFNET,P_CLHTR,P_FNET,P_HTR,P_TOTDFLUC,P_TOTDFLUX,P_TOTUFLUC,P_TOTUFLUX,&
3 & P_TAVEL,PZ,P_TZ,P_TBOUND,PFRAC,P_SEMISS,P_SEMISLW,K_IREFLECT) 
4
5!     Reformatted for F90 by JJMorcrette, ECMWF, 980714
6!     Speed-up by D.Salmond, ECMWF, 9907
7!     Bug-fix by M.J. Iacono, AER, Inc., 9911
8!     Bug-fix by JJMorcrette, ECMWF, 991209 (RAT1, RAT2 initialization)
9!     Speed-up by D. Salmond, ECMWF, 9912
10!     Bug-fix by JJMorcrette, ECMWF, 0005 (extrapolation T<160K)
11!     Speed-up by D. Salmond, ECMWF, 000515
12
13!-* This program calculates the upward fluxes, downward fluxes,
14!   and heating rates for an arbitrary atmosphere.  The input to
15!   this program is the atmospheric profile and all Planck function
16!   information.  First-order "numerical" quadrature is used for the
17!   angle integration, i.e. only one exponential is computed per layer
18!   per g-value per band.  Cloud overlap is treated with a generalized
19!   maximum/random method in which adjacent cloud layers are treated
20!   with maximum overlap, and non-adjacent cloud groups are treated
21!   with random overlap.  For adjacent cloud layers, cloud information
22!   is carried from the previous two layers.
23
24USE PARKIND1  ,ONLY : JPIM     ,JPRB
25USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
26
27USE PARRRTM  , ONLY : JPBAND   ,JPGPT   ,JPLAY
28USE YOERRTAB , ONLY : BPADE
29USE YOERRTWN , ONLY : TOTPLNK  ,DELWAVE
30USE YOERRTFTR, ONLY : NGB
31
32IMPLICIT NONE
33
34INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV
35INTEGER(KIND=JPIM),INTENT(IN)    :: K_ISTART
36INTEGER(KIND=JPIM),INTENT(IN)    :: K_IEND
37INTEGER(KIND=JPIM),INTENT(IN)    :: K_ICLDLYR(JPLAY) ! Cloud indicator
38REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDFRAC(JPLAY) ! Cloud fraction
39REAL(KIND=JPRB)                  :: Z_CLDFRAC(JPLAY) ! Cloud fraction
40REAL(KIND=JPRB)   ,INTENT(IN)    :: P_TAUCLD(JPLAY,JPBAND) ! Spectral optical thickness
41REAL(KIND=JPRB)   ,INTENT(IN)    :: P_ABSS1(JPGPT*JPLAY)
42REAL(KIND=JPRB)   ,INTENT(IN)    :: P_OD(JPGPT,JPLAY)
43REAL(KIND=JPRB)   ,INTENT(IN)    :: P_TAUSF1(JPGPT*JPLAY)
44REAL(KIND=JPRB)                  :: P_CLFNET(0:JPLAY) ! Argument NOT used
45REAL(KIND=JPRB)                  :: P_CLHTR(0:JPLAY) ! Argument NOT used
46REAL(KIND=JPRB)                  :: P_FNET(0:JPLAY) ! Argument NOT used
47REAL(KIND=JPRB)                  :: P_HTR(0:JPLAY) ! Argument NOT used
48REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TOTDFLUC(0:JPLAY)
49REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TOTDFLUX(0:JPLAY)
50REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TOTUFLUC(0:JPLAY)
51REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_TOTUFLUX(0:JPLAY)
52REAL(KIND=JPRB)   ,INTENT(IN)    :: P_TAVEL(JPLAY)
53REAL(KIND=JPRB)                  :: PZ(0:JPLAY) ! Argument NOT used
54REAL(KIND=JPRB)   ,INTENT(IN)    :: P_TZ(0:JPLAY)
55REAL(KIND=JPRB)   ,INTENT(IN)    :: P_TBOUND
56REAL(KIND=JPRB)   ,INTENT(IN)    :: PFRAC(JPGPT,JPLAY)
57REAL(KIND=JPRB)   ,INTENT(IN)    :: P_SEMISS(JPBAND)
58REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_SEMISLW
59INTEGER(KIND=JPIM)               :: K_IREFLECT ! Argument NOT used
60!- from PROFILE             
61!- from SP             
62!- from SURFACE             
63INTEGER(KIND=JPIM) :: INDLAY(JPLAY),INDLEV(0:JPLAY)
64
65REAL(KIND=JPRB) :: Z_BBU1(JPGPT*JPLAY),Z_BBUTOT1(JPGPT*JPLAY)
66REAL(KIND=JPRB) :: Z_TLAYFRAC(JPLAY),Z_TLEVFRAC(0:JPLAY)
67REAL(KIND=JPRB) :: Z_BGLEV(JPGPT)
68!-- DS_000515
69REAL(KIND=JPRB) :: Z_PLVL(JPBAND+1,0:JPLAY),Z_PLAY(JPBAND+1,0:JPLAY),Z_WTNUM(3)
70!-- DS_000515
71REAL(KIND=JPRB) :: Z_ODCLDNW(JPGPT,JPLAY)
72REAL(KIND=JPRB) :: Z_SEMIS(JPGPT),Z_RADUEMIT(JPGPT)
73
74REAL(KIND=JPRB) :: Z_RADCLRU1(JPGPT) ,Z_RADCLRD1(JPGPT)
75REAL(KIND=JPRB) :: Z_RADLU1(JPGPT)   ,Z_RADLD1(JPGPT)
76!-- DS_000515
77REAL(KIND=JPRB) :: Z_TRNCLD(JPLAY,JPBAND+1)
78!-- DS_000515
79REAL(KIND=JPRB) :: Z_ABSCLDNW(JPGPT,JPLAY)
80REAL(KIND=JPRB) :: Z_ATOT1(JPGPT*JPLAY)
81
82REAL(KIND=JPRB) :: Z_SURFEMIS(JPBAND),Z_PLNKEMIT(JPBAND)
83
84! dimension of arrays required for cloud overlap calculations
85
86REAL(KIND=JPRB) :: Z_CLRRADU(jpgpt),Z_CLDRADU(jpgpt),Z_OLDCLD(jpgpt)
87REAL(KIND=JPRB) :: Z_OLDCLR(jpgpt),Z_RAD(jpgpt),Z_FACCLD1(jplay+1),Z_FACCLD2(jplay+1)
88REAL(KIND=JPRB) :: Z_FACCLR1(jplay+1),Z_FACCLR2(jplay+1)
89REAL(KIND=JPRB) :: Z_FACCMB1(jplay+1),Z_FACCMB2(jplay+1)
90REAL(KIND=JPRB) :: Z_FACCLD1D(0:jplay),Z_FACCLD2D(0:jplay),Z_FACCLR1D(0:jplay)
91REAL(KIND=JPRB) :: Z_FACCLR2D(0:jplay),Z_FACCMB1D(0:jplay),Z_FACCMB2D(0:jplay)
92REAL(KIND=JPRB) :: Z_CLRRADD(jpgpt),Z_CLDRADD(jpgpt)
93INTEGER(KIND=JPIM) :: istcld(jplay+1),istcldd(0:jplay)
94!******
95
96!REAL_B :: ZPLVL(JPGPT+1,JPLAY)  ,ZPLAY(JPGPT+1,JPLAY)
97!REAL_B :: ZTRNCLD(JPGPT+1,JPLAY),ZTAUCLD(JPGPT+1,JPLAY)
98
99INTEGER(KIND=JPIM) :: IBAND, ICLDDN, IENT, INDBOUND, INDEX, IPR, I_LAY, I_LEV, I_NBI
100
101REAL(KIND=JPRB) :: Z_BBD, Z_BBDTOT, Z_BGLAY, Z_CLDSRC, Z_DBDTLAY, Z_DBDTLEV,&
102 & Z_DELBGDN, Z_DELBGUP, Z_DRAD1, Z_DRADCL1, Z_FACTOT1, &
103 & Z_FMAX, Z_FMIN, Z_GASSRC, Z_ODSM, Z_PLANKBND, Z_RADCLD, Z_RADD, Z_RADMOD, Z_RAT1, Z_RAT2, Z_SUMPL, &
104 & Z_SUMPLEM, Z_TBNDFRAC, Z_TRNS, Z_TTOT, Z_URAD1, Z_URADCL1, ZEXTAU 
105REAL(KIND=JPRB) :: ZHOOK_HANDLE
106
107
108
109REAL(KIND=JPRB)                  :: CLFNET(0:JPLAY) ! Argument NOT used
110REAL(KIND=JPRB)                  :: CLHTR(0:JPLAY) ! Argument NOT used
111REAL(KIND=JPRB)                  :: FNET(0:JPLAY) ! Argument NOT used
112REAL(KIND=JPRB)                  :: HTR(0:JPLAY) ! Argument NOT used
113
114
115
116!--------------------------------------------------------------------------
117! Input
118!  JPLAY                 ! Maximum number of model layers
119!  JPGPT                 ! Total number of g-point subintervals
120!  JPBAND                ! Number of longwave spectral bands
121!  SECANG                ! Diffusivity angle
122!  WTNUM                 ! Weight for radiance to flux conversion
123!  KLEV                  ! Number of model layers
124!  PAVEL(JPLAY)          ! Mid-layer pressures (hPa)
125!  PZ(0:JPLAY)           ! Interface pressures (hPa)
126!  TAVEL(JPLAY)          ! Mid-layer temperatures (K)
127!  TZ(0:JPLAY)           ! Interface temperatures (K)
128!  TBOUND                ! Surface temperature
129!  CLDFRAC(JPLAY)        ! Layer cloud fraction
130!  TAUCLD(JPLAY,JPBAND)  ! Layer cloud optical thickness
131!  ITR
132!  PFRAC(JPGPT,JPLAY)    ! Planck function fractions
133!  ICLDLYR(JPLAY)        ! Flag for cloudy layers
134!  ICLD                  ! Flag for cloudy column
135!  IREFLECT              ! Flag for specular reflection
136!  SEMISS(JPBAND)        ! Surface spectral emissivity
137!  BPADE                 ! Pade constant
138!  OD                    ! Clear-sky optical thickness
139!  TAUSF1                !
140!  ABSS1                 ! 
141
142!  ABSS(JPGPT*JPLAY)     !
143!  ABSCLD(JPLAY)         !
144!  ATOT(JPGPT*JPLAY)     !
145!  ODCLR(JPGPT,JPLAY)    !
146!  ODCLD(JPBAND,JPLAY)   !
147!  EFCLFR1(JPBAND,JPLAY) ! Effective cloud fraction
148!  RADLU(JPGPT)          ! Upward radiance
149!  URAD                  ! Spectrally summed upward radiance
150!  RADCLRU(JPGPT)        ! Clear-sky upward radiance
151!  CLRURAD               ! Spectrally summed clear-sky upward radiance
152!  RADLD(JPGPT)          ! Downward radiance
153!  DRAD                  ! Spectrally summed downward radiance
154!  RADCLRD(JPGPT)        ! Clear-sky downward radiance
155!  CLRDRAD               ! Spectrally summed clear-sky downward radiance
156
157! Output
158!  TOTUFLUX(0:JPLAY)     ! Upward longwave flux
159!  TOTDFLUX(0:JPLAY)     ! Downward longwave flux
160!  TOTUFLUC(0:JPLAY)     ! Clear-sky upward longwave flux
161!  TOTDFLUC(0:JPLAY)     ! Clear-sky downward longwave flux
162
163! Maximum/Random cloud overlap variables
164! for upward radiaitve transfer
165!  FACCLR2  fraction of clear radiance from previous layer that needs to
166!           be switched to cloudy stream
167!  FACCLR1  fraction of the radiance that had been switched in the previous
168!           layer from cloudy to clear that needs to be switched back to
169!           cloudy in the current layer
170!  FACCLD2  fraction of cloudy radiance from previous layer that needs to
171!           be switched to clear stream
172!           be switched to cloudy stream
173!  FACCLD1  fraction of the radiance that had been switched in the previous
174!           layer from clear to cloudy that needs to be switched back to
175!           clear in the current layer
176! for downward radiaitve transfer
177!  FACCLR2D fraction of clear radiance from previous layer that needs to
178!           be switched to cloudy stream
179!  FACCLR1D fraction of the radiance that had been switched in the previous
180!           layer from cloudy to clear that needs to be switched back to
181!           cloudy in the current layer
182!  FACCLD2D fraction of cloudy radiance from previous layer that needs to
183!           be switched to clear stream
184!           be switched to cloudy stream
185!  FACCLD1D fraction of the radiance that had been switched in the previous
186!           layer from clear to cloudy that needs to be switched back to
187!           clear in the current layer
188
189!--------------------------------------------------------------------------
190
191! CORRECTION PROVISOIRE BUG POTENTIEL MPLFH
192! on initialise le niveau klev+1 de p_cldfrac, tableau surdimensionne
193! a 100 mais apparemment non initialise en klev+1
194Z_CLDFRAC(1:KLEV)=P_CLDFRAC(1:KLEV)
195Z_CLDFRAC(KLEV+1)=0.0_JPRB
196IF (LHOOK) CALL DR_HOOK('RRTM_RTRN1A_140GP',0,ZHOOK_HANDLE)
197Z_WTNUM(1)=0.5_JPRB
198Z_WTNUM(2)=0.0_JPRB
199Z_WTNUM(3)=0.0_JPRB
200
201DO I_LAY = 0, KLEV
202ENDDO
203!-start JJM_000511
204IF (P_TBOUND < 339._JPRB .AND. P_TBOUND >= 160._JPRB ) THEN
205  INDBOUND = P_TBOUND - 159._JPRB
206  Z_TBNDFRAC = P_TBOUND - INT(P_TBOUND)
207ELSEIF (P_TBOUND >= 339._JPRB ) THEN
208  INDBOUND = 180
209  Z_TBNDFRAC = P_TBOUND - 339._JPRB
210ELSEIF (P_TBOUND < 160._JPRB ) THEN
211  INDBOUND = 1
212  Z_TBNDFRAC = P_TBOUND - 160._JPRB
213ENDIF 
214!-end JJM_000511
215 
216DO I_LAY = 0, KLEV
217  P_TOTUFLUC(I_LAY) = 0.0_JPRB
218  P_TOTDFLUC(I_LAY) = 0.0_JPRB
219  P_TOTUFLUX(I_LAY) = 0.0_JPRB
220  P_TOTDFLUX(I_LAY) = 0.0_JPRB
221!-start JJM_000511
222  IF (P_TZ(I_LAY) < 339._JPRB .AND. P_TZ(I_LAY) >= 160._JPRB ) THEN
223    INDLEV(I_LAY) = P_TZ(I_LAY) - 159._JPRB
224    Z_TLEVFRAC(I_LAY) = P_TZ(I_LAY) - INT(P_TZ(I_LAY))
225  ELSEIF (P_TZ(I_LAY) >= 339._JPRB ) THEN
226    INDLEV(I_LAY) = 180
227    Z_TLEVFRAC(I_LAY) = P_TZ(I_LAY) - 339._JPRB
228  ELSEIF (P_TZ(I_LAY) < 160._JPRB ) THEN
229    INDLEV(I_LAY) = 1
230    Z_TLEVFRAC(I_LAY) = P_TZ(I_LAY) - 160._JPRB
231  ENDIF   
232!-end JJM_000511
233ENDDO
234
235!_start_jjm 991209
236DO I_LEV=0,KLEV
237  Z_FACCLD1(I_LEV+1) = 0.0_JPRB
238  Z_FACCLD2(I_LEV+1) = 0.0_JPRB
239  Z_FACCLR1(I_LEV+1) = 0.0_JPRB
240  Z_FACCLR2(I_LEV+1) = 0.0_JPRB
241  Z_FACCMB1(I_LEV+1) = 0.0_JPRB
242  Z_FACCMB2(I_LEV+1) = 0.0_JPRB
243  Z_FACCLD1D(I_LEV) = 0.0_JPRB
244  Z_FACCLD2D(I_LEV) = 0.0_JPRB
245  Z_FACCLR1D(I_LEV) = 0.0_JPRB
246  Z_FACCLR2D(I_LEV) = 0.0_JPRB
247  Z_FACCMB1D(I_LEV) = 0.0_JPRB
248  Z_FACCMB2D(I_LEV) = 0.0_JPRB
249ENDDO 
250
251Z_RAT1 = 0.0_JPRB
252Z_RAT2 = 0.0_JPRB
253
254!_end_jjm 991209
255
256Z_SUMPL   = 0.0_JPRB
257Z_SUMPLEM = 0.0_JPRB
258
259ISTCLD(1) = 1
260ISTCLDD(KLEV) = 1
261
262DO I_LEV = 1, KLEV
263!-- DS_000515
264!-start JJM_000511
265  IF (P_TAVEL(I_LEV) < 339._JPRB .AND. P_TAVEL(I_LEV) >= 160._JPRB ) THEN
266    INDLAY(I_LEV) = P_TAVEL(I_LEV) - 159._JPRB
267    Z_TLAYFRAC(I_LEV) = P_TAVEL(I_LEV) - INT(P_TAVEL(I_LEV))
268  ELSEIF (P_TAVEL(I_LEV) >= 339._JPRB ) THEN
269    INDLAY(I_LEV) = 180
270    Z_TLAYFRAC(I_LEV) = P_TAVEL(I_LEV) - 339._JPRB
271  ELSEIF (P_TAVEL(I_LEV) < 160._JPRB ) THEN
272    INDLAY(I_LEV) = 1
273    Z_TLAYFRAC(I_LEV) = P_TAVEL(I_LEV) - 160._JPRB
274  ENDIF 
275!-end JJM_000511
276ENDDO
277!-- DS_000515
278
279!-- DS_000515
280!OCL SCALAR
281
282DO I_LEV = 1, KLEV
283  IF (K_ICLDLYR(I_LEV) == 1) THEN
284
285!mji   
286    ISTCLD(I_LEV+1) = 0
287    IF (I_LEV  ==  KLEV) THEN
288      Z_FACCLD1(I_LEV+1) = 0.0_JPRB
289      Z_FACCLD2(I_LEV+1) = 0.0_JPRB
290      Z_FACCLR1(I_LEV+1) = 0.0_JPRB
291      Z_FACCLR2(I_LEV+1) = 0.0_JPRB
292!-- DS_000515     
293!SB debug >>
294     Z_FACCMB1(I_LEV+1) =0.0_JPRB
295     Z_FACCMB2(I_LEV+1) =0.0_JPRB
296!SB debug <<
297!mji      ISTCLD(LEV+1) = _ZERO_
298    ELSEIF (Z_CLDFRAC(I_LEV+1)  >=  Z_CLDFRAC(I_LEV)) THEN
299      Z_FACCLD1(I_LEV+1) = 0.0_JPRB
300      Z_FACCLD2(I_LEV+1) = 0.0_JPRB
301      IF (ISTCLD(I_LEV)  ==  1) THEN
302!mji        ISTCLD(LEV+1) = 0
303        Z_FACCLR1(I_LEV+1) = 0.0_JPRB
304!mji       
305        Z_FACCLR2(I_LEV+1) = 0.0_JPRB
306        IF (Z_CLDFRAC(I_LEV) < 1.0_JPRB) THEN
307          Z_FACCLR2(I_LEV+1) = (Z_CLDFRAC(I_LEV+1)-Z_CLDFRAC(I_LEV))/&
308           & (1.0_JPRB-Z_CLDFRAC(I_LEV)) 
309        ENDIF 
310!SB debug >>
311      Z_FACCLR2(I_LEV) = 0.0_JPRB
312      Z_FACCLD2(I_LEV) = 0.0_JPRB
313!SB debug <<
314      ELSE
315        Z_FMAX = MAX(Z_CLDFRAC(I_LEV),Z_CLDFRAC(I_LEV-1))
316!mji
317        IF (Z_CLDFRAC(I_LEV+1)  >  Z_FMAX) THEN
318          Z_FACCLR1(I_LEV+1) = Z_RAT2
319          Z_FACCLR2(I_LEV+1) = (Z_CLDFRAC(I_LEV+1)-Z_FMAX)/(1.0_JPRB-Z_FMAX)
320!mji         
321        ELSEIF (Z_CLDFRAC(I_LEV+1) < Z_FMAX) THEN
322          Z_FACCLR1(I_LEV+1) = (Z_CLDFRAC(I_LEV+1)-Z_CLDFRAC(I_LEV))/&
323           & (Z_CLDFRAC(I_LEV-1)-Z_CLDFRAC(I_LEV)) 
324          Z_FACCLR2(I_LEV+1) = 0.0_JPRB
325!mji
326        ELSE
327          Z_FACCLR1(I_LEV+1) = Z_RAT2 
328          Z_FACCLR2(I_LEV+1) = 0.0_JPRB
329        ENDIF
330      ENDIF
331      IF (Z_FACCLR1(I_LEV+1) > 0.0_JPRB .OR. Z_FACCLR2(I_LEV+1) > 0.0_JPRB) THEN
332        Z_RAT1 = 1.0_JPRB
333        Z_RAT2 = 0.0_JPRB
334!SB debug >>
335!      ENDIF
336      ELSE
337        Z_RAT1 = 0.0_JPRB
338        Z_RAT2 = 0.0_JPRB
339      ENDIF
340!SB debug <<
341    ELSE
342      Z_FACCLR1(I_LEV+1) = 0.0_JPRB
343      Z_FACCLR2(I_LEV+1) = 0.0_JPRB
344      IF (ISTCLD(I_LEV)  ==  1) THEN
345!mji        ISTCLD(LEV+1) = 0
346        Z_FACCLD1(I_LEV+1) = 0.0_JPRB
347        Z_FACCLD2(I_LEV+1) = (Z_CLDFRAC(I_LEV)-Z_CLDFRAC(I_LEV+1))/Z_CLDFRAC(I_LEV)
348!SB debug >>
349        Z_FACCLR2(I_LEV) = 0.0_JPRB
350        Z_FACCLD2(I_LEV) = 0.0_JPRB
351!SB debug <<
352      ELSE
353        Z_FMIN = MIN(Z_CLDFRAC(I_LEV),Z_CLDFRAC(I_LEV-1))
354        IF (Z_CLDFRAC(I_LEV+1)  <=  Z_FMIN) THEN
355          Z_FACCLD1(I_LEV+1) = Z_RAT1
356          Z_FACCLD2(I_LEV+1) = (Z_FMIN-Z_CLDFRAC(I_LEV+1))/Z_FMIN
357        ELSE
358          Z_FACCLD1(I_LEV+1) = (Z_CLDFRAC(I_LEV)-Z_CLDFRAC(I_LEV+1))/&
359           & (Z_CLDFRAC(I_LEV)-Z_FMIN) 
360          Z_FACCLD2(I_LEV+1) = 0.0_JPRB
361        ENDIF
362      ENDIF
363      IF (Z_FACCLD1(I_LEV+1) > 0.0_JPRB .OR. Z_FACCLD2(I_LEV+1) > 0.0_JPRB) THEN
364        Z_RAT1 = 0.0_JPRB
365        Z_RAT2 = 1.0_JPRB
366!SB debug >>
367!      ENDIF
368      ELSE
369        Z_RAT1 = 0.0_JPRB
370        Z_RAT2 = 0.0_JPRB
371      ENDIF
372!SB debug <<
373    ENDIF
374!fcc
375
376!SB debug >>
377!    IF (I_LEV == 1) THEN
378!      Z_FACCMB1(I_LEV+1) = 0.
379!      Z_FACCMB2(I_LEV+1) = Z_FACCLD1(I_LEV+1) * Z_FACCLR2(I_LEV)
380!    ELSE
381!      Z_FACCMB1(I_LEV+1) = Z_FACCLR1(I_LEV+1) * Z_FACCLD2(I_LEV) *Z_CLDFRAC(I_LEV-1)
382!      Z_FACCMB2(I_LEV+1) = Z_FACCLD1(I_LEV+1) * Z_FACCLR2(I_LEV) *&
383!       & (1.0_JPRB - Z_CLDFRAC(I_LEV-1))   
384!    ENDIF
385     if(istcld(i_lev).ne.1.and.i_lev.ne.1) then
386        z_faccmb1(i_lev+1) = max(0.,min(z_cldfrac(i_lev+1)-z_cldfrac(i_lev), &
387               z_cldfrac(i_lev-1)-z_cldfrac(i_lev)))
388        z_faccmb2(i_lev+1) = max(0.,min(z_cldfrac(i_lev)-z_cldfrac(i_lev+1), &
389               z_cldfrac(i_lev)-z_cldfrac(i_lev-1)))
390     endif
391!SB debug <<
392!end fcc
393  ELSE
394!-- DS_000515
395    ISTCLD(I_LEV+1) = 1
396  ENDIF
397ENDDO
398
399!_start_jjm 991209
400Z_RAT1 = 0.0_JPRB
401Z_RAT2 = 0.0_JPRB
402!_end_jjm 991209
403
404!-- DS_000515
405!OCL SCALAR
406
407DO I_LEV = KLEV, 1, -1
408  IF (K_ICLDLYR(I_LEV) == 1) THEN
409!mji
410    ISTCLDD(I_LEV-1) = 0 
411    IF (I_LEV  ==  1) THEN
412      Z_FACCLD1D(I_LEV-1) = 0.0_JPRB
413      Z_FACCLD2D(I_LEV-1) = 0.0_JPRB
414      Z_FACCLR1D(I_LEV-1) = 0.0_JPRB
415      Z_FACCLR2D(I_LEV-1) = 0.0_JPRB
416      Z_FACCMB1D(I_LEV-1) = 0.0_JPRB
417      Z_FACCMB2D(I_LEV-1) = 0.0_JPRB
418!mji      ISTCLDD(LEV-1) = _ZERO_
419    ELSEIF (Z_CLDFRAC(I_LEV-1)  >=  Z_CLDFRAC(I_LEV)) THEN
420      Z_FACCLD1D(I_LEV-1) = 0.0_JPRB
421      Z_FACCLD2D(I_LEV-1) = 0.0_JPRB
422      IF (ISTCLDD(I_LEV)  ==  1) THEN
423!mji        ISTCLDD(LEV-1) = 0
424        Z_FACCLR1D(I_LEV-1) = 0.0_JPRB
425        Z_FACCLR2D(I_LEV-1) = 0.0_JPRB
426        IF (Z_CLDFRAC(I_LEV) < 1.0_JPRB) THEN
427          Z_FACCLR2D(I_LEV-1) = (Z_CLDFRAC(I_LEV-1)-Z_CLDFRAC(I_LEV))/&
428           & (1.0_JPRB-Z_CLDFRAC(I_LEV)) 
429        ENDIF
430!SB debug >>
431       z_facclr2d(i_lev)=0.0_JPRB       
432       z_faccld2d(i_lev)=0.0_JPRB
433!SB debug <<
434      ELSE
435        Z_FMAX = MAX(Z_CLDFRAC(I_LEV),Z_CLDFRAC(I_LEV+1))
436!mji
437        IF (Z_CLDFRAC(I_LEV-1)  >  Z_FMAX) THEN
438          Z_FACCLR1D(I_LEV-1) = Z_RAT2
439          Z_FACCLR2D(I_LEV-1) = (Z_CLDFRAC(I_LEV-1)-Z_FMAX)/(1.0_JPRB-Z_FMAX)
440!mji
441        ELSEIF (Z_CLDFRAC(I_LEV-1) < Z_FMAX) THEN
442          Z_FACCLR1D(I_LEV-1) = (Z_CLDFRAC(I_LEV-1)-Z_CLDFRAC(I_LEV))/&
443           & (Z_CLDFRAC(I_LEV+1)-Z_CLDFRAC(I_LEV)) 
444          Z_FACCLR2D(I_LEV-1) = 0.0_JPRB
445!mji
446        ELSE         
447          Z_FACCLR1D(I_LEV-1) = Z_RAT2
448          Z_FACCLR2D(I_LEV-1) = 0.0_JPRB
449        ENDIF
450      ENDIF
451      IF (Z_FACCLR1D(I_LEV-1) > 0.0_JPRB .OR. Z_FACCLR2D(I_LEV-1) > 0.0_JPRB)THEN
452        Z_RAT1 = 1.0_JPRB
453        Z_RAT2 = 0.0_JPRB
454!SB debug >>
455!      ENDIF
456      else
457        Z_RAT1 = 0.0_JPRB
458        Z_RAT2 = 0.0_JPRB
459      endif       
460!SB debug <<
461    ELSE
462      Z_FACCLR1D(I_LEV-1) = 0.0_JPRB
463      Z_FACCLR2D(I_LEV-1) = 0.0_JPRB
464      IF (ISTCLDD(I_LEV)  ==  1) THEN
465!mji        ISTCLDD(LEV-1) = 0
466        Z_FACCLD1D(I_LEV-1) = 0.0_JPRB
467        Z_FACCLD2D(I_LEV-1) = (Z_CLDFRAC(I_LEV)-Z_CLDFRAC(I_LEV-1))/Z_CLDFRAC(I_LEV)
468!SB debug >>
469        z_facclr2d(i_lev)=0.0_JPRB
470        z_faccld2d(i_lev)=0.0_JPRB
471!SB debug <<
472      ELSE
473        Z_FMIN = MIN(Z_CLDFRAC(I_LEV),Z_CLDFRAC(I_LEV+1))
474        IF (Z_CLDFRAC(I_LEV-1)  <=  Z_FMIN) THEN
475          Z_FACCLD1D(I_LEV-1) = Z_RAT1
476          Z_FACCLD2D(I_LEV-1) = (Z_FMIN-Z_CLDFRAC(I_LEV-1))/Z_FMIN
477        ELSE
478          Z_FACCLD1D(I_LEV-1) = (Z_CLDFRAC(I_LEV)-Z_CLDFRAC(I_LEV-1))/&
479           & (Z_CLDFRAC(I_LEV)-Z_FMIN) 
480          Z_FACCLD2D(I_LEV-1) = 0.0_JPRB
481        ENDIF
482      ENDIF
483      IF (Z_FACCLD1D(I_LEV-1) > 0.0_JPRB .OR. Z_FACCLD2D(I_LEV-1) > 0.0_JPRB)THEN
484        Z_RAT1 = 0.0_JPRB
485        Z_RAT2 = 1.0_JPRB
486!SB debug >>
487!      ENDIF
488      ELSE
489        Z_RAT1 = 0.0_JPRB
490        Z_RAT2 = 0.0_JPRB
491      ENDIF
492!SB debug <<
493    ENDIF
494!SB debug >>
495!    Z_FACCMB1D(I_LEV-1) = Z_FACCLR1D(I_LEV-1) * Z_FACCLD2D(I_LEV) *Z_CLDFRAC(I_LEV+1)
496!    Z_FACCMB2D(I_LEV-1) = Z_FACCLD1D(I_LEV-1) * Z_FACCLR2D(I_LEV) *&
497!     & (1.0_JPRB - Z_CLDFRAC(I_LEV+1)) 
498    if (istcldd(i_lev).ne.1.and.i_lev.ne.1) then
499       z_faccmb1d(i_lev-1) = max(0.,min(z_cldfrac(i_lev+1)-z_cldfrac(i_lev), &
500                            z_cldfrac(i_lev-1)-z_cldfrac(i_lev)))
501       z_faccmb2d(i_lev-1) = max(0.,min(z_cldfrac(i_lev)-z_cldfrac(i_lev+1), &
502                    z_cldfrac(i_lev)-z_cldfrac(i_lev-1)))
503    endif
504!SB debug <<
505  ELSE
506    ISTCLDD(I_LEV-1) = 1
507  ENDIF
508ENDDO
509
510!- Loop over frequency bands.
511
512DO IBAND = K_ISTART, K_IEND
513  Z_DBDTLEV = TOTPLNK(INDBOUND+1,IBAND)-TOTPLNK(INDBOUND,IBAND)
514  Z_PLANKBND = DELWAVE(IBAND) * (TOTPLNK(INDBOUND,IBAND) + Z_TBNDFRAC * Z_DBDTLEV)
515  Z_DBDTLEV = TOTPLNK(INDLEV(0)+1,IBAND) -TOTPLNK(INDLEV(0),IBAND)
516!-- DS_000515
517  Z_PLVL(IBAND,0) = DELWAVE(IBAND)&
518   & * (TOTPLNK(INDLEV(0),IBAND) + Z_TLEVFRAC(0)*Z_DBDTLEV) 
519
520  Z_SURFEMIS(IBAND) = P_SEMISS(IBAND)
521  Z_PLNKEMIT(IBAND) = Z_SURFEMIS(IBAND) * Z_PLANKBND
522  Z_SUMPLEM  = Z_SUMPLEM + Z_PLNKEMIT(IBAND)
523  Z_SUMPL    = Z_SUMPL   + Z_PLANKBND
524!--DS
525ENDDO
526!---
527
528!-- DS_000515
529DO I_LEV = 1, KLEV
530  DO IBAND = K_ISTART, K_IEND
531! print *,'RTRN1A: I_LEV JPLAY IBAND INDLAY',I_LEV,JPLAY,IBAND,INDLAY(I_LEV)
532!----             
533!- Calculate the integrated Planck functions for at the
534!  level and layer temperatures.
535!  Compute cloud transmittance for cloudy layers.
536    Z_DBDTLEV = TOTPLNK(INDLEV(I_LEV)+1,IBAND) - TOTPLNK(INDLEV(I_LEV),IBAND)
537    Z_DBDTLAY = TOTPLNK(INDLAY(I_LEV)+1,IBAND) - TOTPLNK(INDLAY(I_LEV),IBAND)
538!-- DS_000515
539    Z_PLAY(IBAND,I_LEV) = DELWAVE(IBAND)&
540     & *(TOTPLNK(INDLAY(I_LEV),IBAND)+Z_TLAYFRAC(I_LEV)*Z_DBDTLAY) 
541    Z_PLVL(IBAND,I_LEV) = DELWAVE(IBAND)&
542     & *(TOTPLNK(INDLEV(I_LEV),IBAND)+Z_TLEVFRAC(I_LEV)*Z_DBDTLEV) 
543    IF (K_ICLDLYR(I_LEV) > 0) THEN
544      ZEXTAU = MIN( P_TAUCLD(I_LEV,IBAND), 200._JPRB)
545      Z_TRNCLD(I_LEV,IBAND) = EXP( -ZEXTAU )
546    ENDIF
547!-- DS_000515
548  ENDDO
549
550ENDDO
551
552P_SEMISLW = Z_SUMPLEM / Z_SUMPL
553
554!--DS
555!O IPR = 1, JPGPT
556! NBI = NGB(IPR)
557! DO LEV =  1 , KLEV
558!-- DS_000515
559!   ZPLAY(IPR,LEV) = PLAY(LEV,NGB(IPR))
560!   ZPLVL(IPR,LEV) = PLVL(LEV-1,NGB(IPR))
561!   ZTAUCLD(IPR,LEV) = TAUCLD(LEV,NGB(IPR))
562!   ZTRNCLD(IPR,LEV) = TRNCLD(LEV,NGB(IPR))
563!-- DS_000515
564! ENDDO
565!NDDO
566!----     
567
568!- For cloudy layers, set cloud parameters for radiative transfer.
569DO I_LEV = 1, KLEV
570  IF (K_ICLDLYR(I_LEV) > 0) THEN
571    DO IPR = 1, JPGPT
572!--DS         
573!            NBI = NGB(IPR)
574      Z_ODCLDNW(IPR,I_LEV) = P_TAUCLD(I_LEV,NGB(IPR))
575      Z_ABSCLDNW(IPR,I_LEV) = 1.0_JPRB - Z_TRNCLD(I_LEV,NGB(IPR))
576!----           
577!            EFCLFRNW(IPR,LEV) = ABSCLDNW(IPR,LEV) * CLDFRAC(LEV)
578    ENDDO
579  ENDIF
580ENDDO
581
582!- Initialize for radiative transfer.
583DO IPR = 1, JPGPT
584  Z_RADCLRD1(IPR) = 0.0_JPRB
585  Z_RADLD1(IPR)   = 0.0_JPRB
586  I_NBI = NGB(IPR)
587  Z_SEMIS(IPR) = Z_SURFEMIS(I_NBI)
588  Z_RADUEMIT(IPR) = PFRAC(IPR,1) * Z_PLNKEMIT(I_NBI)
589!-- DS_000515
590  Z_BGLEV(IPR) = PFRAC(IPR,KLEV) * Z_PLVL(I_NBI,KLEV)
591ENDDO
592
593!- Downward radiative transfer.
594!  *** DRAD1 holds summed radiance for total sky stream
595!  *** DRADCL1 holds summed radiance for clear sky stream
596
597ICLDDN = 0
598DO I_LEV = KLEV, 1, -1
599  Z_DRAD1   = 0.0_JPRB
600  Z_DRADCL1 = 0.0_JPRB
601
602  IF (K_ICLDLYR(I_LEV) == 1) THEN
603
604!  *** Cloudy layer
605    ICLDDN = 1
606    IENT = JPGPT * (I_LEV-1)
607    DO IPR = 1, JPGPT
608      INDEX = IENT + IPR
609!--DS           
610!            NBI = NGB(IPR)
611      Z_BGLAY = PFRAC(IPR,I_LEV) * Z_PLAY(NGB(IPR),I_LEV)
612!----           
613      Z_DELBGUP     = Z_BGLEV(IPR) - Z_BGLAY
614      Z_BBU1(INDEX) = Z_BGLAY + P_TAUSF1(INDEX) * Z_DELBGUP
615!--DS           
616      Z_BGLEV(IPR) = PFRAC(IPR,I_LEV) * Z_PLVL(NGB(IPR),I_LEV-1)
617!----           
618      Z_DELBGDN = Z_BGLEV(IPR) - Z_BGLAY
619      Z_BBD = Z_BGLAY + P_TAUSF1(INDEX) * Z_DELBGDN
620!- total-sky downward flux         
621      Z_ODSM = P_OD(IPR,I_LEV) + Z_ODCLDNW(IPR,I_LEV)
622      Z_FACTOT1 = Z_ODSM / (BPADE + Z_ODSM)
623      Z_BBUTOT1(INDEX) = Z_BGLAY + Z_FACTOT1 * Z_DELBGUP
624      Z_ATOT1(INDEX) = P_ABSS1(INDEX) + Z_ABSCLDNW(IPR,I_LEV)&
625       & - P_ABSS1(INDEX) * Z_ABSCLDNW(IPR,I_LEV) 
626      Z_BBDTOT = Z_BGLAY + Z_FACTOT1 * Z_DELBGDN
627      Z_GASSRC = Z_BBD * P_ABSS1(INDEX)
628!***
629      IF (ISTCLDD(I_LEV)  ==  1) THEN
630        Z_CLDRADD(IPR) = Z_CLDFRAC(I_LEV) * Z_RADLD1(IPR)
631        Z_CLRRADD(IPR) = Z_RADLD1(IPR) - Z_CLDRADD(IPR)
632        Z_OLDCLD(IPR) = Z_CLDRADD(IPR)
633        Z_OLDCLR(IPR) = Z_CLRRADD(IPR)
634        Z_RAD(IPR) = 0.0_JPRB
635      ENDIF
636      Z_TTOT = 1.0_JPRB - Z_ATOT1(INDEX)
637      Z_CLDSRC = Z_BBDTOT * Z_ATOT1(INDEX)
638     
639! Separate RT equations for clear and cloudy streams     
640      Z_CLDRADD(IPR) = Z_CLDRADD(IPR) * Z_TTOT + Z_CLDFRAC(I_LEV) * Z_CLDSRC
641      Z_CLRRADD(IPR) = Z_CLRRADD(IPR) * (1.0_JPRB-P_ABSS1(INDEX)) +&
642       & (1.0_JPRB - Z_CLDFRAC(I_LEV)) * Z_GASSRC 
643
644!  Total sky downward radiance
645      Z_RADLD1(IPR) = Z_CLDRADD(IPR) + Z_CLRRADD(IPR)
646      Z_DRAD1 = Z_DRAD1 + Z_RADLD1(IPR)
647     
648!  Clear-sky downward radiance         
649      Z_RADCLRD1(IPR) = Z_RADCLRD1(IPR)+(Z_BBD-Z_RADCLRD1(IPR))*P_ABSS1(INDEX)
650      Z_DRADCL1 = Z_DRADCL1 + Z_RADCLRD1(IPR)
651
652!* Code to account for maximum/random overlap:
653!   Performs RT on the radiance most recently switched between clear and
654!   cloudy streams
655      Z_RADMOD = Z_RAD(IPR) * (Z_FACCLR1D(I_LEV-1) * (1.0_JPRB-P_ABSS1(INDEX)) +&
656       & Z_FACCLD1D(I_LEV-1) *  Z_TTOT) - &
657       & Z_FACCMB1D(I_LEV-1) * Z_GASSRC + &
658       & Z_FACCMB2D(I_LEV-1) * Z_CLDSRC 
659       
660!   Computes what the clear and cloudy streams would have been had no
661!   radiance been switched       
662      Z_OLDCLD(IPR) = Z_CLDRADD(IPR) - Z_RADMOD
663      Z_OLDCLR(IPR) = Z_CLRRADD(IPR) + Z_RADMOD
664     
665!   Computes the radiance to be switched between clear and cloudy.     
666      Z_RAD(IPR) = -Z_RADMOD + Z_FACCLR2D(I_LEV-1)*Z_OLDCLR(IPR) -&
667       & Z_FACCLD2D(I_LEV-1)*Z_OLDCLD(IPR) 
668      Z_CLDRADD(IPR) = Z_CLDRADD(IPR) + Z_RAD(IPR)
669      Z_CLRRADD(IPR) = Z_CLRRADD(IPR) - Z_RAD(IPR)
670!***
671
672    ENDDO
673
674  ELSE
675
676!  *** Clear layer
677!  *** DRAD1 holds summed radiance for total sky stream
678!  *** DRADCL1 holds summed radiance for clear sky stream
679
680    IENT = JPGPT * (I_LEV-1)
681    IF (ICLDDN == 1) THEN
682      DO IPR = 1, JPGPT
683        INDEX = IENT + IPR
684!--DS         
685!           NBI = NGB(IPR)
686        Z_BGLAY = PFRAC(IPR,I_LEV) * Z_PLAY(NGB(IPR),I_LEV)
687!----           
688        Z_DELBGUP     = Z_BGLEV(IPR) - Z_BGLAY
689        Z_BBU1(INDEX) = Z_BGLAY + P_TAUSF1(INDEX) * Z_DELBGUP
690!--DS           
691        Z_BGLEV(IPR) = PFRAC(IPR,I_LEV) * Z_PLVL(NGB(IPR),I_LEV-1)
692!----                     
693        Z_DELBGDN = Z_BGLEV(IPR) - Z_BGLAY
694        Z_BBD = Z_BGLAY + P_TAUSF1(INDEX) * Z_DELBGDN
695       
696!- total-sky downward radiance
697        Z_RADLD1(IPR) = Z_RADLD1(IPR)+(Z_BBD-Z_RADLD1(IPR))*P_ABSS1(INDEX)
698        Z_DRAD1 = Z_DRAD1 + Z_RADLD1(IPR)
699       
700!- clear-sky downward radiance
701!-  Set clear sky stream to total sky stream as long as layers
702!-  remain clear.  Streams diverge when a cloud is reached.
703        Z_RADCLRD1(IPR) = Z_RADCLRD1(IPR)+(Z_BBD-Z_RADCLRD1(IPR))*P_ABSS1(INDEX)
704        Z_DRADCL1 = Z_DRADCL1 + Z_RADCLRD1(IPR)
705      ENDDO
706           
707    ELSE
708       
709      DO IPR = 1, JPGPT
710        INDEX = IENT + IPR
711!--DS         
712!           NBI = NGB(IPR)
713        Z_BGLAY = PFRAC(IPR,I_LEV) * Z_PLAY(NGB(IPR),I_LEV)
714!----           
715        Z_DELBGUP     = Z_BGLEV(IPR) - Z_BGLAY
716        Z_BBU1(INDEX) = Z_BGLAY + P_TAUSF1(INDEX) * Z_DELBGUP
717!--DS           
718        Z_BGLEV(IPR) = PFRAC(IPR,I_LEV) * Z_PLVL(NGB(IPR),I_LEV-1)
719!----                     
720        Z_DELBGDN = Z_BGLEV(IPR) - Z_BGLAY
721        Z_BBD = Z_BGLAY + P_TAUSF1(INDEX) * Z_DELBGDN
722!- total-sky downward flux         
723        Z_RADLD1(IPR) = Z_RADLD1(IPR)+(Z_BBD-Z_RADLD1(IPR))*P_ABSS1(INDEX)
724        Z_DRAD1 = Z_DRAD1 + Z_RADLD1(IPR)
725!- clear-sky downward flux         
726!-  Set clear sky stream to total sky stream as long as layers
727!-  remain clear.  Streams diverge when a cloud is reached.
728        Z_RADCLRD1(IPR) = Z_RADLD1(IPR)
729      ENDDO
730      Z_DRADCL1 = Z_DRAD1
731    ENDIF
732   
733  ENDIF
734
735  P_TOTDFLUC(I_LEV-1) = Z_DRADCL1 * Z_WTNUM(1)
736  P_TOTDFLUX(I_LEV-1) = Z_DRAD1   * Z_WTNUM(1)
737
738ENDDO
739
740! Spectral reflectivity and reflectance
741! Includes the contribution of spectrally varying longwave emissivity
742! and reflection from the surface to the upward radiative transfer.
743! Note: Spectral and Lambertian reflections are identical for the one
744! angle flux integration used here.
745
746Z_URAD1   = 0.0_JPRB
747Z_URADCL1 = 0.0_JPRB
748
749!start JJM_000511
750!IF (IREFLECT  ==  0) THEN
751!- Lambertian reflection.
752DO IPR = 1, JPGPT
753! Clear-sky radiance
754!    RADCLD = _TWO_ * (RADCLRD1(IPR) * WTNUM(1) )
755  Z_RADCLD = Z_RADCLRD1(IPR)
756  Z_RADCLRU1(IPR) = Z_RADUEMIT(IPR) + (1.0_JPRB - Z_SEMIS(IPR)) * Z_RADCLD
757  Z_URADCL1 = Z_URADCL1 + Z_RADCLRU1(IPR)
758
759! Total sky radiance
760!    RADD = _TWO_ * (RADLD1(IPR) * WTNUM(1) )
761  Z_RADD = Z_RADLD1(IPR)
762  Z_RADLU1(IPR) = Z_RADUEMIT(IPR) + (1.0_JPRB - Z_SEMIS(IPR)) * Z_RADD
763  Z_URAD1 = Z_URAD1 + Z_RADLU1(IPR)
764ENDDO
765P_TOTUFLUC(0) = Z_URADCL1 * 0.5_JPRB
766P_TOTUFLUX(0) = Z_URAD1 * 0.5_JPRB
767!ELSE
768!!- Specular reflection.
769!  DO IPR = 1, JPGPT
770!    RADCLU = RADUEMIT(IPR)
771!    RADCLRU1(IPR) = RADCLU + (_ONE_ - SEMIS(IPR)) * RADCLRD1(IPR)
772!    URADCL1 = URADCL1 + RADCLRU1(IPR)
773
774!    RADU = RADUEMIT(IPR)
775!    RADLU1(IPR) = RADU + (_ONE_ - SEMIS(IPR)) * RADLD1(IPR)
776!    URAD1 = URAD1 + RADLU1(IPR)
777!  ENDDO
778!  TOTUFLUC(0) = URADCL1 * WTNUM(1)
779!  TOTUFLUX(0) = URAD1   * WTNUM(1)
780!ENDIF
781
782!- Upward radiative transfer.
783!- *** URAD1 holds the summed radiance for total sky stream
784!- *** URADCL1 holds the summed radiance for clear sky stream
785DO I_LEV = 1, KLEV
786  Z_URAD1   = 0.0_JPRB
787  Z_URADCL1 = 0.0_JPRB
788
789! Check flag for cloud in current layer
790  IF (K_ICLDLYR(I_LEV) == 1) THEN
791
792!- *** Cloudy layer
793    IENT = JPGPT * (I_LEV-1)
794    DO IPR = 1, JPGPT
795      INDEX = IENT + IPR
796!- total-sky upward flux         
797      Z_GASSRC = Z_BBU1(INDEX) * P_ABSS1(INDEX)
798
799!- If first cloudy layer in sequence, split up radiance into clear and
800!    cloudy streams depending on cloud fraction
801      IF (ISTCLD(I_LEV)  ==  1) THEN
802        Z_CLDRADU(IPR) = Z_CLDFRAC(I_LEV) * Z_RADLU1(IPR)
803        Z_CLRRADU(IPR) = Z_RADLU1(IPR) - Z_CLDRADU(IPR)
804        Z_OLDCLD(IPR) = Z_CLDRADU(IPR)
805        Z_OLDCLR(IPR) = Z_CLRRADU(IPR)
806        Z_RAD(IPR) = 0.0_JPRB
807      ENDIF
808      Z_TTOT = 1.0_JPRB - Z_ATOT1(INDEX)
809      Z_TRNS = 1.0_JPRB - P_ABSS1(INDEX)
810      Z_CLDSRC = Z_BBUTOT1(INDEX) * Z_ATOT1(INDEX)
811
812!- Separate RT equations for clear and cloudy streams     
813      Z_CLDRADU(IPR) = Z_CLDRADU(IPR) * Z_TTOT + Z_CLDFRAC(I_LEV) * Z_CLDSRC
814      Z_CLRRADU(IPR) = Z_CLRRADU(IPR) * Z_TRNS +(1.0_JPRB - Z_CLDFRAC(I_LEV)) * Z_GASSRC
815!***
816
817!- total sky upward flux
818      Z_RADLU1(IPR) = Z_CLDRADU(IPR) + Z_CLRRADU(IPR)
819      Z_URAD1 = Z_URAD1 + Z_RADLU1(IPR)
820     
821!- clear-sky upward flux
822      Z_RADCLRU1(IPR) = Z_RADCLRU1(IPR) + (Z_BBU1(INDEX)-Z_RADCLRU1(IPR))&
823       & *P_ABSS1(INDEX) 
824      Z_URADCL1 = Z_URADCL1 + Z_RADCLRU1(IPR)
825
826!* Code to account for maximum/random overlap:
827!   Performs RT on the radiance most recently switched between clear and
828!   cloudy streams
829      Z_RADMOD = Z_RAD(IPR) * (Z_FACCLR1(I_LEV+1) * Z_TRNS +&
830       & Z_FACCLD1(I_LEV+1) *  Z_TTOT) - &
831       & Z_FACCMB1(I_LEV+1) * Z_GASSRC + &
832       & Z_FACCMB2(I_LEV+1) * Z_CLDSRC 
833       
834!   Computes what the clear and cloudy streams would have been had no
835!   radiance been switched       
836      Z_OLDCLD(IPR) = Z_CLDRADU(IPR) - Z_RADMOD
837      Z_OLDCLR(IPR) = Z_CLRRADU(IPR) + Z_RADMOD
838     
839!   Computes the radiance to be switched between clear and cloudy.     
840      Z_RAD(IPR) = -Z_RADMOD + Z_FACCLR2(I_LEV+1)*Z_OLDCLR(IPR) -&
841       & Z_FACCLD2(I_LEV+1)*Z_OLDCLD(IPR) 
842      Z_CLDRADU(IPR) = Z_CLDRADU(IPR) + Z_RAD(IPR)
843      Z_CLRRADU(IPR) = Z_CLRRADU(IPR) - Z_RAD(IPR)
844!***
845    ENDDO
846
847  ELSE
848
849!- *** Clear layer
850    IENT = JPGPT * (I_LEV-1)
851    DO IPR = 1, JPGPT
852      INDEX = IENT + IPR
853!- total-sky upward flux         
854      Z_RADLU1(IPR) = Z_RADLU1(IPR)+(Z_BBU1(INDEX)-Z_RADLU1(IPR))*P_ABSS1(INDEX)
855      Z_URAD1 = Z_URAD1 + Z_RADLU1(IPR)
856!- clear-sky upward flux
857!   Upward clear and total sky streams must be separate because surface
858!   reflectance is different for each.
859      Z_RADCLRU1(IPR) = Z_RADCLRU1(IPR)+(Z_BBU1(INDEX)-Z_RADCLRU1(IPR))*P_ABSS1(INDEX)
860      Z_URADCL1 = Z_URADCL1 + Z_RADCLRU1(IPR)
861    ENDDO
862
863  ENDIF
864
865  P_TOTUFLUC(I_LEV) = Z_URADCL1 * Z_WTNUM(1)
866  P_TOTUFLUX(I_LEV) = Z_URAD1   * Z_WTNUM(1)
867
868ENDDO
869
870!* Convert radiances to fluxes and heating rates for total and clear sky.
871! ** NB: moved to calling routine
872!      TOTUFLUC(0) = TOTUFLUC(0) * FLUXFAC
873!      TOTDFLUC(0) = TOTDFLUC(0) * FLUXFAC
874!      TOTUFLUX(0) = TOTUFLUX(0) * FLUXFAC
875!      TOTDFLUX(0) = TOTDFLUX(0) * FLUXFAC
876
877!      CLFNET(0) = (P_TOTUFLUC(0) - P_TOTDFLUC(0))
878!      FNET(0)   = (P_TOTUFLUX(0) - P_TOTDFLUX(0))
879!      DO LEV = 1, KLEV
880!        TOTUFLUC(LEV) = TOTUFLUC(LEV) * FLUXFAC
881!        TOTDFLUC(LEV) = TOTDFLUC(LEV) * FLUXFAC
882!         CLFNET(LEV) =(P_TOTUFLUC(LEV) - P_TOTDFLUC(LEV))
883
884!        TOTUFLUX(LEV) = TOTUFLUX(LEV) * FLUXFAC
885!        TOTDFLUX(LEV) = TOTDFLUX(LEV) * FLUXFAC
886!        FNET(LEV) = (P_TOTUFLUX(LEV) - P_TOTDFLUX(LEV))
887!        L = LEV - 1
888
889!- Calculate Heating Rates.
890!        CLHTR(L)=HEATFAC*(CLFNET(L)-CLFNET(LEV))/(PZ(L)-PZ(LEV))
891!        HTR(L)  =HEATFAC*(FNET(L)  -FNET(LEV))  /(PZ(L)-PZ(LEV))
892!      END DO
893!      CLHTR(KLEV) = 0.0
894!      HTR(KLEV)   = 0.0
895
896
897
898IF (LHOOK) CALL DR_HOOK('RRTM_RTRN1A_140GP',1,ZHOOK_HANDLE)
899END SUBROUTINE RRTM_RTRN1A_140GP
Note: See TracBrowser for help on using the repository browser.