source: LMDZ5/branches/Cold_pool_death/libf/phymar/rrtm_rtrn1a_140gp.F90 @ 3966

Last change on this file since 3966 was 2160, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
File size: 25.5 KB
Line 
1SUBROUTINE RRTM_RTRN1A_140GP (KLEV,ISTART,IEND,ICLDLYR,CLDFRAC,TAUCLD,ABSS1 &
2  &, OD,TAUSF1,CLFNET,CLHTR,FNET,HTR,TOTDFLUC,TOTDFLUX,TOTUFLUC,TOTUFLUX &
3  &, TAVEL,PZ,TZ,TBOUND,PFRAC,SEMISS,SEMISLW,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
24
25#include "tsmbkind.h"
26
27USE PARRRTM  , ONLY : JPBAND   ,JPGPT   ,JPLAY
28USE YOERRTAB , ONLY : BPADE
29USE YOERRTWN , ONLY : TOTPLNK  ,DELWAVE
30USE YOERRTFTR, ONLY : NGB
31
32IMPLICIT NONE
33
34!     DUMMY INTEGER SCALARS
35INTEGER_M :: KLEV
36INTEGER_M :: ISTART
37INTEGER_M :: IEND
38
39INTEGER_M :: ICLDLYR(JPLAY)      ! Cloud indicator
40REAL_B :: CLDFRAC(JPLAY)         ! Cloud fraction
41REAL_B :: TAUCLD(JPLAY,JPBAND)   ! Spectral optical thickness
42REAL_B :: ABSS1 (JPGPT*JPLAY)
43REAL_B :: OD    (JPGPT,JPLAY)
44REAL_B :: TAUSF1(JPGPT*JPLAY)
45REAL_B :: CLFNET  (0:JPLAY)
46REAL_B :: CLHTR   (0:JPLAY)
47REAL_B :: FNET    (0:JPLAY)
48REAL_B :: HTR     (0:JPLAY)
49REAL_B :: TOTDFLUC(0:JPLAY)
50REAL_B :: TOTDFLUX(0:JPLAY)
51REAL_B :: TOTUFLUC(0:JPLAY)
52REAL_B :: TOTUFLUX(0:JPLAY)
53
54!- from PROFILE             
55REAL_B :: TAVEL(JPLAY)
56REAL_B :: PZ(0:JPLAY)
57REAL_B :: TZ(0:JPLAY)
58REAL_B :: TBOUND
59
60!- from SP             
61REAL_B :: PFRAC(JPGPT,JPLAY)
62
63!- from SURFACE             
64REAL_B :: SEMISS(JPBAND)
65REAL_B :: SEMISLW
66INTEGER_M :: IREFLECT
67
68INTEGER_M :: INDLAY(JPLAY),INDLEV(0:JPLAY)
69
70REAL_B :: BBU1(JPGPT*JPLAY),BBUTOT1(JPGPT*JPLAY)
71REAL_B :: TLAYFRAC(JPLAY),TLEVFRAC(0:JPLAY)
72REAL_B :: BGLEV(JPGPT)
73!-- DS_000515
74REAL_B :: PLVL(0:JPLAY,JPBAND+1),PLAY(0:JPLAY,JPBAND+1),WTNUM(3)
75!-- DS_000515
76REAL_B :: ODCLD(JPBAND,JPLAY),EFCLFR1(JPBAND,JPLAY)
77REAL_B :: ODCLDNW(JPGPT,JPLAY)
78REAL_B :: SEMIS(JPGPT),RADUEMIT(JPGPT)
79
80REAL_B :: RADCLRU1(JPGPT) ,RADCLRD1(JPGPT)
81REAL_B :: RADLU1(JPGPT)   ,RADLD1(JPGPT)
82REAL_B :: ABSCLD1(JPBAND,JPLAY)
83!-- DS_000515
84REAL_B :: TRNCLD(JPLAY,JPBAND+1)
85!-- DS_000515
86REAL_B :: ABSCLDNW(JPGPT,JPLAY)
87REAL_B :: ATOT1(JPGPT*JPLAY)
88
89REAL_B :: SURFEMIS(JPBAND),PLNKEMIT(JPBAND)
90
91! dimension of arrays required for cloud overlap calculations
92
93REAL_B :: clrradu(jpgpt),cldradu(jpgpt),oldcld(jpgpt)
94REAL_B :: oldclr(jpgpt),rad(jpgpt),faccld1(jplay+1),faccld2(jplay+1)
95REAL_B :: facclr1(jplay+1),facclr2(jplay+1)
96REAL_B :: faccmb1(jplay+1),faccmb2(jplay+1)
97REAL_B :: faccld1d(0:jplay),faccld2d(0:jplay),facclr1d(0:jplay)
98REAL_B :: facclr2d(0:jplay),faccmb1d(0:jplay),faccmb2d(0:jplay)
99REAL_B :: clrradd(jpgpt),cldradd(jpgpt)
100INTEGER_M :: istcld(jplay+1),istcldd(0:jplay)
101!******
102
103REAL_B :: ZPLVL(JPGPT+1,JPLAY)  ,ZPLAY(JPGPT+1,JPLAY)
104REAL_B :: ZTRNCLD(JPGPT+1,JPLAY),ZTAUCLD(JPGPT+1,JPLAY)
105
106!     LOCAL INTEGER SCALARS
107INTEGER_M :: IBAND, ICLDDN, IENT, INDBOUND, INDEX, IPR, LAY, LEV, NBI
108
109!     LOCAL REAL SCALARS
110REAL_B :: BBD, BBDTOT, BGLAY, CLDSRC, DBDTLAY, DBDTLEV,&
111          &DELBGDN, DELBGUP, DRAD1, DRADCL1, FACTOT1, &
112          &FMAX, FMIN, GASSRC, ODSM, PLANKBND, RADCLD, &
113          &RADCLU, RADD, RADMOD, RADU, RAT1, RAT2, SUMPL, &
114          &SUMPLEM, TBNDFRAC, TRNS, TTOT, URAD1, URADCL1
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! Local
143!  ABSS(JPGPT*JPLAY)     !
144!  ABSCLD(JPLAY)         !
145!  ATOT(JPGPT*JPLAY)     !
146!  ODCLR(JPGPT,JPLAY)    !
147!  ODCLD(JPBAND,JPLAY)   !
148!  EFCLFR1(JPBAND,JPLAY) ! Effective cloud fraction
149!  RADLU(JPGPT)          ! Upward radiance
150!  URAD                  ! Spectrally summed upward radiance
151!  RADCLRU(JPGPT)        ! Clear-sky upward radiance
152!  CLRURAD               ! Spectrally summed clear-sky upward radiance
153!  RADLD(JPGPT)          ! Downward radiance
154!  DRAD                  ! Spectrally summed downward radiance
155!  RADCLRD(JPGPT)        ! Clear-sky downward radiance
156!  CLRDRAD               ! Spectrally summed clear-sky downward radiance
157!
158! Output
159!  TOTUFLUX(0:JPLAY)     ! Upward longwave flux
160!  TOTDFLUX(0:JPLAY)     ! Downward longwave flux
161!  TOTUFLUC(0:JPLAY)     ! Clear-sky upward longwave flux
162!  TOTDFLUC(0:JPLAY)     ! Clear-sky downward longwave flux
163!
164! Maximum/Random cloud overlap variables
165! for upward radiaitve transfer
166!  FACCLR2  fraction of clear radiance from previous layer that needs to
167!           be switched to cloudy stream
168!  FACCLR1  fraction of the radiance that had been switched in the previous
169!           layer from cloudy to clear that needs to be switched back to
170!           cloudy in the current layer
171!  FACCLD2  fraction of cloudy radiance from previous layer that needs to
172!           be switched to clear stream
173!           be switched to cloudy stream
174!  FACCLD1  fraction of the radiance that had been switched in the previous
175!           layer from clear to cloudy that needs to be switched back to
176!           clear in the current layer
177! for downward radiaitve transfer
178!  FACCLR2D fraction of clear radiance from previous layer that needs to
179!           be switched to cloudy stream
180!  FACCLR1D fraction of the radiance that had been switched in the previous
181!           layer from cloudy to clear that needs to be switched back to
182!           cloudy in the current layer
183!  FACCLD2D fraction of cloudy radiance from previous layer that needs to
184!           be switched to clear stream
185!           be switched to cloudy stream
186!  FACCLD1D fraction of the radiance that had been switched in the previous
187!           layer from clear to cloudy that needs to be switched back to
188!           clear in the current layer
189!
190!--------------------------------------------------------------------------
191
192WTNUM(1)=_HALF_
193WTNUM(2)=_ZERO_
194WTNUM(3)=_ZERO_
195
196!-start JJM_000511
197IF (TBOUND < 339._JPRB .AND. TBOUND >= 160._JPRB ) THEN
198  INDBOUND = TBOUND - 159._JPRB
199  TBNDFRAC = TBOUND - INT(TBOUND)
200ELSE IF (TBOUND >= 339._JPRB ) THEN
201  INDBOUND = 180
202  TBNDFRAC = TBOUND - 339._JPRB
203ELSE IF (TBOUND < 160._JPRB ) THEN
204  INDBOUND = 1
205  TBNDFRAC = TBOUND - 160._JPRB
206ENDIF 
207!-end JJM_000511
208 
209DO LAY = 0, KLEV
210  TOTUFLUC(LAY) = _ZERO_
211  TOTDFLUC(LAY) = _ZERO_
212  TOTUFLUX(LAY) = _ZERO_
213  TOTDFLUX(LAY) = _ZERO_
214!-start JJM_000511
215  IF (TZ(LAY) < 339._JPRB .AND. TZ(LAY) >= 160._JPRB ) THEN
216    INDLEV(LAY) = TZ(LAY) - 159._JPRB
217    TLEVFRAC(LAY) = TZ(LAY) - INT(TZ(LAY))
218  ELSE IF (TZ(LAY) >= 339._JPRB ) THEN
219    INDLEV(LAY) = 180
220    TLEVFRAC(LAY) = TZ(LAY) - 339._JPRB
221  ELSE IF (TZ(LAY) < 160._JPRB ) THEN
222    INDLEV(LAY) = 1
223    TLEVFRAC(LAY) = TZ(LAY) - 160._JPRB
224  ENDIF   
225!-end JJM_000511
226ENDDO
227
228!_start_jjm 991209
229DO LEV=0,KLEV
230  FACCLD1(LEV+1) = _ZERO_
231  FACCLD2(LEV+1) = _ZERO_
232  FACCLR1(LEV+1) = _ZERO_
233  FACCLR2(LEV+1) = _ZERO_
234  FACCMB1(LEV+1) = _ZERO_
235  FACCMB2(LEV+1) = _ZERO_
236  FACCLD1D(LEV) = _ZERO_
237  FACCLD2D(LEV) = _ZERO_
238  FACCLR1D(LEV) = _ZERO_
239  FACCLR2D(LEV) = _ZERO_
240  FACCMB1D(LEV) = _ZERO_
241  FACCMB2D(LEV) = _ZERO_
242END DO 
243RAT1 = _ZERO_
244RAT2 = _ZERO_
245!_end_jjm 991209
246
247
248
249SUMPL   = _ZERO_
250SUMPLEM = _ZERO_
251
252ISTCLD(1) = 1
253ISTCLDD(KLEV) = 1
254
255DO LEV = 1, KLEV
256!-- DS_000515
257!-start JJM_000511
258  IF (TAVEL(LEV) < 339._JPRB .AND. TAVEL(LEV) >= 160._JPRB ) THEN
259    INDLAY(LEV) = TAVEL(LEV) - 159._JPRB
260    TLAYFRAC(LEV) = TAVEL(LEV) - INT(TAVEL(LEV))
261  ELSE IF (TAVEL(LEV) >= 339._JPRB ) THEN
262    INDLAY(LEV) = 180
263    TLAYFRAC(LEV) = TAVEL(LEV) - 339._JPRB
264  ELSE IF (TAVEL(LEV) < 160._JPRB ) THEN
265    INDLAY(LEV) = 1
266    TLAYFRAC(LEV) = TAVEL(LEV) - 160._JPRB
267  ENDIF 
268!-end JJM_000511
269END DO
270!-- DS_000515
271
272!-- DS_000515
273!OCL SCALAR
274
275DO LEV = 1, KLEV
276  IF (ICLDLYR(LEV) == 1) THEN
277
278!mji   
279    ISTCLD(LEV+1) = 0
280    IF (LEV  ==  KLEV) THEN
281      FACCLD1(LEV+1) = _ZERO_
282      FACCLD2(LEV+1) = _ZERO_
283      FACCLR1(LEV+1) = _ZERO_
284      FACCLR2(LEV+1) = _ZERO_
285!-- DS_000515     
286!      FACCMB1(LEV+1) = _ZERO_
287!      FACCMB2(LEV+1) = _ZERO_
288!mji      ISTCLD(LEV+1) = _ZERO_
289    ELSEIF (CLDFRAC(LEV+1)  >=  CLDFRAC(LEV)) THEN
290      FACCLD1(LEV+1) = _ZERO_
291      FACCLD2(LEV+1) = _ZERO_
292      IF (ISTCLD(LEV)  ==  1) THEN
293!mji        ISTCLD(LEV+1) = 0
294        FACCLR1(LEV+1) = _ZERO_
295!mji       
296        FACCLR2(LEV+1) = _ZERO_
297        IF (CLDFRAC(LEV) < _ONE_) THEN
298        FACCLR2(LEV+1) = (CLDFRAC(LEV+1)-CLDFRAC(LEV))/&
299         &(_ONE_-CLDFRAC(LEV))
300        END IF   
301      ELSE
302        FMAX = MAX(CLDFRAC(LEV),CLDFRAC(LEV-1))
303!mji
304        IF (CLDFRAC(LEV+1)  >  FMAX) THEN
305          FACCLR1(LEV+1) = RAT2
306          FACCLR2(LEV+1) = (CLDFRAC(LEV+1)-FMAX)/(_ONE_-FMAX)
307!mji         
308        ELSE IF (CLDFRAC(LEV+1) < FMAX) THEN
309          FACCLR1(LEV+1) = (CLDFRAC(LEV+1)-CLDFRAC(LEV))/&
310           &(CLDFRAC(LEV-1)-CLDFRAC(LEV))
311          FACCLR2(LEV+1) = _ZERO_
312!mji
313        ELSE
314          FACCLR1(LEV+1) = RAT2 
315          FACCLR2(LEV+1) = _ZERO_
316        ENDIF
317      ENDIF
318      IF (FACCLR1(LEV+1) > _ZERO_ .OR. FACCLR2(LEV+1) > _ZERO_) THEN
319        RAT1 = _ONE_
320        RAT2 = _ZERO_
321      ENDIF
322    ELSE
323      FACCLR1(LEV+1) = _ZERO_
324      FACCLR2(LEV+1) = _ZERO_
325      IF (ISTCLD(LEV)  ==  1) THEN
326!mji        ISTCLD(LEV+1) = 0
327        FACCLD1(LEV+1) = _ZERO_
328        FACCLD2(LEV+1) = (CLDFRAC(LEV)-CLDFRAC(LEV+1))/CLDFRAC(LEV)
329      ELSE
330        FMIN = MIN(CLDFRAC(LEV),CLDFRAC(LEV-1))
331        IF (CLDFRAC(LEV+1)  <=  FMIN) THEN
332          FACCLD1(LEV+1) = RAT1
333          FACCLD2(LEV+1) = (FMIN-CLDFRAC(LEV+1))/FMIN
334        ELSE
335          FACCLD1(LEV+1) = (CLDFRAC(LEV)-CLDFRAC(LEV+1))/&
336           &(CLDFRAC(LEV)-FMIN)
337          FACCLD2(LEV+1) = _ZERO_
338        ENDIF
339      ENDIF
340      IF (FACCLD1(LEV+1) > _ZERO_ .OR. FACCLD2(LEV+1) > _ZERO_) THEN
341        RAT1 = _ZERO_
342        RAT2 = _ONE_
343      ENDIF
344    ENDIF
345!fcc
346    IF (LEV == 1) THEN
347      FACCMB1(LEV+1) = 0.
348      FACCMB2(LEV+1) = FACCLD1(LEV+1) * FACCLR2(LEV)
349    ELSE
350      FACCMB1(LEV+1) = FACCLR1(LEV+1) * FACCLD2(LEV) *CLDFRAC(LEV-1)
351      FACCMB2(LEV+1) = FACCLD1(LEV+1) * FACCLR2(LEV) *&
352       &(_ONE_ - CLDFRAC(LEV-1))
353    ENDIF
354!end fcc
355  ELSE
356!-- DS_000515
357    ISTCLD(LEV+1) = 1
358  ENDIF
359ENDDO
360
361!_start_jjm 991209
362RAT1 = _ZERO_
363RAT2 = _ZERO_
364!_end_jjm 991209
365
366!-- DS_000515
367!OCL SCALAR
368
369DO LEV = KLEV, 1, -1
370  IF (ICLDLYR(LEV) == 1) THEN
371!mji
372    ISTCLDD(LEV-1) = 0 
373    IF (LEV  ==  1) THEN
374      FACCLD1D(LEV-1) = _ZERO_
375      FACCLD2D(LEV-1) = _ZERO_
376      FACCLR1D(LEV-1) = _ZERO_
377      FACCLR2D(LEV-1) = _ZERO_
378      FACCMB1D(LEV-1) = _ZERO_
379      FACCMB2D(LEV-1) = _ZERO_
380!mji      ISTCLDD(LEV-1) = _ZERO_
381    ELSEIF (CLDFRAC(LEV-1)  >=  CLDFRAC(LEV)) THEN
382      FACCLD1D(LEV-1) = _ZERO_
383      FACCLD2D(LEV-1) = _ZERO_
384      IF (ISTCLDD(LEV)  ==  1) THEN
385!mji        ISTCLDD(LEV-1) = 0
386        FACCLR1D(LEV-1) = _ZERO_
387        FACCLR2D(LEV-1) = _ZERO_
388        IF (CLDFRAC(LEV) < _ONE_) THEN
389          FACCLR2D(LEV-1) = (CLDFRAC(LEV-1)-CLDFRAC(LEV))/&
390           &(_ONE_-CLDFRAC(LEV))
391        END IF
392      ELSE
393        FMAX = MAX(CLDFRAC(LEV),CLDFRAC(LEV+1))
394!mji
395        IF (CLDFRAC(LEV-1)  >  FMAX) THEN
396          FACCLR1D(LEV-1) = RAT2
397          FACCLR2D(LEV-1) = (CLDFRAC(LEV-1)-FMAX)/(_ONE_-FMAX)
398!mji
399        ELSE IF (CLDFRAC(LEV-1) < FMAX) THEN
400          FACCLR1D(LEV-1) = (CLDFRAC(LEV-1)-CLDFRAC(LEV))/&
401           &(CLDFRAC(LEV+1)-CLDFRAC(LEV))
402          FACCLR2D(LEV-1) = _ZERO_
403!mji
404        ELSE         
405          FACCLR1D(LEV-1) = RAT2
406          FACCLR2D(LEV-1) = _ZERO_
407        ENDIF
408      ENDIF
409      IF (FACCLR1D(LEV-1) > _ZERO_ .OR. FACCLR2D(LEV-1) > _ZERO_)THEN
410        RAT1 = _ONE_
411        RAT2 = _ZERO_
412      ENDIF
413    ELSE
414      FACCLR1D(LEV-1) = _ZERO_
415      FACCLR2D(LEV-1) = _ZERO_
416      IF (ISTCLDD(LEV)  ==  1) THEN
417!mji        ISTCLDD(LEV-1) = 0
418        FACCLD1D(LEV-1) = _ZERO_
419        FACCLD2D(LEV-1) = (CLDFRAC(LEV)-CLDFRAC(LEV-1))/CLDFRAC(LEV)
420      ELSE
421        FMIN = MIN(CLDFRAC(LEV),CLDFRAC(LEV+1))
422        IF (CLDFRAC(LEV-1)  <=  FMIN) THEN
423          FACCLD1D(LEV-1) = RAT1
424          FACCLD2D(LEV-1) = (FMIN-CLDFRAC(LEV-1))/FMIN
425        ELSE
426          FACCLD1D(LEV-1) = (CLDFRAC(LEV)-CLDFRAC(LEV-1))/&
427           &(CLDFRAC(LEV)-FMIN)
428          FACCLD2D(LEV-1) = _ZERO_
429        ENDIF
430      ENDIF
431      IF (FACCLD1D(LEV-1) > _ZERO_ .OR. FACCLD2D(LEV-1) > _ZERO_)THEN
432        RAT1 = _ZERO_
433        RAT2 = _ONE_
434      ENDIF
435    ENDIF
436    FACCMB1D(LEV-1) = FACCLR1D(LEV-1) * FACCLD2D(LEV) *CLDFRAC(LEV+1)
437    FACCMB2D(LEV-1) = FACCLD1D(LEV-1) * FACCLR2D(LEV) *&
438     &(_ONE_ - CLDFRAC(LEV+1))
439  ELSE
440    ISTCLDD(LEV-1) = 1
441  ENDIF
442ENDDO
443
444
445!- Loop over frequency bands.
446
447DO IBAND = ISTART, IEND
448  DBDTLEV = TOTPLNK(INDBOUND+1,IBAND)-TOTPLNK(INDBOUND,IBAND)
449  PLANKBND = DELWAVE(IBAND) * (TOTPLNK(INDBOUND,IBAND) + TBNDFRAC * DBDTLEV)
450  DBDTLEV = TOTPLNK(INDLEV(0)+1,IBAND) -TOTPLNK(INDLEV(0),IBAND)
451!-- DS_000515
452  PLVL(0,IBAND) = DELWAVE(IBAND)&
453   &* (TOTPLNK(INDLEV(0),IBAND) + TLEVFRAC(0)*DBDTLEV)
454
455  SURFEMIS(IBAND) = SEMISS(IBAND)
456  PLNKEMIT(IBAND) = SURFEMIS(IBAND) * PLANKBND
457  SUMPLEM  = SUMPLEM + PLNKEMIT(IBAND)
458  SUMPL    = SUMPL   + PLANKBND
459!--DS
460ENDDO
461!---
462
463!-- DS_000515
464DO IBAND = ISTART, IEND
465  DO LEV = 1, KLEV
466!----             
467!- Calculate the integrated Planck functions for at the
468!  level and layer temperatures.
469!  Compute cloud transmittance for cloudy layers.
470    DBDTLEV = TOTPLNK(INDLEV(LEV)+1,IBAND) - TOTPLNK(INDLEV(LEV),IBAND)
471    DBDTLAY = TOTPLNK(INDLAY(LEV)+1,IBAND) - TOTPLNK(INDLAY(LEV),IBAND)
472!-- DS_000515
473    PLAY(LEV,IBAND) = DELWAVE(IBAND)&
474     &*(TOTPLNK(INDLAY(LEV),IBAND)+TLAYFRAC(LEV)*DBDTLAY)
475    PLVL(LEV,IBAND) = DELWAVE(IBAND)&
476     &*(TOTPLNK(INDLEV(LEV),IBAND)+TLEVFRAC(LEV)*DBDTLEV)
477    IF (ICLDLYR(LEV) > 0) THEN
478      TRNCLD(LEV,IBAND) = EXP(-TAUCLD(LEV,IBAND))
479    ENDIF
480!-- DS_000515
481  ENDDO
482
483ENDDO
484
485SEMISLW = SUMPLEM / SUMPL
486
487!--DS
488DO IPR = 1, JPGPT
489  NBI = NGB(IPR)
490  DO LEV =  1 , KLEV
491!-- DS_000515
492    ZPLAY(IPR,LEV) = PLAY(LEV,NBI)
493    ZPLVL(IPR,LEV) = PLVL(LEV-1,NBI)
494    ZTAUCLD(IPR,LEV) = TAUCLD(LEV,NBI)
495    ZTRNCLD(IPR,LEV) = TRNCLD(LEV,NBI)
496!-- DS_000515
497  ENDDO
498ENDDO
499!----     
500
501!- For cloudy layers, set cloud parameters for radiative transfer.
502DO LEV = 1, KLEV
503  IF (ICLDLYR(LEV) > 0) THEN
504    DO IPR = 1, JPGPT
505!--DS         
506!            NBI = NGB(IPR)
507      ODCLDNW(IPR,LEV) = ZTAUCLD(IPR,LEV)
508      ABSCLDNW(IPR,LEV) = _ONE_ - ZTRNCLD(IPR,LEV)
509!----           
510!            EFCLFRNW(IPR,LEV) = ABSCLDNW(IPR,LEV) * CLDFRAC(LEV)
511    ENDDO
512  ENDIF
513ENDDO
514
515!- Initialize for radiative transfer.
516DO IPR = 1, JPGPT
517  RADCLRD1(IPR) = _ZERO_
518  RADLD1(IPR)   = _ZERO_
519  NBI = NGB(IPR)
520  SEMIS(IPR) = SURFEMIS(NBI)
521  RADUEMIT(IPR) = PFRAC(IPR,1) * PLNKEMIT(NBI)
522!-- DS_000515
523  BGLEV(IPR) = PFRAC(IPR,KLEV) * PLVL(KLEV,NBI)
524ENDDO
525
526!- Downward radiative transfer.
527!  *** DRAD1 holds summed radiance for total sky stream
528!  *** DRADCL1 holds summed radiance for clear sky stream
529
530ICLDDN = 0
531DO LEV = KLEV, 1, -1
532  DRAD1   = _ZERO_
533  DRADCL1 = _ZERO_
534
535  IF (ICLDLYR(LEV) == 1) THEN
536
537!  *** Cloudy layer
538    ICLDDN = 1
539    IENT = JPGPT * (LEV-1)
540    DO IPR = 1, JPGPT
541      INDEX = IENT + IPR
542!--DS           
543!            NBI = NGB(IPR)
544      BGLAY = PFRAC(IPR,LEV) * ZPLAY(IPR,LEV)
545!----           
546      DELBGUP     = BGLEV(IPR) - BGLAY
547      BBU1(INDEX) = BGLAY + TAUSF1(INDEX) * DELBGUP
548!--DS           
549      BGLEV(IPR) = PFRAC(IPR,LEV) * ZPLVL(IPR,LEV)
550!----           
551      DELBGDN = BGLEV(IPR) - BGLAY
552      BBD = BGLAY + TAUSF1(INDEX) * DELBGDN
553!- total-sky downward flux         
554      ODSM = OD(IPR,LEV) + ODCLDNW(IPR,LEV)
555      FACTOT1 = ODSM / (BPADE + ODSM)
556      BBUTOT1(INDEX) = BGLAY + FACTOT1 * DELBGUP
557      ATOT1(INDEX) = ABSS1(INDEX) + ABSCLDNW(IPR,LEV)&
558       &- ABSS1(INDEX) * ABSCLDNW(IPR,LEV)
559      BBDTOT = BGLAY + FACTOT1 * DELBGDN
560      GASSRC = BBD * ABSS1(INDEX)
561!***
562      IF (ISTCLDD(LEV)  ==  1) THEN
563        CLDRADD(IPR) = CLDFRAC(LEV) * RADLD1(IPR)
564        CLRRADD(IPR) = RADLD1(IPR) - CLDRADD(IPR)
565        OLDCLD(IPR) = CLDRADD(IPR)
566        OLDCLR(IPR) = CLRRADD(IPR)
567        RAD(IPR) = _ZERO_
568      ENDIF
569      TTOT = _ONE_ - ATOT1(INDEX)
570      CLDSRC = BBDTOT * ATOT1(INDEX)
571     
572! Separate RT equations for clear and cloudy streams     
573      CLDRADD(IPR) = CLDRADD(IPR) * TTOT + CLDFRAC(LEV) * CLDSRC
574      CLRRADD(IPR) = CLRRADD(IPR) * (_ONE_-ABSS1(INDEX)) +&
575       &(_ONE_ - CLDFRAC(LEV)) * GASSRC
576
577!  Total sky downward radiance
578      RADLD1(IPR) = CLDRADD(IPR) + CLRRADD(IPR)
579      DRAD1 = DRAD1 + RADLD1(IPR)
580     
581!  Clear-sky downward radiance         
582      RADCLRD1(IPR) = RADCLRD1(IPR)+(BBD-RADCLRD1(IPR))*ABSS1(INDEX)
583      DRADCL1 = DRADCL1 + RADCLRD1(IPR)
584
585!* Code to account for maximum/random overlap:
586!   Performs RT on the radiance most recently switched between clear and
587!   cloudy streams
588      RADMOD = RAD(IPR) * (FACCLR1D(LEV-1) * (_ONE_-ABSS1(INDEX)) +&
589       &FACCLD1D(LEV-1) *  TTOT) - &
590       &FACCMB1D(LEV-1) * GASSRC + &
591       &FACCMB2D(LEV-1) * CLDSRC
592       
593!   Computes what the clear and cloudy streams would have been had no
594!   radiance been switched       
595      OLDCLD(IPR) = CLDRADD(IPR) - RADMOD
596      OLDCLR(IPR) = CLRRADD(IPR) + RADMOD
597     
598!   Computes the radiance to be switched between clear and cloudy.     
599      RAD(IPR) = -RADMOD + FACCLR2D(LEV-1)*OLDCLR(IPR) -&
600       &FACCLD2D(LEV-1)*OLDCLD(IPR)
601      CLDRADD(IPR) = CLDRADD(IPR) + RAD(IPR)
602      CLRRADD(IPR) = CLRRADD(IPR) - RAD(IPR)
603!***
604
605    ENDDO
606
607  ELSE
608
609!  *** Clear layer
610!  *** DRAD1 holds summed radiance for total sky stream
611!  *** DRADCL1 holds summed radiance for clear sky stream
612
613    IENT = JPGPT * (LEV-1)
614    IF (ICLDDN == 1) THEN
615      DO IPR = 1, JPGPT
616        INDEX = IENT + IPR
617!--DS         
618!           NBI = NGB(IPR)
619        BGLAY = PFRAC(IPR,LEV) * ZPLAY(IPR,LEV)
620!----           
621        DELBGUP     = BGLEV(IPR) - BGLAY
622        BBU1(INDEX) = BGLAY + TAUSF1(INDEX) * DELBGUP
623!--DS           
624        BGLEV(IPR) = PFRAC(IPR,LEV) * ZPLVL(IPR,LEV)
625!----                     
626        DELBGDN = BGLEV(IPR) - BGLAY
627        BBD = BGLAY + TAUSF1(INDEX) * DELBGDN
628       
629!- total-sky downward radiance
630        RADLD1(IPR) = RADLD1(IPR)+(BBD-RADLD1(IPR))*ABSS1(INDEX)
631        DRAD1 = DRAD1 + RADLD1(IPR)
632       
633!- clear-sky downward radiance
634!-  Set clear sky stream to total sky stream as long as layers
635!-  remain clear.  Streams diverge when a cloud is reached.
636        RADCLRD1(IPR) = RADCLRD1(IPR)+(BBD-RADCLRD1(IPR))*ABSS1(INDEX)
637        DRADCL1 = DRADCL1 + RADCLRD1(IPR)
638      ENDDO
639           
640    ELSE
641       
642      DO IPR = 1, JPGPT
643        INDEX = IENT + IPR
644!--DS         
645!           NBI = NGB(IPR)
646        BGLAY = PFRAC(IPR,LEV) * ZPLAY(IPR,LEV)
647!----           
648        DELBGUP     = BGLEV(IPR) - BGLAY
649        BBU1(INDEX) = BGLAY + TAUSF1(INDEX) * DELBGUP
650!--DS           
651        BGLEV(IPR) = PFRAC(IPR,LEV) * ZPLVL(IPR,LEV)
652!----                     
653        DELBGDN = BGLEV(IPR) - BGLAY
654        BBD = BGLAY + TAUSF1(INDEX) * DELBGDN
655!- total-sky downward flux         
656        RADLD1(IPR) = RADLD1(IPR)+(BBD-RADLD1(IPR))*ABSS1(INDEX)
657        DRAD1 = DRAD1 + RADLD1(IPR)
658!- clear-sky downward flux         
659!-  Set clear sky stream to total sky stream as long as layers
660!-  remain clear.  Streams diverge when a cloud is reached.
661        RADCLRD1(IPR) = RADLD1(IPR)
662      ENDDO
663      DRADCL1 = DRAD1
664    ENDIF
665   
666  ENDIF
667
668  TOTDFLUC(LEV-1) = DRADCL1 * WTNUM(1)
669  TOTDFLUX(LEV-1) = DRAD1   * WTNUM(1)
670
671ENDDO
672
673
674
675
676
677! Spectral reflectivity and reflectance
678! Includes the contribution of spectrally varying longwave emissivity
679! and reflection from the surface to the upward radiative transfer.
680! Note: Spectral and Lambertian reflections are identical for the one
681! angle flux integration used here.
682
683URAD1   = _ZERO_
684URADCL1 = _ZERO_
685
686!start JJM_000511
687!IF (IREFLECT  ==  0) THEN
688!- Lambertian reflection.
689  DO IPR = 1, JPGPT
690! Clear-sky radiance
691!    RADCLD = _TWO_ * (RADCLRD1(IPR) * WTNUM(1) )
692    RADCLD = RADCLRD1(IPR)
693    RADCLRU1(IPR) = RADUEMIT(IPR) + (_ONE_ - SEMIS(IPR)) * RADCLD
694    URADCL1 = URADCL1 + RADCLRU1(IPR)
695
696! Total sky radiance
697!    RADD = _TWO_ * (RADLD1(IPR) * WTNUM(1) )
698    RADD = RADLD1(IPR)
699    RADLU1(IPR) = RADUEMIT(IPR) + (_ONE_ - SEMIS(IPR)) * RADD
700    URAD1 = URAD1 + RADLU1(IPR)
701  ENDDO
702  TOTUFLUC(0) = URADCL1 * _HALF_
703  TOTUFLUX(0) = URAD1 * _HALF_
704!ELSE
705!!- Specular reflection.
706!  DO IPR = 1, JPGPT
707!    RADCLU = RADUEMIT(IPR)
708!    RADCLRU1(IPR) = RADCLU + (_ONE_ - SEMIS(IPR)) * RADCLRD1(IPR)
709!    URADCL1 = URADCL1 + RADCLRU1(IPR)
710!
711!    RADU = RADUEMIT(IPR)
712!    RADLU1(IPR) = RADU + (_ONE_ - SEMIS(IPR)) * RADLD1(IPR)
713!    URAD1 = URAD1 + RADLU1(IPR)
714!  ENDDO
715!  TOTUFLUC(0) = URADCL1 * WTNUM(1)
716!  TOTUFLUX(0) = URAD1   * WTNUM(1)
717!ENDIF
718
719!- Upward radiative transfer.
720!- *** URAD1 holds the summed radiance for total sky stream
721!- *** URADCL1 holds the summed radiance for clear sky stream
722DO LEV = 1, KLEV
723  URAD1   = _ZERO_
724  URADCL1 = _ZERO_
725
726! Check flag for cloud in current layer
727  IF (ICLDLYR(LEV) == 1) THEN
728
729!- *** Cloudy layer
730    IENT = JPGPT * (LEV-1)
731    DO IPR = 1, JPGPT
732      INDEX = IENT + IPR
733!- total-sky upward flux         
734      GASSRC = BBU1(INDEX) * ABSS1(INDEX)
735!
736!- If first cloudy layer in sequence, split up radiance into clear and
737!    cloudy streams depending on cloud fraction
738      IF (ISTCLD(LEV)  ==  1) THEN
739        CLDRADU(IPR) = CLDFRAC(LEV) * RADLU1(IPR)
740        CLRRADU(IPR) = RADLU1(IPR) - CLDRADU(IPR)
741        OLDCLD(IPR) = CLDRADU(IPR)
742        OLDCLR(IPR) = CLRRADU(IPR)
743        RAD(IPR) = _ZERO_
744      ENDIF
745      TTOT = _ONE_ - ATOT1(INDEX)
746      TRNS = _ONE_ - ABSS1(INDEX)
747      CLDSRC = BBUTOT1(INDEX) * ATOT1(INDEX)
748!     
749!- Separate RT equations for clear and cloudy streams     
750      CLDRADU(IPR) = CLDRADU(IPR) * TTOT + CLDFRAC(LEV) * CLDSRC
751      CLRRADU(IPR) = CLRRADU(IPR) * TRNS +(_ONE_ - CLDFRAC(LEV)) * GASSRC
752!***
753
754!- total sky upward flux
755      RADLU1(IPR) = CLDRADU(IPR) + CLRRADU(IPR)
756      URAD1 = URAD1 + RADLU1(IPR)
757     
758!- clear-sky upward flux
759      RADCLRU1(IPR) = RADCLRU1(IPR) + (BBU1(INDEX)-RADCLRU1(IPR))&
760       &*ABSS1(INDEX)
761      URADCL1 = URADCL1 + RADCLRU1(IPR)
762
763!* Code to account for maximum/random overlap:
764!   Performs RT on the radiance most recently switched between clear and
765!   cloudy streams
766      RADMOD = RAD(IPR) * (FACCLR1(LEV+1) * TRNS +&
767       &FACCLD1(LEV+1) *  TTOT) - &
768       &FACCMB1(LEV+1) * GASSRC + &
769       &FACCMB2(LEV+1) * CLDSRC
770       
771!   Computes what the clear and cloudy streams would have been had no
772!   radiance been switched       
773      OLDCLD(IPR) = CLDRADU(IPR) - RADMOD
774      OLDCLR(IPR) = CLRRADU(IPR) + RADMOD
775     
776!   Computes the radiance to be switched between clear and cloudy.     
777      RAD(IPR) = -RADMOD + FACCLR2(LEV+1)*OLDCLR(IPR) -&
778       &FACCLD2(LEV+1)*OLDCLD(IPR)
779      CLDRADU(IPR) = CLDRADU(IPR) + RAD(IPR)
780      CLRRADU(IPR) = CLRRADU(IPR) - RAD(IPR)
781!***
782    ENDDO
783
784  ELSE
785
786!- *** Clear layer
787    IENT = JPGPT * (LEV-1)
788    DO IPR = 1, JPGPT
789      INDEX = IENT + IPR
790!- total-sky upward flux         
791      RADLU1(IPR) = RADLU1(IPR)+(BBU1(INDEX)-RADLU1(IPR))*ABSS1(INDEX)
792      URAD1 = URAD1 + RADLU1(IPR)
793!- clear-sky upward flux
794!   Upward clear and total sky streams must be separate because surface
795!   reflectance is different for each.
796      RADCLRU1(IPR) = RADCLRU1(IPR)+(BBU1(INDEX)-RADCLRU1(IPR))*ABSS1(INDEX)
797      URADCL1 = URADCL1 + RADCLRU1(IPR)
798    ENDDO
799
800  ENDIF
801
802  TOTUFLUC(LEV) = URADCL1 * WTNUM(1)
803  TOTUFLUX(LEV) = URAD1   * WTNUM(1)
804
805ENDDO
806
807
808!* Convert radiances to fluxes and heating rates for total and clear sky.
809! ** NB: moved to calling routine
810!      TOTUFLUC(0) = TOTUFLUC(0) * FLUXFAC
811!      TOTDFLUC(0) = TOTDFLUC(0) * FLUXFAC
812!      TOTUFLUX(0) = TOTUFLUX(0) * FLUXFAC
813!      TOTDFLUX(0) = TOTDFLUX(0) * FLUXFAC
814
815!      CLFNET(0) = TOTUFLUC(0) - TOTDFLUC(0)
816!      FNET(0)   = TOTUFLUX(0) - TOTDFLUX(0)
817!      DO LEV = 1, KLEV
818!        TOTUFLUC(LEV) = TOTUFLUC(LEV) * FLUXFAC
819!        TOTDFLUC(LEV) = TOTDFLUC(LEV) * FLUXFAC
820!        CLFNET(LEV) = TOTUFLUC(LEV) - TOTDFLUC(LEV)
821
822!        TOTUFLUX(LEV) = TOTUFLUX(LEV) * FLUXFAC
823!        TOTDFLUX(LEV) = TOTDFLUX(LEV) * FLUXFAC
824!        FNET(LEV) = TOTUFLUX(LEV) - TOTDFLUX(LEV)
825!        L = LEV - 1
826
827!- Calculate Heating Rates.
828!        CLHTR(L)=HEATFAC*(CLFNET(L)-CLFNET(LEV))/(PZ(L)-PZ(LEV))
829!        HTR(L)  =HEATFAC*(FNET(L)  -FNET(LEV))  /(PZ(L)-PZ(LEV))
830!      END DO
831!      CLHTR(KLEV) = 0.0
832!      HTR(KLEV)   = 0.0
833
834RETURN
835END SUBROUTINE RRTM_RTRN1A_140GP
Note: See TracBrowser for help on using the repository browser.