      SUBROUTINE CALHEAT(ISTL_)  
C
C   Subroutine CALHEAT takes the information from the atmospheric boundary
C   file and the wind forcing file and calculates the net heat flux across
C   the water surface boundary.  The heat flux is then used to update the
C   water temperature either in the surface cells, or distributed across
C   the cells in the vertical and into the bottom.  The subroutine has
C   three options these are:
C
C    ISOPT(2)=1: Full surface and internal heat transfer calculation
C                using meteorologic data from input stream.
C                IASWRAD=0:  ADSORB SW SOLR RAD TO ALL LAYERS AND BED
C                IASWRAD=1:  ADSORB SW SOLR RAD TO TO SURFACE LAYER
C    ISOPT(2)=2: Transient equilibrium surface heat transfer calculation
C                using external equilibrium temperature and heat transfer
C                coefficient data from the meteorologic input data.
C    ISOPT(2)=3: Equilibrium surface heat transfer calculation using constant
C                equilibrium temperature and heat transfer coefficients
C                HEQT and TEMO read in through input stream.
C    ISOPT(2)=4: Equilibrium surface heat transfer calculation using algorithm
C                from CE-QUAL-W2.
C
C   The heat flux terms are derived from a paper by Rosati
C   and Miyakoda (1988) entitled "A General Circulation Model for Upper Ocean
C   Simulation".  The heat flux is prescribed by term for the following
C   influxes and outfluxes:
C
C     - Short Wave Incoming Radiation (+)
C     - Net Long Wave Radiation (+/-)
C     - Sensible Heat Flux (convection -)
C     - Latent Heat Flux (evaporation +/-)
C
C   Two formulations of the Latent Heat Flux are provided.  The first is from
C   the Rosati and Miyakoda paper, the second is an alternate formulation by
C   Brady, Graves, and Geyer (1969).  The second formulation was taken from
C   "Hydrodynamics and Transport for Water Quality Modeling" (Martin and
C   McCutcheon, 1999).  The Rosati and Miyakoda formulation will have zero
C   evaporative cooling or heating if wind speed goes to zero.  The Brady,
C   Graves, and Geyer formulation provides for a minimum evaporative cooling
C   under zero wind speed.
C
C
C VARIABLE LIST:
C
C   CLOUDT  = Cloud cover (0 to 10)
C   HCON    = Sensible heat flux (W/m2)
C   HLAT    = Latent heat flux (W/m2)
C   HLWBR   = Net longwave radiation (atmospheric long wave plus back
C             radiation, W/m2)
C   SOLSWRT = Short wave incoming radiation (W/m2)
C   SVPW    = Saturation vapor pressure in mb based upon the water surface
C             temperature
C   TATMT   = Temperature of air above water surface (deg C)
C   TEM     = Water temperature in cell (deg C)
C   VPA     = Vapor pressure of air at near surface air temperature (mb)
C   WINDST  = Wind speed at 10 meters over cell surface (m/s)
C
C MODIFICATION HISTORY:
C
C   Date       Author             Comments
C   ---------- ------------------ ---------------------------------------------
C   06/01/1992 John M. Hamrick    Orignial author
C   11/07/2000 Steven Peene       Cleaned code, provided more detailed
C                                 descriptions, added alternate formulation
C                                 for latent heat flux, separated out
C                                 individual heat flux terms
C   11/01/2005 Paul M. Craig      Added Option 4, to use the Equilibrium Temperature 
C                                 algorithym from CE-QUAL-W2.  Also added the sub-option
C                                 under this option to couple or decouple the bottom temperature
C                                 to the water column temperatures.
C                                 Added the ability to input spatially variable bed temps and
C                                 thermally active bed thicknesses.
C                                 Also cleaned up the code and added structure.
C
C CHANGE RECORD  
C **  SUBROUTINE CALHEAT CALCULATES SURFACE AND INTERNAL HEAT SOURCES  
C **  AND SINKS IN THE HEAT (TEM) TRANSPORT EQUATION  
C  
      USE GLOBAL  
      USE OMP_LIB
      REAL,SAVE,ALLOCATABLE,DIMENSION(:,:)::NETRAD  
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::TBEDTHK
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::HDEP
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::RADBOT
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::FLUXTB
      !REAL,SAVE ::    PTIME
      !REAL,SAVE ::    PMCTOL
      REAL            K_ABOVE
