File: C:\NOAA\NEMS_11731\src\atmos\phys\module_RA_GFDL.F90

1     !-----------------------------------------------------------------------
2     !
3           MODULE MODULE_RA_GFDL
4     !
5     !-----------------------------------------------------------------------
6     !
7     !***  THE RADIATION DRIVERS AND PACKAGES
8     !
9     !---------------------
10     !--- Modifications ---
11     !---------------------
12     !--- 24 Feb 2010 - Ferrier 
13     ! 1) Removed EQUIVALENCE for 2D saved arrays (TABLE1, etc.), made them 
14     !    local in subroutine TABLE, and used results to define 1D saved arrays 
15     !    (T1, etc.) for use in other subroutines (produced upper-level warm bias).
16     ! 2) Removed "mep-test" block in RADFS (comment not correct).
17     ! 3) 
18     !-----------------------------------------------------------------------
19     !
20           USE MODULE_INCLUDE
21           USE MODULE_CONSTANTS,ONLY : CAPPA,CP,EP_2,EPSQ,G,PI
22     !
23           USE MODULE_MP_ETANEW,ONLY : RHgrd,T_ICE,FPVS,QAUT0,XMImax,XMIexp  &
24                                      ,MDImin,MDImax,MASSI,FLARGE1,FLARGE2   &
25                                      ,NLImin,NLImax,GPVS
26     !
27     !-----------------------------------------------------------------------
28     !
29           IMPLICIT NONE
30     !
31     !-----------------------------------------------------------------------
32     !
33           PRIVATE
34     !
35           PUBLIC :: GFDL_INIT,RDTEMP,TIME_MEASURE,GFDL,CAL_MON_DAY,ZENITH
36     !
37     !-----------------------------------------------------------------------
38     !
39           REAL :: DPD=360./365.                                             &
40                  ,RLAG=14.8125
41     !
42     !-----------------------------------------------------------------------
43     !-----------------------------------------------------------------------
44     !***  FOR GFDL RADIATION 
45     !-----------------------------------------------------------------------
46     !-----------------------------------------------------------------------
47     !
48           INTEGER,PARAMETER :: NL=81
49           INTEGER,PARAMETER :: NBLY=15
50           REAL,PARAMETER :: DEGRAD=3.1415926/180.
51           REAL,PARAMETER :: RTHRESH=1.E-15,RTD=1./DEGRAD
52     
53           INTEGER, SAVE, DIMENSION(3)     :: LTOP
54           REAL   , SAVE, DIMENSION(37,NL) :: XDUO3N,XDO3N2,XDO3N3,XDO3N4
55           REAL   , SAVE, DIMENSION(NL)    :: PRGFDL
56           REAL   , SAVE                   :: AB15WD,SKO2D,SKC1R,SKO3R
57     
58           REAL   , SAVE :: SOURCE(28,NBLY), DSRCE(28,NBLY)
59     
60           REAL   ,SAVE, DIMENSION(5040):: T1,T2,T4,EM1V,EM1VW,EM3V
61           REAL   ,SAVE                 :: R1,RSIN1,RCOS1,RCOS2
62     ! Created by CO2 initialization
63           REAL,   SAVE, ALLOCATABLE, DIMENSION(:,:) :: CO251,CDT51,CDT58,C2D51,&
64                                                C2D58,CO258
65           REAL,   SAVE, ALLOCATABLE, DIMENSION(:)   :: STEMP,GTEMP,CO231,CO238, &
66                                                C2D31,C2D38,CDT31,CDT38, &
67                                                CO271,CO278,C2D71,C2D78, &
68                                                CDT71,CDT78
69           REAL,   SAVE, ALLOCATABLE, DIMENSION(:)   :: CO2M51,CO2M58,CDTM51,CDTM58, &
70                                                C2DM51,C2DM58
71     ! Used by CO2 initialization
72           REAL   ,SAVE, DIMENSION(109) :: PA, XA, CA, ETA, SEXPV
73           REAL   ,SAVE, DIMENSION(109,109) :: TRANSA
74           REAL   ,SAVE  :: CORE,UEXP,SEXP
75     
76           REAL,SAVE,DIMENSION(4) :: PTOPC
77     !
78     !--- Used for Gaussian look up tables
79     !
80           REAL, PRIVATE,PARAMETER :: XSDmax=3.1, DXSD=.01
81           INTEGER, PRIVATE,PARAMETER :: NXSD=XSDmax/DXSD
82           REAL, DIMENSION(NXSD),PRIVATE,SAVE :: AXSD
83           REAL, PRIVATE :: RSQR
84           LOGICAL, PRIVATE,SAVE :: SDprint=.FALSE.
85     !
86     !--- Important parameters for cloud properties - see extensive comments in
87     !    DO 580 loop within subroutine RADTN 
88     !
89           REAL, PARAMETER ::  &
90          &   TRAD_ice=0.5*T_ice      & !--- Very tunable parameter
91          &,  TRADK_ice=TRAD_ice+273.15   & !--- Very tunable parameter
92          &,  ABSCOEF_W=800.          & !--- Very tunable parameter
93          &,  ABSCOEF_I=500.          & !--- Very tunable parameter
94          &,  SECANG=-1.66            & !--- Very tunable parameter
95          &,  CLDCOEF_LW=1.5          & !--- Enhance LW cloud depths
96          &,  ABSCOEF_LW=SECANG*CLDCOEF_LW  & !--- Final factor for cloud emissivities
97          &,  Qconv=0.1e-3            & !--- Very tunable parameter
98          &,  CTauCW=ABSCOEF_W*Qconv  &
99          &,  CTauCI=ABSCOEF_I*Qconv
100     !
101     !-----------------------------------------------------------------------
102     !  Assign co2 and trace gases amount (units are parts/part by volumn)
103     !
104        REAL,PARAMETER :: co2=300.e-6
105     !
106     
107     
108     !-----------------------------------------------------------------------
109     !
110           CONTAINS
111     !
112     !-----------------------------------------------------------------------
113     !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
114     !-----------------------------------------------------------------------
115           SUBROUTINE RDTEMP(NTSD,DT,JULDAY,JULYR,IHRST,GLAT,GLON            &
116                            ,CZEN,CZMEAN,T,RSWTT,RLWTT                       &
117                            ,IDS,IDE,JDS,JDE,LM                              &
118                            ,IMS,IME,JMS,JME                                 &
119                            ,ITS,ITE,JTS,JTE                                 &
120                            ,ITS_B1,ITE_B1,JTS_B1,JTE_B1)
121     !***********************************************************************
122     !$$$  SUBPROGRAM DOCUMENTATION BLOCK
123     !                .      .    .     
124     ! SUBPROGRAM:    RDTEMP      RADIATIVE TEMPERATURE CHANGE
125     !   PRGRMMR: BLACK           ORG: W/NP22     DATE: 93-12-29
126     !     
127     ! ABSTRACT:
128     !     RDTEMP APPLIES THE TEMPERATURE TENDENCIES DUE TO
129     !     RADIATION AT ALL LAYERS AT EACH ADJUSTMENT TIME STEP
130     !     
131     ! PROGRAM HISTORY LOG:
132     !   87-09-??  BLACK      - ORIGINATOR
133     !   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
134     !   95-11-20  ABELES     - PARALLEL OPTIMIZATION
135     !   98-10-30  BLACK      - MODIFIED FOR DISTRIBUTED MEMORY
136     !   02-06-07  BLACK      - WRF CODING STANDARDS
137     !   02-09-09  WOLFE      - CONVERTING TO GLOBAL INDEXING
138     !   06-08-29  BLACK      - INTO NMMB
139     !     
140     ! USAGE: CALL RDTEMP FROM SUBROUTINE PHY_RUN
141     !  
142     ! ATTRIBUTES:
143     !   LANGUAGE: FORTRAN 90
144     !   MACHINE : IBM 
145     !$$$  
146     !-----------------------------------------------------------------------
147     !
148           IMPLICIT NONE
149     !
150     !-----------------------------------------------------------------------
151     !
152           INTEGER,INTENT(IN) :: IDE,IDS,IME,IMS,ITE,ITE_B1,ITS,ITS_B1       &
153                                ,JDE,JDS,JME,JMS,JTE,JTE_B1,JTS,JTS_B1       &
154                                ,LM
155     !
156           INTEGER,INTENT(IN) :: IHRST,JULDAY,JULYR,NTSD
157     !
158           REAL,INTENT(IN) :: DT
159     !
160           REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CZMEAN,GLAT,GLON 
161     !
162           REAL,DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(IN) :: RLWTT,RSWTT
163     !
164           REAL,DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(INOUT) :: T
165     !
166           REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: CZEN
167     !
168     !-----------------------------------------------------------------------
169     !***  LOCAL VARIABLES
170     !-----------------------------------------------------------------------
171     !
172           INTEGER :: I,J,JDAY,JMONTH,K
173     !
174           INTEGER,DIMENSION(3) :: IDAT
175     !
176           REAL :: DAYI,HOUR,TIMES,TTNDKL
177     !
178           REAL,DIMENSION(IMS:IME,JMS:JME) :: CZEN2,XLAT2,XLON2
179     !
180           REAL,DIMENSION(ITS:ITE,JTS:JTE) :: FACTR
181     !
182           REAL :: DEGRAD=3.1415926/180.
183           real :: xlat1,xlon1
184     !
185     !-----------------------------------------------------------------------
186     !-----------------------------------------------------------------------
187     !
188     !***  GET CURRENT VALUE OF COS(ZENITH ANGLE)
189     !
190           TIMES=NTSD*DT
191     !
192           DO J=JTS,JTE
193           DO I=ITS,ITE
194             XLAT2(I,J)=GLAT(I,J)
195             XLON2(I,J)=GLON(I,J)
196           ENDDO
197           ENDDO
198     !
199           CALL CAL_MON_DAY(JULDAY,JULYR,JMONTH,JDAY)
200     
201           IDAT(1)=JMONTH
202           IDAT(2)=JDAY
203           IDAT(3)=JULYR
204     !
205           CALL ZENITH(TIMES,DAYI,HOUR,IDAT,IHRST,XLON2,XLAT2,CZEN2          &
206          &           ,ITS,ITE,JTS,JTE                                       &
207          &           ,IDS,IDE,JDS,JDE,1,LM+1                                &
208          &           ,IMS,IME,JMS,JME,1,LM+1                                &
209          &           ,ITS,ITE,JTS,JTE,1,LM)
210     !
211           DO J=JTS,JTE
212           DO I=ITS,ITE
213             CZEN(I,J)=CZEN2(I,J)
214             IF(CZMEAN(I,J)>0.)THEN 
215               FACTR(I,J)=CZEN(I,J)/CZMEAN(I,J)
216             ELSE
217               FACTR(I,J)=0.
218             ENDIF
219           ENDDO
220           ENDDO
221     !
222           DO K=1,LM
223             DO J=JTS_B1,JTE_B1
224             DO I=ITS_B1,ITE_B1
225               TTNDKL=RSWTT(I,J,K)*FACTR(I,J)+RLWTT(I,J,K)
226               T(I,J,K)=T(I,J,K)+TTNDKL*DT
227             ENDDO
228             ENDDO
229           ENDDO
230     !-----------------------------------------------------------------------
231     !
232           END SUBROUTINE RDTEMP
233     !
234     !-----------------------------------------------------------------------
235     !-----------------------------------------------------------------------
236     !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
237     !-----------------------------------------------------------------------
238           SUBROUTINE TIME_MEASURE(START_YEAR,START_MONTH,START_DAY          &
239                                  ,START_HOUR,START_MINUTE,START_SECOND      &
240                                  ,NTIMESTEP,DT                              &
241                                  ,JULDAY,JULYR,JULIAN,XTIME)
242     !-----------------------------------------------------------------------
243     !
244           IMPLICIT NONE
245     !
246     !-----------------------------------------------------------------------
247     !
248           INTEGER,INTENT(IN) :: START_YEAR,START_MONTH,START_DAY,START_HOUR &
249          &                     ,START_MINUTE,START_SECOND,NTIMESTEP
250     !
251           REAL,INTENT(IN) :: DT
252     !
253           INTEGER,INTENT(OUT) :: JULDAY,JULYR
254     !
255           REAL,INTENT(OUT) :: JULIAN,XTIME
256     !
257     !-----------------------------------------------------------------------
258     !***  LOCAL VARIABLES
259     !-----------------------------------------------------------------------
260     !
261           INTEGER :: N
262     !
263           INTEGER,DIMENSION(12),SAVE :: MONTH=(/31,28,31,30,31,30           &
264                                                ,31,31,30,31,30,31/)
265     !
266           REAL :: SUM
267     !
268     !-----------------------------------------------------------------------
269     !***********************************************************************
270     !-----------------------------------------------------------------------
271     !
272           JULYR=START_YEAR
273     !
274           IF(MOD(START_YEAR,4)==0)THEN
275             MONTH(2)=29
276           ENDIF
277     !
278                 JULDAY=0
279             julcount: DO N=1,12
280                     IF(N==START_MONTH)EXIT julcount
281                     JULDAY=JULDAY+MONTH(N)
282             ENDDO julcount
283     !
284           JULDAY=JULDAY+START_DAY  ! The day of the year the forecast begins
285                                    ! 12Z 2 January --> 2
286     !
287           SUM=(START_HOUR+(START_MINUTE+START_SECOND/60.)/60.)/24.
288           JULIAN=JULDAY-(1.-SUM)  ! The exact day the forecast begins
289                                   ! 12Z 2 January --> 1.5
290     !
291           XTIME=NTIMESTEP*DT/60. ! Minutes since start of forecast
292     !
293     !-----------------------------------------------------------------------
294     !
295           END SUBROUTINE TIME_MEASURE
296     !
297     !---------------------------------------------------------------------
298     !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
299     !---------------------------------------------------------------------
300           SUBROUTINE GFDL_INIT(EMISS,SFULL,SHALF,PPTOP,                     &
301          &                     JULYR,MONTH,IDAY,GMT,                        &
302          &                     CO2TF,                                       &
303          &                     IDS, IDE, JDS, JDE, KDS, KDE,                &
304          &                     IMS, IME, JMS, JME, KMS, KME,                &
305          &                     ITS, ITE, JTS, JTE, KTS, KTE              )
306     !-----------------------------------------------------------------------
307           IMPLICIT NONE
308     !-----------------------------------------------------------------------
309           INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                     &
310          &                     ,IMS,IME,JMS,JME,KMS,KME                     &
311          &                     ,ITS,ITE,JTS,JTE,KTS,KTE
312           INTEGER,INTENT(IN) :: JULYR,MONTH,IDAY,CO2TF
313           REAL,INTENT(IN) :: GMT,PPTOP
314           REAL,DIMENSION(KMS:KME),INTENT(IN) :: SFULL, SHALF
315           REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: EMISS
316     !
317           INTEGER :: I,IHRST,J,N
318           REAL :: PCLD,XSD,PI,SQR2PI
319           REAL :: SSLP=1013.25
320           REAL, PARAMETER :: PTOP_HI=150.,PTOP_MID=350.,PTOP_LO=642.,       &
321          &                   PLBTM=105000.
322     !-----------------------------------------------------------------------
323     !***********************************************************************
324     !-----------------------------------------------------------------------
325     !
326     !--- In case Ferrier microphysics is not called, initialize lookup tables for 
327     !    saturation vapor pressures (only FPVS is used in radiation, which calculates
328     !    vapor pressure w/r/t water for T>=0C and w/r/t ice for T<0C).
329     !
330             CALL GPVS
331     !
332     !***  INITIALIZE DIAGNOSTIC LOW,MIDDLE,HIGH CLOUD LAYER PRESSURE LIMITS.
333     !
334           LTOP(1)=0
335           LTOP(2)=0
336           LTOP(3)=0
337     !
338           DO N=1,KTE
339             PCLD=(SSLP-PPTOP*10.)*SHALF(N)+PPTOP*10.
340             IF(PCLD>=PTOP_LO)LTOP(1)=N
341             IF(PCLD>=PTOP_MID)LTOP(2)=N
342             IF(PCLD>=PTOP_HI)LTOP(3)=N
343     !       PRINT *,N,PCLD,SHALF(N),PSTAR,PPTOP
344           ENDDO
345     !***  
346     !***  ASSIGN THE PRESSURES FOR CLOUD DOMAIN BOUNDARIES
347     !***
348           PTOPC(1)=PLBTM
349           PTOPC(2)=PTOP_LO*100.
350           PTOPC(3)=PTOP_MID*100.
351           PTOPC(4)=PTOP_HI*100.
352     !
353     !***  USE CALL TO CONRAD FOR DIRECT READ OF CO2 FUNCTIONS
354     !***  OTHERWISE CALL CO2O3.
355     !
356           IF(CO2TF==1)THEN
357             CALL CO2O3(SFULL,SHALF,PPTOP,KME-KMS,KME-KMS+1,KME-KMS+2)
358           ELSE
359             CALL CONRAD(KDS,KDE,KMS,KME,KTS,KTE)
360           ENDIF
361     !
362           CALL O3CLIM
363     !	write(0,*) 'call TABLE'
364           CALL TABLE
365     !	write(0,*) 'return call TABLE'
366           IHRST=NINT(GMT)
367           CALL SOLARD(IHRST,IDAY,MONTH,JULYR)
368     !
369     !***  FOR NOW, GFDL RADIATION ASSUMES EMISSIVITY = 1.0
370     !
371           DO J=JTS,JTE
372           DO I=ITS,ITE
373             EMISS(I,J) = 1.0
374           ENDDO
375           ENDDO
376     !
377     !---  Calculate the area under the Gaussian curve at the start of the
378     !---  model run and build the look up table AXSD
379     !
380           PI=ACOS(-1.)
381           SQR2PI=SQRT(2.*PI)
382           RSQR=1./SQR2PI
383           DO I=1,NXSD
384             XSD=REAL(I)*DXSD
385             AXSD(I)=GAUSIN(XSD)
386             if (SDprint) print *,'I, XSD, AXSD =',I,XSD,AXSD(I)
387           ENDDO
388     !! !***
389     !! !***  MESO STANDARD DEVIATION OF EK AND MAHRT'S CLOUD COVER ALOGRITHM
390     !! !***
391            if (SDprint) print *,                                            &
392          & 'RHgrd,T_ICE,NLImin,NLImax,FLARGE1,FLARGE2,MDImin,MDImax=',&
393          &  RHgrd,T_ICE,NLImin,NLImax,FLARGE1,FLARGE2,MDImin,MDImax
394     !
395     !-----------------------------------------------------------------------
396           END SUBROUTINE GFDL_INIT
397     !-----------------------------------------------------------------------
398     !
399     !
400     !-----------------------------------------------------------------------
401           SUBROUTINE GFDL(DT,THRATEN,THRATENLW,THRATENSW,CLDFRA,PI3D        &
402          &                ,XLAND,P8W,T                                      &
403          &                ,QV,QW,QI,QS                                      & 
404          &                ,TSK2D,GLW,RSWIN,GSW,RSWINC                       &
405          &                ,RSWTOA,RLWTOA,CZMEAN                             & 
406          &                ,GLAT,GLON,HTOP,HBOT,ALBEDO,CUPPT                 &
407          &                ,SNOW,G,GMT                                       &
408     !BSF => for NAMX changes, pass in surface emissivity (SFCEMS) [different for snow]
409          &                ,NSTEPRA,NPHS,ITIMESTEP                           &
410          &                ,XTIME,JULIAN                                     &
411          &                ,JULYR,JULDAY,GFDL_LW,GFDL_SW                     &
412          &                ,CFRACL,CFRACM,CFRACH                             &
413          &                ,ACFRST,NCFRST,ACFRCV,NCFRCV                      &
414          &                ,IDS,IDE,JDS,JDE,KDS,KDE                          &
415          &                ,IMS,IME,JMS,JME,KMS,KME                          &
416          &                ,ITS,ITE,JTS,JTE,KTS,KTE)
417     !-----------------------------------------------------------------------
418           IMPLICIT NONE
419     !-----------------------------------------------------------------------
420           INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                     &
421          &                     ,IMS,IME,JMS,JME,KMS,KME                     &
422          &                     ,ITS,ITE,JTS,JTE,KTS,KTE,ITIMESTEP           &
423          &                     ,NPHS,NSTEPRA
424      
425           INTEGER,INTENT(IN) :: julyr,julday   
426           INTEGER,INTENT(INOUT),DIMENSION(ims:ime,jms:jme) :: NCFRST        &
427                                                              ,NCFRCV
428           REAL,INTENT(IN) :: DT,GMT,G,XTIME,JULIAN
429     
430           REAL,INTENT(INOUT),DIMENSION(ims:ime, jms:jme, kts:kte):: CLDFRA
431           REAL,INTENT(INOUT),DIMENSION(ims:ime, jms:jme, kts:kte):: THRATEN &
432                                                        ,THRATENLW,THRATENSW
433           REAL,INTENT(IN),DIMENSION(ims:ime, jms:jme, kms:kme)::   p8w
434           REAL,INTENT(IN),DIMENSION(ims:ime, jms:jme, kts:kte)::     t,     &
435          &                                                          qs,     &
436          &                                                          qv,     &
437          &                                                          qw,     &
438          &                                                        PI3D
439           REAL, INTENT(IN), DIMENSION(ims:ime, jms:jme):: ALBEDO,SNOW,      &
440          &                                                TSK2D,XLAND
441           REAL, INTENT(IN), DIMENSION(ims:ime, jms:jme):: GLAT,GLON
442           REAL, INTENT(INOUT), DIMENSION(ims:ime, jms:jme):: HTOP,HBOT,CUPPT
443           REAL, INTENT(INOUT), DIMENSION(ims:ime, jms:jme):: RSWTOA,        &
444          &                                                   RLWTOA,        &
445          &                                                   ACFRST,        &
446          &                                                   ACFRCV
447           REAL,INTENT(INOUT),DIMENSION(ims:ime, jms:jme):: GLW,GSW
448           REAL,INTENT(OUT),DIMENSION(ims:ime, jms:jme):: CZMEAN             &
449          &                                           ,RSWIN,RSWINC        &
450          &                                           ,CFRACL,CFRACM,CFRACH
451           LOGICAL, INTENT(IN) :: gfdl_lw,gfdl_sw
452           REAL, OPTIONAL, INTENT(IN), DIMENSION(ims:ime, jms:jme, kts:kte):: QI
453     
454           REAL, DIMENSION(ims:ime, jms:jme, kms:kme):: PFLIP,QIFLIP,QFLIP,  &
455          &                                             QWFLIP
456           REAL, DIMENSION(its:ite, jts:jte, kts:kte)::TENDS,TENDL
457           REAL, DIMENSION(ims:ime, jms:jme):: CUTOP,CUBOT
458           INTEGER :: IDAT(3),Jmonth,Jday,IHOUR
459           INTEGER :: I,J,K,IHRST
460     ! begin debugging radiation
461           integer :: imd,jmd
462           real :: FSWrat
463     ! end debugging radiation
464     !-----------------------------------------------------------------------
465     !***********************************************************************
466     !-----------------------------------------------------------------------
467     
468           IF(GFDL_LW.AND.GFDL_SW )GO TO 100
469     !
470           DO K=KTS,KTE
471             DO J=JTS,JTE
472               DO I=ITS,ITE
473                 CLDFRA(I,J,K)=0.
474               ENDDO
475             ENDDO
476           ENDDO
477     !
478     !- Note that the effects of rain are ignored in this radiation package (BSF 2005-01-25)
479     !
480           DO K=KTS,KTE
481             DO J=JTS,JTE
482             DO I=ITS,ITE
483               QFLIP (I,J,K)=MAX(0.,QV(I,J,K)/(1.+QV(I,J,K)))
484               QWFLIP(I,J,K)=MAX(QW(I,J,K),0.)      !Modified
485     ! Note that QIFLIP will contain QS+QI if both are passed in, otherwise just QS 
486     !     Eta MP now outputs QS instead of QI (JD 2006-05-12)
487               QIFLIP(I,J,K)=MAX(QS(I,J,K),0.)
488               IF(PRESENT(QI))QIFLIP(I,J,K)=QIFLIP(I,J,K)+QI(I,J,K)
489     !
490     !***  USE MONOTONIC HYDROSTATIC PRESSURE INTERPOLATED TO MID-LEVEL
491     !
492               PFLIP(I,J,K)=0.5*(P8W(I,J,K)+P8W(I,J,K+1))
493             ENDDO
494             ENDDO
495           ENDDO
496     !
497           DO J=JTS,JTE
498           DO I=ITS,ITE
499             CUBOT(I,J)=KTE+1-HBOT(I,J)
500             CUTOP(I,J)=KTE+1-HTOP(I,J)
501           ENDDO
502           ENDDO
503     !
504           CALL CAL_MON_DAY(JULDAY,JULYR,JMONTH,JDAY)     
505     !
506           IDAT(1)=JMONTH
507           IDAT(2)=JDAY
508           IDAT(3)=JULYR
509     
510           IHRST  =NINT(GMT)
511           IHOUR  =MOD((IHRST+NINT(XTIME/60.0)),24)
512           CALL SOLARD(IHOUR,JDAY,JMONTH,JULYR)
513     
514     !-----------------------------------------------------------------------
515           CALL RADTN (DT,T,QFLIP,QWFLIP,QIFLIP,                             &
516          &            PFLIP,P8W,XLAND,TSK2D,                                &
517          &            GLAT,GLON,CUTOP,CUBOT,ALBEDO,CUPPT,                   &
518          &            ACFRCV,NCFRCV,ACFRST,NCFRST,                          &
519          &            SNOW,GLW,GSW,RSWIN,RSWINC,                            &
520     !BSF => for NAMX changes, pass in surface emissivity (SFCEMS) [different for snow]
521          &            IDAT,IHRST,XTIME,JULIAN,                              &
522          &            NSTEPRA,NSTEPRA,NPHS,ITIMESTEP,                       &
523          &            TENDS,TENDL,CLDFRA,RSWTOA,RLWTOA,CZMEAN,              &
524          &            CFRACL,CFRACM,CFRACH,                                 &
525          &            IDS,IDE,JDS,JDE,KDS,KDE,                              &
526          &            IMS,IME,JMS,JME,KMS,KME,                              &
527          &            ITS,ITE,JTS,JTE,KTS,KTE                              )
528     !-----------------------------------------------------------------------
529     ! begin debugging radiation
530     !     imd=(ims+ime)/2
531     !     jmd=(jms+jme)/2
532     !     FSWrat=0.
533     !     if (RSWIN(imd,jmd) .gt. 0.)   &
534     !        FSWrat=(RSWIN(imd,jmd)-GSW(imd,jmd))/RSWIN(imd,jmd)
535     !     write(6,"(2a,2i5,5f9.2,f8.4,i3,2f8.4)") & 
536     !       '{rad4 imd,jmd,GSW,RSWIN,RSWOUT=RSWIN-GSW,RSWINC,GLW,' &
537     !      ,'ACFRCV,NCFRCV,ALBEDO,RSWOUT/RSWIN = '   &
538     !      ,imd,jmd, GSW(imd,jmd),RSWIN(imd,jmd)  &
539     !      ,RSWIN(imd,jmd)-GSW(imd,jmd),RSWINC(imd,jmd),GLW(imd,jmd) &
540     !      ,ACFRCV(imd,jmd),NCFRCV(imd,jmd),ALBEDO(imd,jmd),FSWrat
541     ! end debugging radiation
542     !
543     !--- Need to save LW & SW tendencies since radiation calculates both and this block
544     !    is skipped when GFDL SW is called, both only if GFDL LW is also called
545     !    
546     
547     
548            
549           IF(GFDL_LW)THEN
550             DO K = KTS,KTE
551               DO J=JTS,JTE
552               DO I=ITS,ITE
553                 THRATENLW(I,J,K)=TENDL(I,J,K)/PI3D(I,J,K)
554                 THRATENSW(I,J,K)=TENDS(I,J,K)/PI3D(I,J,K)
555                 THRATEN(I,J,K)  =THRATEN(I,J,K) + THRATENLW(I,J,K)
556               ENDDO
557               ENDDO
558             ENDDO
559           ENDIF
560     !
561     !*** THIS ASSUMES THAT LONGWAVE IS CALLED FIRST IN THE RADIATION_DRIVER.
562     !    Only gets executed if a different LW scheme (not GFDL) is called
563     !
564           IF(GFDL_SW)THEN
565             DO K=KTS,KTE
566               DO J=JTS,JTE
567               DO I=ITS,ITE
568                 THRATENSW(I,J,K)=TENDS(I,J,K)/PI3D(I,J,K)
569               ENDDO
570               ENDDO
571             ENDDO
572           ENDIF
573     !
574     !***  RESET ACCUMULATED CONVECTIVE CLOUD TOP/BOT AND CONVECTIVE PRECIP
575     !***  FOR NEXT INTERVAL BETWEEN RADIATION CALLS
576     !
577           DO J=JTS,JTE
578           DO I=ITS,ITE
579             CUPPT(I,J)=0.
580           ENDDO
581           ENDDO
582     !
583       100 IF(GFDL_SW)THEN
584             DO K=KTS,KTE
585               DO J=JTS,JTE
586               DO I=ITS,ITE
587                 THRATEN(I,J,K)=THRATEN(I,J,K)+THRATENSW(I,J,K)
588               ENDDO
589               ENDDO
590             ENDDO
591           ENDIF
592     
593     
594     !
595       END SUBROUTINE GFDL
596     !
597     !-----------------------------------------------------------------------
598           SUBROUTINE RADTN(DT,T,Q,QCW,QICE,                                 &
599          &                 PFLIP,P8W,XLAND,TSK2D,                           &
600          &                 GLAT,GLON,CUTOP,CUBOT,ALB,CUPPT,                 &
601          &                 ACFRCV,NCFRCV,ACFRST,NCFRST,                     &
602          &                 SNO,GLW,GSW,RSWIN,RSWINC,                        &
603     !BSF => for NAMX changes, pass in surface emissivity (SFCEMS) [different for snow]
604          &                 IDAT,IHRST,XTIME,JULIAN,                         &
605          &                 NRADS,NRADL,NPHS,NTSD,                           &
606          &                 TENDS,TENDL,CLDFRA,RSWTOA,RLWTOA,CZMEAN,         &
607          &                 CFRACL,CFRACM,CFRACH,                            &
608          &                 ids,ide, jds,jde, kds,kde,                       &
609          &                 ims,ime, jms,jme, kms,kme,                       &
610          &                 its,ite, jts,jte, kts,kte                       )
611     !-----------------------------------------------------------------------
612           IMPLICIT NONE
613     !-----------------------------------------------------------------------
614     
615     ! GLAT : geodetic latitude in radians of the mass points on the computational grid.
616     
617     ! CZEN : instantaneous cosine of the solar zenith angle.
618     
619     ! CUTOP : (REAL) model layer number that is highest in the atmosphere
620     !        in which convective cloud occurred since the previous call to the
621     !        radiation driver.
622     
623     ! CUBOT : (REAL) model layer number that is lowest in the atmosphere
624     !        in which convective cloud occurred since the previous call to the
625     !        radiation driver.
626     
627     ! ALB  : is no longer used in the operational radiation.  Prior to 24 July 2001
628     !        ALB was the climatological albedo that was modified within RADTN to
629     !        account for vegetation fraction and snow.
630     !
631     ! ALB  : reintroduced as the dynamic albedo from LSM
632     
633     ! CUPPT: accumulated convective precipitation (meters) since the
634     !        last call to the radiation.
635     
636     ! TSK2D : skin temperature
637     
638     ! IHE and IHW are relative location indices needed to locate neighboring
639     !       points on the Eta's Arakawa E grid since arrays are indexed locally on
640     !       each MPI task rather than globally.  IHE refers to the adjacent grid
641     !       point (a V point) to the east of the mass point being considered.  IHW
642     !       is the adjacent grid point to the west of the given mass point.
643     
644     ! IRAD is a relic from older code that is no longer needed.
645     
646     ! ACFRCV : sum of the convective cloud fractions that were computed
647     !          during each call to the radiation between calls to the subroutines that
648     !          do the forecast output.
649     
650     ! NCFRCV : the total number of times in which the convective cloud
651     !          fraction was computed to be greater than zero in the radiation between
652     !          calls to the output routines.  In the post-processor, ACFRCV is divided
653     !          by NCFRCV to yield an average convective cloud fraction.
654     
655     !          ACFRST and NCFRST are the analogs for stratiform cloud cover.
656     
657     !          VEGFRC is the fraction of the gridbox with vegetation.
658     
659     !          LVL holds the number of model layers that lie below the ground surface
660     !          at each point.  Clearly for sigma coordinates LVL is zero everywhere.
661     
662     ! CTHK  :  an assumed maximum thickness of stratiform clouds currently set
663     !          to 20000 Pascals.  I think this is relevant for computing "low",
664     !          "middle", and "high" cloud fractions which are post-processed but which
665     !          do not feed back into the integration.
666     
667     ! IDAT  : a 3-element integer array holding the month, day, and year,
668     !        respectively, of the date for the start time of the free forecast.
669     
670     ! ABCFF : holds coefficients for various absorption bands.  You can see
671     !         where they are set in GFDLRD.F.
672     
673     ! LTOP  : a 3-element integer array holding the model layer that is at or
674     !         immediately below the specified pressure levels for the tops 
675     !         of "high" (15000 Pa), "middle" (35000 Pa), and "low" (64200 Pa) 
676     !         stratiform clouds.  These are for the diagnostic cloud layers 
677     !         needed in the output but not in the integration.
678     
679     ! NRADS : integer number of fundamental timesteps (our smallest
680     !         timestep, i.e., the one for inertial gravity wave adjustment) 
681     !         between updates of the shortwave tendencies.  
682     
683     ! NRADL : integer number of fundamental timesteps between updates of
684     !         the longwave tendencies.  
685     
686     ! NTSD   : integer counter of the fundamental timesteps that have
687     !         elapsed since the start of the forecast.
688     
689     ! GLW : incoming longwave radiation at the surface
690     ! GSW : NET (down minus up, or incoming minus outgoing) all-sky shortwave radiation at the surface
691     ! RSWIN  : total (clear + cloudy sky) incoming (downward) solar radiation at the surface
692     ! RSWINC : clear sky incoming (downward) solar radiation at the surface
693     
694     ! TENDS,TENDL : shortwave,longwave (respectively) temperature tendency
695     
696     ! CLDFRA : 3D cloud fraction
697     
698     ! RSWTOA, RLWTOA : outgoing shortwave, longwave (respectively) fluxes at top of atmosphere
699     
700     ! CZMEAN : time-average cosine of the zenith angle
701     
702     ! CFRACL,CFRACM,CFRACH : low, middle, & high (diagnosed) cloud fractions
703     
704     ! XTIME : time since simulation start (minutes)
705                                                                                                                                                   
706     ! JULIAN: Day of year (0.0 at 00Z Jan 1st)
707     
708     !**********************************************************************
709     !****************************** NOTE **********************************
710     !**********************************************************************
711     !*** DUE TO THE RESETTING OF CONVECTIVE PRECIP AND CONVECTIVE CLOUD
712     !*** TOPS AND BOTTOMS, SHORTWAVE MUST NOT BE CALLED LESS FREQUENTLY
713     !*** THAN LONGWAVE.
714     !**********************************************************************
715     !****************************** NOTE **********************************
716     !**********************************************************************
717     !-----------------------------------------------------------------------
718           INTEGER, INTENT(IN)        :: ids,ide, jds,jde, kds,kde ,         &
719          &                              ims,ime, jms,jme, kms,kme ,         &
720          &                              its,ite, jts,jte, kts,kte
721           INTEGER, INTENT(IN)        :: NRADS,NRADL,NTSD,NPHS 
722           REAL   , INTENT(IN)        :: DT,XTIME,JULIAN
723           INTEGER, INTENT(IN), DIMENSION(3) :: IDAT
724     !-----------------------------------------------------------------------
725           INTEGER            :: LM1,LP1,LM
726           INTEGER, INTENT(IN)               :: IHRST
727     !
728           REAL, PARAMETER :: EPSQ1=1.E-5,EPSQ=1.E-12,EPSO3=1.E-10,H0=0.     &
729          &, H1=1.,HALF=.5,T0C=273.15,CUPRATE=24.*1000.,HPINC=HALF*1.E1      &
730     !------------------------ For Clouds ----------------------------------
731          &, CLFRmin=0.01, TAUCmax=4.161                                     &
732     !--- Parameters used for new cloud cover scheme
733          &, XSDmin=-XSDmax, DXSD1=-DXSD, STSDM=0.01, CVSDM=.04              &
734          &, DXSD2=HALF*DXSD, DXSD2N=-DXSD2, PCLDY=0.25
735     !
736           INTEGER, PARAMETER :: NB=12,KSMUD=0
737           INTEGER,PARAMETER :: K15=SELECTED_REAL_KIND(15)
738           REAL (KIND=K15) :: DDX,EEX,PROD
739     !-----------------------------------------------------------------------
740           LOGICAL :: SHORT,LONG
741           LOGICAL :: BITX,BITY,BITZ,BITW,BIT1,BIT2,BITC,BITCP1,BITSP1
742           LOGICAL, SAVE :: CNCLD=.TRUE.
743           LOGICAL :: NEW_CLOUD
744     !-----------------------------------------------------------------------
745           REAL, INTENT(IN), DIMENSION(ims:ime,jms:jme) :: XLAND,TSK2D
746           REAL, INTENT(IN), DIMENSION(ims:ime, jms:jme, kms:kme):: P8W      &
747          &                                                        ,PFLIP    &
748          &                                                        ,Q,QCW    &
749          &                                                        ,QICE     &
750          &                                                        ,T
751     
752           REAL, INTENT(OUT), DIMENSION(ims:ime, jms:jme):: GLW,GSW,CZMEAN   &
753          &                                                ,RSWIN,RSWINC     &
754          &                                                ,CFRACL,CFRACM    &
755          &                                                ,CFRACH
756           REAL, INTENT(OUT),DIMENSION(ims:ime,jms:jme,kts:kte) :: CLDFRA
757     
758     !     REAL,   INTENT(IN), DIMENSION(kms:kme)   :: ETAD
759     !     REAL,   INTENT(IN), DIMENSION(kms:kme)   :: AETA
760     !-----------------------------------------------------------------------
761           REAL, INTENT(IN), DIMENSION(ims:ime,jms:jme) :: CUTOP,CUBOT,CUPPT
762           REAL,   INTENT(IN   ), DIMENSION(ims:ime,jms:jme)  :: ALB,SNO
763     !BSF => for NAMX changes, pass in surface emissivity (SFCEMS) [different for snow]
764           REAL,   INTENT(IN   ), DIMENSION(ims:ime,jms:jme)  :: GLAT,GLON
765     !-----------------------------------------------------------------------
766           REAL,   DIMENSION(ims:ime,jms:jme)  :: CZEN
767           INTEGER, DIMENSION(its:ite, jts:jte):: LMH
768     !-----------------------------------------------------------------------
769     !     INTEGER,INTENT(IN), DIMENSION(jms:jme) :: IHE,IHW
770     !-----------------------------------------------------------------------
771           REAL,   INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: ACFRCV,ACFRST &
772                                                               ,RSWTOA,RLWTOA
773           INTEGER,INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: NCFRCV,NCFRST
774     !-----------------------------------------------------------------------
775           REAL,   INTENT(INOUT),DIMENSION(its:ite,jts:jte,kts:kte) :: TENDL,&
776          &                                                            TENDS
777     !-----------------------------------------------------------------------
778           REAL :: CTHK(3)
779           DATA CTHK/20000.0,20000.0,20000.0/
780     
781           REAL,DIMENSION(10),SAVE :: CC,PPT
782     !-----------------------------------------------------------------------
783           REAL,SAVE :: ABCFF(NB)
784           INTEGER,DIMENSION(its:ite,jts:jte) :: LVL
785           REAL,   DIMENSION(its:ite, jts:jte):: PDSL,FNE,FSE,TL
786           REAL,   DIMENSION(  0:kte)  :: CLDAMT
787           REAL,   DIMENSION(its:ite,3):: CLDCFR
788           INTEGER,   DIMENSION(its:ite,3):: MBOT,MTOP
789           REAL,   DIMENSION(its:ite)  :: PSFC,TSKN,ALBEDO,XLAT,COSZ,        &
790          &                               SLMSK,FLWUP,                       &
791          &                               FSWDN,FSWUP,FSWDNS,FSWUPS,FLWDNS,  &
792          &                               FLWUPS,FSWDNSC
793     
794           REAL,   DIMENSION(its:ite,kts:kte) :: PMID,TMID
795           REAL,   DIMENSION(its:ite,kts:kte) :: QMID,THMID,OZN,POZN
796           REAL,   DIMENSION(its:ite,jts:jte) :: TOT 
797     
798           REAL,   DIMENSION(its:ite,kts:kte+1) :: PINT,EMIS,CAMT,TAUcld
799           INTEGER,DIMENSION(its:ite,kts:kte+1) :: KBTM,KTOP
800           INTEGER,DIMENSION(its:ite)   :: NCLDS,KCLD 
801           REAL,   DIMENSION(its:ite)   :: TAUDAR
802           REAL,   DIMENSION(its:ite,NB,kts:kte+1) ::RRCL,TTCL
803     
804           REAL,   DIMENSION(its:ite,kts:kte):: CSMID,CCMID,QWMID,QIMID
805           REAL,SAVE :: P400=40000.
806           INTEGER,SAVE :: NFILE=14
807     
808     !-----------------------------------------------------------------------
809           REAL    :: CLSTP,TIME,DAYI,HOUR,ADDL,RANG
810           REAL    :: TIMES,EXNER,APES,SNOFAC,CCLIMIT,CLIMIT,P1,P2,CC1,CC2
811           REAL    :: PMOD,CLFR1,CTAU,WV,ARG,CLDMAX
812           REAL    :: CL1,CL2,CR1,DPCL,QSUM,PRS1,PRS2,DELP,TCLD,DD,EE,AA,FF
813           REAL    :: BB,GG,FCTR,PDSLIJ,CFRAVG,SNOMM
814           REAL    :: THICK,CONVPRATE,CLFR,ESAT,QSAT,RHUM,QCLD
815           REAL    :: RHtot,RRHO,FLARGE,FSMALL,DSNOW,SDM,QPCLDY,DIFCLD
816           REAL    :: TauC,CTauL,CTauS,  CFSmax,CFCmax
817           INTEGER :: I,J,MYJS,MYJE,MYIS,MYIE,NTSPH,NRADPP,ITIMSW,ITIMLW,    &
818          &           JD,II
819           INTEGER :: L,N,LML,LVLIJ,IR,KNTLYR,LL,NC,L400,NMOD,LTROP,IWKL
820           INTEGER :: LCNVB,LCNVT
821           INTEGER :: NLVL,MALVL,LLTOP,LLBOT,KBT2,KTH1,KBT1,KTH2,KTOP1
822           INTEGER :: NBAND,NCLD,LBASE,NKTP,NBTM,KS
823           INTEGER :: INDEXS,IXSD
824           DATA    CC/0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/
825           DATA    PPT/0.,.14,.31,.70,1.6,3.4,7.7,17.,38.,85./
826           DATA ABCFF/2*4.0E-5,.002,.035,.377,1.95,9.40,44.6,190.,989.,      &
827          &           2706.,39011./
828     ! begin debugging radiation
829           integer :: imd,jmd, Jndx
830           real :: FSWrat
831           imd=(ims+ime)/2
832           jmd=(jms+jme)/2
833     ! end debugging radiation
834     !
835     !=======================================================================
836     !
837           MYJS=jts
838           MYJE=jte
839           MYIS=its
840           MYIE=ite
841           LM=kte
842           LM1=LM-1
843           LP1=LM+1
844     !
845           DO J=JTS,JTE
846           DO I=ITS,ITE
847             LMH(I,J)=KME-1
848             LVL(I,J)=0
849           ENDDO
850           ENDDO
851     !**********************************************************************
852     !***  THE FOLLOWING CODE IS EXECUTED EACH TIME THE RADIATION IS CALLED.
853     !**********************************************************************
854     !----------------------CONVECTION--------------------------------------
855     !  NRADPP IS THE NUMBER OF TIME STEPS TO ACCUMULATE CONVECTIVE PRECIP
856     !     FOR RADIATION
857     !   NOTE: THIS WILL NOT WORK IF NRADS AND NRADL ARE DIFFERENT UNLESS
858     !         THEY ARE INTEGER MULTIPLES OF EACH OTHER
859     !  CLSTP IS THE NUMBER OF HOURS OF THE ACCUMULATION PERIOD
860     !
861           NTSPH=NINT(3600./DT)
862           NRADPP=MIN(NRADS,NRADL)
863           CLSTP=1.0*NRADPP/NTSPH
864           CONVPRATE=CUPRATE/CLSTP
865     !----------------------CONVECTION--------------------------------------
866     !***
867     !***  STATE WHETHER THE SHORT OR LONGWAVE COMPUTATIONS ARE TO BE DONE.
868     !***
869           SHORT=.TRUE. 
870           LONG=.TRUE. 
871           ITIMSW=0
872           ITIMLW=0
873           IF(SHORT)ITIMSW=1
874           IF(LONG) ITIMLW=1
875     !***
876     !***  FIND THE MEAN COSINE OF THE SOLAR ZENITH ANGLE 
877     !***  BETWEEN THE CURRENT TIME AND THE NEXT TIME RADIATION IS
878     !***  CALLED.  ONLY AVERAGE IF THE SUN IS ABOVE THE HORIZON.
879     !***
880           TIME=XTIME*60.
881     !-----------------------------------------------------------------------
882           CALL ZENITH(TIME,DAYI,HOUR,IDAT,IHRST,GLON,GLAT,CZEN,             &
883          &            ITS,ITE,JTS,JTE,                                      &
884          &            ids,ide, jds,jde, kds,kde,                            &
885          &            ims,ime, jms,jme, kms,kme,                            &
886          &            its,ite, jts,jte, kts,kte                         )
887     !-----------------------------------------------------------------------
888     !
889           ADDL=0.
890           IF(MOD(IDAT(3),4).EQ.0)ADDL=1.
891           RANG=2.*PI*(DAYI-RLAG)/(365.+ADDL)
892           RSIN1=SIN(RANG)
893           RCOS1=COS(RANG)
894           RCOS2=COS(2.*RANG)
895     !
896     !-----------------------------------------------------------------------
897           IF(SHORT)THEN
898             DO J=MYJS,MYJE
899             DO I=MYIS,MYIE
900               CZMEAN(I,J)=0.
901               TOT(I,J)=0.
902             ENDDO
903             ENDDO
904     !
905             DO II=0,NRADS,NPHS
906               TIMES=XTIME*60.+II*DT
907               CALL ZENITH(TIMES,DAYI,HOUR,IDAT,IHRST,GLON,GLAT,CZEN,        &
908          &                ITS,ITE,JTS,JTE,                                  &
909          &                ids,ide, jds,jde, kds,kde,                        &
910          &                ims,ime, jms,jme, kms,kme,                        &
911          &                its,ite, jts,jte, kts,kte                       )
912     
913               DO J=MYJS,MYJE
914               DO I=MYIS,MYIE
915                 IF(CZEN(I,J).GT.0.)THEN
916                   CZMEAN(I,J)=CZMEAN(I,J)+CZEN(I,J)
917                   TOT(I,J)=TOT(I,J)+1.
918                 ENDIF
919               ENDDO
920               ENDDO
921             ENDDO
922             DO J=MYJS,MYJE
923             DO I=MYIS,MYIE
924               IF(TOT(I,J).GT.0.)CZMEAN(I,J)=CZMEAN(I,J)/TOT(I,J)
925             ENDDO
926             ENDDO
927           ENDIF
928     !
929     !***  Do not modify pressure for ozone concentrations below the top layer
930     !***
931           DO L=2,LM
932           DO I=MYIS,MYIE
933             POZN(I,L)=H1
934           ENDDO
935           ENDDO
936     !-----------------------------------------------------------------------
937     !
938     !***********************************************************************
939     !***  THIS IS THE BEGINNING OF THE PRIMARY LOOP THROUGH THE DOMAIN
940     !***********************************************************************
941     !                        *********************
942                              DO 700 J = MYJS, MYJE
943     !                        *********************
944           DO 125 L=1,LM
945           DO I=MYIS,MYIE
946             TMID(I,L)=T(I,J,1)
947             QMID(I,L)=EPSQ
948             QWMID(I,L)=0.
949             QIMID(I,L)=0.
950             CSMID(I,L)=0.
951             CCMID(I,L)=0.
952             OZN(I,L)=EPSO3
953             TENDS(I,J,L)=0.
954             TENDL(I,J,L)=0.
955           ENDDO
956       125 CONTINUE
957     !
958           DO 140 N=1,3
959           DO I=MYIS,MYIE
960             CLDCFR(I,N)=0.
961             MTOP(I,N)=0
962             MBOT(I,N)=0
963           ENDDO
964       140 CONTINUE
965     !***
966     !***  FILL IN WORKING ARRAYS WHERE VALUES AT L=LM ARE THOSE THAT
967     !***  ARE ACTUALLY AT ETA LEVEL L=LMH.
968     !***
969           DO 200 I=MYIS,MYIE
970           LML=LMH(I,J)
971           LVLIJ=LVL(I,J)
972     !
973           DO L=1,LML
974             PMID(I,L+LVLIJ)=PFLIP(I,J,L)
975             PINT(I,L+LVLIJ+1)=P8W(I,J,L+1)
976             EXNER=(1.E5/PMID(I,L+LVLIJ))**CAPPA
977             TMID(I,L+LVLIJ)=T(I,J,L)
978             THMID(I,L+LVLIJ)=T(I,J,L)*EXNER
979             QMID(I,L+LVLIJ)=MAX(EPSQ, Q(I,J,L))
980     !--- Note that rain is ignored, only effects from cloud water and 
981     !    ice (cloud ice + snow) are considered
982             QWMID(I,L+LVLIJ)=QCW(I,J,L)
983             QIMID(I,L+LVLIJ)=QICE(I,J,L)
984           ENDDO
985     !***
986     !***  FILL IN ARTIFICIAL VALUES ABOVE THE TOP OF THE DOMAIN.
987     !***  PRESSURE DEPTHS OF THESE LAYERS IS 1 HPA.
988     !***  TEMPERATURES ABOVE ARE ALREADY ISOTHERMAL WITH (TRUE) LAYER 1.
989     !***
990           IF(LVLIJ.GT.0)THEN
991             KNTLYR=0
992     !
993             DO L=LVLIJ,1,-1
994               KNTLYR=KNTLYR+1
995               PMID(I,L)=P8W(I,J,1)-REAL(2*KNTLYR-1)*HPINC
996               PINT(I,L+1)=PMID(I,L)+HPINC
997               EXNER=(1.E5/PMID(I,L))**CAPPA
998               THMID(I,L)=TMID(I,L)*EXNER
999             ENDDO
1000           ENDIF
1001     !
1002           IF(LVLIJ.EQ.0) THEN
1003              PINT(I,1)=P8W(I,J,1)
1004           ELSE
1005              PINT(I,1)=PMID(I,1)-HPINC
1006           ENDIF
1007       200 CONTINUE
1008     !***
1009     !***  FILL IN THE SURFACE PRESSURE, SKIN TEMPERATURE, GEODETIC LATITUDE,
1010     !***  ZENITH ANGLE, SEA MASK, AND ALBEDO.  THE SKIN TEMPERATURE IS
1011     !***  NEGATIVE OVER WATER.
1012     !***
1013           DO 250 I=MYIS,MYIE
1014           PSFC(I)=P8W(I,J,KME)
1015           APES=(PSFC(I)*1.E-5)**CAPPA
1016           IF((XLAND(I,J)-1.5).GT.0.)THEN
1017             TSKN(I)=-TSK2D(I,J)
1018           ELSE
1019             TSKN(I)=TSK2D(I,J)
1020           ENDIF
1021     
1022           SLMSK(I)=XLAND(I,J)-1.
1023     !
1024     !     SNO(I,J)=AMAX1(SNO(I,J),0.)
1025     !BSF => for NAMX changes, pass in surface emissivity (SFCEMS) [different for snow]
1026           SNOMM=AMAX1(SNO(I,J),0.)
1027           SNOFAC=AMIN1(SNOMM/0.02, 1.0)
1028     !!!!  ALBEDO(I)=ALB(I,J)+(1.0-0.01*VEGFRC(I,J))*SNOFAC*(SNOALB-ALB(I,J))
1029           ALBEDO(I)=ALB(I,J)
1030     !
1031           XLAT(I)=GLAT(I,J)*RTD
1032           COSZ(I)=CZMEAN(I,J)
1033       250 CONTINUE
1034     !-----------------------------------------------------------------------
1035     !---  COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION  (Ferrier, Nov '04)
1036     !
1037     !--- Assumes Gaussian-distributed probability density functions (PDFs) for
1038     !    total relative humidity (RHtot) within the grid for convective and
1039     !    grid-scale cloud processes.  The standard deviation of RHtot is assumed
1040     !    to be larger for convective clouds than grid-scale (stratiform) clouds.
1041     !-----------------------------------------------------------------------
1042     !
1043           DO I=MYIS,MYIE
1044             LML=LMH(I,J)
1045             LVLIJ=LVL(I,J)
1046             DO 255 L=1,LML
1047                 LL=L+LVLIJ
1048                 WV=QMID(I,LL)/(1.-QMID(I,LL))       !--- Water vapor mixing ratio
1049                 QCLD=QWMID(I,LL)+QIMID(I,LL)        !--- Total cloud water + ice mixing ratio
1050                 IF (QCLD .LE. EPSQ) GO TO 255       !--- Skip if no condensate is present
1051                 CLFR=H0
1052                 WV=QMID(I,LL)/(1.-QMID(I,LL))       !--- Water vapor mixing ratio
1053                    
1054         !
1055         !--- Saturation vapor pressure w/r/t water ( >=0C ) or ice ( <0C )
1056         !
1057                 ESAT=1000.*FPVS(TMID(I,LL))         !--- Saturation vapor pressure (Pa)
1058                 QSAT=EP_2*ESAT/(PMID(I,LL)-ESAT)    !--- Saturation mixing ratio
1059                 RHUM=WV/QSAT                        !--- Relative humidity
1060         !
1061         !--- Revised cloud cover parameterization (temporarily ignore rain)
1062         !
1063                 RHtot=(WV+QCLD)/QSAT                !--- Total relative humidity
1064     !!    !
1065     !!    !--- QOVRCST is the amount of cloud condensate associated with full
1066     !!    !    overcast, PCLDY is an arbitrary factor for partial cloudiness
1067     !!    !
1068     !!            TCLD=TMID(I,LL)-T0C                 !--- Air temp in deg C
1069     !!            RRHO=(R_D*TMID(I,LL)*(1.+EP_1*QMID(I,LL)))/PMID(I,LL)
1070     !!            IF (TCLD .GE. 0.) THEN
1071     !!               QOVRCST(I,LL)=QAUT0*RRHO
1072     !!            ELSE
1073     !!               IF (TCLD.GE.-8. .AND. TCLD.LE.-3.) THEN
1074     !!                  FLARGE=FLARGE1
1075     !!               ELSE
1076     !!                  FLARGE=FLARGE2
1077     !!               ENDIF
1078     !!               FSMALL=(1.-FLARGE)/FLARGE
1079     !!               DSNOW=XMImax*EXP(XMIexp*TCLD)
1080     !!               INDEXS=MAX(MDImin, MIN(MDImax, INT(DSNOW)))
1081     !!               QOVRCST(I,LL)=NLImax*( FSMALL*MASSI(MDImin)              &
1082     !!     &                               +MASSI(INDEXS) )*RRHO
1083     !!            ENDIF                 !--- End IF (TCLD .GE. 0.)
1084     !!            QOVRCST(I,LL)=PCLDY*QOVRCST(I,LL)
1085                 LCNVT=NINT(CUTOP(I,J))+LVLIJ
1086                 LCNVT=MIN(LM,LCNVT)
1087                 LCNVB=NINT(CUBOT(I,J))+LVLIJ
1088                 LCNVB=MIN(LM,LCNVB)
1089                 IF (LL.GE.LCNVT .AND. LL.LE.LCNVB) THEN
1090                    SDM=CVSDM
1091                 ELSE
1092                    SDM=STSDM
1093                 ENDIF
1094                 ARG=(RHtot-RHgrd)/SDM
1095                 IF (ARG.LE.DXSD2 .AND. ARG.GE.DXSD2N) THEN
1096                    CLFR=HALF
1097                 ELSE IF (ARG .GT. DXSD2) THEN
1098                    IF (ARG .GE. XSDmax) THEN
1099                       CLFR=H1
1100                    ELSE
1101                       IXSD=INT(ARG/DXSD+HALF)
1102                       IXSD=MIN(NXSD, MAX(IXSD,1))
1103                       CLFR=HALF+AXSD(IXSD)
1104                    ENDIF              !--- End IF (ARG .GE. XSDmax)
1105                 ELSE
1106                    IF (ARG .LE. XSDmin) THEN
1107                       CLFR=H0
1108                    ELSE
1109                       IXSD=INT(ARG/DXSD1+HALF)
1110                       IXSD=MIN(NXSD, MAX(IXSD,1))
1111                       CLFR=HALF-AXSD(IXSD)
1112                       IF (CLFR .LT. CLFRmin) CLFR=H0
1113                    ENDIF        !--- End IF (ARG .LE. XSDmin) 
1114                 ENDIF           !--- IF (ARG.LE.DXSD2 .AND. ARG.GE.DXSD2N)
1115                 CSMID(I,LL)=CLFR
1116     !!  !
1117     !!  !--- Here the condensate is adjusted to be only over the cloudy area
1118     !!  !
1119     !!            IF (CLFR.GT.0. .AND. QCLD.LE.0.) THEN
1120     !!  !
1121     !!  !--- Put in modest amounts of cloud water & cloud ice for partially cloudy grids
1122     !!  !
1123     !!               QPCLDY=MIN(.01*QSAT, QOVRCST(I,LL))
1124     !!               IF (TCLD .GE. H0) THEN
1125     !!                  QWMID(I,LL)=QPCLDY
1126     !!               ELSE
1127     !!                  QIMID(I,LL)=QPCLDY
1128     !!               ENDIF
1129     !!            ENDIF          !--- End IF (CLFR.GT.0. .AND. QCLD.LE.0.) 
1130     255       CONTINUE         !--- End DO L=1,LML
1131           ENDDO                !--- End DO I=MYIS,MYIE
1132     !
1133     !***********************************************************************
1134     !******************  END OF GRID-SCALE CLOUD FRACTIONS  ****************
1135     !
1136     !---  COMPUTE CONVECTIVE CLOUD COVER FOR RADIATION 
1137     !
1138     !--- The parameterization of Slingo (1987, QJRMS, Table 1, p. 904) is 
1139     !    used for convective cloud fraction as a function of precipitation 
1140     !    rate.  Cloud fractions have been increased by 20% for each rainrate
1141     !    interval so that shallow, nonprecipitating convection is ascribed a
1142     !    constant cloud fraction of 0.1  (Ferrier, Feb '02).
1143     !***********************************************************************
1144     !
1145           IF (CNCLD) THEN
1146             DO I=MYIS,MYIE
1147     !
1148     !***  CLOUD TOPS AND BOTTOMS COME FROM CUCNVC
1149     !     Convective clouds need to be at least 2 model layers thick
1150     !
1151               IF (CUBOT(I,J)-CUTOP(I,J) .GT. 1.0) THEN
1152      !--- Compute convective cloud fractions if appropriate  (Ferrier, Feb '02)
1153                 CLFR=CC(1)
1154                 PMOD=CUPPT(I,J)*CONVPRATE
1155                 IF (PMOD .GT. PPT(1)) THEN
1156                   DO NC=1,10
1157                     IF(PMOD.GT.PPT(NC)) NMOD=NC
1158                   ENDDO
1159                   IF (NMOD .GE. 10) THEN
1160                     CLFR=CC(10)
1161                   ELSE
1162                     CC1=CC(NMOD)
1163                     CC2=CC(NMOD+1)
1164                     P1=PPT(NMOD)
1165                     P2=PPT(NMOD+1)
1166                     CLFR=CC1+(CC2-CC1)*(PMOD-P1)/(P2-P1)
1167                   ENDIF      !--- End IF (NMOD .GE. 10) ...
1168                   CLFR=MIN(H1, CLFR)
1169                 ENDIF        !--- End IF (PMOD .GT. PPT(1)) ...
1170       !
1171       !***  ADD LVL TO BE CONSISTENT WITH OTHER WORKING ARRAYS
1172       !
1173                 LVLIJ=LVL(I,J)
1174                 LCNVT=NINT(CUTOP(I,J))+LVLIJ
1175                 LCNVT=MIN(LM,LCNVT)
1176                 LCNVB=NINT(CUBOT(I,J))+LVLIJ
1177                 LCNVB=MIN(LM,LCNVB)
1178     !! !
1179     !! !---- For debugging
1180     !! !
1181     !!      WRITE(6,"(2(A,I3),2(A,I2),2(A,F5.2),2(A,I2),A,F6.4)") 
1182     !!     & ' J=',J,' I=',I,' LCNVB=',LCNVB,' LCNVT=',LCNVT
1183     !!     &, ' CUBOT=',CUBOT(I,J),' CUTOP=',CUTOP(I,J)
1184     !!     &,' LVL=',LVLIJ,' LMH=',LMH(I,J),' CCMID=',CLFR
1185     !! !
1186        !
1187        !--- Build in small amounts of subgrid-scale convective condensate 
1188        !    (simple assumptions), but only if the convective cloud fraction 
1189        !    exceeds that of the grid-scale cloud fraction
1190        !
1191                 DO LL=LCNVT,LCNVB
1192                   ARG=MAX(H0, H1-CSMID(I,LL))
1193                   CCMID(I,LL)=MIN(ARG,CLFR)
1194                 ENDDO           !--- End DO LL=LCNVT,LCNVB
1195               ENDIF             !--- IF (CUBOT(I,J)-CUTOP(I,J) .GT. 1.0) ...
1196             ENDDO               !--- End DO I loop
1197           ENDIF                 !--- End IF (CNCLD) ...
1198     !
1199     !*********************************************************************
1200     !***************  END OF CONVECTIVE CLOUD FRACTIONS  *****************
1201     !*********************************************************************
1202     !***
1203     !***  DETERMINE THE FRACTIONAL CLOUD COVERAGE FOR HIGH, MID
1204     !***  AND LOW OF CLOUDS FROM THE CLOUD COVERAGE AT EACH LEVEL
1205     !***
1206     !***  NOTE: THIS IS FOR DIAGNOSTICS ONLY!!!
1207     !***
1208     !***
1209            DO 500 I=MYIS,MYIE
1210     !!
1211            DO L=0,LM
1212              CLDAMT(L)=0.
1213            ENDDO
1214     !!  
1215     !!***  NOW GOES LOW, MIDDLE, HIGH
1216     !!
1217            DO 480 NLVL=1,3
1218            CLDMAX=0.
1219            MALVL=LM
1220            LLTOP=LM+1-LTOP(NLVL)+LVL(I,J)
1221     !!***
1222     !!***  GO TO THE NEXT CLOUD LAYER IF THE TOP OF THE CLOUD-TYPE IN
1223     !!***  QUESTION IS BELOW GROUND OR IS IN THE LOWEST LAYER ABOVE GROUND.
1224     !!***
1225            IF(LLTOP.GE.LM)GO TO 480
1226     !!
1227            IF(NLVL.GT.1)THEN
1228              LLBOT=LM+1-LTOP(NLVL-1)-1+LVL(I,J)
1229              LLBOT=MIN(LLBOT,LM1)
1230            ELSE
1231              LLBOT=LM1
1232            ENDIF
1233     !!
1234            DO 435 L=LLTOP,LLBOT
1235            CLDAMT(L)=AMAX1(CSMID(I,L),CCMID(I,L))
1236            IF(CLDAMT(L).GT.CLDMAX)THEN
1237              MALVL=L
1238              CLDMAX=CLDAMT(L)
1239            ENDIF
1240        435 CONTINUE
1241     !!*********************************************************************
1242     !! NOW, CALCULATE THE TOTAL CLOUD FRACTION IN THIS PRESSURE DOMAIN
1243     !! USING THE METHOD DEVELOPED BY Y.H., K.A.C. AND A.K. (NOV., 1992).
1244     !! IN THIS METHOD, IT IS ASSUMED THAT SEPERATED CLOUD LAYERS ARE
1245     !! RADOMLY OVERLAPPED AND ADJACENT CLOUD LAYERS ARE MAXIMUM OVERLAPPED.
1246     !! VERTICAL LOCATION OF EACH TYPE OF CLOUD IS DETERMINED BY THE THICKEST
1247     !! CONTINUING CLOUD LAYERS IN THE DOMAIN.
1248     !!*********************************************************************
1249            CL1=0.0
1250            CL2=0.0
1251            KBT1=LLBOT
1252            KBT2=LLBOT
1253            KTH1=0
1254            KTH2=0
1255     !!
1256            DO 450 LL=LLTOP,LLBOT
1257            L=LLBOT-LL+LLTOP
1258            BIT1=.FALSE.
1259            CR1=CLDAMT(L)
1260            BITX=(PINT(I,L).GE.PTOPC(NLVL+1)).AND.                           &
1261           &     (PINT(I,L).LT.PTOPC(NLVL)).AND.                             &
1262           &     (CLDAMT(L).GT.0.0)
1263            BIT1=BIT1.OR.BITX
1264            IF(.NOT.BIT1)GO TO 450
1265     !!***
1266     !!***  BITY=T: FIRST CLOUD LAYER; BITZ=T:CONSECUTIVE CLOUD LAYER
1267     !!***  NOTE:  WE ASSUME THAT THE THICKNESS OF EACH CLOUD LAYER IN THE
1268     !!***         DOMAIN IS LESS THAN 200 MB TO AVOID TOO MUCH COOLING OR
1269     !!***         HEATING. SO WE SET CTHK(NLVL)=200*E2. BUT THIS LIMIT MAY
1270     !!***         WORK WELL FOR CONVECTIVE CLOUDS. MODIFICATION MAY BE
1271     !!***         NEEDED IN THE FUTURE.
1272     !!***
1273            BITY=BITX.AND.(KTH2.LE.0)
1274            BITZ=BITX.AND.(KTH2.GT.0)
1275     !!
1276            IF(BITY)THEN
1277              KBT2=L
1278              KTH2=1
1279            ENDIF
1280     !!
1281            IF(BITZ)THEN
1282              KTOP1=KBT2-KTH2+1
1283              DPCL=PMID(I,KBT2)-PMID(I,KTOP1)
1284              IF(DPCL.LT.CTHK(NLVL))THEN
1285                KTH2=KTH2+1
1286              ELSE
1287                KBT2=KBT2-1
1288              ENDIF
1289            ENDIF
1290            IF(BITX)CL2=AMAX1(CL2,CR1)
1291     !!***
1292     !!*** AT THE DOMAIN BOUNDARY OR SEPARATED CLD LAYERS, RANDOM OVERLAP.
1293     !!*** CHOOSE THE THICKEST OR THE LARGEST FRACTION AMT AS THE CLD
1294     !!*** LAYER IN THAT DOMAIN.
1295     !!***
1296            BIT2=.FALSE.
1297            BITY=BITX.AND.(CLDAMT(L-1).LE.0.0.OR. &
1298                 PINT(I,L-1).LT.PTOPC(NLVL+1))
1299            BITZ=BITY.AND.CL1.GT.0.0
1300            BITW=BITY.AND.CL1.LE.0.0
1301            BIT2=BIT2.OR.BITY
1302            IF(.NOT.BIT2)GO TO 450
1303     !!
1304            IF(BITZ)THEN
1305              KBT1=INT((CL1*KBT1+CL2*KBT2)/(CL1+CL2))
1306              KTH1=INT((CL1*KTH1+CL2*KTH2)/(CL1+CL2))+1
1307              CL1=CL1+CL2-CL1*CL2
1308            ENDIF
1309     !!
1310            IF(BITW)THEN
1311              KBT1=KBT2
1312              KTH1=KTH2
1313              CL1=CL2
1314            ENDIF
1315     !!
1316            IF(BITY)THEN
1317              KBT2=LLBOT
1318              KTH2=0
1319              CL2=0.0
1320            ENDIF
1321        450 CONTINUE
1322     !
1323            CLDCFR(I,NLVL)=AMIN1(1.0,CL1)
1324            MTOP(I,NLVL)=MIN(KBT1,KBT1-KTH1+1)
1325            MBOT(I,NLVL)=KBT1
1326        480 CONTINUE
1327        500 CONTINUE
1328     
1329     !***
1330     !***  SET THE UN-NEEDED TAUDAR TO ONE
1331     !***
1332           DO I=MYIS,MYIE
1333             TAUDAR(I)=1.0
1334           ENDDO
1335     !----------------------------------------------------------------------
1336     ! NOW, CALCULATE THE CLOUD RADIATIVE PROPERTIES AFTER DAVIS (1982),
1337     ! HARSHVARDHAN ET AL (1987) AND Y.H., K.A.C. AND A.K. (1993).
1338     ! 
1339     ! UPDATE: THE FOLLOWING PARTS ARE MODIFIED, AFTER Y.T.H. (1994), TO 
1340     !         CALCULATE THE RADIATIVE PROPERTIES OF CLOUDS ON EACH MODEL
1341     !         LAYER. BOTH CONVECTIVE AND STRATIFORM CLOUDS ARE USED
1342     !         IN THIS CALCULATIONS.
1343     !
1344     !                                     QINGYUN ZHAO   95-3-22
1345     !
1346     !----------------------------------------------------------------------
1347     !
1348     !***
1349     !*** INITIALIZE ARRAYS FOR USES LATER
1350     !***
1351     
1352           DO 600 I=MYIS,MYIE
1353           LML=LMH(I,J)
1354           LVLIJ=LVL(I,J)
1355     !
1356     !***
1357     !*** NOTE: LAYER=1 IS THE SURFACE, AND LAYER=2 IS THE FIRST CLOUD
1358     !***       LAYER ABOVE THE SURFACE AND SO ON.
1359     !***
1360           EMIS(I,1)=1.0
1361           KTOP(I,1)=LP1
1362           KBTM(I,1)=LP1
1363           CAMT(I,1)=1.0
1364           TAUCLD(I,1)=0.0
1365           KCLD(I)=2
1366     !
1367           DO NBAND=1,NB
1368             RRCL(I,NBAND,1)=0.0
1369             TTCL(I,NBAND,1)=1.0
1370           ENDDO
1371     !
1372           DO 510 L=2,LP1
1373           CAMT(I,L)=0.0
1374           TAUCLD(I,L)=0.0
1375           KTOP(I,L)=1
1376           KBTM(I,L)=1
1377           EMIS(I,L)=0.0
1378     !
1379           DO NBAND=1,NB
1380             RRCL(I,NBAND,L)=0.0
1381             TTCL(I,NBAND,L)=1.0
1382           ENDDO
1383       510 CONTINUE
1384     
1385     !### End changes so far
1386     !***
1387     !*** NOW CALCULATE THE AMOUNT, TOP, BOTTOM AND TYPE OF EACH CLOUD LAYER
1388     !*** CLOUD TYPE=1: STRATIFORM CLOUD
1389     !***       TYPE=2: CONVECTIVE CLOUD
1390     !*** WHEN BOTH CONVECTIVE AND STRATIFORM CLOUDS EXIST AT THE SAME POINT,
1391     !*** SELECT CONVECTIVE CLOUD WITH THE HIGHER CLOUD FRACTION.
1392     !*** CLOUD LAYERS ARE SEPARATED BY TOTAL ABSENCE OF CLOUDINESS.
1393     !*** NOTE: THERE IS ONLY ONE CONVECTIVE CLOUD LAYER IN ONE COLUMN.
1394     !*** KTOP AND KBTM ARE THE TOP AND BOTTOM OF EACH CLOUD LAYER IN TERMS
1395     !*** OF MODEL LEVEL.
1396     !***
1397           NEW_CLOUD=.TRUE.
1398     !
1399           DO L=2,LML
1400             LL=LML-L+1+LVLIJ                        !-- Model layer
1401             CLFR=MAX(CCMID(I,LL),CSMID(I,LL))       !-- Cloud fraction in layer
1402             CLFR1=MAX(CCMID(I,LL+1),CSMID(I,LL+1))  !-- Cloud fraction in lower layer
1403     !-------------------
1404             IF (CLFR .GE. CLFRMIN) THEN
1405     !--- Cloud present at level
1406               IF (NEW_CLOUD) THEN
1407     !--- New cloud layer
1408                 IF(L==2.AND.CLFR1>=CLFRmin)THEN
1409                   KBTM(I,KCLD(I))=LL+1
1410                   CAMT(I,KCLD(I))=CLFR1
1411                 ELSE
1412                   KBTM(I,KCLD(I))=LL
1413                   CAMT(I,KCLD(I))=CLFR
1414                 ENDIF
1415                 NEW_CLOUD=.FALSE.
1416               ELSE
1417     !--- Existing cloud layer
1418                 CAMT(I,KCLD(I))=AMAX1(CAMT(I,KCLD(I)), CLFR)
1419               ENDIF        ! End IF (NEW_CLOUD .EQ. 0) ...
1420             ELSE IF (CLFR1 .GE. CLFRMIN) THEN
1421     !--- Cloud is not present at level but did exist at lower level, then ...
1422               IF (L .EQ. 2) THEN
1423     !--- For the case of ground fog
1424                 KBTM(I,KCLD(I))=LL+1
1425                 CAMT(I,KCLD(I))=CLFR1
1426               ENDIF
1427               KTOP(I,KCLD(I))=LL+1
1428               NEW_CLOUD=.TRUE.
1429               KCLD(I)=KCLD(I)+1
1430               CAMT(I,KCLD(I))=0.0
1431             ENDIF
1432     !-------------------
1433           ENDDO      !--- End DO L loop
1434     !***
1435     !*** THE REAL NUMBER OF CLOUD LAYERS IS (THE FIRST IS THE GROUND;
1436     !*** THE LAST IS THE SKY):
1437     !***
1438           NCLDS(I)=KCLD(I)-2
1439           NCLD=NCLDS(I)
1440     !***
1441     !***  NOW CALCULATE CLOUD RADIATIVE PROPERTIES
1442     !***
1443           IF(NCLD.GE.1)THEN
1444     !***
1445     !*** NOTE: THE FOLLOWING CALCULATIONS, THE UNIT FOR PRESSURE IS MB!!!
1446     !***
1447             DO 580 NC=2,NCLD+1
1448     !
1449             TauC=0.    !--- Total optical depth for each cloud layer (solar & longwave)
1450             QSUM=0.0
1451             NKTP=LP1
1452             NBTM=0
1453             BITX=CAMT(I,NC).GE.CLFRMIN
1454             NKTP=MIN(NKTP,KTOP(I,NC))
1455             NBTM=MAX(NBTM,KBTM(I,NC))
1456     !
1457             DO LL=NKTP,NBTM
1458               IF(LL.GE.KTOP(I,NC).AND.LL.LE.KBTM(I,NC).AND.BITX)THEN
1459                 PRS1=PINT(I,LL)*0.01
1460                 PRS2=PINT(I,LL+1)*0.01
1461                 DELP=PRS2-PRS1
1462                 TCLD=TMID(I,LL)-T0C
1463                 QSUM=QSUM+QMID(I,LL)*DELP*(PRS1+PRS2)                       &     
1464          &           /(120.1612*SQRT(TMID(I,LL)))
1465     !
1466     !***********************************************************************
1467     !****  IMPORTANT NOTES concerning input cloud optical properties  ******
1468     !***********************************************************************
1469     !
1470     !--- The simple optical depth parameterization from eq. (1) of Harshvardhan
1471     !    et al. (1989, JAS, p. 1924; hereafter referred to as HRCD by authorship)
1472     !    is used for convective cloud properties with some simple changes.
1473     !
1474     !--- The optical depth Tau is Tau=CTau*DELP, where values of CTau are
1475     !    described below.
1476     !      1) CTau=0.08*(Qc/Q0) for cloud water mixing ratio (Qc), where
1477     !         Q0 is assumed to be the threshold mixing ratio for "thick anvils",
1478     !         as noted in the 2nd paragraph after eq. (1) in Harshvardhan et al.
1479     !         (1989).  A value of Q0=0.1 g/kg is assumed based on experience w/
1480     !         cloud observations, and it is intended only to be a crude scaling
1481     !         factor for "order of magnitude" effects.  The functional dependence
1482     !         on mixing ratio is based on Stephens (1978, JAS, p. 2124, eq. 7).
1483     !         Result: CTau=800.*Qc => note that the "800." factor is referred to
1484     !         as an absorption coefficient
1485     !      2) For an assumed value of Q0=1 g/kg for "thick anvils", then 
1486     !         CTau=80.*Qc, or an absorption coefficient that is an order of 
1487     !         magnitude less.
1488     !      => ABSCOEF_W can vary from 100. to 1000. !!
1489     !      3) From p. 3105 of Dudhia (1989), values of 
1490     !         0.14 (m**2/g) * 1000 (g/kg) / 9.81 (m/s**2) = 14.27 /Pa
1491     !         => 14.27 (/Pa) * 100 (Pa/mb) = 1427 /mb
1492     !      4) From Dudhia's SW radiation, ABSCOEF_W ~ 1000.  after units conversion
1493     !      5) Again from p. 3105 of Dudhia (1989), he notes that ice absorption 
1494     !         coefficients are roughly half those of cloud water, it was decided
1495     !         to keep this simple and assume half that of water.
1496     !      => ABSCOEF_I=0.5*ABSCOEF_W
1497     !
1498     !--- For convection, the following is assumed:
1499     !      1) A characteristic water/ice mixing ratio (Qconv)
1500     !      2) A temperature threshold for water or ice (TRAD_ice)
1501     !
1502     !-----------------------------------------------------------------------
1503     !
1504                 CTau=0.
1505     !-- For crude estimation of convective cloud optical depths
1506                 IF (CCMID(I,LL) .GE. CLFRmin) THEN
1507                   IF (TCLD .GE. TRAD_ice) THEN
1508                     CTau=CTauCW            !--- Convective cloud water
1509                   ELSE
1510                     CTau=CTauCI            !--- Convective ice
1511                   ENDIF
1512     !              CTau=CTau*CCMID(I,LL)    !--- Reduce by convective cloud fraction
1513                 ENDIF
1514     !
1515     !-- For crude estimation of grid-scale cloud optical depths
1516     !
1517     !--   => The following 2 lines were intended to reduce cloud optical depths further 
1518     !        than what's parameterized in the NAM and what's theoretically justified
1519     !            CTau=CTau+CSMID(I,LL)*   &
1520     !     &           ( ABSCOEF_W*QWMID(I,LL)+ABSCOEF_I*QIMID(I,LL) )
1521                 CTau=CTau+ABSCOEF_W*QWMID(I,LL)+ABSCOEF_I*QIMID(I,LL)
1522                 TauC=TauC+DELP*CTau
1523                 TAUcld(I,NC)=TauC
1524               ENDIF      !--- End IF(LL.GE.KTOP(I,NC) ....
1525             ENDDO        !--- End DO LL
1526     !
1527     !!!!    IF(BITX)EMIS(I,NC)=1.0-EXP(-0.75*TauC)
1528             IF(BITX)EMIS(I,NC)=1.0-EXP(ABSCOEF_LW*TauC)
1529             IF(QSUM.GE.EPSQ1)THEN
1530     !
1531               DO 570 NBAND=1,NB
1532               IF(BITX)THEN
1533                 PROD=ABCFF(NBAND)*QSUM
1534                 DDX=TauC/(TauC+PROD)
1535                 EEX=1.0-DDX
1536                 IF(ABS(EEX).GE.1.E-8)THEN
1537                   DD=DDX
1538                   EE=EEX
1539                   FF=1.0-DD*0.85
1540                   AA=MIN(50.0,SQRT(3.0*EE*FF)*TauC)
1541                   AA=EXP(-AA)
1542                   BB=FF/EE
1543                   GG=SQRT(BB)
1544                   DD=(GG+1.0)*(GG+1.0)-(GG-1.0)*(GG-1.0)*AA*AA
1545                   RRCL(I,NBAND,NC)=MAX(0.1E-5,(BB-1.0)*(1.0-AA*AA)/DD)
1546                   TTCL(I,NBAND,NC)=AMAX1(0.1E-5,4.0*GG*AA/DD)
1547                 ENDIF
1548               ENDIF
1549       570     CONTINUE
1550             ENDIF
1551       580   CONTINUE
1552     !
1553           ENDIF
1554     !
1555       600 CONTINUE
1556     !*********************************************************************
1557     !******************  COMPUTE OZONE AT MIDLAYERS  *********************
1558     !*********************************************************************
1559     !
1560     !***  MODIFY PRESSURE AT THE TOP MODEL LAYER TO ACCOUNT FOR THE TOTAL
1561     !***  OZONE FROM MODEL TOP (PINT_1) TO THE TOP OF THE ATMOSPHERE (0 MB)
1562     !
1563           DO I=MYIS,MYIE
1564             FCTR=PINT(I,2)/(PINT(I,2)-PINT(I,1))
1565             POZN(I,1)=FCTR*(PMID(I,1)-PINT(I,1))
1566           ENDDO
1567     !
1568           CALL OZON2D(LM,POZN,XLAT,OZN,                                &
1569                       MYIS,MYIE,                                       &
1570                       ids,ide, jds,jde, kds,kde,                       &
1571                       ims,ime, jms,jme, kms,kme,                       &
1572                       its,ite, jts,jte, kts,kte                        )
1573     !
1574     !***  
1575     !***  NOW THE VARIABLES REQUIRED BY RADFS HAVE BEEN CALCULATED.
1576     !***
1577     !----------------------------------------------------------------------
1578     !***
1579     !***  CALL THE GFDL RADIATION DRIVER
1580     !***
1581     !***
1582           Jndx=J
1583           CALL RADFS &
1584          &     (PSFC,PMID,PINT,QMID,TMID,OZN,TSKN,SLMSK,ALBEDO,XLAT         &
1585     !BSF => for NAMX changes, pass in surface emissivity (SFCEMS) [different for snow]
1586          &,     TAUcld,CAMT,KTOP,KBTM,NCLDS,EMIS,RRCL,TTCL                  &
1587          &,     COSZ,TAUDAR,1                                               &
1588          &,     1,0                                                         &
1589          &,     ITIMSW,ITIMLW                                               &
1590          &,     TENDS(ITS:ITE,J,KTS:KTE),TENDL(ITS:ITE,J,KTS:KTE)           &
1591          &,     FLWUP,FSWUP,FSWDN,FSWDNS,FSWUPS,FLWDNS,FLWUPS,FSWDNSC       &
1592          &,     ids,ide, jds,jde, kds,kde                                   &
1593          &,     ims,ime, jms,jme, kms,kme                                   &
1594     ! begin debugging radiation
1595          &,     its,ite, jts,jte, kts,kte                                   &
1596          &,     imd,jmd, Jndx                                       )
1597     ! end debugging radiation
1598     !----------------------------------------------------------------------
1599           IF(LONG)THEN
1600     !
1601     !--  All fluxes in W/m**2
1602     !--- GLW    => downward longwave at the surface (formerly RLWIN) 
1603     !--- RLWTOA => outgoing longwave at the top of the atmosphere
1604     !-- Note:  RLWOUT & SIGT4 have been removed because they are no longer being used!
1605     !
1606             DO I=MYIS,MYIE
1607               GLW(I,J)=FLWDNS(I)
1608               RLWTOA(I,J)=FLWUP(I)
1609             ENDDO
1610           ENDIF
1611     !
1612           IF(SHORT)THEN
1613     !
1614     !--  All fluxes in W/m**2
1615     !--- GSW    => NET shortwave at the surface 
1616     !--- RSWIN  => incoming shortwave at the surface (all sky)
1617     !--- RSWINC => clear-sky incoming shortwave at the surface
1618     !--- RSWTOA => outgoing (reflected) shortwave at the top of the atmosphere 
1619     !
1620             DO I=MYIS,MYIE
1621               GSW(I,J)=FSWDNS(I)-FSWUPS(I)
1622               RSWIN(I,J) =FSWDNS(I)
1623               RSWINC(I,J)=FSWDNSC(I)
1624               RSWTOA(I,J)=FSWUP(I)
1625             ENDDO
1626           ENDIF
1627     !
1628     !***  ARRAYS ACFRST AND ACFRCV ACCUMULATE AVERAGE STRATIFORM AND
1629     !***  CONVECTIVE CLOUD FRACTIONS, RESPECTIVELY. 
1630     !***  ACCUMLATE THESE VARIABLES ONLY ONCE PER RADIATION CALL.
1631     !
1632     !***  ASSUME RANDOM OVERLAP BETWEEN LOW, MIDDLE, & HIGH LAYERS.
1633     !
1634     !***  UPDATE NEW 3D CLOUD FRACTION (CLDFRA)
1635     !
1636           DO I=MYIS,MYIE
1637             CFRACL(I,J)=CLDCFR(I,1)
1638             CFRACM(I,J)=CLDCFR(I,2)
1639             CFRACH(I,J)=CLDCFR(I,3)
1640             IF(CNCLD)THEN
1641               CFSmax=0.   !-- Maximum cloud fraction (stratiform component)
1642               CFCmax=0.   !-- Maximum cloud fraction (convective component)
1643               DO L=1,LMH(I,J)
1644                 LL=L+LVL(I,J)
1645                 CFSmax=MAX(CFSmax, CSMID(I,LL) )
1646                 CFCmax=MAX(CFCmax, CCMID(I,LL) )
1647               ENDDO
1648               ACFRST(I,J)=ACFRST(I,J)+CFSmax
1649               NCFRST(I,J)=NCFRST(I,J)+1
1650               ACFRCV(I,J)=ACFRCV(I,J)+CFCmax
1651               NCFRCV(I,J)=NCFRCV(I,J)+1
1652             ELSE
1653       !--- Count only locations with grid-scale cloudiness, ignore convective clouds
1654       !    (option not used, but if so set to the total cloud fraction)
1655               CFRAVG=1.-(1.-CFRACL(I,J))*(1.-CFRACM(I,J))*(1.-CFRACH(I,J))
1656               ACFRST(I,J)=ACFRST(I,J)+CFRAVG
1657               NCFRST(I,J)=NCFRST(I,J)+1
1658             ENDIF
1659     !--- Flip 3D cloud fractions in the vertical and save time
1660             LML=LMH(I,J)
1661             DO L=1,LML
1662               CLDFRA(I,J,L)=MAX(CCMID(I,L),CSMID(I,L))
1663             ENDDO
1664           ENDDO      !-- I index
1665     !***
1666     !***  THIS ROW IS FINISHED. GO TO NEXT
1667     !***
1668     !                        *********************
1669       700                          CONTINUE
1670     !                        *********************
1671     !----------------------------------------------------------------------
1672     !***
1673     !***  CALLS TO RADIATION THIS TIME STEP ARE COMPLETE.
1674     !***
1675     !----------------------------------------------------------------------
1676     ! begin debugging radiation
1677     !     FSWrat=0.
1678     !     if (RSWIN(imd,jmd) .gt. 0.)  &
1679     !        FSWrat=(RSWIN(imd,jmd)-GSW(imd,jmd))/RSWIN(imd,jmd)
1680     !     write(6,"(2a,2i5,7f9.2)") &
1681     !       '{rad3 imd,jmd,GSW,RSWIN,RSWOUT=RSWIN-GSW,RSWINC,GLW,' &
1682     !      ,'ALBEDO,RSWOUT/RSWIN = '&
1683     !      ,imd,jmd, GSW(imd,jmd),RSWIN(imd,jmd)  &
1684     !      ,RSWIN(imd,jmd)-GSW(imd,jmd),RSWINC(imd,jmd),GLW(imd,jmd) &
1685     !      ,ALB(imd,jmd),FSWrat
1686     ! end debugging radiation
1687     !----------------------------------------------------------------------
1688     !
1689     !--- Need to save LW & SW tendencies since radiation calculates both and this block
1690     
1691           END SUBROUTINE RADTN
1692     
1693     !----------------------------------------------------------------------
1694     
1695           REAL FUNCTION GAUSIN(xsd)
1696           REAL, PARAMETER :: crit=1.e-3
1697           REAL A1,A2,RN,B1,B2,B3,SUM,xsd
1698     !
1699     !  This function calculate area under the Gaussian curve between mean
1700     !  and xsd # of standard deviation (03/22/2004  Hsin-mu Lin)
1701     !
1702           a1=xsd*RSQR
1703           a2=exp(-0.5*xsd**2)
1704           rn=1.
1705           b1=1.
1706           b2=1.
1707           b3=1.
1708           sum=1.
1709           do while (b2 .gt. crit)
1710              rn=rn+1.
1711              b2=xsd**2/(2.*rn-1.)
1712              b3=b1*b2
1713              sum=sum+b3
1714              b1=b3
1715           enddo
1716           GAUSIN=a1*a2*sum
1717           RETURN
1718           END FUNCTION GAUSIN
1719     
1720     !----------------------------------------------------------------------
1721     
1722           SUBROUTINE ZENITH(TIMES,DAYI,HOUR,IDAT,IHRST,GLON,GLAT,CZEN,     &
1723                             MYIS,MYIE,MYJS,MYJE,                           &
1724                             IDS,IDE, JDS,JDE, KDS,KDE,                     &
1725                             IMS,IME, JMS,JME, KMS,KME,                     &
1726                             ITS,ITE, JTS,JTE, KTS,KTE)
1727     !----------------------------------------------------------------------
1728           IMPLICIT NONE
1729     !----------------------------------------------------------------------
1730           INTEGER, INTENT(IN)        :: IDS,IDE, JDS,JDE, KDS,KDE ,        &
1731                                         IMS,IME, JMS,JME, KMS,KME ,        &
1732                                         ITS,ITE, JTS,JTE, KTS,KTE
1733           INTEGER, INTENT(IN)        :: MYJS,MYJE,MYIS,MYIE
1734     
1735           REAL,    INTENT(IN)        :: TIMES
1736           REAL,    INTENT(OUT)       :: HOUR,DAYI
1737           INTEGER, INTENT(IN)        :: IHRST
1738     
1739           INTEGER, INTENT(IN), DIMENSION(3) :: IDAT 
1740           REAL,    INTENT(IN), DIMENSION(IMS:IME,JMS:JME) :: GLAT,GLON
1741           REAL,    INTENT(OUT), DIMENSION(IMS:IME,JMS:JME) :: CZEN
1742     
1743           REAL,    PARAMETER :: GSTC1=24110.54841,GSTC2=8640184.812866,    &
1744                                 GSTC3=9.3104E-2,GSTC4=-6.2E-6,             &
1745                                 PI=3.1415926,PI2=2.*PI,PIH=0.5*PI,         &
1746                                 DEG2RD=3.1415926/180.,OBLIQ=23.440*DEG2RD, &
1747                                 ZEROJD=2451545.0
1748     
1749           REAL    :: DAY,YFCTR,ADDDAY,STARTYR,DATJUL,DIFJD,SLONM,   &
1750                      ANOM,SLON,DEC,RA,DATJ0,TU,STIM0,SIDTIM,HRANG
1751           REAL    :: HRLCL,SINALT
1752           INTEGER :: KMNTH,KNT,IDIFYR,J,I
1753           LOGICAL :: LEAP
1754     !-----------------------------------------------------------------------
1755     !-----------------------------------------------------------------------
1756           INTEGER :: MONTH (12)
1757     !-----------------------------------------------------------------------
1758           DATA MONTH/31,28,31,30,31,30,31,31,30,31,30,31/
1759     !***********************************************************************
1760     !     SAVE MONTH
1761     
1762           DAY=0.
1763           LEAP=.FALSE.
1764           IF(MOD(IDAT(3),4).EQ.0)THEN
1765             MONTH(2)=29
1766             LEAP=.TRUE.
1767           ENDIF
1768           IF(IDAT(1).GT.1)THEN
1769             KMNTH=IDAT(1)-1
1770             DO 10 KNT=1,KMNTH
1771             DAY=DAY+REAL(MONTH(KNT))
1772        10   CONTINUE
1773           ENDIF
1774     !***
1775     !***  CALCULATE EXACT NUMBER OF DAYS FROM BEGINNING OF YEAR TO
1776     !***  FORECAST TIME OF INTEREST 
1777     !***
1778           DAY=DAY+REAL(IDAT(2)-1)+(REAL(IHRST)+TIMES/3600.)/24.
1779           DAYI=REAL(INT(DAY)+1)
1780           HOUR=(DAY-DAYI+1.)*24.
1781           YFCTR=2000.-IDAT(3)
1782     !-----------------------------------------------------------------------
1783     !***
1784     !***  FIND CELESTIAL LONGITUDE OF THE SUN THEN THE SOLAR DECLINATION AND
1785     !***  RIGHT ASCENSION.
1786     !***
1787     !-----------------------------------------------------------------------
1788           IDIFYR=IDAT(3)-2000
1789     !***
1790     !***  FIND JULIAN DATE OF START OF THE RELEVANT YEAR
1791     !***  ADDING IN LEAP DAYS AS NEEDED
1792     !***
1793           IF(IDIFYR.LT.0)THEN
1794             ADDDAY=REAL(IDIFYR/4)
1795           ELSE
1796             ADDDAY=REAL((IDIFYR+3)/4)
1797           ENDIF
1798           STARTYR=ZEROJD+IDIFYR*365.+ADDDAY-0.5
1799     !***
1800     !***  THE JULIAN DATE OF THE TIME IN QUESTION
1801     !***
1802           DATJUL=STARTYR+DAY
1803     !
1804     !***  DIFFERENCE OF ACTUAL JULIAN DATE FROM JULIAN DATE
1805     !***  AT 00H 1 January 2000
1806     !
1807           DIFJD=DATJUL-ZEROJD
1808     !
1809     !***  MEAN GEOMETRIC LONGITUDE OF THE SUN
1810     !
1811           SLONM=(280.460+0.9856474*DIFJD)*DEG2RD+YFCTR*PI2
1812     !
1813     !***  THE MEAN ANOMOLY
1814     !
1815           ANOM=(357.528+0.9856003*DIFJD)*DEG2RD
1816     !
1817     !***  APPARENT GEOMETRIC LONGITUDE OF THE SUN
1818     !
1819           SLON=SLONM+(1.915*SIN(ANOM)+0.020*SIN(2.*ANOM))*DEG2RD
1820           IF(SLON.GT.PI2)SLON=SLON-PI2
1821     !
1822     !***  DECLINATION AND RIGHT ASCENSION
1823     ! 
1824           DEC=ASIN(SIN(SLON)*SIN(OBLIQ))
1825           RA=ACOS(COS(SLON)/COS(DEC))
1826           IF(SLON.GT.PI)RA=PI2-RA
1827     !***
1828     !***  FIND THE GREENWICH SIDEREAL TIME THEN THE LOCAL SOLAR
1829     !***  HOUR ANGLE.
1830     !***
1831           DATJ0=STARTYR+DAYI-1.
1832           TU=(DATJ0-2451545.)/36525.
1833           STIM0=GSTC1+TU*(GSTC2+GSTC3*TU+GSTC4*TU*TU)
1834           SIDTIM=STIM0/3600.+YFCTR*24.+1.00273791*HOUR
1835           SIDTIM=SIDTIM*15.*DEG2RD
1836           IF(SIDTIM.LT.0.)SIDTIM=SIDTIM+PI2
1837           IF(SIDTIM.GT.PI2)SIDTIM=SIDTIM-PI2
1838           HRANG=SIDTIM-RA
1839     !
1840           DO 100 J=MYJS,MYJE
1841           DO 100 I=MYIS,MYIE
1842           HRLCL=HRANG+GLON(I,J)+PI2
1843     !***
1844     !***  THE ZENITH ANGLE IS THE COMPLEMENT OF THE ALTITUDE THUS THE
1845     !***  COSINE OF THE ZENITH ANGLE EQUALS THE SINE OF THE ALTITUDE.
1846     !***
1847           SINALT=SIN(DEC)*SIN(GLAT(I,J))+COS(DEC)*COS(HRLCL)* &
1848            COS(GLAT(I,J))
1849           IF(SINALT.LT.0.)SINALT=0.
1850           CZEN(I,J)=SINALT
1851       100 CONTINUE
1852     !***
1853     !***  IF THE FORECAST IS IN A DIFFERENT YEAR THAN THE START TIME,
1854     !***  RESET DAYI TO THE PROPER DAY OF THE NEW YEAR (IT MUST NOT BE
1855     !***  RESET BEFORE THE SOLAR ZENITH ANGLE IS COMPUTED).
1856     !***
1857           IF(DAYI.GT.365.)THEN
1858             IF(.NOT.LEAP)THEN
1859               DAYI=DAYI-365.
1860             ELSEIF(LEAP.AND.DAYI.GT.366.)THEN
1861               DAYI=DAYI-366.
1862             ENDIF
1863           ENDIF
1864     !
1865           END SUBROUTINE ZENITH
1866     !-----------------------------------------------------------------------
1867     
1868       SUBROUTINE OZON2D (LK,POZN,XLAT,QO3,                                &
1869                          MYIS,MYIE,                                       &
1870                          ids,ide, jds,jde, kds,kde,                       &
1871                          ims,ime, jms,jme, kms,kme,                       &
1872                          its,ite, jts,jte, kts,kte                        )
1873     !----------------------------------------------------------------------
1874      IMPLICIT NONE
1875     !----------------------------------------------------------------------
1876           INTEGER, INTENT(IN)        :: ids,ide, jds,jde, kds,kde ,      &
1877                                         ims,ime, jms,jme, kms,kme ,      &
1878                                         its,ite, jts,jte, kts,kte  
1879           INTEGER, INTENT(IN)        :: LK,MYIS,MYIE
1880           REAL,    INTENT(IN), DIMENSION(its:ite,kts:kte) :: POZN
1881           REAL,    INTENT(IN), DIMENSION(its:ite)  :: XLAT
1882           REAL,    INTENT(INOUT), DIMENSION(its:ite,kts:kte) :: QO3
1883     !----------------------------------------------------------------------
1884           INTEGER, PARAMETER ::  NL=81,NLP1=NL+1,LNGTH=37*NL
1885     !----------------------------------------------------------------------
1886     !----------------------------------------------------------------------
1887           INTEGER,DIMENSION(its:ite)    :: IARG,JJROW
1888           REAL,   DIMENSION(its:ite)    :: TTHAN
1889           REAL,   DIMENSION(its:ite,NL) :: QO3O3
1890     
1891           INTEGER :: I,K,NUMITR,ILOG,IT,NHALF
1892           REAL    :: TH2,DO3V,DO3VP,APHI,APLO
1893     !----------------------------------------------------------------------
1894     !
1895           DO I=ITS,ITE
1896             IARG(I)=0.
1897           ENDDO
1898     !
1899           DO I=MYIS,MYIE
1900             TH2=0.2*XLAT(I)
1901             JJROW(I)=19.001-TH2
1902             TTHAN(I)=(19-JJROW(I))-TH2
1903             IARG(I)=MIN(JJROW(I),36)
1904           ENDDO
1905     !
1906     !***  SEASONAL AND SPATIAL INTERPOLATION DONE BELOW.
1907     !
1908           DO K=1,NL
1909           DO I=MYIS,MYIE
1910             DO3V=XDUO3N(JJROW(I),K)+RSIN1*XDO3N2(JJROW(I),K)  &
1911                        +RCOS1*XDO3N3(JJROW(I),K)  &
1912                        +RCOS2*XDO3N4(JJROW(I),K)
1913             DO3VP=XDUO3N(IARG(I)+1,K)+RSIN1*XDO3N2(IARG(I)+1,K) &
1914                         +RCOS1*XDO3N3(IARG(I)+1,K) &
1915                         +RCOS2*XDO3N4(IARG(I)+1,K)
1916     !
1917     !***  NOW LATITUDINAL INTERPOLATION
1918     !***  AND CONVERT O3 INTO MASS MIXING RATIO (ORIG DATA MPY BY 1.E4)
1919     ! 
1920             QO3O3(I,K)=1.E-4*(DO3V+TTHAN(I)*(DO3VP-DO3V))
1921           ENDDO
1922           ENDDO
1923     !***
1924     !***  VERTICAL INTERPOLATION FOR EACH GRIDPOINT (LINEAR IN LN P)
1925     !***
1926           NUMITR=0
1927           ILOG=NL
1928        20 CONTINUE
1929           ILOG=(ILOG+1)/2
1930             IF(ILOG.EQ.1)GO TO 25
1931             NUMITR=NUMITR+1
1932             GO TO 20
1933        25 CONTINUE
1934     !
1935           DO 60 K=1,LK
1936     !
1937           NHALF=(NL+1)/2
1938           DO I=MYIS,MYIE
1939             JJROW(I)=NHALF
1940           ENDDO
1941     !
1942           DO 40 IT=1,NUMITR
1943           NHALF=(NHALF+1)/2
1944           DO I=MYIS,MYIE
1945             IF(POZN(I,K).LT.PRGFDL(JJROW(I)-1))THEN
1946               JJROW(I)=JJROW(I)-NHALF
1947             ELSEIF(POZN(I,K).GE.PRGFDL(JJROW(I)))THEN
1948               JJROW(I)=JJROW(I)+NHALF
1949             ENDIF
1950             JJROW(I)=MIN(JJROW(I),NL)
1951             JJROW(I)=MAX(JJROW(I),2)
1952           ENDDO
1953        40 CONTINUE
1954     !
1955           DO 50 I=MYIS,MYIE
1956           IF(POZN(I,K).LT.PRGFDL(1))THEN
1957             QO3(I,K)=QO3O3(I,1)
1958           ELSE IF(POZN(I,K).GT.PRGFDL(NL))THEN
1959             QO3(I,K)=QO3O3(I,NL)
1960           ELSE
1961             APLO=ALOG(PRGFDL(JJROW(I)-1))
1962             APHI=ALOG(PRGFDL(JJROW(I)))
1963             QO3(I,K)=QO3O3(I,JJROW(I))+(ALOG(POZN(I,K))-APHI)/ &
1964                        (APLO-APHI)* &
1965                        (QO3O3(I,JJROW(I)-1)-QO3O3(I,JJROW(I)))
1966           ENDIF
1967        50 CONTINUE
1968     !
1969        60 CONTINUE
1970     
1971       END SUBROUTINE OZON2D
1972     !-----------------------------------------------------------------------
1973     
1974           SUBROUTINE O3INT(PHALF,DDUO3N,DDO3N2,DDO3N3,DDO3N4, &
1975                      ids,ide, jds,jde, kds,kde,            &
1976                      ims,ime, jms,jme, kms,kme,            &
1977                      its,ite, jts,jte, kts,kte             )
1978     !----------------------------------------------------------------------
1979      IMPLICIT NONE
1980     !----------------------------------------------------------------------
1981           INTEGER, INTENT(IN)        :: ids,ide, jds,jde, kds,kde ,      &
1982                                         ims,ime, jms,jme, kms,kme ,      &
1983                                         its,ite, jts,jte, kts,kte
1984     
1985     !$$$  SUBPROGRAM DOCUMENTATION BLOCK
1986     !                .      .    .                                       .
1987     ! SUBPROGRAM:    O3INT       COMPUTE ZONAL MEAN OZONE FOR ETA LYRS
1988     !   PRGMMR: KENNETH CAMPANA  ORG: W/NMC23    DATE: 89-07-07
1989     !           MICHAEL BALDWIN  ORG: W/NMC22    DATE: 92-06-08
1990     !
1991     ! ABSTRACT: THIS CODE WRITTEN AT GFDL...
1992     !   CALCULATES SEASONAL ZONAL MEAN OZONE,EVERY 5 DEG OF LATITUDE,
1993     !   FOR CURRENT MODEL VERTICAL COORDINATE. OUTPUT DATA IN G/G * 1.E4
1994     !   CODE IS CALLED ONLY ONCE.
1995     !
1996     ! PROGRAM HISTORY LOG:
1997     !   84-01-01  FELS AND SCHWARZKOPF,GFDL.
1998     !   89-07-07  K. CAMPANA - ADAPTED STAND-ALONE CODE FOR IN-LINE USE.
1999     !   92-06-08  M. BALDWIN - UPDATE TO RUN IN ETA MODEL
2000     !
2001     ! USAGE:    CALL O3INT(O3,SIGL) OLD
2002     !   INPUT ARGUMENT LIST:
2003     !     PHALF    - MID LAYER PRESSURE (K=LM+1 IS MODEL SURFACE)
2004     !   OUTPUT ARGUMENT LIST:
2005     !     DDUO3N   - ZONAL MEAN OZONE DATA IN ALL MODEL LAYERS (G/G*1.E4)
2006     !     DDO3N2     DIMENSIONED(L,N),WHERE L(=37) IS LATITUDE BETWEEN
2007     !     DDO3N3     N AND S POLES,N=NUM OF VERTICAL LYRS(K=1 IS TOP LYR)
2008     !     DDO3N4     AND SEASON-WIN,SPR,SUM,FALL.
2009     !        IN COMMON
2010     !
2011     !   OUTPUT FILES:
2012     !     OUTPUT   - PRINT FILE.
2013     !
2014     ! ATTRIBUTES:
2015     !   LANGUAGE: FORTRAN 200.
2016     !
2017     !$$$
2018     !....     PROGRAM O3INT FROM DAN SCHWARZKOPF-GETS ZONAL MEAN O3
2019     !..    OUTPUT O3 IS WINTER,SPRING,SUMMER,FALL (NORTHERN HEMISPHERE)
2020     !-----------------------------------------------------------------------
2021     !-----------------------------------------------------------------------
2022     !     *********************************************************
2023           INTEGER :: N,NP,NP2,NM1
2024     !     *********************************************************
2025     !-----------------------------------------------------------------------
2026     !***
2027     !***  SEASONAL CLIMATOLOGIES OF O3 (OBTAINED FROM A PREVIOUSLY RUN
2028     !***  CODE WHICH INTERPOLATES O3 TO USER VERTICAL COORDINATE).
2029     !***  DEFINED AS 5 DEG LAT MEANS N.P.->S.P.
2030     !***
2031           REAL, INTENT(OUT), DIMENSION(37,kte):: DDUO3N,DDO3N2,DDO3N3,DDO3N4
2032     
2033     !     *********************************************************
2034           REAL,DIMENSION(33) :: P2
2035           REAL,DIMENSION(37) :: PH2
2036           REAL,DIMENSION(45) :: PH1
2037           REAL,DIMENSION(48) :: P1
2038           REAL,DIMENSION(81) :: P,RSTD,RDATA
2039           REAL,DIMENSION(82) :: QI,PH
2040           REAL,DIMENSION(19,kts:kte) :: DDUO3(19,kts:kte)
2041           REAL,DIMENSION(10,41) :: RO31,RO32
2042           REAL,DIMENSION(19,41) :: DUO3N
2043           REAL,DIMENSION(19) :: TEMPN
2044           REAL,DIMENSION(10,9) :: O3HI2
2045           REAL,DIMENSION(10,25) :: O3HI
2046           REAL,DIMENSION(10,16) :: O3LO1,O3LO2,O3LO3,O3LO4,O3HI1
2047           REAL,DIMENSION(10,40) :: RO3M
2048           REAL,DIMENSION(10,41) :: RO3
2049           REAL,DIMENSION(37,kts:kte) :: O35DEG
2050           REAL,DIMENSION(kts:kte) :: RBAR
2051           REAL,DIMENSION(kts:kte+1) :: PHALF
2052     
2053           INTEGER :: NKK,NK,NKP,K,L,NCASE,ITAPE,IPLACE,NKMM,NKM,KI,KK,KQ,JJ,KEN
2054           REAL :: O3RD,O3TOT,O3DU
2055     
2056           EQUIVALENCE (O3HI1(1,1),O3HI(1,1)),(O3HI2(1,1),O3HI(1,17))
2057           EQUIVALENCE (PH1(1),PH(1)),(PH2(1),PH(46))
2058           EQUIVALENCE (P1(1),P(1)),(P2(1),P(49))
2059           DATA PH1/      0., &
2060                0.1027246E-04, 0.1239831E-04, 0.1491845E-04, 0.1788053E-04, &
2061                0.2135032E-04, 0.2540162E-04, 0.3011718E-04, 0.3558949E-04, &
2062                0.4192172E-04, 0.4922875E-04, 0.5763817E-04, 0.6729146E-04, &
2063                0.7834518E-04, 0.9097232E-04, 0.1053635E-03, 0.1217288E-03, &
2064                0.1402989E-03, 0.1613270E-03, 0.1850904E-03, 0.2119495E-03, &
2065                0.2423836E-03, 0.2768980E-03, 0.3160017E-03, 0.3602623E-03, &
2066                0.4103126E-03, 0.4668569E-03, 0.5306792E-03, 0.6026516E-03, &
2067                0.6839018E-03, 0.7759249E-03, 0.8803303E-03, 0.9987843E-03, &
2068                0.1133178E-02, 0.1285955E-02, 0.1460360E-02, 0.1660001E-02, &
2069                0.1888764E-02, 0.2151165E-02, 0.2452466E-02, 0.2798806E-02, &
2070                0.3197345E-02, 0.3656456E-02, 0.4185934E-02, 0.4797257E-02/
2071           DATA PH2/ &
2072                0.5503893E-02, 0.6321654E-02, 0.7269144E-02, 0.8368272E-02, &
2073                0.9644873E-02, 0.1112946E-01, 0.1285810E-01, 0.1487354E-01, &
2074                0.1722643E-01, 0.1997696E-01, 0.2319670E-01, 0.2697093E-01, &
2075                0.3140135E-01, 0.3660952E-01, 0.4274090E-01, 0.4996992E-01, &
2076                0.5848471E-01, 0.6847525E-01, 0.8017242E-01, 0.9386772E-01, &
2077                0.1099026E+00, 0.1286765E+00, 0.1506574E+00, 0.1763932E+00, &
2078                0.2065253E+00, 0.2415209E+00, 0.2814823E+00, 0.3266369E+00, &
2079                0.3774861E+00, 0.4345638E+00, 0.4984375E+00, 0.5697097E+00, &
2080                0.6490189E+00, 0.7370409E+00, 0.8344896E+00, 0.9421190E+00, &
2081                0.1000000E+01/
2082           DATA P1/ &
2083                0.9300000E-05, 0.1129521E-04, 0.1360915E-04, 0.1635370E-04, &
2084                0.1954990E-04, 0.2331653E-04, 0.2767314E-04, 0.3277707E-04, &
2085                0.3864321E-04, 0.4547839E-04, 0.5328839E-04, 0.6234301E-04, &
2086                0.7263268E-04, 0.8450696E-04, 0.9793231E-04, 0.1133587E-03, &
2087                0.1307170E-03, 0.1505832E-03, 0.1728373E-03, 0.1982122E-03, &
2088                0.2266389E-03, 0.2592220E-03, 0.2957792E-03, 0.3376068E-03, &
2089                0.3844381E-03, 0.4379281E-03, 0.4976965E-03, 0.5658476E-03, &
2090                0.6418494E-03, 0.7287094E-03, 0.8261995E-03, 0.9380076E-03, &
2091                0.1063498E-02, 0.1207423E-02, 0.1369594E-02, 0.1557141E-02, &
2092                0.1769657E-02, 0.2015887E-02, 0.2295520E-02, 0.2620143E-02, &
2093                0.2989651E-02, 0.3419469E-02, 0.3909867E-02, 0.4481491E-02, &
2094                0.5135272E-02, 0.5898971E-02, 0.6774619E-02, 0.7799763E-02/
2095           DATA P2/ &
2096                0.8978218E-02, 0.1036103E-01, 0.1195488E-01, 0.1382957E-01, &
2097                0.1599631E-01, 0.1855114E-01, 0.2151235E-01, 0.2501293E-01, &
2098                0.2908220E-01, 0.3390544E-01, 0.3952926E-01, 0.4621349E-01, &
2099                0.5403168E-01, 0.6330472E-01, 0.7406807E-01, 0.8677983E-01, &
2100                0.1015345E+00, 0.1189603E+00, 0.1391863E+00, 0.1630739E+00, &
2101                0.1908004E+00, 0.2235461E+00, 0.2609410E+00, 0.3036404E+00, &
2102                0.3513750E+00, 0.4055375E+00, 0.4656677E+00, 0.5335132E+00, &
2103                0.6083618E+00, 0.6923932E+00, 0.7845676E+00, 0.8875882E+00, &
2104                0.1000000E+01/
2105           DATA O3HI1/ &
2106            .55,.50,.45,.45,.40,.35,.35,.30,.30,.30, &
2107            .55,.51,.46,.47,.42,.38,.37,.36,.35,.35, &
2108            .55,.53,.48,.49,.44,.42,.41,.40,.38,.38, &
2109            .60,.55,.52,.52,.50,.47,.46,.44,.42,.41, &
2110            .65,.60,.55,.56,.53,.52,.50,.48,.45,.45, &
2111            .75,.65,.60,.60,.55,.55,.55,.50,.48,.47, &
2112            .80,.75,.75,.75,.70,.70,.65,.63,.60,.60, &
2113            .90,.85,.85,.80,.80,.75,.75,.74,.72,.71, &
2114            1.10,1.05,1.00,.90,.90,.90,.85,.83,.80,.80, &
2115            1.40,1.30,1.25,1.25,1.25,1.20,1.15,1.10,1.05,1.00, &
2116            1.7,1.7,1.6,1.6,1.6,1.6,1.6,1.6,1.5,1.5, &
2117            2.1,2.0,1.9,1.9,1.9,1.8,1.8,1.8,1.7,1.7, &
2118            2.4,2.3,2.2,2.2,2.2,2.1,2.1,2.1,2.0,2.0, &
2119            2.7,2.5,2.5,2.5,2.5,2.5,2.4,2.4,2.3,2.3, &
2120            2.9,2.8,2.7,2.7,2.7,2.7,2.7,2.7,2.6,2.6, &
2121            3.1,3.1,3.0,3.0,3.0,3.0,3.0,3.0,2.9,2.8/
2122           DATA O3HI2/ &
2123            3.3,3.4,3.4,3.6,3.7,3.9,4.0,4.1,4.0,3.8, &
2124            3.6,3.8,3.9,4.2,4.7,5.3,5.6,5.7,5.5,5.2, &
2125            4.1,4.3,4.7,5.2,6.0,6.7,7.0,6.8,6.4,6.2, &
2126            5.4,5.7,6.0,6.6,7.3,8.0,8.4,7.7,7.1,6.7, &
2127            6.7,6.8,7.0,7.6,8.3,10.0,9.6,8.2,7.5,7.2, &
2128            9.2,9.3,9.4,9.6,10.3,10.6,10.0,8.5,7.7,7.3, &
2129            12.6,12.1,12.0,12.1,11.7,11.0,10.0,8.6,7.8,7.4, &
2130            14.2,13.5,13.1,12.8,11.9,10.9,9.8,8.5,7.8,7.5, &
2131            14.3,14.0,13.4,12.7,11.6,10.6,9.3,8.4,7.6,7.3/
2132           DATA O3LO1/ &
2133            14.9,14.2,13.3,12.5,11.2,10.3,9.5,8.6,7.5,7.4, &
2134            14.5,14.1,13.0,11.8,10.5,9.8,9.2,7.9,7.4,7.4, &
2135            11.8,11.5,10.9,10.5,9.9,9.6,8.9,7.5,7.2,7.2, &
2136            7.3,7.7,7.8,8.4,8.4,8.5,7.9,7.4,7.1,7.1, &
2137            4.1,4.4,5.3,6.6,6.9,7.5,7.4,7.2,7.0,6.9, &
2138            1.8,1.9,2.5,3.3,4.5,5.8,6.3,6.3,6.4,6.1, &
2139            0.4,0.5,0.8,1.2,2.7,3.6,4.6,4.7,5.0,5.2, &
2140            .10,.15,.20,.50,1.4,2.1,3.0,3.2,3.5,3.9, &
2141            .07,.10,.12,.30,1.0,1.4,1.8,1.9,2.3,2.5, &
2142            .06,.08,.10,.15,.60,.80,1.4,1.5,1.5,1.6, &
2143            .05,.05,.06,.09,.20,.40,.70,.80,.90,.90, &
2144            .05,.05,.06,.08,.10,.13,.20,.25,.30,.40, &
2145            .05,.05,.05,.06,.07,.07,.08,.09,.10,.13, &
2146            .05,.05,.05,.05,.06,.06,.06,.06,.07,.07, &
2147            .05,.05,.05,.05,.05,.05,.05,.06,.06,.06, &
2148            .04,.04,.04,.04,.04,.04,.04,.05,.05,.05/
2149           DATA O3LO2/ &
2150            14.8,14.2,13.8,12.2,11.0,9.8,8.5,7.8,7.4,6.9, &
2151            13.2,13.0,12.5,11.3,10.4,9.0,7.8,7.5,7.0,6.6, &
2152            10.6,10.6,10.7,10.1,9.4,8.6,7.5,7.0,6.5,6.1, &
2153            7.0,7.3,7.5,7.5,7.5,7.3,6.7,6.4,6.0,5.8, &
2154            3.8,4.0,4.7,5.0,5.2,5.9,5.8,5.6,5.5,5.5, &
2155            1.4,1.6,2.4,3.0,3.7,4.1,4.6,4.8,5.1,5.0, &
2156            .40,.50,.90,1.2,2.0,2.7,3.2,3.6,4.3,4.1, &
2157            .07,.10,.20,.30,.80,1.4,2.1,2.4,2.7,3.0, &
2158            .06,.07,.09,.15,.30,.70,1.2,1.4,1.6,2.0, &
2159            .05,.05,.06,.12,.15,.30,.60,.70,.80,.80, &
2160            .04,.05,.06,.08,.09,.15,.30,.40,.40,.40, &
2161            .04,.04,.05,.055,.06,.09,.12,.13,.15,.15, &
2162            .03,.03,.045,.052,.055,.06,.07,.07,.06,.07, &
2163            .03,.03,.04,.051,.052,.052,.06,.06,.05,.05, &
2164            .02,.02,.03,.05,.05,.05,.04,.04,.04,.04, &
2165            .02,.02,.02,.04,.04,.04,.03,.03,.03,.03/
2166           DATA O3LO3/ &
2167            14.5,14.0,13.5,11.3,11.0,10.0,9.0,8.3,7.5,7.3, &
2168            13.5,13.2,12.5,11.1,10.4,9.7,8.2,7.8,7.4,6.8, &
2169            10.8,10.9,11.0,10.4,10.0,9.6,7.9,7.5,7.0,6.7, &
2170            7.3,7.5,7.8,8.5,9.0,8.5,7.7,7.4,6.9,6.5, &
2171            4.1,4.5,5.3,6.2,7.3,7.7,7.3,7.0,6.6,6.4, &
2172            1.8,2.0,2.2,3.8,4.3,5.6,6.2,6.2,6.4,6.2, &
2173            .30,.50,.60,1.5,2.8,3.7,4.5,4.7,5.5,5.6, &
2174            .09,.10,.15,.60,1.2,2.1,3.0,3.5,4.0,4.3, &
2175            .06,.08,.10,.30,.60,1.1,1.9,2.2,2.9,3.0, &
2176            .04,.05,.06,.15,.45,.60,1.1,1.3,1.6,1.8, &
2177            .04,.04,.04,.08,.20,.30,.55,.60,.75,.90, &
2178            .04,.04,.04,.05,.06,.10,.12,.15,.20,.25, &
2179            .04,.04,.03,.04,.05,.06,.07,.07,.07,.08, &
2180            .03,.03,.04,.05,.05,.05,.05,.05,.05,.05, &
2181            .03,.03,.03,.04,.04,.04,.05,.05,.04,.04, &
2182            .02,.02,.02,.04,.04,.04,.04,.04,.03,.03/
2183           DATA O3LO4/ &
2184            14.2,13.8,13.2,12.5,11.7,10.5,8.6,7.8,7.5,6.6, &
2185            12.5,12.4,12.2,11.7,10.8,9.8,7.8,7.2,6.5,6.1, &
2186            10.6,10.5,10.4,10.1,9.6,9.0,7.1,6.8,6.1,5.9, &
2187            7.0,7.4,7.9,7.8,7.6,7.3,6.2,6.1,5.8,5.6, &
2188            4.2,4.6,5.1,5.6,5.9,5.9,5.9,5.8,5.6,5.3, &
2189            2.1,2.3,2.6,2.9,3.5,4.3,4.8,4.9,5.1,5.1, &
2190            0.7,0.8,1.0,1.5,2.0,2.8,3.5,3.6,3.7,4.0, &
2191            .15,.20,.40,.50,.60,1.4,2.1,2.2,2.3,2.5, &
2192            .08,.10,.15,.25,.30,.90,1.2,1.3,1.4,1.6, &
2193            .07,.08,.10,.14,.20,.50,.70,.90,.90,.80, &
2194            .05,.06,.08,.12,.14,.20,.35,.40,.60,.50, &
2195            .05,.05,.08,.09,.09,.09,.11,.12,.15,.18, &
2196            .04,.05,.06,.07,.07,.08,.08,.08,.08,.08, &
2197            .04,.04,.05,.07,.07,.07,.07,.07,.06,.05, &
2198            .02,.02,.04,.05,.05,.05,.05,.05,.04,.04, &
2199            .02,.02,.03,.04,.04,.04,.04,.04,.03,.03/
2200     
2201           N=kte;NP=N+1;NP2=N+2;NM1=N-1
2202     
2203           NKK=41
2204           NK=81
2205           NKP=NK+1
2206           DO 24 K=1,NP
2207        24 PHALF(K)=PHALF(K)*0.01*1.0E+03
2208           DO 25 K=1,NK
2209           PH(K)=PH(K)*1013250.
2210        25 P(K)=P(K)*1013250.
2211           PH(NKP)=PH(NKP)*1013250.
2212     !***LOAD ARRAYS RO31,RO32,AS IN DICKS PGM.
2213           DO 1010 K=1,25
2214           DO 1010 L=1,10
2215             RO31(L,K)=O3HI(L,K)
2216             RO32(L,K)=O3HI(L,K)
2217     1010  CONTINUE
2218     !
2219           DO 3000 NCASE=1,4
2220           ITAPE=NCASE+50
2221           IPLACE=2
2222           IF (NCASE.EQ.2) IPLACE=4
2223           IF (NCASE.EQ.3) IPLACE=1
2224           IF (NCASE.EQ.4) IPLACE=3
2225     !***NCASE=1: SPRING (IN N.H.)
2226     !***NCASE=2: FALL   (IN N.H.)
2227     !***NCASE=3: WINTER (IN N.H.)
2228     !***NCASE=4: SUMMER (IN N.H.)
2229           IF (NCASE.EQ.1.OR.NCASE.EQ.2) THEN
2230              DO 1011 K=26,41
2231              DO 1011 L=1,10
2232                RO31(L,K)=O3LO1(L,K-25)
2233                RO32(L,K)=O3LO2(L,K-25)
2234     1011     CONTINUE
2235           ENDIF
2236           IF (NCASE.EQ.3.OR.NCASE.EQ.4) THEN
2237              DO 1031 K=26,41
2238              DO 1031 L=1,10
2239                RO31(L,K)=O3LO3(L,K-25)
2240                RO32(L,K)=O3LO4(L,K-25)
2241     1031     CONTINUE
2242           ENDIF
2243           DO 30 KK=1,NKK
2244           DO 31 L=1,10
2245           DUO3N(L,KK)=RO31(11-L,KK)
2246        31 DUO3N(L+9,KK)=RO32(L,KK)
2247           DUO3N(10,KK)=.5*(RO31(1,KK)+RO32(1,KK))
2248        30 CONTINUE
2249     !***FOR NCASE=2 OR NCASE=4,REVERSE LATITUDE ARRANGEMENT OF CORR. SEASON
2250           IF (NCASE.EQ.2.OR.NCASE.EQ.4) THEN
2251              DO 1024 KK=1,NKK
2252              DO 1025 L=1,19
2253                TEMPN(L)=DUO3N(20-L,KK)
2254     1025     CONTINUE
2255              DO 1026 L=1,19
2256                DUO3N(L,KK)=TEMPN(L)
2257     1026     CONTINUE
2258     1024     CONTINUE
2259           ENDIF
2260     !***DUO3N NOW IS O3 PROFILE FOR APPROPRIATE SEASON,AT STD. PRESSURE
2261     !      LEVELS
2262     !KAC  WRITE (6,800) DUO3N
2263     !***BEGIN LATITUDE (10 DEG) LOOP
2264           DO 33 L=1,19
2265           DO 22 KK=1,NKK
2266        22 RSTD(KK)=DUO3N(L,KK)
2267           NKM=NK-1
2268           NKMM=NK-3
2269     !     BESSELS HALF-POINT INTERPOLATION FORMULA
2270           DO 60 K=4,NKMM,2
2271           KI=K/2
2272        60 RDATA(K)=.5*(RSTD(KI)+RSTD(KI+1))-(RSTD(KI+2)-RSTD(KI+1)-RSTD(KI)+ &
2273           RSTD(KI-1))/16.
2274           RDATA(2)=.5*(RSTD(2)+RSTD(1))
2275           RDATA(NKM)=.5*(RSTD(NKK)+RSTD(NKK-1))
2276     !     PUT UNCHANGED DATA INTO NEW ARRAY
2277           DO 61 K=1,NK,2
2278           KQ=(K+1)/2
2279        61 RDATA(K)=RSTD(KQ)
2280     !---NOTE TO NMC: THIS WRITE IS COMMENTED OUT TO REDUCE PRINTOUT
2281     !     WRITE (6,798) RDATA
2282     !     CALCULATE LAYER-MEAN OZONE MIXING RATIO FOR EACH MODEL LEVEL
2283           DO 99 KK=1,N
2284           RBAR(KK)=0.
2285     !     LOOP TO CALCULATE SUMS TO GET LAYER OZONE MEAN
2286           DO 98 K=1,NK
2287           IF(PH(K+1).LT.PHALF(KK)) GO TO 98
2288           IF(PH(K).GT.PHALF(KK+1)) GO TO 98
2289           IF(PH(K+1).LT.PHALF(KK+1).AND.PH(K).LT.PHALF(KK)) RBAR(KK)=RBAR(KK &
2290           )+RDATA(K)*(PH(K+1)-PHALF(KK))
2291           IF(PH(K+1).LT.PHALF(KK+1).AND.PH(K).GE.PHALF(KK)) RBAR(KK)=RBAR(KK &
2292           )+RDATA(K)*(PH(K+1)-PH(K))
2293           IF(PH(K+1).GT.PHALF(KK+1).AND.PH(K).GT.PHALF(KK)) RBAR(KK)=RBAR(KK &
2294           )+RDATA(K)*(PHALF(KK+1)-PH(K))
2295        98 CONTINUE
2296           RBAR(KK)=RBAR(KK)/(PHALF(KK+1)-PHALF(KK))
2297           IF(RBAR(KK).GT..0000) GO TO 99
2298     !     CODE TO COVER CASE WHEN MODEL RESOLUTION IS SO FINE THAT NO VALUE
2299     !     OF P(K) IN THE OZONE DATA ARRAY FALLS BETWEEN PHALF(KK+1) AND
2300     !     PHALF(KK).   PROCEDURE IS TO SIMPLY GRAB THE NEAREST VALUE FROM
2301     !     RDATA
2302           DO 29 K=1,NK
2303           IF(PH(K).LT.PHALF(KK).AND.PH(K+1).GE.PHALF(KK+1)) RBAR(KK)=RDATA(K)
2304        29 CONTINUE
2305        99 CONTINUE
2306     !     CALCULATE TOTAL OZONE
2307           O3RD=0.
2308           DO 89 KK=1,80
2309        89 O3RD=O3RD+RDATA(KK)*(PH(KK+1)-PH(KK))
2310           O3RD=O3RD+RDATA(81)*(P(81)-PH(81))
2311           O3RD=O3RD/980.
2312           O3TOT=0.
2313           DO 88 KK=1,N
2314        88 O3TOT=O3TOT+RBAR(KK)*(PHALF(KK+1)-PHALF(KK))
2315           O3TOT=O3TOT/980.
2316     !     UNITS ARE MICROGRAMS/CM**2
2317           O3DU=O3TOT/2.144
2318     !     O3DU UNITS ARE DOBSON UNITS (10**-3 ATM-CM)
2319     !--NOTE TO NMC: THIS IS COMMENTED OUT TO SAVE PRINTOUT
2320     !     WRITE (6,796) O3RD,O3TOT,O3DU
2321           DO 23 KK=1,N
2322        23 DDUO3(L,KK)=RBAR(KK)*.01
2323        33 CONTINUE
2324     !***END OF LATITUDE LOOP
2325     !
2326     !***CREATE 5 DEG OZONE QUANTITIES BY LINEAR INTERPOLATION OF
2327     !      10 DEG VALUES
2328           DO 1060 KK=1,N
2329             DO 1061 L=1,19
2330               O35DEG(2*L-1,KK)=DDUO3(L,KK)
2331     1061    CONTINUE
2332             DO 1062 L=1,18
2333               O35DEG(2*L,KK)=0.5*(DDUO3(L,KK)+DDUO3(L+1,KK))
2334     1062    CONTINUE
2335     1060  CONTINUE
2336     !***OUTPUT TO UNIT (ITAPE) THE OZONE VALUES FOR LATER USE
2337     !O222  ***************************************************
2338     !C          WRITE (66) O35DEG
2339           IF (IPLACE.EQ.1) THEN
2340           DO 302 JJ=1,37
2341            DO 302 KEN=1,N
2342             DDUO3N(JJ,KEN) = O35DEG(JJ,KEN)
2343       302 CONTINUE
2344           ELSE IF (IPLACE.EQ.2) THEN
2345           DO 312 JJ=1,37
2346            DO 312 KEN=1,N
2347             DDO3N2(JJ,KEN) = O35DEG(JJ,KEN)
2348       312 CONTINUE
2349           ELSE IF (IPLACE.EQ.3) THEN
2350           DO 322 JJ=1,37
2351            DO 322 KEN=1,N
2352             DDO3N3(JJ,KEN) = O35DEG(JJ,KEN)
2353       322 CONTINUE
2354           ELSE IF (IPLACE.EQ.4) THEN
2355           DO 332 JJ=1,37
2356            DO 332 KEN=1,N
2357             DDO3N4(JJ,KEN) = O35DEG(JJ,KEN)
2358       332 CONTINUE
2359           END IF
2360     !O222  ***************************************************
2361     3000  CONTINUE
2362     !***END OF LOOP OVER CASES
2363           RETURN
2364        1  FORMAT(10F4.2)
2365         2 FORMAT(10X,E14.7,1X,E14.7,1X,E14.7,1X,E14.7,1X)
2366        3  FORMAT(10E12.5)
2367       797 FORMAT(10F7.2)
2368       799 FORMAT(19F6.4)
2369       800 FORMAT(19F6.2)
2370       102 FORMAT(' O3 IPLACE=',I4)
2371      1033 FORMAT(19F6.5)
2372       101 FORMAT(5X,1H*,F6.5,1H,,F6.5,1H,,F6.5,1H,,F6.5,1H,,F6.5,1H,,F6.5, &
2373           1H,,F6.5,1H,,F6.5,1H,,F6.5,1H,)
2374           
2375           END SUBROUTINE O3INT
2376     !----------------------------------------------------------------
2377     
2378       SUBROUTINE CLO89(CLDFAC,CAMT,NCLDS,KBTM,KTOP                  &
2379           ,          ids,ide, jds,jde, kds,kde                      &
2380           ,          ims,ime, jms,jme, kms,kme                      &
2381           ,          its,ite, jts,jte, kts,kte                      )
2382     !----------------------------------------------------------------------
2383      IMPLICIT NONE
2384     !----------------------------------------------------------------------
2385           INTEGER, INTENT(IN)        :: ids,ide, jds,jde, kds,kde ,      &
2386                                         ims,ime, jms,jme, kms,kme ,      &
2387                                         its,ite, jts,jte, kts,kte
2388     !----------------------------------------------------------------------
2389     
2390     !     ************************************************************
2391     !     *                                                          *
2392     !     * THIS SUBROUTINE WAS MODIFIED TO BE USED IN THE ETA MODEL *
2393     !     *                                                          *
2394     !     *                            Q. ZHAO    95-3-22            *
2395     !     *                                                          *
2396     !     ************************************************************
2397     
2398           REAL,    INTENT(OUT),DIMENSION(its:ite,kts:kte+1,kts:kte+1) :: CLDFAC
2399           REAL,    INTENT(IN), DIMENSION(its:ite,kts:kte+1) :: CAMT
2400           INTEGER, INTENT(IN), DIMENSION(its:ite,kts:kte+1) :: KBTM,KTOP
2401           INTEGER, INTENT(IN), DIMENSION(its:ite)           :: NCLDS
2402     
2403           REAL,    DIMENSION(kts:kte+1,kts:kte+1,64) :: CLDIPT
2404           REAL,    DIMENSION(kts:kte+1) :: CLDROW
2405           INTEGER:: IQ,ITOP,I,J,JTOP,IR,IP,K1,K2,KB,K,KP,KT,NC
2406           REAL   :: XCLD
2407     
2408           INTEGER :: L,LP1,LP2,LP3,LM1,LM2,LM3,MYIS,MYIE
2409     
2410         !  DIMENSION CLDIPT(LP1,LP1, 64 )
2411         !  DIMENSION NCLDS(IDIM1:IDIM2),KTOP(IDIM1:IDIM2,LP1), &
2412         !            KBTM(IDIM1:IDIM2,LP1)
2413         !  DIMENSION CLDROW(LP1)
2414         !  DIMENSION CAMT(IDIM1:IDIM2,LP1),CLDFAC(IDIM1:IDIM2,LP1,LP1)
2415     
2416           L=kte
2417           LP1=L+1;  LP2=L+2;  LP3=L+3
2418           LM1=L-1;  LM2=L-2;  LM3=L-3
2419           MYIS=its; MYIE=ite
2420     
2421     !
2422           DO 1 IQ=MYIS,MYIE,64
2423           ITOP=IQ+63
2424           IF(ITOP.GT.MYIE) ITOP=MYIE
2425           JTOP=ITOP-IQ+1
2426           DO 11 IP=1,JTOP
2427           IR=IQ+IP-1
2428           IF (NCLDS(IR).EQ.0) THEN
2429             DO 25 J=1,LP1
2430             DO 25 I=1,LP1
2431             CLDIPT(I,J,IP)=1.
2432     25      CONTINUE
2433           ENDIF
2434           IF (NCLDS(IR).GE.1) THEN
2435               XCLD=1.-CAMT(IR,2)
2436                K1=KTOP(IR,2)+1
2437                K2=KBTM(IR,2)
2438               DO 27 J=1,LP1
2439                   CLDROW(J)=1.
2440     27        CONTINUE
2441               DO 29 J=1,K2
2442                   CLDROW(J)=XCLD
2443     29        CONTINUE
2444               KB=MAX(K1,K2+1)
2445               DO 33 K=KB,LP1
2446               DO 33 KP=1,LP1
2447                    CLDIPT(KP,K,IP)=CLDROW(KP)
2448     33        CONTINUE
2449               DO 37 J=1,LP1
2450                   CLDROW(J)=1.
2451     37        CONTINUE
2452               DO 39 J=K1,LP1
2453                   CLDROW(J)=XCLD
2454     39        CONTINUE
2455               KT=MIN(K1-1,K2)
2456               DO 43 K=1,KT
2457               DO 43 KP=1,LP1
2458                   CLDIPT(KP,K,IP)=CLDROW(KP)
2459     43        CONTINUE
2460               IF(K2+1.LE.K1-1) THEN
2461                 DO 31 J=K2+1,K1-1
2462                 DO 31 I=1,LP1
2463                     CLDIPT(I,J,IP)=1.
2464     31          CONTINUE
2465               ELSE IF(K1.LE.K2) THEN
2466                 DO 32 J=K1,K2
2467                 DO 32 I=1,LP1
2468                     CLDIPT(I,J,IP)=XCLD
2469     32          CONTINUE
2470               ENDIF
2471           ENDIF
2472     
2473           IF (NCLDS(IR).GE.2) THEN
2474             DO 21 NC=2,NCLDS(IR)
2475               XCLD=1.-CAMT(IR,NC+1)
2476                K1=KTOP(IR,NC+1)+1
2477                K2=KBTM(IR,NC+1)
2478               DO 47 J=1,LP1
2479                   CLDROW(J)=1.
2480     47        CONTINUE
2481               DO 49 J=1,K2
2482                   CLDROW(J)=XCLD
2483     49        CONTINUE
2484               KB=MAX(K1,K2+1)
2485               DO 53 K=KB,LP1
2486               DO 53 KP=1,LP1
2487                    CLDIPT(KP,K,IP)=CLDIPT(KP,K,IP)*CLDROW(KP)
2488     53        CONTINUE
2489               DO 57 J=1,LP1
2490                   CLDROW(J)=1.
2491     57        CONTINUE
2492               DO 59 J=K1,LP1
2493                   CLDROW(J)=XCLD
2494     59        CONTINUE
2495               KT=MIN(K1-1,K2)
2496               DO 63 K=1,KT
2497               DO 63 KP=1,LP1
2498                   CLDIPT(KP,K,IP)=CLDIPT(KP,K,IP)*CLDROW(KP)
2499     63        CONTINUE
2500               IF(K1.LE.K2) THEN
2501                 DO 52 J=K1,K2
2502                 DO 52 I=1,LP1
2503                     CLDIPT(I,J,IP)=CLDIPT(I,J,IP)*XCLD
2504     52          CONTINUE
2505               ENDIF
2506     21        CONTINUE
2507           ENDIF
2508     11    CONTINUE
2509           DO 71 J=1,LP1
2510           DO 71 I=1,LP1
2511           DO 71 IP=1,JTOP
2512           IR=IQ+IP-1
2513           CLDFAC(IR,I,J)=CLDIPT(I,J,IP)
2514     71    CONTINUE
2515     1     CONTINUE
2516     
2517       END SUBROUTINE CLO89
2518     !----------------------------------------------------------------
2519           SUBROUTINE LWR88(HEATRA,GRNFLX,TOPFLX,                         &
2520                            PRESS,TEMP,RH2O,QO3,CLDFAC,                   &
2521                            CAMT,NCLDS,KTOP,KBTM,                         &
2522                            BO3RND,AO3RND, &
2523                            APCM,BPCM,ATPCM,BTPCM,ACOMB,BCOMB,BETACM,     &
2524                            ZERO,ONE,H18E3,P0INV,H6P08108,DIFFCTR,        &
2525                            GINV,H3M4,BETINW,RATH2OMW,GP0INV,P0,P0XZP8,   &
2526                            P0XZP2,H3M3,H1M3,H1M2,H25E2,B0,B2,B1,B3,HAF,  &
2527                            TEN,HP1,FOUR,HM1EZ,                           &
2528                            RADCON,QUARTR,TWO,                            &
2529                            HM6666M2,HMP66667,HMP5, HP166666,H41666M2,    &
2530                            RADCON1,H16E1, H28E1,H44194M2,H1P41819,       &
2531                            ids,ide, jds,jde, kds,kde,                    &
2532                            ims,ime, jms,jme, kms,kme,                    &
2533                            its,ite, jts,jte, kts,kte ,Jndx               )
2534     !---------------------------------------------------------------------
2535      IMPLICIT NONE
2536     !----------------------------------------------------------------------
2537     
2538           INTEGER, INTENT(IN)        :: ids,ide, jds,jde, kds,kde ,      &
2539                                         ims,ime, jms,jme, kms,kme ,      &
2540                                         its,ite, jts,jte, kts,kte , Jndx
2541           REAL,    INTENT(IN)        :: ZERO,ONE,H18E3,P0INV,H6P08108,DIFFCTR
2542           REAL,    INTENT(IN)        :: GINV,H3M4,BETINW,RATH2OMW,GP0INV
2543           REAL,    INTENT(IN)        :: P0XZP8,P0XZP2,H3M3,P0,H1M3
2544           REAL,    INTENT(IN)        :: H1M2,H25E2,B0,B1,B2,B3,HAF
2545           REAL,    INTENT(IN)        :: TEN,HP1,FOUR,HM1EZ         
2546           REAL,    INTENT(IN)        :: RADCON,QUARTR,TWO
2547           REAL,    INTENT(IN)        :: HM6666M2,HMP66667,HMP5, HP166666,H41666M2
2548           REAL,    INTENT(IN) :: RADCON1,H16E1, H28E1,H44194M2,H1P41819
2549     !----------------------------------------------------------------------
2550           REAL, INTENT(IN), DIMENSION(3) :: BO3RND,AO3RND
2551           REAL,INTENT(IN),DIMENSION(NBLY) :: APCM,BPCM,ATPCM,BTPCM,ACOMB, &
2552                                              BCOMB,BETACM
2553     
2554           REAL,    INTENT(IN),DIMENSION(its:ite,kts:kte+1,kts:kte+1) :: CLDFAC
2555           REAL,    INTENT(IN), DIMENSION(its:ite,kts:kte+1) :: CAMT
2556           INTEGER, INTENT(IN), DIMENSION(its:ite,kts:kte+1) :: KBTM,KTOP
2557           INTEGER, INTENT(IN), DIMENSION(its:ite)           :: NCLDS
2558          
2559           REAL,    INTENT(IN), DIMENSION(its:ite,kts:kte+1) :: PRESS,TEMP
2560           REAL,    INTENT(IN), DIMENSION(its:ite,kts:kte)   :: RH2O,QO3
2561           REAL,    INTENT(OUT), DIMENSION(its:ite,kts:kte)   :: HEATRA
2562           REAL,    INTENT(OUT), DIMENSION(its:ite)           :: GRNFLX,TOPFLX
2563     
2564           REAL,    DIMENSION(its:ite,kts:kte+1) :: TEXPSL,TOTPHI,TOTO3,CNTVAL,&
2565                                                    TPHIO3,TOTVO2,TSTDAV,TDAV, & 
2566                                                    VSUM3,CO2R1,D2CD21,DCO2D1, &
2567                                                    CO2R2,D2CD22,DCO2D2,CO2SP1,&
2568                                                    CO2SP2,CO2R,DCO2DT,D2CDT2, &
2569                                                    TLSQU,DIFT
2570           REAL,    DIMENSION(its:ite,kts:kte)   :: DELP2,DELP,CO2NBL,&
2571                                                    QH2O,VV,VAR1,VAR2,VAR3,VAR4
2572           REAL,    DIMENSION(its:ite,kts:kte+1) :: P,T
2573           REAL,    DIMENSION(its:ite,kts:kte)   :: CO2MR,CO2MD,CO2M2D
2574           REAL,    DIMENSION(its:ite,kts:kte*2+1):: EMPL
2575     
2576           REAL,    DIMENSION(its:ite)           :: EMX1,EMX2,VSUM1,VSUM2,A1,A2 
2577           REAL,    DIMENSION(its:ite,kts:kte+1,kts:kte+1) :: CO21
2578     !
2579     !
2580     !****COMPUTE FLUX PRESSURES (P) AND DIFFERENCES (DELP2,DELP)
2581     !****COMPUTE FLUX LEVEL TEMPERATURES (T) AND CONTINUUM TEMPERATURE
2582     !    CORRECTIONS (TEXPSL)
2583         
2584           INTEGER :: K, I,KP, KK
2585           INTEGER :: L,LP1,LP2,LP3,LM1,LM2,LM3,MYIS,MYIE,LLP1,LL
2586     
2587           L=kte
2588           LP1=L+1;  LP2=L+2;  LP3=L+3; LLP1 = 2*L + 1
2589           LM1=L-1;  LM2=L-2;  LM3=L-3; LL = 2*L
2590           MYIS=its; MYIE=ite
2591     
2592     
2593           DO 103 K=2,L
2594           DO 103 I=MYIS,MYIE
2595           P(I,K)=HAF*(PRESS(I,K-1)+PRESS(I,K))
2596           T(I,K)=HAF*(TEMP(I,K-1)+TEMP(I,K))
2597     103   CONTINUE
2598           DO 105 I=MYIS,MYIE
2599           P(I,1)=ZERO
2600           P(I,LP1)=PRESS(I,LP1)
2601           T(I,1)=TEMP(I,1)
2602           T(I,LP1)=TEMP(I,LP1)
2603     105   CONTINUE
2604           DO 107 K=1,L
2605           DO 107 I=MYIS,MYIE
2606           DELP2(I,K)=P(I,K+1)-P(I,K)
2607           DELP(I,K)=ONE/DELP2(I,K)
2608     107   CONTINUE
2609     !****COMPUTE ARGUMENT FOR CONT.TEMP.COEFF.
2610     !    (THIS IS 1800.(1./TEMP-1./296.))
2611           DO 125 K=1,LP1
2612           DO 125 I=MYIS,MYIE
2613           TEXPSL(I,K)=H18E3/TEMP(I,K)-H6P08108
2614     !...THEN TAKE EXPONENTIAL
2615           TEXPSL(I,K)=EXP(TEXPSL(I,K))
2616     125   CONTINUE
2617     !***COMPUTE OPTICAL PATHS FOR H2O AND O3, USING THE DIFFUSIVITY
2618     !   APPROXIMATION FOR THE ANGULAR INTEGRATION (1.66). OBTAIN THE
2619     !   UNWEIGHTED VALUES(VAR1,VAR3) AND THE WEIGHTED VALUES(VAR2,VAR4).
2620     !   THE QUANTITIES H3M4(.0003) AND H3M3(.003) APPEARING IN THE VAR2 AND
2621     !   VAR4 EXPRESSIONS ARE THE APPROXIMATE VOIGT CORRECTIONS FOR H2O AND
2622     !   O3,RESPECTIVELY.
2623     !
2624           DO 131 K=1,L
2625           DO 131 I=MYIS,MYIE
2626           QH2O(I,K)=RH2O(I,K)*DIFFCTR
2627     !---VV IS THE LAYER-MEAN PRESSURE (IN ATM),WHICH IS NOT THE SAME AS
2628     !   THE LEVEL PRESSURE (PRESS)
2629           VV(I,K)=HAF*(P(I,K+1)+P(I,K))*P0INV
2630           VAR1(I,K)=DELP2(I,K)*QH2O(I,K)*GINV
2631           VAR3(I,K)=DELP2(I,K)*QO3(I,K)*DIFFCTR*GINV
2632           VAR2(I,K)=VAR1(I,K)*(VV(I,K)+H3M4)
2633           VAR4(I,K)=VAR3(I,K)*(VV(I,K)+H3M3)
2634     !  COMPUTE OPTICAL PATH FOR THE H2O CONTINUUM, USING ROBERTS COEFFS.
2635     !  (BETINW),AND TEMP. CORRECTION (TEXPSL). THE DIFFUSIVITY FACTOR
2636     !  (WHICH CANCELS OUT IN THIS EXPRESSION) IS ASSUMED TO BE 1.66. THE
2637     !  USE OF THE DIFFUSIVITY FACTOR HAS BEEN SHOWN TO BE A SIGNIFICANT
2638     !  SOURCE OF ERROR IN THE CONTINUUM CALCS.,BUT THE TIME PENALTY OF
2639     !  AN ANGULAR INTEGRATION IS SEVERE.
2640     !
2641           CNTVAL(I,K)=TEXPSL(I,K)*RH2O(I,K)*VAR2(I,K)*BETINW/ &
2642                        (RH2O(I,K)+RATH2OMW)
2643     131   CONTINUE
2644     !   COMPUTE SUMMED OPTICAL PATHS FOR H2O,O3 AND CONTINUUM
2645           DO 201 I=MYIS,MYIE
2646           TOTPHI(I,1)=ZERO
2647           TOTO3(I,1)=ZERO
2648           TPHIO3(I,1)=ZERO
2649           TOTVO2(I,1)=ZERO
2650     201   CONTINUE
2651           DO 203 K=2,LP1
2652           DO 203 I=MYIS,MYIE
2653           TOTPHI(I,K)=TOTPHI(I,K-1)+VAR2(I,K-1)
2654           TOTO3(I,K)=TOTO3(I,K-1)+VAR3(I,K-1)
2655           TPHIO3(I,K)=TPHIO3(I,K-1)+VAR4(I,K-1)
2656           TOTVO2(I,K)=TOTVO2(I,K-1)+CNTVAL(I,K-1)
2657     203   CONTINUE
2658     !---EMX1 IS THE ADDITIONAL PRESSURE-SCALED MASS FROM PRESS(L) TO
2659     !   P(L). IT IS USED IN NEARBY LAYER AND EMISS CALCULATIONS.
2660     !---EMX2 IS THE ADDITIONAL PRESSURE-SCALED MASS FROM PRESS(L) TO
2661     !   P(LP1). IT IS USED IN CALCULATIONS BETWEEN FLUX LEVELS L AND LP1.
2662     !
2663           DO 801 I=MYIS,MYIE
2664           EMX1(I)=QH2O(I,L)*PRESS(I,L)*(PRESS(I,L)-P(I,L))*GP0INV
2665           EMX2(I)=QH2O(I,L)*PRESS(I,L)*(P(I,LP1)-PRESS(I,L))*GP0INV
2666     801   CONTINUE
2667     !---EMPL IS THE PRESSURE SCALED MASS FROM P(K) TO PRESS(K) (INDEX 2-LP1)
2668     !   OR TO PRESS(K+1) (INDEX LP2-LL)
2669           DO 811 K=1,L
2670           DO 811 I=MYIS,MYIE
2671           EMPL(I,K+1)=QH2O(I,K)*P(I,K+1)*(P(I,K+1)-PRESS(I,K))*GP0INV
2672     811   CONTINUE
2673           DO 812 K=1,LM1
2674           DO 812 I=MYIS,MYIE
2675           EMPL(I,LP2+K-1)=QH2O(I,K+1)*P(I,K+1)*(PRESS(I,K+1)-P(I,K+1)) &
2676                          *GP0INV
2677     812   CONTINUE
2678           DO 821 I=MYIS,MYIE
2679           EMPL(I,1)=VAR2(I,L)
2680           EMPL(I,LLP1)=EMPL(I,LL)
2681     821   CONTINUE
2682     !***COMPUTE WEIGHTED TEMPERATURE (TDAV) AND PRESSURE (TSTDAV) INTEGRALS
2683     !   FOR USE IN OBTAINING TEMP. DIFFERENCE BET. SOUNDING AND STD.
2684     !   TEMP. SOUNDING (DIFT)
2685           DO 161 I=MYIS,MYIE
2686           TSTDAV(I,1)=ZERO
2687           TDAV(I,1)=ZERO
2688     161   CONTINUE
2689           DO 162 K=1,LP1
2690           DO 162 I=MYIS,MYIE
2691           VSUM3(I,K)=TEMP(I,K)-STEMP(K)
2692     162   CONTINUE
2693           DO 163 K=1,L
2694           DO 165 I=MYIS,MYIE
2695           VSUM2(I)=GTEMP(K)*DELP2(I,K)
2696           VSUM1(I)=VSUM2(I)*VSUM3(I,K)
2697           TSTDAV(I,K+1)=TSTDAV(I,K)+VSUM2(I)
2698           TDAV(I,K+1)=TDAV(I,K)+VSUM1(I)
2699     165   CONTINUE
2700     163   CONTINUE
2701     !
2702     !****EVALUATE COEFFICIENTS FOR CO2 PRESSURE INTERPOLATION (A1,A2)
2703           DO 171 I=MYIS,MYIE
2704           A1(I)=(PRESS(I,LP1)-P0XZP8)/P0XZP2
2705           A2(I)=(P0-PRESS(I,LP1))/P0XZP2
2706     171   CONTINUE
2707     !***PERFORM CO2 PRESSURE INTERPOLATION ON ALL INPUTTED TRANSMISSION
2708     !   FUNCTIONS AND TEMP. DERIVATIVES
2709     !---SUCCESSIVELY COMPUTING CO2R,DCO2DT AND D2CDT2 IS DONE TO SAVE
2710     !   STORAGE (AT A SLIGHT LOSS IN COMPUTATION TIME)
2711           DO 184 K=1,LP1
2712           DO 184 I=MYIS,MYIE
2713             CO2R1(I,K)=A1(I)*CO231(K)+A2(I)*CO238(K)
2714             D2CD21(I,K)=H1M3*(A1(I)*C2D31(K)+A2(I)*C2D38(K))
2715             DCO2D1(I,K)=H1M2*(A1(I)*CDT31(K)+A2(I)*CDT38(K))
2716             CO2R2(I,K)=A1(I)*CO271(K)+A2(I)*CO278(K)
2717             D2CD22(I,K)=H1M3*(A1(I)*C2D71(K)+A2(I)*C2D78(K))
2718             DCO2D2(I,K)=H1M2*(A1(I)*CDT71(K)+A2(I)*CDT78(K))
2719     184   CONTINUE
2720           DO 190 K=1,L
2721           DO 190 I=MYIS,MYIE
2722             CO2MR(I,K)=A1(I)*CO2M51(K)+A2(I)*CO2M58(K)
2723             CO2MD(I,K)=H1M2*(A1(I)*CDTM51(K)+A2(I)*CDTM58(K))
2724             CO2M2D(I,K)=H1M3*(A1(I)*C2DM51(K)+A2(I)*C2DM58(K))
2725     190   CONTINUE
2726     !***COMPUTE CO2 TEMPERATURE INTERPOLATIONS FOR ALL BANDS,USING DIFT
2727     !
2728     !   THE CASE WHERE K=1 IS HANDLED FIRST. WE ARE NOW REPLACING
2729     !   3-DIMENSIONAL ARRAYS BY 2-D ARRAYS, TO SAVE SPACE. THUS THIS
2730     !   CALCULATION IS FOR (I,KP,1)
2731           DO 211 KP=2,LP1
2732           DO 211 I=MYIS,MYIE
2733           DIFT(I,KP)=TDAV(I,KP)/TSTDAV(I,KP)
2734     211   CONTINUE
2735           DO 212 I=MYIS,MYIE
2736           CO21(I,1,1)=1.0
2737           CO2SP1(I,1)=1.0
2738           CO2SP2(I,1)=1.0
2739     212   CONTINUE
2740           DO 215 KP=2,LP1
2741           DO 215 I=MYIS,MYIE
2742     !---CALCULATIONS FOR KP>1 FOR K=1
2743           CO2R(I,KP)=A1(I)*CO251(KP,1)+A2(I)*CO258(KP,1)
2744           DCO2DT(I,KP)=H1M2*(A1(I)*CDT51(KP,1)+A2(I)*CDT58(KP,1))
2745           D2CDT2(I,KP)=H1M3*(A1(I)*C2D51(KP,1)+A2(I)*C2D58(KP,1))
2746           CO21(I,KP,1)=CO2R(I,KP)+DIFT(I,KP)*(DCO2DT(I,KP)+ &
2747                        HAF*DIFT(I,KP)*D2CDT2(I,KP))
2748     !---CALCULATIONS FOR (EFFECTIVELY) KP=1,K>KP. THESE USE THE
2749     !   SAME VALUE OF DIFT DUE TO SYMMETRY
2750           CO2R(I,KP)=A1(I)*CO251(1,KP)+A2(I)*CO258(1,KP)
2751           DCO2DT(I,KP)=H1M2*(A1(I)*CDT51(1,KP)+A2(I)*CDT58(1,KP))
2752           D2CDT2(I,KP)=H1M3*(A1(I)*C2D51(1,KP)+A2(I)*C2D58(1,KP))
2753           CO21(I,1,KP)=CO2R(I,KP)+DIFT(I,KP)*(DCO2DT(I,KP)+ &
2754                        HAF*DIFT(I,KP)*D2CDT2(I,KP))
2755     215   CONTINUE
2756     !   THE TRANSMISSION FUNCTIONS USED IN SPA88 MAY BE COMPUTED NOW.
2757     !---(IN THE 250 LOOP,DIFT REALLY SHOULD BE (I,1,K), BUT DIFT IS
2758     !    INVARIANT WITH RESPECT TO K,KP,AND SO (I,1,K)=(I,K,1))
2759           DO 250 K=2,LP1
2760           DO 250 I=MYIS,MYIE
2761           CO2SP1(I,K)=CO2R1(I,K)+DIFT(I,K)*(DCO2D1(I,K)+HAF*DIFT(I,K)* &
2762            D2CD21(I,K))
2763           CO2SP2(I,K)=CO2R2(I,K)+DIFT(I,K)*(DCO2D2(I,K)+HAF*DIFT(I,K)* &
2764            D2CD22(I,K))
2765     250   CONTINUE
2766     !
2767     !   NEXT THE CASE WHEN K=2...L
2768           DO 220 K=2,L
2769           DO 222 KP=K+1,LP1
2770           DO 222 I=MYIS,MYIE
2771           DIFT(I,KP)=(TDAV(I,KP)-TDAV(I,K))/ &
2772                         (TSTDAV(I,KP)-TSTDAV(I,K))
2773           CO2R(I,KP)=A1(I)*CO251(KP,K)+A2(I)*CO258(KP,K)
2774           DCO2DT(I,KP)=H1M2*(A1(I)*CDT51(KP,K)+A2(I)*CDT58(KP,K))
2775           D2CDT2(I,KP)=H1M3*(A1(I)*C2D51(KP,K)+A2(I)*C2D58(KP,K))
2776           CO21(I,KP,K)=CO2R(I,KP)+DIFT(I,KP)*(DCO2DT(I,KP)+ &
2777                        HAF*DIFT(I,KP)*D2CDT2(I,KP))
2778           CO2R(I,KP)=A1(I)*CO251(K,KP)+A2(I)*CO258(K,KP)
2779           DCO2DT(I,KP)=H1M2*(A1(I)*CDT51(K,KP)+A2(I)*CDT58(K,KP))
2780           D2CDT2(I,KP)=H1M3*(A1(I)*C2D51(K,KP)+A2(I)*C2D58(K,KP))
2781           CO21(I,K,KP)=CO2R(I,KP)+DIFT(I,KP)*(DCO2DT(I,KP)+ &
2782                        HAF*DIFT(I,KP)*D2CDT2(I,KP))
2783     222   CONTINUE
2784     220   CONTINUE
2785     !   FINALLY THE CASE WHEN K=KP,K=2..LP1
2786           DO 206 K=2,LP1
2787           DO 206 I=MYIS,MYIE
2788           DIFT(I,K)=HAF*(VSUM3(I,K)+VSUM3(I,K-1))
2789           CO2R(I,K)=A1(I)*CO251(K,K)+A2(I)*CO258(K,K)
2790           DCO2DT(I,K)=H1M2*(A1(I)*CDT51(K,K)+A2(I)*CDT58(K,K))
2791           D2CDT2(I,K)=H1M3*(A1(I)*C2D51(K,K)+A2(I)*C2D58(K,K))
2792           CO21(I,K,K)=CO2R(I,K)+DIFT(I,K)*(DCO2DT(I,K)+ &
2793                        HAF*DIFT(I,K)*D2CDT2(I,K))
2794     206   CONTINUE
2795     !--- WE AREN'T DOING NBL TFS ON THE 100 CM-1 BANDS .
2796           DO 260 K=1,L
2797           DO 260 I=MYIS,MYIE
2798           CO2NBL(I,K)=CO2MR(I,K)+VSUM3(I,K)*(CO2MD(I,K)+HAF* &
2799            VSUM3(I,K)*CO2M2D(I,K))
2800     260   CONTINUE
2801     !***COMPUTE TEMP. COEFFICIENT BASED ON T(K) (SEE REF.2)
2802           DO 264 K=1,LP1
2803           DO 264 I=MYIS,MYIE
2804           IF (T(I,K).LE.H25E2) THEN
2805              TLSQU(I,K)=B0+(T(I,K)-H25E2)* &
2806                                 (B1+(T(I,K)-H25E2)* &
2807                              (B2+B3*(T(I,K)-H25E2)))
2808           ELSE
2809              TLSQU(I,K)=B0
2810           ENDIF
2811     264   CONTINUE
2812     !***APPLY TO ALL CO2 TFS
2813           DO 280 K=1,LP1
2814           DO 282 KP=1,LP1
2815           DO 282 I=MYIS,MYIE
2816           CO21(I,KP,K)=CO21(I,KP,K)*(ONE-TLSQU(I,KP))+TLSQU(I,KP)
2817     282   CONTINUE
2818     280   CONTINUE
2819           DO 284 K=1,LP1
2820           DO 286 I=MYIS,MYIE
2821           CO2SP1(I,K)=CO2SP1(I,K)*(ONE-TLSQU(I,1))+TLSQU(I,1)
2822           CO2SP2(I,K)=CO2SP2(I,K)*(ONE-TLSQU(I,1))+TLSQU(I,1)
2823     286   CONTINUE
2824     284   CONTINUE
2825           DO 288 K=1,L
2826           DO 290 I=MYIS,MYIE
2827           CO2NBL(I,K)=CO2NBL(I,K)*(ONE-TLSQU(I,K))+TLSQU(I,K)
2828     290   CONTINUE
2829     288   CONTINUE
2830     
2831           CALL FST88(HEATRA,GRNFLX,TOPFLX, &
2832                      QH2O,PRESS,P,DELP,DELP2,TEMP,T, &
2833                      CLDFAC,NCLDS,KTOP,KBTM,CAMT, &
2834                      CO21,CO2NBL,CO2SP1,CO2SP2, &
2835                      VAR1,VAR2,VAR3,VAR4,CNTVAL, &
2836                      TOTO3,TPHIO3,TOTPHI,TOTVO2, &
2837                      EMX1,EMX2,EMPL, &
2838                      BO3RND,AO3RND, &
2839                      APCM,BPCM,ATPCM,BTPCM,ACOMB,BCOMB,BETACM,     &
2840                      TEN,HP1,HAF,ONE,FOUR,HM1EZ,       &
2841                      RADCON,QUARTR,TWO,  &
2842                      HM6666M2,HMP66667,HMP5, &
2843                      HP166666,H41666M2,RADCON1, &
2844                      H16E1, H28E1, H25E2, H44194M2,H1P41819, &
2845                      ids,ide, jds,jde, kds,kde,                    &
2846                      ims,ime, jms,jme, kms,kme,                    &
2847                      its,ite, jts,jte, kts,kte,Jndx                 )
2848     
2849       END SUBROUTINE LWR88
2850     !---------------------------------------------------------------------
2851       SUBROUTINE FST88(HEATRA,GRNFLX,TOPFLX, &
2852                            QH2O,PRESS,P,DELP,DELP2,TEMP,T, &
2853                            CLDFAC,NCLDS,KTOP,KBTM,CAMT, &
2854                            CO21,CO2NBL,CO2SP1,CO2SP2, &
2855                            VAR1,VAR2,VAR3,VAR4,CNTVAL, &
2856                            TOTO3,TPHIO3,TOTPHI,TOTVO2, &
2857                            EMX1,EMX2,EMPL, &
2858                            BO3RND,AO3RND, &
2859                            APCM,BPCM,ATPCM,BTPCM,ACOMB,BCOMB,BETACM,     &
2860                            TEN,HP1,HAF,ONE,FOUR,HM1EZ,       &
2861                            RADCON,QUARTR,TWO, &
2862                            HM6666M2,HMP66667,HMP5, &
2863                            HP166666,H41666M2,RADCON1, &
2864                            H16E1, H28E1, H25E2, H44194M2,H1P41819, &
2865                            ids,ide, jds,jde, kds,kde,                    &
2866                            ims,ime, jms,jme, kms,kme,                    &
2867                            its,ite, jts,jte, kts,kte,Jndx                )
2868     !---------------------------------------------------------------------
2869      IMPLICIT NONE
2870     !----------------------------------------------------------------------
2871     
2872           INTEGER, INTENT(IN)        :: ids,ide, jds,jde, kds,kde ,      &
2873                                         ims,ime, jms,jme, kms,kme ,      &
2874                                         its,ite, jts,jte, kts,kte ,Jndx
2875     
2876           REAL,    INTENT(IN)        :: TEN,HP1,HAF,ONE,FOUR,HM1EZ
2877           REAL,    INTENT(IN)        :: RADCON,QUARTR,TWO
2878           REAL,    INTENT(IN)        :: HM6666M2,HMP66667,HMP5
2879           REAL,    INTENT(IN)        :: HP166666,H41666M2,RADCON1,H16E1, H28E1 
2880           REAL,    INTENT(IN)        :: H25E2,H44194M2,H1P41819
2881     
2882           REAL,INTENT(IN),DIMENSION(NBLY) :: APCM,BPCM,ATPCM,BTPCM,ACOMB, &
2883                                              BCOMB,BETACM
2884     
2885           REAL, INTENT(IN), DIMENSION(its:ite,kts:kte*2+1) :: EMPL
2886           REAL, INTENT(IN), DIMENSION(its:ite,kts:kte+1) :: TOTO3,TPHIO3,TOTPHI,CNTVAL,&
2887                                                             CO2SP1,CO2SP2   
2888     
2889           REAL,    INTENT(IN),DIMENSION(its:ite,kts:kte+1,kts:kte+1) :: CLDFAC
2890           REAL,    INTENT(IN), DIMENSION(its:ite,kts:kte+1) :: CAMT,TOTVO2
2891           INTEGER, INTENT(IN), DIMENSION(its:ite,kts:kte+1) :: KBTM,KTOP
2892           INTEGER, INTENT(IN), DIMENSION(its:ite)           :: NCLDS
2893           REAL,    INTENT(IN), DIMENSION(its:ite,kts:kte)   :: QH2O
2894           REAL,    INTENT(IN), DIMENSION(its:ite,kts:kte+1) :: PRESS,TEMP
2895           REAL,    INTENT(OUT), DIMENSION(its:ite,kts:kte)  :: HEATRA
2896           REAL,    INTENT(OUT), DIMENSION(its:ite)          :: GRNFLX,TOPFLX
2897           REAL,    INTENT(IN), DIMENSION(its:ite,kts:kte+1) :: P,T
2898           REAL,    INTENT(INOUT), DIMENSION(its:ite,kts:kte+1,kts:kte+1) :: CO21
2899           REAL,    INTENT(IN), DIMENSION(its:ite,kts:kte)   :: CO2NBL,DELP2, &
2900                                                                DELP,&
2901                                                    VAR1,VAR2,VAR3,VAR4
2902           REAL, INTENT(IN), DIMENSION(3) :: BO3RND,AO3RND
2903           REAL, INTENT(IN), DIMENSION(its:ite)   :: EMX1,EMX2
2904           
2905           REAL, DIMENSION(its:ite,kts:kte*2+1) :: TPL,EMD,ALP,C,CSUB,CSUB2
2906           REAL, DIMENSION(its:ite,kts:kte*2+1) :: C2
2907           INTEGER, DIMENSION(its:ite,kts:kte+1) :: IXO
2908           REAL, DIMENSION(its:ite,kts:kte+1) :: VTMP3,FXO,DT,FXOE2,DTE2, &
2909                                                 SS1,CSOUR,TC,OSS,CSS,DTC,SS2,&
2910                                                 AVEPHI,E1CTS1,E1FLX,  &
2911                                                 E1CTW1,DSORC,EMISS,FAC1,&
2912                                                 TO3SP,OVER1D,CNTTAU,TOTEVV,&
2913                                                 CO2SP,FLX,AVMO3, &
2914                                                 AVPHO3,AVVO2,CONT1D,TO31D,EMISDG,&
2915                                                 DELPR1
2916           REAL, DIMENSION(its:ite,kts:kte+1) :: EMISSB,DELPR2,CONTDG,TO3DG,HEATEM,&
2917                                                 VSUM1,FLXNET,Z1
2918     
2919           REAL, DIMENSION(its:ite,kts:kte+1,NBLY) :: SORC
2920           REAL, DIMENSION(its:ite,kts:kte)   :: E1CTS2,E1CTW2,TO3SPC,RLOG,EXCTS,&
2921                                                 CTSO3,CTS
2922           REAL, DIMENSION(its:ite)   :: GXCTS,FLX1E1
2923           REAL, DIMENSION(its:ite)   :: PTOP,PBOT,FTOP,FBOT,DELPTC
2924           REAL, DIMENSION(its:ite,2) :: FXOSP,DTSP,EMSPEC
2925           INTEGER :: K, I,KP,LLM2,J1,J3,KMAX,KMIN,KCLDS,ICNT,LLM1
2926           INTEGER :: L,LP1,LP2,LP3,LM1,LM2,LM3,MYIS,MYIE,LLP1,LL,KK,KLEN
2927     
2928           L=kte
2929           LP1=L+1;  LP2=L+2;  LP3=L+3; LLP1 = 2*L + 1
2930           LM1=L-1;  LM2=L-2;  LM3=L-3; LL = 2*L
2931           LLM2 = LL-2; LLM1=LL-1
2932           MYIS=its; MYIE=ite
2933     
2934     !
2935           DO 101 K=1,LP1
2936           DO 101 I=MYIS,MYIE
2937     !---TEMP. INDICES FOR E1,SOURCE
2938           VTMP3(I,K)=AINT(TEMP(I,K)*HP1)
2939           FXO(I,K)=VTMP3(I,K)-9.
2940           DT(I,K)=TEMP(I,K)-TEN*VTMP3(I,K)
2941     !---INTEGER INDEX FOR SOURCE (USED IMMEDIATELY)
2942           IXO(I,K)=FXO(I,K)
2943     101   CONTINUE
2944           DO 103 k=1,L
2945           DO 103 I=MYIS,MYIE
2946     !---TEMP. INDICES FOR E2 (KP=1 LAYER NOT USED IN FLUX CALCULATIONS)
2947           VTMP3(I,K)=AINT(T(I,K+1)*HP1)
2948           FXOE2(I,K)=VTMP3(I,K)-9.
2949           DTE2(I,K)=T(I,K+1)-TEN*VTMP3(I,K)
2950     103   CONTINUE
2951     !---SPECIAL CASE TO HANDLE KP=LP1 LAYER AND SPECIAL E2 CALCS.
2952           DO 105 I=MYIS,MYIE
2953           FXOE2(I,LP1)=FXO(I,L)
2954           DTE2(I,LP1)=DT(I,L)
2955           FXOSP(I,1)=FXOE2(I,LM1)
2956           FXOSP(I,2)=FXO(I,LM1)
2957           DTSP(I,1)=DTE2(I,LM1)
2958           DTSP(I,2)=DT(I,LM1)
2959     105   CONTINUE
2960     !
2961     !---SOURCE FUNCTION FOR COMBINED BAND 1
2962           DO 4114 I=MYIS,MYIE
2963           DO 4114 K=1,LP1
2964             VTMP3(I,K)=SOURCE(IXO(I,K),1)
2965             DSORC(I,K)=DSRCE(IXO(I,K),1)
2966     4114   CONTINUE
2967           DO 4112 K=1,LP1
2968           DO 4112 I=MYIS,MYIE
2969           SORC(I,K,1)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
2970     4112   CONTINUE
2971     !---SOURCE FUNCTION FOR COMBINED BAND 2
2972           DO 4214 I=MYIS,MYIE
2973           DO 4214 K=1,LP1
2974             VTMP3(I,K)=SOURCE(IXO(I,K),2)
2975             DSORC(I,K)=DSRCE(IXO(I,K),2)
2976     4214   CONTINUE
2977           DO 4212 K=1,LP1
2978           DO 4212 I=MYIS,MYIE
2979           SORC(I,K,2)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
2980     4212   CONTINUE
2981     !---SOURCE FUNCTION FOR COMBINED BAND 3
2982           DO 4314 I=MYIS,MYIE
2983           DO 4314 K=1,LP1
2984             VTMP3(I,K)=SOURCE(IXO(I,K),3)
2985             DSORC(I,K)=DSRCE(IXO(I,K),3)
2986     4314   CONTINUE
2987           DO 4312 K=1,LP1
2988           DO 4312 I=MYIS,MYIE
2989           SORC(I,K,3)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
2990     4312   CONTINUE
2991     !---SOURCE FUNCTION FOR COMBINED BAND 4
2992           DO 4414 I=MYIS,MYIE
2993           DO 4414 K=1,LP1
2994             VTMP3(I,K)=SOURCE(IXO(I,K),4)
2995             DSORC(I,K)=DSRCE(IXO(I,K),4)
2996     4414   CONTINUE
2997           DO 4412 K=1,LP1
2998           DO 4412 I=MYIS,MYIE
2999           SORC(I,K,4)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
3000     4412   CONTINUE
3001     !---SOURCE FUNCTION FOR COMBINED BAND 5
3002           DO 4514 I=MYIS,MYIE
3003           DO 4514 K=1,LP1
3004             VTMP3(I,K)=SOURCE(IXO(I,K),5)
3005             DSORC(I,K)=DSRCE(IXO(I,K),5)
3006     4514   CONTINUE
3007           DO 4512 K=1,LP1
3008           DO 4512 I=MYIS,MYIE
3009           SORC(I,K,5)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
3010     4512   CONTINUE
3011     !---SOURCE FUNCTION FOR COMBINED BAND 6
3012           DO 4614 I=MYIS,MYIE
3013           DO 4614 K=1,LP1
3014             VTMP3(I,K)=SOURCE(IXO(I,K),6)
3015             DSORC(I,K)=DSRCE(IXO(I,K),6)
3016     4614   CONTINUE
3017           DO 4612 K=1,LP1
3018           DO 4612 I=MYIS,MYIE
3019           SORC(I,K,6)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
3020     4612   CONTINUE
3021     !---SOURCE FUNCTION FOR COMBINED BAND 7
3022           DO 4714 I=MYIS,MYIE
3023           DO 4714 K=1,LP1
3024             VTMP3(I,K)=SOURCE(IXO(I,K),7)
3025             DSORC(I,K)=DSRCE(IXO(I,K),7)
3026     4714   CONTINUE
3027           DO 4712 K=1,LP1
3028           DO 4712 I=MYIS,MYIE
3029           SORC(I,K,7)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
3030     4712   CONTINUE
3031     !---SOURCE FUNCTION FOR COMBINED BAND 8
3032           DO 4814 I=MYIS,MYIE
3033           DO 4814 K=1,LP1
3034             VTMP3(I,K)=SOURCE(IXO(I,K),8)
3035             DSORC(I,K)=DSRCE(IXO(I,K),8)
3036     4814   CONTINUE
3037           DO 4812 K=1,LP1
3038           DO 4812 I=MYIS,MYIE
3039           SORC(I,K,8)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
3040     4812   CONTINUE
3041     !---SOURCE FUNCTION FOR BAND 9 (560-670 CM-1)
3042           DO 4914 I=MYIS,MYIE
3043           DO 4914 K=1,LP1
3044             VTMP3(I,K)=SOURCE(IXO(I,K),9)
3045             DSORC(I,K)=DSRCE(IXO(I,K),9)
3046     4914   CONTINUE
3047           DO 4912 K=1,LP1
3048           DO 4912 I=MYIS,MYIE
3049           SORC(I,K,9)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
3050     4912   CONTINUE
3051     !---SOURCE FUNCTION FOR BAND 10 (670-800 CM-1)
3052           DO 5014 I=MYIS,MYIE
3053           DO 5014 K=1,LP1
3054             VTMP3(I,K)=SOURCE(IXO(I,K),10)
3055             DSORC(I,K)=DSRCE(IXO(I,K),10)
3056     5014  CONTINUE
3057           DO 5012 K=1,LP1
3058           DO 5012 I=MYIS,MYIE
3059           SORC(I,K,10)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
3060     5012   CONTINUE
3061     !---SOURCE FUNCTION FOR BAND 11 (800-900 CM-1)
3062           DO 5114 I=MYIS,MYIE
3063           DO 5114 K=1,LP1
3064             VTMP3(I,K)=SOURCE(IXO(I,K),11)
3065             DSORC(I,K)=DSRCE(IXO(I,K),11)
3066     5114   CONTINUE
3067           DO 5112 K=1,LP1
3068           DO 5112 I=MYIS,MYIE
3069           SORC(I,K,11)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
3070     5112   CONTINUE
3071     !---SOURCE FUNCTION FOR BAND 12 (900-990 CM-1)
3072           DO 5214 I=MYIS,MYIE
3073           DO 5214 K=1,LP1
3074             VTMP3(I,K)=SOURCE(IXO(I,K),12)
3075             DSORC(I,K)=DSRCE(IXO(I,K),12)
3076     5214   CONTINUE
3077           DO 5212 K=1,LP1
3078           DO 5212 I=MYIS,MYIE
3079           SORC(I,K,12)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
3080     5212   CONTINUE
3081     !---SOURCE FUNCTION FOR BAND 13 (990-1070 CM-1)
3082           DO 5314 I=MYIS,MYIE
3083           DO 5314 K=1,LP1
3084             VTMP3(I,K)=SOURCE(IXO(I,K),13)
3085             DSORC(I,K)=DSRCE(IXO(I,K),13)
3086     5314   CONTINUE
3087           DO 5312 K=1,LP1
3088           DO 5312 I=MYIS,MYIE
3089           SORC(I,K,13)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
3090     5312   CONTINUE
3091     !---SOURCE FUNCTION FOR BAND 14 (1070-1200 CM-1)
3092           DO 5414 I=MYIS,MYIE
3093           DO 5414 K=1,LP1
3094             VTMP3(I,K)=SOURCE(IXO(I,K),14)
3095             DSORC(I,K)=DSRCE(IXO(I,K),14)
3096     5414   CONTINUE
3097           DO 5412 K=1,LP1
3098           DO 5412 I=MYIS,MYIE
3099           SORC(I,K,14)=VTMP3(I,K)+DT(I,K)*DSORC(I,K)
3100     5412   CONTINUE
3101     !
3102     !---OBTAIN SPECIAL SOURCE FUNCTIONS FOR THE 15 UM BAND (CSOUR)
3103     !   AND THE WINDOW REGION (SS1)
3104           DO 131 K=1,LP1
3105           DO 131 I=MYIS,MYIE
3106           SS1(I,K)=SORC(I,K,11)+SORC(I,K,12)+SORC(I,K,14)
3107     131   CONTINUE
3108           DO 143 K=1,LP1
3109           DO 143 I=MYIS,MYIE
3110           CSOUR(I,K)=SORC(I,K,9)+SORC(I,K,10)
3111     143   CONTINUE
3112     !
3113     !---COMPUTE TEMP**4 (TC) AND VERTICAL TEMPERATURE DIFFERENCES
3114     !   (OSS,CSS,SS2,DTC). ALL THESE WILL BE USED LATER IN FLUX COMPUTA-
3115     !   TIONS.
3116     !
3117           DO 901 K=1,LP1
3118           DO 901 I=MYIS,MYIE
3119           TC(I,K)=TEMP(I,K)*TEMP(I,K)*TEMP(I,K)*TEMP(I,K)
3120     901   CONTINUE
3121           DO 903 K=1,L
3122           DO 903 I=MYIS,MYIE
3123           OSS(I,K+1)=SORC(I,K+1,13)-SORC(I,K,13)
3124           CSS(I,K+1)=CSOUR(I,K+1)-CSOUR(I,K)
3125           DTC(I,K+1)=TC(I,K+1)-TC(I,K)
3126           SS2(I,K+1)=SS1(I,K+1)-SS1(I,K)
3127     903   CONTINUE
3128     !
3129     !
3130     !---THE FOLLOWIMG IS A DRASTIC REWRITE OF THE RADIATION CODE TO
3131     !    (LARGELY) ELIMINATE THREE-DIMENSIONAL ARRAYS. THE CODE WORKS
3132     !    ON THE FOLLOWING PRINCIPLES:
3133     !
3134     !          LET K = FIXED FLUX LEVEL, KP = VARYING FLUX LEVEL
3135     !          THEN FLUX(K)=SUM OVER KP : (DELTAB(KP)*TAU(KP,K))
3136     !               OVER ALL KP'S, FROM 1 TO LP1.
3137     !
3138     !          WE CAN BREAK DOWN THE CALCULATIONS FOR ALL K'S AS FOLLOWS:
3139     !
3140     !          FOR ALL K'S K=1 TO LP1:
3141     !              FLUX(K)=SUM OVER KP : (DELTAB(KP)*TAU(KP,K))  (1)
3142     !                      OVER ALL KP'S, FROM K+1 TO LP1
3143     !          AND
3144     !              FOR KP FROM K+1 TO LP1:
3145     !                 FLUX(KP) = DELTAB(K)*TAU(K,KP)              (2)
3146     !
3147     !          NOW IF TAU(K,KP)=TAU(KP,K) (SYMMETRICAL ARRAYS)
3148     !          WE CAN COMPUTE A 1-DIMENSIONAL ARRAY TAU1D(KP) FROM
3149     !          K+1 TO LP1, EACH TIME K IS INCREMENTED.
3150     !          EQUATIONS (1) AND (2) THEN BECOME:
3151     !
3152     !             TAU1D(KP) = (VALUES FOR TAU(KP,K) AT THE PARTICULAR K)
3153     !             FLUX(K) = SUM OVER KP : (DELTAB(KP)*TAU1D(KP))   (3)
3154     !             FLUX(KP) = DELTAB(K)*TAU1D(KP)                   (4)
3155     !
3156     !         THE TERMS FOR TAU (K,K) AND OTHER SPECIAL TERMS (FOR
3157     !         NEARBY LAYERS) MUST, OF COURSE, BE HANDLED SEPARATELY, AND
3158     !         WITH CARE.
3159     !
3160     !      COMPUTE "UPPER TRIANGLE" TRANSMISSION FUNCTIONS FOR
3161     !      THE 9.6 UM BAND (TO3SP) AND THE 15 UM BAND (OVER1D). ALSO,
3162     !      THE
3163     !      STAGE 1...COMPUTE O3 ,OVER TRANSMISSION FCTNS AND AVEPHI
3164     !---DO K=1 CALCULATION (FROM FLUX LAYER KK TO THE TOP) SEPARATELY
3165     !   AS VECTORIZATION IS IMPROVED,AND OZONE CTS TRANSMISSIVITY
3166     !   MAY BE EXTRACTED HERE.
3167           DO 3021 K=1,L
3168           DO 3021 I=MYIS,MYIE
3169           AVEPHI(I,K)=TOTPHI(I,K+1)
3170     3021  CONTINUE
3171     !---IN ORDER TO PROPERLY EVALUATE EMISS INTEGRATED OVER THE (LP1)
3172     !   LAYER, A SPECIAL EVALUATION OF EMISS IS DONE. THIS REQUIRES
3173     !   A SPECIAL COMPUTATION OF AVEPHI, AND IT IS STORED IN THE
3174     !   (OTHERWISE VACANT) LP1'TH POSITION
3175     !
3176           DO 803 I=MYIS,MYIE
3177           AVEPHI(I,LP1)=AVEPHI(I,LM1)+EMX1(I)
3178     803   CONTINUE
3179     !   COMPUTE FLUXES FOR K=1
3180           CALL E1E290(E1CTS1,E1CTS2,E1FLX,E1CTW1,E1CTW2,EMISS, &
3181                       FXO,DT,FXOE2,DTE2,AVEPHI,TEMP,T,         &
3182                       H16E1,TEN,HP1,H28E1,HAF,                 &
3183                       ids,ide, jds,jde, kds,kde,               &
3184                       ims,ime, jms,jme, kms,kme,               &
3185                       its,ite, jts,jte, kts,kte                )
3186     
3187           DO 302 K=1,L
3188           DO 302 I=MYIS,MYIE
3189           FAC1(I,K)=BO3RND(2)*TPHIO3(I,K+1)/TOTO3(I,K+1)
3190           TO3SPC(I,K)=HAF*(FAC1(I,K)* &
3191               (SQRT(ONE+(FOUR*AO3RND(2)*TOTO3(I,K+1))/FAC1(I,K))-ONE))
3192     !   FOR K=1, TO3SP IS USED INSTEAD OF TO31D (THEY ARE EQUAL IN THIS
3193     !   CASE); TO3SP IS PASSED TO SPA90, WHILE TO31D IS A WORK-ARRAY.
3194           TO3SP(I,K)=EXP(HM1EZ*(TO3SPC(I,K)+SKO3R*TOTVO2(I,K+1)))
3195           OVER1D(I,K)=EXP(HM1EZ*(SQRT(AB15WD*TOTPHI(I,K+1))+ &
3196                       SKC1R*TOTVO2(I,K+1)))
3197     !---BECAUSE ALL CONTINUUM TRANSMISSIVITIES ARE OBTAINED FROM THE
3198     !  2-D QUANTITY CNTTAU (AND ITS RECIPROCAL TOTEVV) WE STORE BOTH
3199     !  OF THESE HERE. FOR K=1, CONT1D EQUALS CNTTAU
3200           CNTTAU(I,K)=EXP(HM1EZ*TOTVO2(I,K+1))
3201           TOTEVV(I,K)=1./CNTTAU(I,K)
3202     302   CONTINUE
3203           DO 3022 K=1,L
3204           DO 3022 I=MYIS,MYIE
3205           CO2SP(I,K+1)=OVER1D(I,K)*CO21(I,1,K+1)
3206     3022  CONTINUE
3207           DO 3023 K=1,L
3208           DO 3023 I=MYIS,MYIE
3209           CO21(I,K+1,1)=CO21(I,K+1,1)*OVER1D(I,K)
3210     3023  CONTINUE
3211     !---RLOG IS THE NBL AMOUNT FOR THE 15 UM BAND CALCULATION
3212           DO 1808 I=MYIS,MYIE
3213           RLOG(I,1)=OVER1D(I,1)*CO2NBL(I,1)
3214     1808  CONTINUE
3215     !---THE TERMS WHEN KP=1 FOR ALL K ARE THE PHOTON EXCHANGE WITH
3216     !   THE TOP OF THE ATMOSPHERE, AND ARE OBTAINED DIFFERENTLY THAN
3217     !   THE OTHER CALCULATIONS
3218           DO 305 K=2,LP1
3219           DO 305 I=MYIS,MYIE
3220           FLX(I,K)= (TC(I,1)*E1FLX(I,K) &
3221                     +SS1(I,1)*CNTTAU(I,K-1) &
3222                     +SORC(I,1,13)*TO3SP(I,K-1) &
3223                     +CSOUR(I,1)*CO2SP(I,K)) &
3224                     *CLDFAC(I,1,K)
3225     
3226     	if (I .eq. 50 .and. K .eq. 10) then
3227     !	write(0,*) 'E1FLX(I,K), FLX(I,K): ', E1FLX(I,K), FLX(I,K)
3228     	endif
3229     305   CONTINUE
3230           DO 307 I=MYIS,MYIE
3231           FLX(I,1)= TC(I,1)*E1FLX(I,1)+SS1(I,1)+SORC(I,1,13) &
3232                     +CSOUR(I,1)
3233     307   CONTINUE
3234     !---THE KP TERMS FOR K=1...
3235           DO 303 KP=2,LP1
3236           DO 303 I=MYIS,MYIE
3237           FLX(I,1)=FLX(I,1)+(OSS(I,KP)*TO3SP(I,KP-1) &
3238                             +SS2(I,KP)*CNTTAU(I,KP-1) &
3239                             +CSS(I,KP)*CO21(I,KP,1) &
3240                             +DTC(I,KP)*EMISS(I,KP-1))*CLDFAC(I,KP,1)
3241     303   CONTINUE
3242     !          SUBROUTINE SPA88 IS CALLED TO OBTAIN EXACT CTS FOR WATER
3243     !     CO2 AND O3, AND APPROXIMATE CTS CO2 AND O3 CALCULATIONS.
3244     !
3245           CALL SPA88(EXCTS,CTSO3,GXCTS,SORC,CSOUR, &
3246                      CLDFAC,TEMP,PRESS,VAR1,VAR2, &
3247                      P,DELP,DELP2,TOTVO2,TO3SP,TO3SPC, &
3248                      CO2SP1,CO2SP2,CO2SP,              &
3249                      APCM,BPCM,ATPCM,BTPCM,ACOMB,BCOMB,BETACM,     &
3250                      H25E2,ONE,H44194M2,H1P41819,HAF,HM1EZ,TWO,    &
3251                      RADCON,                                 &
3252                      ids,ide, jds,jde, kds,kde,                    &
3253                      ims,ime, jms,jme, kms,kme,                    &
3254                      its,ite, jts,jte, kts,kte                     )
3255     
3256     !
3257     !    THIS SECTION COMPUTES THE EMISSIVITY CTS HEATING RATES FOR 2
3258     !    EMISSIVITY BANDS: THE 0-160,1200-2200 CM-1 BAND AND THE 800-
3259     !    990,1070-1200 CM-1 BAND. THE REMAINING CTS COMTRIBUTIONS ARE
3260     !    CONTAINED IN CTSO3, COMPUTED IN SPA88.
3261     !
3262           DO 998 I=MYIS,MYIE
3263           VTMP3(I,1)=1.
3264     998   CONTINUE
3265           DO 999 K=1,L
3266           DO 999 I=MYIS,MYIE
3267           VTMP3(I,K+1)=CNTTAU(I,K)*CLDFAC(I,K+1,1)
3268     999   CONTINUE
3269           DO 1001 K=1,L
3270           DO 1001 I=MYIS,MYIE
3271           CTS(I,K)=RADCON*DELP(I,K)*(TC(I,K)* &
3272                (E1CTW2(I,K)*CLDFAC(I,K+1,1)-E1CTW1(I,K)*CLDFAC(I,K,1)) + &
3273                 SS1(I,K)*(VTMP3(I,K+1)-VTMP3(I,K)))
3274     1001  CONTINUE
3275     !
3276           DO 1011 K=1,L
3277           DO 1011 I=MYIS,MYIE
3278           VTMP3(I,K)=TC(I,K)*(CLDFAC(I,K,1)*(E1CTS1(I,K)-E1CTW1(I,K)) - &
3279                             CLDFAC(I,K+1,1)*(E1CTS2(I,K)-E1CTW2(I,K)))
3280     1011  CONTINUE
3281           DO 1012 I=MYIS,MYIE
3282           FLX1E1(I)=TC(I,LP1)*CLDFAC(I,LP1,1)* &
3283                     (E1CTS1(I,LP1)-E1CTW1(I,LP1))
3284     1012  CONTINUE
3285           DO 1014 K=1,L
3286           DO 1013 I=MYIS,MYIE
3287           FLX1E1(I)=FLX1E1(I)+VTMP3(I,K)
3288     1013  CONTINUE
3289     1014  CONTINUE
3290     !
3291     !---NOW REPEAT FLUX CALCULATIONS FOR THE K=2..LM1  CASES.
3292     !   CALCULATIONS FOR FLUX LEVEL L AND LP1 ARE DONE SEPARATELY, AS ALL
3293     !   EMISSIVITY AND CO2 CALCULATIONS ARE SPECIAL CASES OR NEARBY LAYERS.
3294     !
3295           DO 321 K=2,LM1
3296           KLEN=K
3297     !
3298           DO 3218 KK=1,LP1-K
3299           DO 3218 I=MYIS,MYIE
3300           AVEPHI(I,KK+K-1)=TOTPHI(I,KK+K)-TOTPHI(I,K)
3301     3218  CONTINUE
3302           DO 1803 I=MYIS,MYIE
3303           AVEPHI(I,LP1)=AVEPHI(I,LM1)+EMX1(I)
3304     1803   CONTINUE
3305     !---COMPUTE EMISSIVITY FLUXES (E2) FOR THIS CASE. NOTE THAT
3306     !   WE HAVE OMITTED THE NEARBY LATER CASE (EMISS(I,K,K)) AS WELL
3307     !   AS ALL CASES WITH K=L OR LP1. BUT THESE CASES HAVE ALWAYS
3308     !   BEEN HANDLED AS SPECIAL CASES, SO WE MAY AS WELL COMPUTE
3309     !    THEIR FLUXES SEPARASTELY.
3310     !
3311           CALL E290(EMISSB,EMISS,AVEPHI,KLEN,FXOE2,DTE2,  &
3312                            H16E1,HP1,H28E1,HAF,TEN,       &
3313                            ids,ide, jds,jde, kds,kde,     &
3314                            ims,ime, jms,jme, kms,kme,     &
3315                            its,ite, jts,jte, kts,kte      )
3316     
3317           DO 322 KK=1,LP1-K
3318           DO 322 I=MYIS,MYIE
3319           AVMO3(I,KK+K-1)=TOTO3(I,KK+K)-TOTO3(I,K)
3320           AVPHO3(I,KK+K-1)=TPHIO3(I,KK+K)-TPHIO3(I,K)
3321           AVVO2(I,KK+K-1)=TOTVO2(I,KK+K)-TOTVO2(I,K)
3322           CONT1D(I,KK+K-1)=CNTTAU(I,KK+K-1)*TOTEVV(I,K-1)
3323     322   CONTINUE
3324     !
3325           DO 3221 KK=1,LP1-K
3326           DO 3221 I=MYIS,MYIE
3327           FAC1(I,K+KK-1)=BO3RND(2)*AVPHO3(I,K+KK-1)/AVMO3(I,K+KK-1)
3328           VTMP3(I,K+KK-1)=HAF*(FAC1(I,K+KK-1)* &
3329             (SQRT(ONE+(FOUR*AO3RND(2)*AVMO3(I,K+KK-1))/ &
3330              FAC1(I,K+KK-1))-ONE))
3331           TO31D(I,K+KK-1)=EXP(HM1EZ*(VTMP3(I,K+KK-1) &
3332                              +SKO3R*AVVO2(I,K+KK-1)))
3333           OVER1D(I,K+KK-1)=EXP(HM1EZ*(SQRT(AB15WD*AVEPHI(I,K+KK-1))+ &
3334                       SKC1R*AVVO2(I,K+KK-1)))
3335           CO21(I,K+KK,K)=OVER1D(I,K+KK-1)*CO21(I,K+KK,K)
3336     3221  CONTINUE
3337           DO 3223 KP=K+1,LP1
3338           DO 3223 I=MYIS,MYIE
3339           CO21(I,K,KP)=OVER1D(I,KP-1)*CO21(I,K,KP)
3340     3223  CONTINUE
3341     !---RLOG IS THE NBL AMOUNT FOR THE 15 UM BAND CALCULATION
3342           DO 1804 I=MYIS,MYIE
3343           RLOG(I,K)=OVER1D(I,K)*CO2NBL(I,K)
3344     1804  CONTINUE
3345     !---THE KP TERMS FOR ARBIRRARY K..
3346           DO 3423 KP=K+1,LP1
3347           DO 3423 I=MYIS,MYIE
3348     	if (I .eq. 50 .and. K .eq. 10) then
3349     !	write(0,*) 'I,L, FLX(I,K) (c): ', I,K, FLX(I,K)
3350     	endif
3351           FLX(I,K)=FLX(I,K)+(OSS(I,KP)*TO31D(I,KP-1) &
3352                             +SS2(I,KP)*CONT1D(I,KP-1) &
3353                             +CSS(I,KP)*CO21(I,KP,K) &
3354                             +DTC(I,KP)*EMISS(I,KP-1))*CLDFAC(I,KP,K)
3355     	if (I .eq. 50 .and. K .eq. 10) then
3356     !	write(0,*) 'I,K, FLX(I,K) (d): ', I,K, FLX(I,K)
3357     	endif
3358     3423  CONTINUE
3359           DO 3425 KP=K+1,LP1
3360           DO 3425 I=MYIS,MYIE
3361           FLX(I,KP)=FLX(I,KP)+(OSS(I,K)*TO31D(I,KP-1) &
3362                              +SS2(I,K)*CONT1D(I,KP-1) &
3363                              +CSS(I,K)*CO21(I,K,KP) &
3364                              +DTC(I,K)*EMISSB(I,KP-1))*CLDFAC(I,K,KP)
3365     3425  CONTINUE
3366     321   CONTINUE
3367     !
3368           DO 821 I=MYIS,MYIE
3369           TPL(I,1)=TEMP(I,L)
3370           TPL(I,LP1)=HAF*(T(I,LP1)+TEMP(I,L))
3371           TPL(I,LLP1)=HAF*(T(I,L)+TEMP(I,L))
3372     821   CONTINUE
3373           DO 823 K=2,L
3374           DO 823 I=MYIS,MYIE
3375           TPL(I,K)=T(I,K)
3376           TPL(I,K+L)=T(I,K)
3377     823   CONTINUE
3378     !
3379     !---E2 FUNCTIONS ARE REQUIRED IN THE NBL CALCULATIONS FOR 2 CASES,
3380     !   DENOTED (IN OLD CODE) AS (L,LP1) AND (LP1,LP1)
3381           DO 833 I=MYIS,MYIE
3382           AVEPHI(I,1)=VAR2(I,L)
3383           AVEPHI(I,2)=VAR2(I,L)+EMPL(I,L)
3384     833   CONTINUE
3385           CALL E2SPEC(EMISS,AVEPHI,FXOSP,DTSP,                          &
3386                           H16E1,TEN,H28E1,HP1,                          &
3387                           ids,ide, jds,jde, kds,kde,                    &
3388                           ims,ime, jms,jme, kms,kme,                    &
3389                           its,ite, jts,jte, kts,kte                     )
3390     
3391     !
3392     !     CALL E3V88 FOR NBL H2O TRANSMISSIVITIES
3393                CALL E3V88(EMD,TPL,EMPL, &
3394                           TEN,HP1,H28E1,H16E1,  &
3395                           ids,ide, jds,jde, kds,kde,                    &
3396                           ims,ime, jms,jme, kms,kme,                    &
3397                           its,ite, jts,jte, kts,kte                     )
3398     !
3399     !   COMPUTE NEARBY LAYER AND SPECIAL-CASE TRANSMISSIVITIES FOR EMISS
3400     !    USING METHODS FOR H2O GIVEN IN REF. (4)
3401           DO 851 K=2,L
3402           DO 851 I=MYIS,MYIE
3403     	if (I .eq. 50 .and. K .eq. 10) then
3404     !	write(0,*) 'EMD(I,K), EMD(I,K+L): ', EMD(I,K), EMD(I,K+L)
3405     	endif
3406           EMISDG(I,K)=EMD(I,K+L)+EMD(I,K)
3407     851   CONTINUE
3408     !
3409     !   NOTE THAT EMX1/2 (PRESSURE SCALED PATHS) ARE NOW COMPUTED IN
3410     !   LWR88
3411           DO 861 I=MYIS,MYIE
3412           EMSPEC(I,1)=(EMD(I,1)*EMPL(I,1)-EMD(I,LP1)*EMPL(I,LP1))/ &
3413            EMX1(I) + QUARTR*(EMISS(I,1)+EMISS(I,2))
3414           EMISDG(I,LP1)=TWO*EMD(I,LP1)
3415           EMSPEC(I,2)=TWO*(EMD(I,1)*EMPL(I,1)-EMD(I,LLP1)*EMPL(I,LLP1))/ &
3416            EMX2(I)
3417     861   CONTINUE
3418           DO 331 I=MYIS,MYIE
3419           FAC1(I,L)=BO3RND(2)*VAR4(I,L)/VAR3(I,L)
3420           VTMP3(I,L)=HAF*(FAC1(I,L)* &
3421               (SQRT(ONE+(FOUR*AO3RND(2)*VAR3(I,L))/FAC1(I,L))-ONE))
3422           TO31D(I,L)=EXP(HM1EZ*(VTMP3(I,L)+SKO3R*CNTVAL(I,L)))
3423           OVER1D(I,L)=EXP(HM1EZ*(SQRT(AB15WD*VAR2(I,L))+ &
3424                       SKC1R*CNTVAL(I,L)))
3425           CONT1D(I,L)=CNTTAU(I,L)*TOTEVV(I,LM1)
3426           RLOG(I,L)=OVER1D(I,L)*CO2NBL(I,L)
3427     331   CONTINUE
3428           DO 618 K=1,L
3429           DO 618 I=MYIS,MYIE
3430           RLOG(I,K)=LOG(RLOG(I,K))
3431     618   CONTINUE
3432           DO 601 K=1,LM1
3433           DO 601 I=MYIS,MYIE
3434           DELPR1(I,K+1)=DELP(I,K+1)*(PRESS(I,K+1)-P(I,K+1))
3435           ALP(I,LP1+K-1)=-SQRT(DELPR1(I,K+1))*RLOG(I,K+1)
3436     601   CONTINUE
3437           DO 603 K=1,L
3438           DO 603 I=MYIS,MYIE
3439           DELPR2(I,K+1)=DELP(I,K)*(P(I,K+1)-PRESS(I,K))
3440           ALP(I,K)=-SQRT(DELPR2(I,K+1))*RLOG(I,K)
3441     603   CONTINUE
3442           DO 625 I=MYIS,MYIE
3443           ALP(I,LL)=-RLOG(I,L)
3444           ALP(I,LLP1)=-RLOG(I,L)*SQRT(DELP(I,L)*(P(I,LP1)-PRESS(I,LM1)))
3445     625   CONTINUE
3446     !        THE FIRST COMPUTATION IS FOR THE 15 UM BAND,WITH THE
3447     !     FOR THE COMBINED H2O AND CO2 TRANSMISSION FUNCTION.
3448     !
3449     !       PERFORM NBL COMPUTATIONS FOR THE 15 UM BAND
3450     !***THE STATEMENT FUNCTION SF IN PREV. VERSIONS IS NOW EXPLICITLY
3451     !   EVALUATED.
3452           DO 631 K=1,LLP1
3453           DO 631 I=MYIS,MYIE
3454           C(I,K)=ALP(I,K)*(HMP66667+ALP(I,K)*(QUARTR+ALP(I,K)*HM6666M2))
3455     631   CONTINUE
3456           DO 641 I=MYIS,MYIE
3457           CO21(I,LP1,LP1)=ONE+C(I,L)
3458           CO21(I,LP1,L)=ONE+(DELP2(I,L)*C(I,LL)-(PRESS(I,L)-P(I,L))* &
3459            C(I,LLM1))/(P(I,LP1)-PRESS(I,L))
3460           CO21(I,L,LP1)=ONE+((P(I,LP1)-PRESS(I,LM1))*C(I,LLP1)- &
3461            (P(I,LP1)-PRESS(I,L))*C(I,L))/(PRESS(I,L)-PRESS(I,LM1))
3462     641   CONTINUE
3463           DO 643 K=2,L
3464           DO 643 I=MYIS,MYIE
3465           CO21(I,K,K)=ONE+HAF*(C(I,LM1+K)+C(I,K-1))
3466     643   CONTINUE
3467     !
3468     !    COMPUTE NEARBY-LAYER TRANSMISSIVITIES FOR THE O3 BAND AND FOR THE
3469     !    ONE-BAND CONTINUUM BAND (TO3 AND EMISS2). THE SF2 FUNCTION IS
3470     !    USED. THE METHOD IS THE SAME AS DESCRIBED FOR CO2 IN REF (4).
3471           DO 651 K=1,LM1
3472           DO 651 I=MYIS,MYIE
3473           CSUB(I,K+1)=CNTVAL(I,K+1)*DELPR1(I,K+1)
3474           CSUB(I,LP1+K-1)=CNTVAL(I,K)*DELPR2(I,K+1)
3475     651   CONTINUE
3476     !---THE SF2 FUNCTION IN PREV. VERSIONS IS NOW EXPLICITLY EVALUATED
3477           DO 655 K=1,LLM2
3478           DO 655 I=MYIS,MYIE
3479           CSUB2(I,K+1)=SKO3R*CSUB(I,K+1)
3480           C(I,K+1)=CSUB(I,K+1)*(HMP5+CSUB(I,K+1)* &
3481                     (HP166666-CSUB(I,K+1)*H41666M2))
3482           C2(I,K+1)=CSUB2(I,K+1)*(HMP5+CSUB2(I,K+1)* &
3483                      (HP166666-CSUB2(I,K+1)*H41666M2))
3484     655   CONTINUE
3485           DO 661 I=MYIS,MYIE
3486           CONTDG(I,LP1)=1.+C(I,LLM1)
3487           TO3DG(I,LP1)=1.+C2(I,LLM1)
3488     661   CONTINUE
3489           DO 663 K=2,L
3490           DO 663 I=MYIS,MYIE
3491           CONTDG(I,K)=ONE+HAF*(C(I,K)+C(I,LM1+K))
3492           TO3DG(I,K)=ONE+HAF*(C2(I,K)+C2(I,LM1+K))
3493     663   CONTINUE
3494     !---NOW OBTAIN FLUXES
3495     !
3496     !    FOR THE DIAGONAL TERMS...
3497           DO 871 K=2,LP1
3498           DO 871 I=MYIS,MYIE
3499           FLX(I,K)=FLX(I,K)+(DTC(I,K)*EMISDG(I,K) &
3500                            +SS2(I,K)*CONTDG(I,K) &
3501                            +OSS(I,K)*TO3DG(I,K) &
3502                            +CSS(I,K)*CO21(I,K,K))*CLDFAC(I,K,K)
3503     871   CONTINUE
3504     !     FOR THE TWO OFF-DIAGONAL TERMS...
3505           DO 873 I=MYIS,MYIE
3506           FLX(I,L)=FLX(I,L)+(CSS(I,LP1)*CO21(I,LP1,L) &
3507                             +DTC(I,LP1)*EMSPEC(I,2) &
3508                             +OSS(I,LP1)*TO31D(I,L) &
3509                             +SS2(I,LP1)*CONT1D(I,L))*CLDFAC(I,LP1,L)
3510           FLX(I,LP1)=FLX(I,LP1)+(CSS(I,L)*CO21(I,L,LP1) &
3511                                 +OSS(I,L)*TO31D(I,L) &
3512                                 +SS2(I,L)*CONT1D(I,L) &
3513                                 +DTC(I,L)*EMSPEC(I,1))*CLDFAC(I,L,LP1)
3514     873   CONTINUE
3515     !
3516     !     FINAL SECTION OBTAINS EMISSIVITY HEATING RATES,
3517     !     TOTAL HEATING RATES AND THE FLUX AT THE GROUND
3518     !
3519     !     .....CALCULATE THE EMISSIVITY HEATING RATES
3520           DO 1101 K=1,L
3521           DO 1101 I=MYIS,MYIE
3522           HEATEM(I,K)=RADCON*(FLX(I,K+1)-FLX(I,K))*DELP(I,K)
3523     
3524     1101  CONTINUE
3525     !     .....CALCULATE THE TOTAL HEATING RATES
3526           DO 1103 K=1,L
3527           DO 1103 I=MYIS,MYIE
3528           HEATRA(I,K)=HEATEM(I,K)-CTS(I,K)-CTSO3(I,K)+EXCTS(I,K)
3529     1103  CONTINUE
3530     !     .....CALCULATE THE FLUX AT EACH FLUX LEVEL USING THE FLUX AT THE
3531     !    TOP (FLX1E1+GXCTS) AND THE INTEGRAL OF THE HEATING RATES (VSUM1)
3532           DO 1111 K=1,L
3533           DO 1111 I=MYIS,MYIE
3534           VSUM1(I,K)=HEATRA(I,K)*DELP2(I,K)*RADCON1
3535     1111  CONTINUE
3536           DO 1115 I=MYIS,MYIE
3537           TOPFLX(I)=FLX1E1(I)+GXCTS(I)
3538           FLXNET(I,1)=TOPFLX(I)
3539     1115  CONTINUE
3540     !---ONLY THE SURFACE VALUE OF FLUX (GRNFLX) IS NEEDED UNLESS
3541     !    THE THICK CLOUD SECTION IS INVOKED.
3542           DO 1123 K=2,LP1
3543           DO 1123 I=MYIS,MYIE
3544           FLXNET(I,K)=FLXNET(I,K-1)+VSUM1(I,K-1)
3545     1123  CONTINUE
3546           DO 1125 I=MYIS,MYIE
3547           GRNFLX(I)=FLXNET(I,LP1)
3548     	if (I .eq. 50) then
3549     !	write(0,*) 'FLXNET(I,LP1), GRNFLX(I): ', FLXNET(I,LP1), GRNFLX(I)
3550     	endif
3551     1125  CONTINUE
3552     !
3553     !     THIS IS THE THICK CLOUD SECTION.OPTIONALLY,IF THICK CLOUD
3554     !     FLUXES ARE TO BE "CONVECTIVELY ADJUSTED",IE,DF/DP IS CONSTANT,
3555     !     FOR CLOUDY PART OF GRID POINT, THE FOLLOWING CODE IS EXECUTED.
3556     !***FIRST,COUNT THE NUMBER OF CLOUDS ALONG THE LAT. ROW. SKIP THE
3557     !   ENTIRE THICK CLOUD COMPUTATION OF THERE ARE NO CLOUDS.
3558           ICNT=0
3559           DO 1301 I=MYIS,MYIE
3560           ICNT=ICNT+NCLDS(I)
3561     1301  CONTINUE
3562           IF (ICNT.EQ.0) GO TO 6999
3563     !---FIND THE MAXIMUM NUMBER OF CLOUDS IN THE LATITUDE ROW
3564           KCLDS=NCLDS(MYIS)
3565           DO 2106 I=MYIS,MYIE
3566           KCLDS=MAX(NCLDS(I),KCLDS)
3567     2106  CONTINUE
3568     !
3569     !
3570     !***OBTAIN THE PRESSURES AND FLUXES OF THE TOP AND BOTTOM OF
3571     !   THE NC'TH CLOUD (IT IS ASSUMED THAT ALL KTOP AND KBTM'S HAVE
3572     !   BEEN DEFINED!).
3573           DO 1361 KK=1,KCLDS
3574           KMIN=LP1
3575           KMAX=0
3576           DO 1362 I=MYIS,MYIE
3577             J1=KTOP(I,KK+1)
3578     !       IF (J1.EQ.1) GO TO 1362
3579             J3=KBTM(I,KK+1)
3580             IF (J3.GT.J1) THEN
3581               PTOP(I)=P(I,J1)
3582               PBOT(I)=P(I,J3+1)
3583               FTOP(I)=FLXNET(I,J1)
3584               FBOT(I)=FLXNET(I,J3+1)
3585     !***OBTAIN THE "FLUX DERIVATIVE" DF/DP (DELPTC)
3586               DELPTC(I)=(FTOP(I)-FBOT(I))/(PTOP(I)-PBOT(I))
3587               KMIN=MIN(KMIN,J1)
3588               KMAX=MAX(KMAX,J3)
3589             ENDIF
3590     1362  CONTINUE
3591           KMIN=KMIN+1
3592     !***CALCULATE THE TOT. FLUX CHG. FROM THE TOP OF THE CLOUD, FOR
3593     !   ALL LEVELS.
3594           DO 1365 K=KMIN,KMAX
3595           DO 1363 I=MYIS,MYIE
3596             IF(KTOP(I,KK+1).LT.K .AND. K.LE.KBTM(I,KK+1)) THEN
3597               Z1(I,K)=(P(I,K)-PTOP(I))*DELPTC(I)+FTOP(I)
3598               FLXNET(I,K)=Z1(I,K)
3599             ENDIF
3600     1363  CONTINUE
3601     1365  CONTINUE
3602     1361  CONTINUE
3603     !***USING THIS FLUX CHG. IN THE CLOUDY PART OF THE GRID BOX, OBTAIN
3604     !   THE NEW FLUXES, WEIGHTING THE CLEAR AND CLOUDY FLUXES:AGAIN, ONLY
3605     !    THE FLUXES IN THICK-CLOUD LEVELS WILL EVENTUALLY BE USED.
3606     !     DO 6051 K=1,LP1
3607     !     DO 6051 I=MYIS,MYIE
3608     !     FLXNET(I,K)=FLXNET(I,K)*(ONE-CAMT(I,NC)) +
3609     !    1            Z1(I,K)*CAMT(I,NC)
3610     !051  CONTINUE
3611     !***MERGE FLXTHK INTO FLXNET FOR APPROPRIATE LEVELS.
3612     !     DO 1401 K=1,LP1
3613     !     DO 1401 I=MYIS,MYIE
3614     !     IF (K.GT.ITOP(I) .AND. K.LE.IBOT(I)
3615     !    1  .AND.  (NC-1).LE.NCLDS(I))  THEN
3616     !          FLXNET(I,K)=FLXTHK(I,K)
3617     !     ENDIF
3618     !401  CONTINUE
3619     !
3620     !******END OF CLOUD LOOP*****
3621     6001  CONTINUE
3622     6999  CONTINUE
3623     !***THE FINAL STEP IS TO RECOMPUTE THE HEATING RATES BASED ON THE
3624     !   REVISED FLUXES:
3625           DO 6101 K=1,L
3626           DO 6101 I=MYIS,MYIE
3627           HEATRA(I,K)=RADCON*(FLXNET(I,K+1)-FLXNET(I,K))*DELP(I,K)
3628     6101  CONTINUE
3629     !     THE THICK CLOUD SECTION ENDS HERE.
3630     
3631       END SUBROUTINE FST88
3632     
3633     !----------------------------------------------------------------------
3634     
3635       SUBROUTINE E1E290(G1,G2,G3,G4,G5,EMISS,FXOE1,DTE1,FXOE2,DTE2,      &
3636                            AVEPHI,TEMP,T,                                &
3637                            H16E1,TEN,HP1,H28E1,HAF,                      &
3638                            ids,ide, jds,jde, kds,kde,                    &
3639                            ims,ime, jms,jme, kms,kme,                    &
3640                            its,ite, jts,jte, kts,kte                     )
3641     !---------------------------------------------------------------------
3642      IMPLICIT NONE
3643     !----------------------------------------------------------------------
3644           INTEGER, INTENT(IN)        :: ids,ide, jds,jde, kds,kde ,      &
3645                                         ims,ime, jms,jme, kms,kme ,      &
3646                                         its,ite, jts,jte, kts,kte
3647           REAL,INTENT(IN) :: H16E1,TEN,HP1,H28E1,HAF
3648     
3649           REAL,INTENT(OUT),DIMENSION(its:ite,kts:kte+1) :: G1,G4,G3,EMISS
3650           REAL,INTENT(IN),DIMENSION(its:ite,kts:kte+1) :: FXOE1,DTE1,FXOE2,DTE2
3651           REAL,INTENT(IN),DIMENSION(its:ite,kts:kte+1) :: AVEPHI,TEMP,T
3652           REAL,INTENT(OUT),DIMENSION(its:ite,kts:kte)   :: G2,G5
3653     
3654           REAL,DIMENSION(its:ite,kts:kte+1) :: TMP3,DU,FYO,WW1,WW2
3655           INTEGER,DIMENSION(its:ite,kts:kte*3+2)   :: IT1
3656           INTEGER,DIMENSION(its:ite,kts:kte+1) :: IVAL
3657     
3658           INTEGER :: K, I,KP,LLM2,J1,J3,KMAX,KMIN,KCLDS,ICNT,LLM1
3659           INTEGER :: L,LP1,LP2,LP3,LM1,LM2,LM3,MYIS,MYIE,LLP1,LL,KK,KLEN
3660     
3661           L=kte
3662           LP1=L+1;  LP2=L+2;  LP3=L+3; LLP1 = 2*L + 1
3663           LM1=L-1;  LM2=L-2;  LM3=L-3; LL = 2*L
3664           LLM2 = LL-2; LLM1=LL-1
3665           MYIS=its; MYIE=ite
3666     
3667     !---FIRST WE OBTAIN THE EMISSIVITIES AS A FUNCTION OF TEMPERATURE
3668     !   (INDEX FXO) AND WATER AMOUNT (INDEX FYO). THIS PART OF THE CODE
3669     !   THUS GENERATES THE E2 FUNCTION. THE FXO INDICES HAVE BEEN
3670     !   OBTAINED IN FST88, FOR CONVENIENCE.
3671     !
3672     !---THIS SUBROUTINE EVALUATES THE K=1 CASE ONLY--
3673     !
3674     !---THIS LOOP REPLACES LOOPS GOING FROMI=1,IMAX AND KP=2,LP1 PLUS
3675     !   THE SPECIAL CASE FOR THE LP1TH LAYER.
3676     
3677           DO 1322 K=1,LP1
3678           DO 1322 I=MYIS,MYIE
3679           TMP3(I,K)=LOG10(AVEPHI(I,K))+H16E1
3680           FYO(I,K)=AINT(TMP3(I,K)*TEN)
3681           DU(I,K)=TMP3(I,K)-HP1*FYO(I,K)
3682           FYO(I,K)=H28E1*FYO(I,K)
3683           IVAL(I,K)=FYO(I,K)+FXOE2(I,K)
3684           EMISS(I,K)=T1(IVAL(I,K))+DU(I,K)*T2(IVAL(I,K)) &
3685                                   +DTE2(I,K)*T4(IVAL(I,K))
3686     1322  CONTINUE
3687     !
3688     !---THE SPECIAL CASE EMISS(I,L) (LAYER KP) IS OBTAINED NOW
3689     !   BY AVERAGING THE VALUES FOR L AND LP1:
3690           DO 1344 I=MYIS,MYIE
3691           EMISS(I,L)=HAF*(EMISS(I,L)+EMISS(I,LP1))
3692     1344  CONTINUE
3693     !
3694     !   CALCULATIONS FOR THE KP=1 LAYER ARE NOT PERFORMED, AS
3695     !   THE RADIATION CODE ASSUMES THAT THE TOP FLUX LAYER (ABOVE THE
3696     !   TOP DATA LEVEL) IS ISOTHERMAL, AND HENCE CONTRIBUTES NOTHING
3697     !   TO THE FLUXES AT OTHER LEVELS.
3698     !
3699     !***THE FOLLOWING IS THE CALCULATION FOR THE E1 FUNCTION, FORMERLY
3700     !    DONE IN SUBROUTINE E1V88. THE MOVE TO E1E288 IS DUE TO THE
3701     !    SAVINGS IN OBTAINING INDEX VALUES (THE TEMP. INDICES HAVE
3702     !    BEEN OBTAINED IN FST88, WHILE THE U-INDICES ARE OBTAINED
3703     !    IN THE E2 CALCS.,WITH K=1).
3704     !
3705     !
3706     !   FOR TERMS INVOLVING TOP LAYER, DU IS NOT KNOWN; IN FACT, WE
3707     !   USE INDEX 2 TO REPERSENT INDEX 1 IN PREV. CODE. THIS MEANS THAT
3708     !    THE IT1 INDEX 1 AND LLP1 HAS TO BE CALCULATED SEPARATELY. THE
3709     !   INDEX LLP2 GIVES THE SAME VALUE AS 1; IT CAN BE OMITTED.
3710           DO 208 I=MYIS,MYIE
3711           IT1(I,1)=FXOE1(I,1)
3712           WW1(I,1)=TEN-DTE1(I,1)
3713           WW2(I,1)=HP1
3714     208   CONTINUE
3715           DO 209 K=1,L
3716           DO 209 I=MYIS,MYIE
3717           IT1(I,K+1)=FYO(I,K)+FXOE1(I,K+1)
3718           IT1(I,LP2+K-1)=FYO(I,K)+FXOE1(I,K)
3719           WW1(I,K+1)=TEN-DTE1(I,K+1)
3720           WW2(I,K+1)=HP1-DU(I,K)
3721     209   CONTINUE
3722           DO 211 KP=1,L
3723           DO 211 I=MYIS,MYIE
3724           IT1(I,KP+LLP1)=FYO(I,KP)+FXOE1(I,1)
3725     211   CONTINUE
3726     !
3727     !
3728     !  G3(I,1) HAS THE SAME VALUES AS G1 (AND DID ALL ALONG)
3729           DO 230 I=MYIS,MYIE
3730           G1(I,1)=WW1(I,1)*WW2(I,1)*EM1V(IT1(I,1))+ &
3731                   WW2(I,1)*DTE1(I,1)*EM1V(IT1(I,1)+1)
3732           G3(I,1)=G1(I,1)
3733     230   CONTINUE
3734           DO 240 K=1,L
3735           DO 240 I=MYIS,MYIE
3736           G1(I,K+1)=WW1(I,K+1)*WW2(I,K+1)*EM1V(IT1(I,K+1))+ &
3737                   WW2(I,K+1)*DTE1(I,K+1)*EM1V(IT1(I,K+1)+1)+ &
3738                   WW1(I,K+1)*DU(I,K)*EM1V(IT1(I,K+1)+28)+ &
3739                   DTE1(I,K+1)*DU(I,K)*EM1V(IT1(I,K+1)+29)
3740           G2(I,K)=WW1(I,K)*WW2(I,K+1)*EM1V(IT1(I,K+LP2-1))+ &
3741                   WW2(I,K+1)*DTE1(I,K)*EM1V(IT1(I,K+LP2-1)+1)+ &
3742                   WW1(I,K)*DU(I,K)*EM1V(IT1(I,K+LP2-1)+28)+ &
3743                   DTE1(I,K)*DU(I,K)*EM1V(IT1(I,K+LP2-1)+29)
3744     240   CONTINUE
3745           DO 241 KP=2,LP1
3746           DO 241 I=MYIS,MYIE
3747           G3(I,KP)=WW1(I,1)*WW2(I,KP)*EM1V(IT1(I,LL+KP))+ &
3748                   WW2(I,KP)*DTE1(I,1)*EM1V(IT1(I,LL+KP)+1)+ &
3749                   WW1(I,1)*DU(I,KP-1)*EM1V(IT1(I,LL+KP)+28)+ &
3750                   DTE1(I,1)*DU(I,KP-1)*EM1V(IT1(I,LL+KP)+29)
3751     241   CONTINUE
3752     !
3753           DO 244 I=MYIS,MYIE
3754           G4(I,1)=WW1(I,1)*WW2(I,1)*EM1VW(IT1(I,1))+ &
3755                   WW2(I,1)*DTE1(I,1)*EM1VW(IT1(I,1)+1)
3756     244   CONTINUE
3757           DO 242 K=1,L
3758           DO 242 I=MYIS,MYIE
3759           G4(I,K+1)=WW1(I,K+1)*WW2(I,K+1)*EM1VW(IT1(I,K+1))+ &
3760                   WW2(I,K+1)*DTE1(I,K+1)*EM1VW(IT1(I,K+1)+1)+ &
3761                   WW1(I,K+1)*DU(I,K)*EM1VW(IT1(I,K+1)+28)+ &
3762                   DTE1(I,K+1)*DU(I,K)*EM1VW(IT1(I,K+1)+29)
3763           G5(I,K)=WW1(I,K)*WW2(I,K+1)*EM1VW(IT1(I,K+LP2-1))+ &
3764                   WW2(I,K+1)*DTE1(I,K)*EM1VW(IT1(I,K+LP2-1)+1)+ &
3765                   WW1(I,K)*DU(I,K)*EM1VW(IT1(I,K+LP2-1)+28)+ &
3766                   DTE1(I,K)*DU(I,K)*EM1VW(IT1(I,K+LP2-1)+29)
3767     242   CONTINUE
3768     !
3769       END SUBROUTINE E1E290
3770     
3771     !----------------------------------------------------------------------
3772     
3773      SUBROUTINE SPA88(EXCTS,CTSO3,GXCTS,SORC,CSOUR,                      &
3774                            CLDFAC,TEMP,PRESS,VAR1,VAR2,                  &
3775                            P,DELP,DELP2,TOTVO2,TO3SP,TO3SPC,             &
3776                            CO2SP1,CO2SP2,CO2SP,                          &
3777                            APCM,BPCM,ATPCM,BTPCM,ACOMB,BCOMB,BETACM,     &
3778                            H25E2,ONE,H44194M2,H1P41819,HAF,HM1EZ,TWO,    &
3779                            RADCON,                                 &
3780                            ids,ide, jds,jde, kds,kde,                    &
3781                            ims,ime, jms,jme, kms,kme,                    &
3782                            its,ite, jts,jte, kts,kte                     )
3783     !---------------------------------------------------------------------
3784      IMPLICIT NONE
3785     !----------------------------------------------------------------------
3786           INTEGER, INTENT(IN)        :: ids,ide, jds,jde, kds,kde ,      &
3787                                         ims,ime, jms,jme, kms,kme ,      &
3788                                         its,ite, jts,jte, kts,kte
3789     
3790           REAL,INTENT(IN) :: H25E2,ONE,H44194M2,H1P41819,HAF,HM1EZ,TWO, &
3791                              RADCON
3792     
3793           REAL,INTENT(IN),DIMENSION(its:ite,kts:kte+1) :: CSOUR
3794           REAL,INTENT(OUT),DIMENSION(its:ite,kts:kte)  :: CTSO3
3795           REAL,INTENT(OUT),DIMENSION(its:ite,kts:kte)  :: EXCTS
3796           REAL,INTENT(OUT),DIMENSION(its:ite)          :: GXCTS
3797           REAL,INTENT(IN),DIMENSION(its:ite,kts:kte+1,NBLY) :: SORC
3798           REAL,INTENT(IN),DIMENSION(its:ite,kts:kte+1,kts:kte+1) :: CLDFAC
3799           REAL,INTENT(IN), DIMENSION(its:ite,kts:kte+1) :: PRESS,TEMP
3800     
3801           REAL,INTENT(IN),DIMENSION(its:ite,kts:kte) :: VAR1,VAR2 
3802           REAL,INTENT(IN),DIMENSION(its:ite,kts:kte+1) :: P
3803           REAL,INTENT(IN),DIMENSION(its:ite,kts:kte)   :: DELP,DELP2,TO3SPC
3804           REAL,INTENT(IN),DIMENSION(its:ite,kts:kte+1) ::TOTVO2,TO3SP,CO2SP1,&
3805                                                          CO2SP2,CO2SP
3806           REAL,INTENT(IN),DIMENSION(NBLY) :: APCM,BPCM,ATPCM,BTPCM,ACOMB, &
3807                                              BCOMB,BETACM
3808     
3809           REAL,DIMENSION(its:ite,kts:kte+1) ::CTMP,CTMP2,CTMP3
3810           REAL,DIMENSION(its:ite,kts:kte)   ::X,Y,FAC1,FAC2,F,FF,AG,AGG, &
3811                                               PHITMP,PSITMP,TOPM,TOPPHI,TT
3812     
3813           INTEGER :: K, I,KP,LLM2,J1,J3,KMAX,KMIN,KCLDS,ICNT,LLM1
3814           INTEGER :: L,LP1,LP2,LP3,LM1,LM2,LM3,MYIS,MYIE,LLP1,LL,KK,KLEN
3815     
3816           L=kte
3817           LP1=L+1;  LP2=L+2;  LP3=L+3; LLP1 = 2*L + 1
3818