INTEGER I,J,kUp,kUpC,kTop

      INTEGER I, J, K, bi, bj

        DO J = 1-OLy, sNy+OLy

           IF ( hFacC(I,J,K,bi,bj) .NE. 0. _d 0 ) 

     &          kTopC(I,J,bi,bj) = K

        IF ( kLev .EQ. kTopC(I,J,bi,bj) ) THEN

        IF ( kLev .EQ. kTopC(I,J,bi,bj) ) THEN

      INTEGER I,J,kUp,kUpC,kTop

     &           - shelfIceHeatFlux(I,J,bi,bj)

     &            shelfIceFreshWaterFlux(I,J,bi,bj)

     &            * salt(I,J,K,bi,bj) * convertEmP2rUnit

     &            shelfIceFreshWaterFlux(I,J,bi,bj)

         DO J = 1, sNy

           K    = kTopC(I,J,bi,bj)

           pLoc = ABS(R_shelfIce(I,J,bi,bj))

            tLoc = theta(I,J,K,bi,bj)

            sLoc = max(salt(I,J,K,bi,bj), 0. _d 0)

            shelfIceHeatFlux(I,J,bi,bj) = 

            shelfIceFreshWaterFlux(I,J,bi,bj) = 

     &           - shelfIceHeatFlux(I,J,bi,bj)

     &           shelfIceFreshWaterFlux(I,J,bi,bj)

      INTEGER I,J,K

         DO J = 1, sNy

           K    = kTopC(I,J,bi,bj)

           pLoc = ABS(R_shelfIce(I,J,bi,bj))

            tLoc = theta(I,J,K,bi,bj)

            sLoc = max(salt(I,J,K,bi,bj), 0. _d 0)

            shelfIceHeatFlux(I,J,bi,bj) = 

            shelfIceFreshWaterFlux(I,J,bi,bj) = 

     &           - shelfIceHeatFlux(I,J,bi,bj)

      INTEGER I,J

      DO J=1,aNy

        fld_atm(I,J)=0.

      DO J=1,MIN(aNy,oNy)

        fld_atm(I,J)=fld_ocn(I,J)

      INTEGER I,J

      DO J=1,oNy

        fld_ocn(I,J)=0.

      DO J=1,MIN(aNy,oNy)

        fld_ocn(I,J)=fld_atm(I,J)

          PPviscAr(I,J,K,bi,bj) = MAX(PPviscMin,PPviscAr(I,J,K,bi,bj))

          PPdiffKr(I,J,K,bi,bj) = MAX(PPdiffMin,PPdiffKr(I,J,K,bi,bj))

       DO J=jMin,jMax

         PPviscAr(I,J,K,bi,bj) = PPviscAr(I,J,K,bi,bj)

     &        * maskC(I,J,K,bi,bj)

         PPdiffKr(I,J,K,bi,bj) = PPdiffKr(I,J,K,bi,bj)

     &        * maskC(I,J,K,bi,bj)

      INTEGER I, J, K

       DO J=jMin,jMax

         IF ( RiNumber(I,J) .LT. RiLimit ) THEN

          denom = 1.0 + PPalpha*RiNumber(I,J)

         PPviscAr(I,J,K,bi,bj) = MAX(PPviscTmp,viscAr)

         PPdiffKr(I,J,K,bi,bj) = MAX(PPviscAr(I,J,K,bi,bj)/denom,

        DO J=jMin,jMax

      DO J= jMin, jMax

        tempu= .5*(  uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)

     &            - (uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )

        tempv= .5*(  vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)

     &            - (vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )

     &       (rhoKm1(I,J) - rhoK(I,J))*recip_drC(K)

        RiNumber(I,J) = buoyFreq/max(RiFlux,epsilon)

      DO J=jMin,jMax

        RiTmp(I,J) = RiNumber(I,J)

      DO J=1-Oly+1,sNy+Oly-1

        RiNumber(I,J) = p5*RiTmp(I,J) 

     &       + p125*RiTmp(I-1,J) + p125*RiTmp(I+1,J) 

     &       + p125*RiTmp(I,J-1) + p125*RiTmp(I,J+1) 

      INTEGER I,J,Km1

      INTEGER I, J

      INTEGER I, J

      INTEGER I, J

      DO J=1,NGP

        CLTOP(J,myThid)=CLTOP(J,mythid)*CLOUDC(J,myThid)

        DO J=1,sNy

          I2 = I+(J-1)*sNx

          DO J=1,NGP

            tmpVar(J) = EnPrec(J,myThid)

     &                *(PRECNV(J,myThid)+PRECLS(J,myThid))

      INTEGER I,J,K, I2,Katm

       DO J=1,NGP

         dTsurf(J,1,myThid) = fMask1(J,1,myThid)*dTsurf(J,1,myThid)

     &                      + fMask1(J,2,myThid)*dTsurf(J,2,myThid)

     &                      + fMask1(J,3,myThid)*dTsurf(J,3,myThid)

         dTsurf(J,1,myThid) = ABS(dTsurf(J,1,myThid))

        DO J=1,NGP

          DQLSC(J,K) = DQLSC(J,K)*dpFac(J,K)

          DTLSC(J,K) = DTLSC(J,K)*dpFac(J,K)

      DO J=1,NGP

        DO K=2,kGrd(J)

          PRECLS(J) = PRECLS(J)-PFACT*DQLSC(J,K)

      INTEGER J, K

      DO J=1,NGP

        DTLSC(J,1) = 0.

        DQLSC(J,1) = 0.

        PRECLS(J)  = 0.

        PSA2(J)    = PSA(J)*PSA(J)

        DO J=1,NGP

          RHREF = RHLSC+DRHLSC*(SIG2/PSA2(J) - 1. _d 0)

          DQLSC(J,K) = MIN(0. _d 0,(RHREF*QSAT(J,K)-QA(J,K)))*RTLSC

          DQLSC(J,K) = MAX(-DQMAX, DQLSC(J,K) )

          DTLSC(J,K) = -TFACT*DQLSC(J,K)

      DO J=1,NGP

        TSKIN(J) = Tsurf(J) + dTskin(J)

        TSFC(J)=273.16 _d 0 + dTskin(J)

      DO J=1,NGP

        FSLAND=1.+MIN(DTHETA,MAX(-DTHETA,TSKIN(J)-T1(J)))*RDTH

        CDENVV(J)=CHL*DENVV(J)*FSLAND

       DO J=1,NGP

        Fstb0 = 1.+MIN(DTHETA,MAX(-DTHETA,TSFC(J) -T1(J)))*RDTH

        Shf0(J) = CHL*DENVV(J)*Fstb0

        dTstb = ( DTHETA+dTstab-ABS(TSKIN(J)-T1(J)) )/dTstab

        dShf(J) = CHL*DENVV(J)*dFstb

       DO J=1,NGP

        EVAP(J) = CDENVV(J)*SWAV(J)*(QSAT0(J,1)-Q0(J))

        Evp0(J) =   Shf0(J)*SWAV(J)*(QSAT0(J,2)-Q0(J))

        dEvp(J) = CDENVV(J)*SWAV(J)*dEvp(J)

     &            + dShf(J)*SWAV(J)*(QSAT0(J,1)-Q0(J))

       DO J=1,NGP

        EVAP(J) = CDENVV(J)*SWAV(J)*(QSAT0(J,1)-Q0(J))

        Evp0(J) = CDENVV(J)*SWAV(J)*(QSAT0(J,2)-Q0(J))

        dEvp(J) = CDENVV(J)*SWAV(J)*dEvp(J)

       DO J=1,NGP

        EVAP(J) = CDENVV(J)*SWAV(J)*MAX(0. _d 0,QSAT0(J,1)-Q0(J))

        Evp0(J) =   Shf0(J)*SWAV(J)*MAX(0. _d 0,QSAT0(J,2)-Q0(J))

        dEvp(J) = CDENVV(J)*SWAV(J)*dEvp(J)

     &            + dShf(J)*SWAV(J)*MAX(0. _d 0,QSAT0(J,1)-Q0(J))

       DO J=1,NGP

        EVAP(J) = CDENVV(J)*SWAV(J)*MAX(0. _d 0,QSAT0(J,1)-Q0(J))

        Evp0(J) = CDENVV(J)*SWAV(J)*MAX(0. _d 0,QSAT0(J,2)-Q0(J))

        dEvp(J) = CDENVV(J)*SWAV(J)*dEvp(J)

       DO J=1,NGP

        SHF(J)  = CDENVV(J)*CP*(TSKIN(J)-T0(J))

        Shf0(J) =   Shf0(J)*CP*(TSFC(J) -T0(J))

        dShf(J) = CDENVV(J)*CP

     &            + dShf(J)*CP*(TSKIN(J)-T0(J))

        dShf(J) = MAX( dShf(J), 0. _d 0 )

        dEvp(J) = MAX( dEvp(J), -dShf(J)/ALHC )

       DO J=1,NGP

        SHF(J)  = CDENVV(J)*CP*(TSKIN(J)-T0(J))

        Shf0(J) = CDENVV(J)*CP*(TSFC(J) -T0(J))

        dShf(J) = CDENVV(J)*CP

      DO J=1,NGP

        TS2     = TSFC(J)*TSFC(J)

        Slr0(J) = SBC*TS2*TS2

        TS2     = TSKIN(J)*TSKIN(J)

        SLRU(J) = SBC*TS2*TS2

        dSlr(J)  = 4. _d 0 *SBC*TS2*TSKIN(J)

      DO J=1,NGP

        sFlx(J,0)= ( SSR(J) + SLRD(J) - EMISloc*Slr0(J) )

     &           - ( Shf0(J) + ALHC*Evp0(J) )

        sFlx(J,1)= ( SSR(J) + SLRD(J) - EMISloc*SLRU(J) )

     &           - ( SHF(J)+ ALHC*EVAP(J) )

        sFlx(J,2)=                    - EMISloc*dSlr(J)

     &           - ( dShf(J) + ALHC*dEvp(J) )

      INTEGER J

      INTEGER K, J

      DO J=1,NGP

        NL1(J)=kGrd(J)-1

      DO J=1,NGP

        CLOUDC(J)=0.

        QCLOUD(J)=0.

        ICLTOP(J)=NLEV+1

        FREFL(J,1)=0.

       DO J=1,NGP

         ALBTOP(J,K)=0.

       DO J=1,NGP

        ICLTOP(J)= kGrd(J)

        DO K=NL1(J),2,-1

         IF (RH(J,K).GT.RHCL1.AND.QA(J,K).GT.QACL1) THEN

            CLOUDC(J)=MAX(CLOUDC(J),RH(J,K)-RHCL1)

            IF (QA(J,K).GT.QACL3) ICLTOP(J)=K

       DO J=1,NGP

        IF (kGrd(J).NE.0)

     &  QCLOUD(J)= MAX( QA(J,kGrd(J)), QA(J,NL1(J)) )

        CLOUDC(J)=MIN(1. _d 0,CLOUDC(J)*RRCL)

        IF (CLOUDC(J).GT.0.0) THEN

          CLOUDC(J)=CLOUDC(J)*MIN(1. _d 0,QCLOUD(J)*RQCL)

          ALBTOP(J,ICLTOP(J))=ALBCL*CLOUDC(J)

          ICLTOP(J)=NLEV+1

      DO J=1,NGP

        ACLOUD(J)=CLOUDC(J)*(ABSCL1+ABSCL2*QCLOUD(J))

      DO J=1,NGP

        DELTAP=ZENIT(J)*DSIG(1)*dpFac(J,1)

        TAU2(J,1,1)=EXP(-DELTAP*ABSDRY)

      DO J=1,NGP

       DO K=2,NL1(J)

         ABS1=ABSDRY+ABSAER*(SIG(K)/PSA(J))**2

         DELTAP=ZENIT(J)*DSIG(K)*dpFac(J,K)

         IF (K.EQ.ICLTOP(J)) THEN

           TAU2(J,K,1)=EXP(-DELTAP*

     &                 (ABS1+ABSWV1*QA(J,K)+2.*ACLOUD(J)))

         ELSE IF (K.GT.ICLTOP(J)) THEN

           TAU2(J,K,1)=EXP(-DELTAP*

     &                 (ABS1+ABSWV1*QA(J,K)+ACLOUD(J)))

           TAU2(J,K,1)=EXP(-DELTAP*(ABS1+ABSWV1*QA(J,K)))

      DO J=1,NGP

       K = kGrd(J)

       ABS1=ABSDRY+ABSAER*(SIG(K)/PSA(J))**2

        DELTAP=ZENIT(J)*DSIG(K)*dpFac(J,K)

        TAU2(J,K,1)=EXP(-DELTAP*(ABS1+ABSWV1*QA(J,K)))

      DO J=1,NGP

       DO K=2,kGrd(J)

         DELTAP=ZENIT(J)*DSIG(K)*dpFac(J,K)

         TAU2(J,K,2)=EXP(-DELTAP*ABSWV2*QA(J,K))

      DO J=1,NGP

        FTOP(J)  =FSOL(J)

        FLUX(J,1)=FSOL(J)*FBAND1-OZUPP(J)

        FLUX(J,2)=FSOL(J)*FBAND2

        STRATC(J)=STRATZ(J)*PSA(J)

      DO J=1,NGP

        DFABS(J,1)=FLUX(J,1)

        FLUX (J,1)=TAU2(J,1,1)*(FLUX(J,1)-OZONE(J)*PSA(J))

        DFABS(J,1)=DFABS(J,1)-FLUX(J,1)

      DO J=1,NGP

       DO K=2,kGrd(J)

         FREFL(J,K)=FLUX(J,1)*ALBTOP(J,K)

         FLUX (J,1)=FLUX(J,1)-FREFL(J,K)

         DFABS(J,K)=FLUX(J,1)

         FLUX (J,1)=TAU2(J,K,1)*FLUX(J,1)

         DFABS(J,K)=DFABS(J,K)-FLUX(J,1)

      DO J=1,NGP

       DO K=2,kGrd(J)

         DFABS(J,K)=DFABS(J,K)+FLUX(J,2)

         FLUX (J,2)=TAU2(J,K,2)*FLUX(J,2)

         DFABS(J,K)=DFABS(J,K)-FLUX(J,2)

      DO J=1,NGP

        FSFC(J,1)=FLUX(J,1)*(1.-ALB(J,1))+FLUX(J,2)

        FSFC(J,2)=FLUX(J,1)*(1.-ALB(J,2))+FLUX(J,2)

        FSFC(J,3)=FLUX(J,1)*(1.-ALB(J,3))+FLUX(J,2)

        FSFC(J,0)=FLUX(J,1)+FLUX(J,2)

        FLUX(J,1)=FLUX(J,1)*ALB(J,0)

        FSFC(J,0)=FSFC(J,0)-FLUX(J,1)

      DO J=1,NGP

       DO K=kGrd(J),1,-1

         DFABS(J,K)=DFABS(J,K)+FLUX(J,1)