C
      IF(.NOT.ALLOCATED(NETRAD))THEN
        ALLOCATE(NETRAD(LCM,KCM))
        ALLOCATE(TBEDTHK(LCM))
        ALLOCATE(HDEP(LCM))
        ALLOCATE(RADBOT(LCM))
        ALLOCATE(FLUXTB(LCM))
	  RADBOT=0.0 !SCJ
	  FLUXTB=0.0 !SCJ

        ! *** Ininitialze Heat Exchange Terms
        IF (PARTID ==0) PRINT *,'CALHEAT: INITIALIZING'
        CALL HEAT_EXCHANGE

        !PMCTOL=0.1
        NETRAD=0.
        HDEP=0.
        TBEDTHK=0.
        IF((ISTOPT(2).EQ.1.AND.IASWRAD.EQ.0).OR.ISTOPT(2).EQ.4)THEN 
          IF(DABEDT.GT.0.)THEN
            IF(PARTID==MASTER_TASK) PRINT *,'CALHEAT: SET CONSTANT THICKNESS TO:',DABEDT
            DO L=2,LA
              TBEDTHK(L)=DABEDT
            ENDDO
          ELSE
            ! *** READ IN THE SPATIALLY VARYING INIT T AND BED THICKNESS (DABEDT) 
            IF (PARTID==MASTER_TASK) PRINT *,'CALHEAT: READ IN SPATIALLY VARYING INIT T AND
     & BED THICKNESS: TEMB.INP'
            DO L=2,LA
              TBEDTHK(L)=ABS(DABEDT)
              IF(ISCI(2).EQ.0)TEMB(L)=ABS(TBEDIT)
            ENDDO
            OPEN(1001,FILE='TEMB.INP',ERR=1000,STATUS='OLD')  
            DO IS=1,4  
              READ(1001,*)  
            ENDDO  
            DO L1 = 2, LA
              READ(1001,*,END=1000) I,J,T1,T2
              L=LIJ(I,J)
              IF(ISCI(2).EQ.0)TEMB(L)=T1
              TBEDTHK(L)=T2
            ENDDO  
 1000       CLOSE(1001)  
          ENDIF
        ENDIF

        IF(ISTOPT(2).EQ.1.OR.ISTOPT(2).EQ.4)THEN 
          DO L=2,LA
            PSHADE(L)=1.0
          ENDDO
          IF(USESHADE)THEN
            ! *** READ IN THE SPATIALLY VARYING INIT T AND BED THICKNESS (DABEDT) 
            IF (PARTID == 0) PRINT *,
     *        'CALHEAT: READ IN SPATIALLY VARYING SHADE: PSHADE.INP'
            OPEN(1001,FILE='PSHADE.INP',ERR=1010,STATUS='OLD')  
            DO IS=1,4  
              READ(1001,*)  
            ENDDO  
            DO L1 = 2, LA
              READ(1001,*,END=1010) I,J,T1
              L=LIJ(I,J)
              PSHADE(L)=T1
            ENDDO  
 1010       CLOSE(1001)  
          ELSE
            IF (PARTID ==0) PRINT *,
     &      'CALHEAT: SET CONST SHADE TO: 1.0 (NO SHADE)'
          ENDIF
        ENDIF

        IF(DEBUG)THEN
          PRINT *,'CALHEAT: Bed Temp(L=2):', TEMB(2)
          OPEN(77,file='calheat.dia',status='unknown')
          CLOSE(77,status='DELETE')
          OPEN(77,file='calheat.dia',status='NEW')
          WRITE(77,998)'TIMEDAY','SRON','ET','TD_C','TA_C','TDEW_F',
     &                  'TAIR_F','FW'
  998     FORMAT(A11,8A9)
        ENDIF
      ENDIF
C  
      ! *** DSLLC BEGIN CHANGE
CPMC      DELT=DT2  
CPMC      IF(ISTL_.EQ.2)THEN  
CPMC        DELT=DT  
CPMC      ENDIF  
CPMC      DELT=DT2  
      IF(ISTL_.EQ.2)THEN  
        IF(ISDYNSTP.EQ.0)THEN  
          DELT=DT  
        ELSE  
          DELT=DTDYN  
        END IF  
      ELSE
        DELT=DT2 
      ENDIF  
      ! *** DSLLC END CHANGE

      ! *** OVERWRITE THE INPUT SOLAR RAD WITH A COMPUTED ONE
      IF(COMPUTESOLRAD)THEN
        CALL SHORT_WAVE_RADIATION(WINDST(2),RHA(2),TATMT(2),CLOUDT(2),
     &                            PATMT(2),SRO,SRON)
        DO L=2,LA
          ! *** USE COMPUTED SRO
          SOLSWRT(L)=SRON
        ENDDO
      ENDIF

      IF(USESHADE)THEN
        DO L=2,LA
          ! *** APPLY PSHADE FACTORS
          SOLSWRT(L)=SOLSWRT(L)*PSHADE(L)
        ENDDO
      ENDIF

      IF(ISTOPT(2).EQ.1)THEN  
        ! *** FULL HEAT BALANCE WITH ATMOSPHERIC LINKAGE

        ! *** SET UP MIN DEPTH
        DO L=2,LA  
          HDEP(L)=MAX(HP(L),0.)  
        ENDDO


        ! *** NET SHORTWAVE SOLAR RADIATION
        IF(IASWRAD.EQ.0.)THEN  
          ! *** ADSORB SW SOLR RAD TO ALL LAYERS AND BED  

        ! NET HEAT FLUX = RSN+RAN-RB-RE-RC  (WATT/M2)
!$OMP PARALLEL PRIVATE(SVPW,CLDFAC,RAN,FW,RE,RC,RB,RSN,
!$OMP& UBED,VBED,USPD,TMPVAL,TFAST,TFAST1,C1,TSLOW,TSLOW1)

          ! *** SURFACE LAYER
          TFAST=SWRATNF*(Z(KC)-1.)
          TFAST1=SWRATNF*(Z(KC-1)-1.)
          TSLOW=SWRATNS*(Z(KC)-1.)
          TSLOW1=SWRATNS*(Z(KC-1)-1.)
        IF (NASER == -1 ) THEN ! USE DeepThunder provided values of longwave radiation
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
          DO L=2,LA  
            SVPW=(10.**((0.7859+0.03477*TEM(L,KC))/  
     &          (1.+0.00412*TEM(L,KC))))  
           
           ! *** Evaporation
            FW=9.2+0.46*WINDST(L)**2;
           !Saturated water vapor pres - actual water vapor press (eqn 2.3.10 in Ji 2008)  
            RE=FW*(SVPW-VPA(L)); 
            ! *** Conduction
            RC=0.47*FW*(TEM(L,KC)-TATMT(L)) 
            NETRAD(L,KC)=LWDNB(L) - LWUPB(L) -RE-RC
          ENDDO
!$OMP END DO NOWAIT
        ELSE   ! use existing implementation based on Swinbank (1963) formula
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
          DO L=2,LA  
            SVPW=(10.**((0.7859+0.03477*TEM(L,KC))/  
     &          (1.+0.00412*TEM(L,KC))))  

            ! *** Net atmospheric radiation (Diffusive)
            CLDFAC=1.0+0.0017*CLOUDT(L)**2
            RAN=0.51531538e-12*(TATMT(L)+273.15)**6*CLDFAC
            ! *** Evaporation
            FW=9.2+0.46*WINDST(L)**2;
           !Saturated water vapor pres - actual water vapor press (eqn 2.3.10 in Ji 2008)  
            RE=FW*(SVPW-VPA(L)); 
            ! *** Conduction
            RC=0.47*FW*(TEM(L,KC)-TATMT(L)) 
            ! *** Longwave back radiation
            ! ***  5.443E-8 = 5.67E-8 * 0.97
            RB=5.443E-8*(TEM(L,KC)+273.15)**4 ! upward radiation flux
            NETRAD(L,KC)=RAN-RB-RE-RC
          ENDDO
!$OMP END DO NOWAIT
        END IF        
          IF(FSWRATF.LT.1.)THEN
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA  
              RSN=SOLSWRT(L)*  
     &           (    FSWRATF*(EXP(TFAST*HDEP(L))-EXP(TFAST1*HDEP(L)))  
     &          +(1.-FSWRATF)*(EXP(TSLOW*HDEP(L))-EXP(TSLOW1*HDEP(L))))    
              NETRAD(L,KC)=NETRAD(L,KC)+RSN
            ENDDO
!$OMP END DO NOWAIT
          ELSE
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA  
              RSN=SOLSWRT(L)*(1.-EXP(TFAST1*HDEP(L)))  
              NETRAD(L,KC)=NETRAD(L,KC)+RSN
            ENDDO
!$OMP END DO NOWAIT
          ENDIF

          ! *** ALL REMAINING LAYERS
          IF(KC.GT.1)THEN  
            DO K=1,KS  
              TFAST=SWRATNF*(Z(K)-1.)
              TFAST1=SWRATNF*(Z(K-1)-1.)
              C2=DELT*DZIC(K)*0.2393E-6
              IF(FSWRATF.LT.1.)THEN
                TSLOW=SWRATNS*(Z(K)-1.)
                TSLOW1=SWRATNS*(Z(K-1)-1.)
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
                DO L=2,LA  
                  RSN=SOLSWRT(L)*   
     &             (    FSWRATF*(EXP(TFAST*HDEP(L))-EXP(TFAST1*HDEP(L)))  
     &           +(1.-FSWRATF)*(EXP(TSLOW*HDEP(L))-EXP(TSLOW1*HDEP(L))))    
                  NETRAD(L,K)=RSN
                ENDDO
!$OMP END DO NOWAIT  
              ELSE
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
                DO L=2,LA  
                  RSN=SOLSWRT(L)*
     &               (EXP(TFAST*HDEP(L))-EXP(TFAST1*HDEP(L)))
                  NETRAD(L,K)=RSN
                ENDDO
!$OMP END DO NOWAIT  
              ENDIF
            ENDDO  
          ENDIF  

          ! *** Distribute heat flux to the bed for each grid cell.
          TFAST=SWRATNF*(Z(0)-1.)
          IF(FSWRATF.LT.1.)THEN
            TSLOW=SWRATNS*(Z(0)-1.)
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA
              LE=LEAST(L)
              LN=LNC(L)
              UBED=0.5*( U(L,1)+U(LE,1) )  
              VBED=0.5*( V(L,1)+V(LN,1) )  
              USPD=SQRT( UBED*UBED+VBED*VBED )  
              TMPVAL=(HTBED1*USPD+HTBED2)*(TEM(L,1)-TEMB(L))  
              NETRAD(L,1)=NETRAD(L,1)-TMPVAL/0.2393E-6

              ! *** UPDATE BOTTOM
              IF(TBEDIT.GT.0.)THEN
                TEMB(L)=TEMB(L) + (TMPVAL + 0.2393E-6*SOLSWRT(L)*
     &               (FSWRATF *EXP(TFAST*HDEP(L))
     &           +(1.-FSWRATF)*EXP(TSLOW*HDEP(L))))*DELT/TBEDTHK(L)  
              ENDIF
            ENDDO  
!$OMP END DO NOWAIT
          ELSE
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA
              LE=LEAST(L)
              LN=LNC(L)
              UBED=0.5*( U(L,1)+U(LE,1) )  
              VBED=0.5*( V(L,1)+V(LN,1) )  
              USPD=SQRT( UBED*UBED+VBED*VBED )  
              TMPVAL=(HTBED1*USPD+HTBED2)*(TEM(L,1)-TEMB(L))  
              NETRAD(L,1)=NETRAD(L,1)-TMPVAL/0.2393E-6

              ! *** UPDATE BOTTOM
              IF(TBEDIT.GT.0.)THEN
                TEMB(L)=TEMB(L) + (TMPVAL + 0.2393E-6*SOLSWRT(L)*
     &                  FSWRATF*EXP(TFAST*HDEP(L)))*DELT/TBEDTHK(L)  
              ENDIF
            ENDDO  
!$OMP END DO NOWAIT
          ENDIF

          ! *** NOW FINALIZE THE TEMPERATURE 
          DO K=1,KC  
            ! *** RHO = 1000.0  Density (kg / m^3)
            ! *** CP  = 4179.0  Specific Heat (J / kg / degC)
            ! *** 0.2393E-6 = 1/RHO/CP
            C1=DELT*DZIC(K)*0.2393E-6
!$OMP DO SCHEDULE(STATIC,CHUNKSIZE)
            DO L=2,LA  
              TEM(L,K)=TEM(L,K)+HPI(L)*C1*NETRAD(L,K)
              IF (TEM(L,K) < 0.85) TEM(L,K) = 0.85   ! include this hard code until development of full ice-model  
            ENDDO  
!$OMP END DO NOWAIT
          ENDDO  
!$OMP END PARALLEL

          IF(ISDRY.GT.0.AND.ISTOPT(2).EQ.1)THEN  
            DO L=2,LA  
              IF(IMASKDRY(L).EQ.1.) TEMB(L)=TATMT(L)  
            ENDDO  
          ENDIF  
  
        ELSE   ! IF(IASWRAD.EQ.1)THEN  

          ! NET HEAT FLUX = RSN+RAN-RB-RE-RC  (WATT/M2)
          C1=DELT*DZIC(KC)*0.2393E-6
          IF (NASER == -1) THEN  ! Deep Thunder computed longwave radiation
            DO L=2,LA  
              SVPW=(10.**((0.7859+0.03477*TEM(L,KC))/  
     &            (1.+0.00412*TEM(L,KC))))  

              ! *** Evaporation
               FW=9.2+0.46*WINDST(L)**2;
               RE=FW*(SVPW-VPA(L));     !VPA fcn(TEMP,AirTemp,RelativeHumidity)

               ! *** Conduction
               RC=0.47*FW*(TEM(L,KC)-TATMT(L))

               NETRAD(L,KC)=LWDNB(L)-LWUPB(L)-RE-RC + SOLSWRT(L)    ! *** ADSORB SW SOLR RAD TO TO SURFACE LAYER  
               ! *** NOW FINALIZE THE TEMPERATURE 
               TEM(L,KC)=TEM(L,KC)+HPI(L)*C1*NETRAD(L,KC)  
             ENDDO  
           ELSE   ! Use existing implementationfor LW rad  based on Swinbank (1963)
             DO L=2,LA  
               SVPW=(10.**((0.7859+0.03477*TEM(L,KC))/  
     &            (1.+0.00412*TEM(L,KC))))  

               ! *** Net atmospheric radiation (Diffusive)
               CLDFAC=1.0+0.0017*CLOUDT(L)**2
               ! ** .5153153831e-12 = 1000.0/3600.0*9.37E-6*2.0411E-7*0.97
               RAN=0.51531538e-12*(TATMT(L)+273.15)**6*CLDFAC
 
               ! *** Evaporation
               FW=9.2+0.46*WINDST(L)**2;
               RE=FW*(SVPW-VPA(L));     !VPA fcn(TEMP,AirTemp,RelativeHumidity)

               ! *** Conduction
               RC=0.47*FW*(TEM(L,KC)-TATMT(L))

               ! *** Longwave back radiation
               ! ***  5.443E-8 = 5.67E-8 * 0.97
               RB=5.443E-8*(TEM(L,KC)+273.15)**4

               NETRAD(L,KC)=RAN-RB-RE-RC + SOLSWRT(L)    ! *** ADSORB SW SOLR RAD TO TO SURFACE LAYER  
               ! *** NOW FINALIZE THE TEMPERATURE 
               TEM(L,KC)=TEM(L,KC)+HPI(L)*C1*NETRAD(L,KC)  
             ENDDO  
           END IF
       ENDIF  
  
      ELSEIF(ISTOPT(2).EQ.2)THEN  
        ! *** IMPLEMENT EXTERNALLY SPECIFIED EQUILIBRIUM TEMPERATURE FORMULATION  
        TMPKC=DELT/DZC(KC)  
        DO ND=1,NDM  
          LF=2+(ND-1)*LDM  
          LL=LF+LDM-1  
          DO L=LF,LL  
            TEM(L,KC)=TEM(L,KC)-TMPKC*CLOUDT(L)*HPI(L)*(TEM(L,KC)
     &          -TATMT(L))  
!            TEM(L,KC)=TEM(L,KC)-TMPKC*HEQT*HPI(L)*(TEM(L,KC)
!     &          -TATMT(L))  !SCJ changed
          ENDDO  
        ENDDO  

      ELSEIF(ISTOPT(2).EQ.3)THEN  
        ! *** IMPLEMENT CONSTANT COEFFICIENT EQUILIBRIUM TEMPERATURE FORMULATION  
        DTHEQT=DELT*HEQT*FLOAT(KC)  
!        DO ND=1,NDM  
!          LF=2+(ND-1)*LDM  
!          LL=LF+LDM-1  
!          DO L=LF,LL  
!            TEM(L,KC)=TEM(L,KC)-DTHEQT*HPI(L)*(TEM(L,KC)-TEMO)  
!          ENDDO  
!        ENDDO 
        DO L=2,LA  
          TEM(L,KC)=TEM(L,KC)-DTHEQT*HPI(L)*(TEM(L,KC)-TEMO)  
        ENDDO

      ELSEIF(ISTOPT(2).EQ.4)THEN  
        ! *** IMPLEMENT W2 EQUILIBRIUM TEMPERATURE FORMULATION 
        IF(.NOT.COMPUTESOLRAD)THEN
          ! *** MUST MAKE AT LEAST ONE CALL TO THIS TO INITIALIZE VARIABLES
          CALL SHORT_WAVE_RADIATION(WINDST(2),RHA(2),TATMT(2),CLOUDT(2),
     &                            PATMT(2),SRO,SRON)
        ENDIF
        ! *** SWRATNF - Background/Clear Water Extinction Coefficient
        ! *** SWRATNS - Light Extinction for TSS (1/m per g/m3)
        ! *** FSWRATF - Fraction of Solar Rad Absobed in the Surface Layer
        ! *** HTBED2  - Bottom Heat Exchange Coefficient (W/m2/s)
        PSHADE_OLD=-1.0
        DO L=2,LA
          IF(PSHADE_OLD.NE.PSHADE(L))THEN
            IF(SOLSWRT(L).gt.0.01.OR.PSHADE_OLD.LT.-.99)then
              CALL EQUILIBRIUM_TEMPERATURE(SOLSWRT(L),ET,CSHE) 
            ENDIF
            PSHADE_OLD=PSHADE(L)
          ENDIF

          ! *** SURFACE HEAT FLUX
          THICK  =HP(L)*DZC(KC)  
          TFLUX = CSHE*(ET-TEM(L,KC))/THICK*DELT
          TEM(L,KC) = TEM(L,KC)+TFLUX 
          ! *** BEGIN PMC
          ! *** TEMPORARY FIX UNTIL BUILD IN ICE SUB-MODEL INTO THE HEAT SUB-MODEL (COOK INLET)
          IF(ISTRAN(1)>0)THEN
            IF( TEM(L,KC)<-1.3 )THEN
              TEM(L,KC) = -1.3*(SAL(L,KC)/35.)
            ENDIF
          ELSE
            IF( TEM(L,KC)<0.1 )THEN
              TEM(L,KC) = 0.1
            ENDIF
          ENDIF
          ! *** END PMC

          ! *** BOTTOM HEAT FLUX
          THICK = HP(L)*DZC(1)  
          TFLUX = HTBED2*(TEMB(L)-TEM(L,1))*DELT
          TEM(L,1) = TEM(L,1)+TFLUX/THICK 
          FLUXTB(L)=TFLUX
        ENDDO

        ! *** Distribute Solar Radiation Across Water Column
        IF(SOLSWRT(2).GT.0.1)THEN
          DO L=2,LA
            K=KC
            IF(ISTRAN(6).GT.0.OR.ISTRAN(7).GT.0)THEN
              TSSS_ABOVE=SNDT(L,K)+SEDT(L,K)
              TSS_ABOVE=TSSS_ABOVE
            ELSE
              TSS_ABOVE=0.0
            ENDIF
            IF(ISTRAN(8).GT.0)THEN
              ! *** Water Quality is Active so account for Chlorophyll and POM
              ! *** If using WQ then use the WQ Coefficients
              WQCHLS_ABOVE=WQCHL(L,K)
              WQCHL_ABOVE=WQCHLS_ABOVE
              POMS_ABOVE=WQV(L,K,4)+WQV(L,K,5) 
              POM_ABOVE =POMS_ABOVE
              GAMMA = WQKEB(1) + WQKETSS*TSS_ABOVE   +
     &                           WQKECHL*WQCHL_ABOVE +
     &                           WQKEPOM*POM_ABOVE
            ELSE
              GAMMA = SWRATNF + SWRATNS*TSS_ABOVE
            ENDIF

            TOP=GAMMA*HP(L)*(Z(K-1)-1.)
            EXPTOP=EXP(TOP)
            K_ABOVE=1.

            ! *** ENSURE AT LEAST THE FSWRATF FRACTION OF SRO IS ATTENUATED
            IF((1.-EXPTOP).GT.FSWRATF)THEN
              SRON=SOLSWRT(L)*EXPTOP
            ELSE
              SRON=SOLSWRT(L)*(1.0-FSWRATF)
            ENDIF
            DO K = KS,1,-1
              ! *** Net Extinction Coefficient
              K_ABOVE=K_ABOVE+1.
              IF(ISTRAN(6).GT.0.OR.ISTRAN(7).GT.0)THEN
                TSSS_ABOVE=TSSS_ABOVE+SNDT(L,K)+SEDT(L,K)
                TSS_ABOVE=TSSS_ABOVE/K_ABOVE
              ENDIF
              IF(ISTRAN(8).GT.0)THEN
                ! *** Water Quality is Active so account for Chlorophyll
                ! *** If using WQ then use the WQ Coefficients
                POMS_ABOVE=POMS_ABOVE+WQV(L,K,4)+WQV(L,K,5) 
                POM_ABOVE=POMS_ABOVE/K_ABOVE

                WQCHLS_ABOVE=WQCHLS_ABOVE+WQCHL(L,K)
                WQCHL_ABOVE=WQCHLS_ABOVE/K_ABOVE
                GAMMA = WQKEB(1) + WQKETSS*TSS_ABOVE   +
     &                             WQKECHL*WQCHL_ABOVE +
     &                             WQKEPOM*POM_ABOVE
              ELSE
                GAMMA = SWRATNF + SWRATNS*TSS_ABOVE
              ENDIF

              BOT=GAMMA*HP(L)*(Z(K-1)-1.)

              ! *** Compute Net Energy
              EXPBOT=EXP(BOT)
              NETRAD(L,K)=SRON*(EXPTOP-EXPBOT)  
              TOP=BOT
              EXPTOP=EXPBOT
            ENDDO  
            RADBOT(L)=EXPBOT*SRON
          ENDDO

          ! *** NOW FINALIZE THE TEMPERATURE 
          DO K=1,KS
            ! *** RHO = 1000.0  Density (kg / m^3)  
            ! *** CP  = 4179.0  Specific Heat (J / kg / degC)
            ! *** 0.2393E-6 = 1/RHO/CP --> Conversion from Watts
            C1=DELT*DZIC(K)*0.2393E-6
            DO L=2,LA  
              TEMO=TEM(L,K)+HPI(L)*C1*NETRAD(L,K)
              !IF(ABS(TEM(L,K)-TEMO).GT.PMCTOL)THEN
              !  IPMC=0
              !ENDIF
              TEM(L,K)=TEMO
            ENDDO  
          ENDDO  
          
          !DO L=2,LA   ! *** PMC PMC PMC
          !  DO K=1,KS
          !    IF(HP(L).LT.1.0.AND.(TEM(L,K)-PMCTOL).GT.TEM(L,K+1)) THEN
          !      IPMC = 0
          !    ENDIF
          !  ENDDO
          !ENDDO

        ENDIF
        ! *** UPDATE BOTTOM
        IF(TBEDIT.GT.0.)THEN
          IF(SOLSWRT(2).GT.0.01.AND.BETAF.GT.0.)THEN
            DO L=2,LA
              TEMB(L)=TEMB(L) + 
     &                (0.2393E-6*RADBOT(L)*DELT*BETAF - FLUXTB(L))
     &                      /TBEDTHK(L)
            ENDDO
          ELSE
            DO L=2,LA
              TEMB(L)=TEMB(L) - FLUXTB(L)/TBEDTHK(L)
            ENDDO
          ENDIF
        ENDIF            
      ENDIF  

      ! *** APPLY DRY CELL CORRECTIONS
      IF(ISDRY.GT.0)THEN  
        DO K=1,KC  
          DO L=2,LA  
            IF(.NOT.LMASKDRY(L))THEN
              TEM(L,K)=TATMT(L)  
              ! *** BEGIN PMC
              ! *** TEMPORARY FIX UNTIL BUILD IN ICE SUB-MODEL INTO THE HEAT SUB-MODEL (COOK INLET)
              IF(ISTRAN(1)>0)THEN
                IF( TEM(L,K)<-1.3 )THEN
                  TEM(L,K) = -1.3*(SAL(L,K)/35.)
                ENDIF
              ELSE
                IF( TEM(L,K)<0.1 )THEN
                  TEM(L,K) = 0.1
                ENDIF
              ENDIF
              ! *** END PMC
            ENDIF
          ENDDO  
        ENDDO  
      ENDIF  

  600 FORMAT(4I5,2E12.4)  

      RETURN  
      END  

************************************************************************
**           S U B R O U T I N E   H E A T  E X C H A N G E           **
**                                                                    **
**                      FROM CE-QUAL-W2 (VER 3.1)                     **
**                                                                    **
************************************************************************

      SUBROUTINE HEAT_EXCHANGE

      USE GLOBAL
      
******* Tupe declaration

        REAL      JDAY
        INTEGER*4 IDAY

******* Allocate/Dimension declaration

        Real ::  MPS_TO_MPH, W_M2_TO_BTU_FT2_DAY 
        Real ::  FLUX_BR_TO_FLUX_SI, BTU_FT2_DAY_TO_W_M2 
        Real ::  REFL 
        Real ::  AFW,BFW,CFW, BCONV
        Real ::  TDEW_F, TAIR_F, WIND_2M
        Real ::  TIMENEXT


******* Data declaration

        DATA MPS_TO_MPH          /2.23714/, 
     .       W_M2_TO_BTU_FT2_DAY /7.60796/,
     .       FLUX_BR_TO_FLUX_SI  /0.23659/
        DATA BTU_FT2_DAY_TO_W_M2 /0.1314/
        DATA AFW /9.2/, BFW /0.46/, CFW /2.0/, BCONV /1.520411/
        DATA REFL/0.06/  ! *** SHORTWAVE RADIATION RELFECTED OFF SURFACE  

******* Function declaration

        DEG_F(X) = X*1.8+32.0
        DEG_C(X) = (X-32.0)*5.0/9.0

        TIMENEXT=TIMEDAY+1./24.
      RETURN

************************************************************************
**                 S H O R T  W A V E  R A D I A T I O N              **
************************************************************************

      ENTRY SHORT_WAVE_RADIATION(WSPD,TD,TAIR,CLD,ATMPR,SRO,SRON)

******* Input Conversions
        IF(TD.LT.1.1.AND.IRELH(NASER).EQ.1)THEN
          ! *** TD IS RELATIVE HUMIDITY.  CONVERT TO DEW POINT
          ! *** Jensen et al. (1990) ASCE Manual No. 70 (see pages 176 & 177)

          ! *** Ambient vapor pressure in kPa 
          VaporP=TD* 0.611*EXP(17.27*TAIR/(TAIR+237.3)) 
 
          ! *** Compute dewpoint temperature (Tdew) in C
          Tdew=(116.9+237.3*LOG(VaporP))/(16.78-LOG(VaporP))

          ! *** Compute wet bulb temperature (Tw)
          !PSR=ATMPR/10.  ! ***  Convert millibars to kPa
          !Gamma=0.00066*PSR
          !Delta=4098.*VaporP/(Tdew+237.3)**2
          !Tw=((Gamma*TAIR)+(Delta*Tdew))/(Gamma+Delta)
        ELSE
          TDEW=TD
        ENDIF

        CLD10 = CLD*10.  ! *** CONVERT FROM 0-1 (EFDC) TO 0-10 (W2)
        IF(WINDH<=0.0)WINDH = 10.0     ! *** FIXING THE WIND SPEED MEASUREMENT HEIGHT

******* British units

        TDEW_F   = DEG_F(TDEW)
        TAIR_F   = DEG_F(TAIR)
        WIND_MPH = WSPD*MPS_TO_MPH
        WIND_2M  = WIND_MPH*(LOG(2.0/0.003)/LOG(WINDH/0.003))

******* Shortwave Radiation

        STANDARD = 15.0*INT(DS_LONG/15.0)

        ! *** Day of the Year
        THOUR    = (TIMEDAY-INT(TIMEDAY))*24.0
        IDAY     =  TIMEDAY-INT(TIMEDAY/365.)*365.
        IDAY     =  IDAY+INT(INT(TIMEDAY/365.)/4.)
        JDAY     = REAL(IDAY)
        PMC1     = (2.*PI*(JDAY-1.))/365.
        EQTNEW   =  0.170*SIN(4.*PI*(JDAY-80.)/373.)-
     &              0.129*SIN(2.*PI*(JDAY-8.)/355.)
        H     = 0.2618*(THOUR-(DS_LONG-STANDARD)*0.066667+EQTNEW-12.0)
        DECL =  0.006918-0.399912*COS(PMC1)  +0.070257*SIN(PMC1)-
     &                   0.006758*COS(2*PMC1)+0.000907*SIN(2*PMC1)-
     &                   0.002697*COS(3*PMC1)+0.001480*SIN(3*PMC1)
        SINAL = SIN(DS_LAT*.017453)*SIN(DECL)+COS(DS_LAT*.017453)
     .          *COS(DECL)*COS(H)
        A0    = 57.2958*ASIN(SINAL)
        IF (A0.GT.0.0) THEN
          SRO  = 2.044*A0+0.1296*A0**2-1.941E-3*A0**3+7.591E-6*A0**4
          SRO  = (1.0-0.0065*CLD10**2)*SRO*24.0
          SRO  = SRO*BTU_FT2_DAY_TO_W_M2
          SRON = SRO*(1.0-REFL) ! *** Adjust for surface reflection
        ELSE
          SRO  = 0.0
          SRON = 0.0
        END IF
      RETURN


************************************************************************
**             E Q U I L I B R I U M  T E M P E R A T U R E           **
************************************************************************

      ENTRY EQUILIBRIUM_TEMPERATURE(SRON,ET,CSHE)

******* British units
        
        ! *** SRON Should already be adjusted for Shading & Reflection
        SRO_BR   = SRON*W_M2_TO_BTU_FT2_DAY 

******* Equilibrium temperature and heat exchange coefficient
        
        ET    = TDEW_F
        TSTAR = (ET+TDEW_F)*0.5
        BETA  = 0.255-(8.5E-3*TSTAR)+(2.04E-4*TSTAR*TSTAR)
        FW    = W_M2_TO_BTU_FT2_DAY*AFW+BCONV*BFW*WIND_2M**CFW
        CSHE  = 15.7+(0.26+BETA)*FW
        RA    = 3.1872E-08*(TAIR_F+459.67)**4
        ETP   = (SRO_BR+RA-1801.0)/CSHE+(CSHE-15.7)
     .            *(0.26*TAIR_F+BETA*TDEW_F)/(CSHE*(0.26+BETA))
        J     = 0
        DO WHILE (ABS(ETP-ET).GT.0.05.AND.J.LT.10)
          ET    = ETP
          TSTAR = (ET+TDEW_F)*0.5
          BETA  = 0.255-(8.5E-3*TSTAR)+(2.04E-4*TSTAR*TSTAR)
          CSHE  = 15.7+(0.26+BETA)*FW
          ETP   = (SRO_BR+RA-1801.0)/CSHE+(CSHE-15.7)
     .            *(0.26*TAIR_F+BETA*TDEW_F)/(CSHE*(0.26+BETA))
          J     = J+1
        END DO

******* SI units

        ! *** RHO = 1000.0  Density (kg / m^3)
        ! *** CP  = 4179.0  Specific Heat (J / kg / degC)
        ! *** 0.2393E-6 = 1/RHO/CP
        ET   = DEG_C(ET)
        CSHE = CSHE*FLUX_BR_TO_FLUX_SI*0.2393E-6

        IF(DEBUG)THEN
          IF(TIMEDAY.GT.TIMENEXT)THEN
            TD_C=DEG_C(TDEW_F)
            TA_C=DEG_C(TAIR_F)
            WRITE(77,999)TIMEDAY,SRON,ET,TD_C,TA_C,TDEW_F,TAIR_F,FW,CSHE
  999       FORMAT(F11.4,8F9.2,E12.4)
            TIMENEXT = TIMENEXT+1./24.
          ENDIF
        ENDIF
      RETURN

      END

