#define c_str(s) ('s\0')
      subroutine depth_init (npt, zin)
c     ------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_data.h'
      include 'comm_new.h'
      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch
      logical sink_region
      common /bathymetry/sink_region,sr_depth, alon_sr, blon_sr, alat_sr, blat_sr
      dimension zin(1)

      if     (initb .eq. 0) then
c.....Constant Depth (Flat Bottom)
         dep_max = zin(nz+1)
         call afill(npt, dept, dep_max)
         depmin = dep_max
      elseif (initb .eq. 1) then
c.....Read Bathymetry Data From a File

         call odb_open(idf_dp, fbdep(1:n_dep), 0)

         call data_on_model_grid(idf_dp, lret, 'bath') 

         call read_zt (idf_dp, lret, npt, 1, 1, 'bath', dept)

         call odb_close(idf_dp)
      elseif (initb .eq. 2) then
         iseed = 10001
         dept(1) = 50 + zin(nz+1)*ran(iseed)
         do i = 2, npt
            dept(i) = 50 +  zin(nz+1)*ran(iseed)
         enddo

         do i = 1, npt
            dept(i) = max(dept(i),zin(1))
         enddo

      elseif (initb .eq. 3) then
c        beta ramp
         
         RAD    = 3.14159265/180.
         rlat = 0.5*(alat + blat)
         fz    = sin(rlat*rad)
         betap = cos(rlat*rad)

         ys = ym(1)
         yn = ym(nyp)
         yl = yn-ys
         
         const = betap*(blat-alat)*rad/fz/yl
        
         do j = 1, nyp
            tp(j) = zin(nz+1)/(1. + const*(ym(j)-ys))
         enddo
         print*,'south,north depths:',tp(1),tp(nyp)

         do k = 1, npt
            j = (iox(k)-1)/nxp + 1
            dept(k) = tp(j)
         enddo

      elseif (initb .eq. 4) then
c        ridge

         depmin = (zin(initb-2) + zin(initb-1))/2.
         depmax = zin(nz+1)
         x1 = xm(1)
         xn = xm(nxp)
         xs = xn-x1
         do i = 1, nxp
            tp(i) = depmin + (xm(i) - x1)*(depmax-depmin)/xs
         enddo
         y1 = ym(1)
         yn = ym(nyp)
         ys = yn-y1
         do j = 1, nyp
            tp(j+nxp) = depmin + (ym(j) - y1)*(depmax-depmin)/ys
         enddo

         do k = 1, npt
            j = (iox(k)-1)/nxp + 1
            i = iox(k) - (j-1)*nxp
            dept(k) = depmax
            if (i.ge.i_ridge_min.and.i.le.i_ridge_max) 
c     *                       dept(k) = (tp(i)+tp(j+nxp))/2.
     *                       dept(k) = tp(j+nxp)
         enddo

      elseif (initb .eq. 5) then
c       isolated sea-mount
        H0 = 4500
        H1 = 4050
        xl = (xm(nxp)-xm(1))/15.2
        x_m = (xm(nxp)+xm(1))/2.
        yl = (ym(nyp)-ym(1))/11.52
        y_m = (ym(nyp)+ym(1))/2.
        do k = 1, npt
           j = (iox(k)-1)/nxp + 1
           i = iox(k) - (j-1)*nxp
           x = (xm(i)-x_m)/xl
           y = (ym(j)-y_m)/yl
           dept(k) = h0 - h1*exp(-(x**2+y**2))
        enddo
      elseif (initb .eq. 6) then
         xw = xm(1)
         xe = xm(nxp)
         xl = xe-xw

         
         do j = 1, nxp
            tp(j) = 4000 - 1000*(xm(j) - xw)/xl
         enddo

         do k = 1, npt
            j = (iox(k)-1)/nxp + 1
            i = iox(k) - (j-1)*nxp
            dept(k) = tp(i)
         enddo
      elseif (initb .eq. 7) then
         ys = ym(1)
         yn = ym(nyp)
         yl = yn-ys
         
         do j = 1, nyp
            tp(j) = 4000 - 1000*(ym(j) - ys)/yl
         enddo

         do k = 1, npt
            j = (iox(k)-1)/nxp + 1
            dept(k) = tp(j)
         enddo
      elseif (initb .eq. 8) then
         xw = xm(1)
         xe = xm(nxp)
         xl = xe-xw
         
         do j = 1, nxp
            tp(j) = 3000 + 1000*(xm(j) - xw)/xl
         enddo

         do k = 1, npt
            j = (iox(k)-1)/nxp + 1
            i = iox(k) - (j-1)*nxp
            dept(k) = tp(i)
         enddo
      elseif (initb .eq. 9) then
         xw = xm(nxp/4)
         xe = xm(3*nxp/4)
         xl = xe-xw
         
         do j = 1,nxp/4-1
            tp(j) = 4000
         enddo
         do j = nxp/4,nxp/2
            tp(j) = 4000 - 1000*(xm(j) - xw)/xl
         enddo
         do j = nxp/2+1,3*nxp/4
            tp(j) = 4000 + 1000*(xm(j) - xe)/xl
         enddo
         do j = 3*nxp/4 + 1, nxp
            tp(j) = 4000 
         enddo

         do k = 1, npt
            j = (iox(k)-1)/nxp + 1
            i = iox(k) - (j-1)*nxp
            dept(k) = tp(i)
         enddo
      endif

      dept(1) = max(dep_min,dept(1))
      do i = 2, npt
         dept(i) = max(dep_min,dept(i))
      enddo

      dep_min = dept(1)
      dep_max = dept(1)
      do i = 2, npt
         dep_max = max(dep_max, dept(i))
         dep_min = min(dep_min, dept(i))
      enddo

      if (sink_region) then
         do k = 1, npt
            j = (iox(k)-1)/nxp + 1
            i = iox(k) - (j-1)*nxp
            xlon = xm(i)
            ylon = ym(j)
            if (xlon.ge.alon_sr.and.xlon.le.blon_sr.and.
     *                       ylon.ge.alat_sr.and.ylon.le.blat_sr) then
                dept(k) = max(dept(k),sr_depth)
            endif
         enddo
      endif

      return
      end

c     ------------------------------------------------
      subroutine gradH_init (nx,ny,npt, zin, gradH, dept,
     *     saxpk,saypk,isk,tp,dx2i,dy2i,ycos,inx,iny,lxx,lyy)
c     ------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      include 'comm_para.h'

      dimension isk(npt),inx(nx),iny(ny)
      dimension tp(npt), saxpk(npt),saypk(npt), dx2i(1),dy2i(1), 
     *          ycos(1), gradH(npt,2), dept(npt), lxx(1), lyy(1)

      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc

      if     (initg .eq. 0) then
c.....Constant Depth (Flat Bottom)
         call afill(npt, gradH, 0.)
      elseif (initg .eq. 1) then
c.....Read Grad(H)/H Data From a File

         call odb_open(idf_grx, fbgrad(1:n_grad)//'.x', 0)
         call odb_open(idf_gry, fbgrad(1:n_grad)//'.y', 0)

         call data_on_model_grid(idf_grx, lret, 'gradx') 
         call data_on_model_grid(idf_gry, lret, 'grady') 

         call read_zt (idf_grx, lret, npt, 1, 1, 'gradx', gradH(1,1))
         call read_zt (idf_gry, lret, npt, 1, 1, 'grady', gradH(1,2))

         call odb_close(idf_grx)
         call odb_close(idf_gry)
      elseif (initg .eq. 2) then
         call grad (nx,ny,npt,dept,dx2i,dy2i,ycos,
     *              inx,iny,isk,gradH(1,1),gradH(1,2),nbx,nby,lxx,lyy)
         do i = 1, npt
            depi = 1./dept(i)
            gradH(i,1) = depi* gradH(i,1)
            gradH(i,2) = depi* gradH(i,2)
            glimx = saxpk(i)*grad_fac
            glimy = saypk(i)*grad_fac
            gradH(i,1) = max(-glimx,min(glimx,gradH(i,1)))
            gradH(i,2) = max(-glimy,min(glimy,gradH(i,2)))
         enddo 
      endif
         
      return
      end

c     ------------------------------------------------
      subroutine Kmap_init (npt,nz,nxp,nyp)
c     ------------------------------------------------
      include 'comm_para.h'
      include 'comm_diff.h'
      include 'comm_data.h'
      include 'diffiso.h'

      real kmapmax,kmapmin

      call afill(npt, Kaspect, 1.)

      if     (initkmap .eq. 0) then
c.....Constant diffusion coefficient
         call afill(npt, Kmap, 1.)
      elseif (initkmap .eq. 1) then
c.....Read Diffusion Data From a File
         call odb_open(idf_dif, fbdif(1:n_dif), 0)
         call data_on_model_grid(idf_dif, lret, 'diff') 
         call read_zt (idf_dif, lret, npt, 1, 1, 'diff', kmap)

         call odb_close(idf_dif)
      elseif (initkmap .eq. 2) then
         emy_max = emy(1)
         do i = 1, npt
            emy_max = max(emy_max,emy(i))
         enddo
         kmapmax = 0
         kmapmin = 1e5
         do i = 1, npt
            Kmap(i)    = emy_max/emy(i)
            kmapmax = max(kmapmax,Kmap(i))
            kmapmin = min(kmapmin,Kmap(i))
         enddo
         print*,'max,min diffusion coef',kmapmax,kmapmin
      elseif (initkmap .eq. 3) then
         emy_max = emy(1)
         do i = 1, npt
            emy_max = max(emy_max,emy(i))
         enddo
         kmapmax = 0
         kmapmin = 1e5
         do i = 1, npt
            Kmap(i)    = emy_max/emy(i)
            Kaspect(i) = emy(i)/emx(i)
            kmapmax = max(kmapmax,Kmap(i))
            kmapmin = min(kmapmin,Kmap(i))
         enddo
         print*,'max,min diffusion factor',kmapmax,kmapmin
      elseif (initkmap .eq. 4) then
         call afill(npt, Kmap, 0.)
         call afill(npt, Kaspect, 0.)
         do i = 1, npt
            j = (iox(i)-1)/nxp + 1
            if (j.eq.1) then
               Kmap(i)    = cint_fac
               Kaspect(i) = bint_fac
            endif
         enddo
      elseif (initkmap .eq. 5) then
         call Kmap_compute(npt,nz,nxp,nyp)
      endif

      return
      end

c     --------------------------------------------------------------------
      subroutine clim_init(npt,nstart,h0,sigma,dzin,hmf,
     *                     hclf,tclf,sclf,dclf,psif,tpf,nsponge,lsponge)
c     --------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      include 'comm_data.h'
      dimension h0(1),hmf(npt,nz),hclf(npt,nz,1),dclf(npt,nz),dzin(nz+1),
     *          tclf(npt,nz,1),sclf(npt,nz,1),tpf(npt,1),ind3(3),sigma(nz)
     *          ,lsponge(nsponge), psif(npt,1)
      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch
c
      do k = nsig+1, nz
         do i = 1, npt
            hclf(i,k,1) = hmf(i,k)
         enddo
      enddo

      if     (icl_h .eq. 0) then
         do k = 1, nsig
            do i = 1, npt
               hclf(i,k,1) = hmf(i,k)
            enddo
         enddo

      elseif (icl_h .eq. 1) then
c........H: sigma structure between Mixed Layer & Thermocline - Static:
         call odb_open(idf_hcl, fbhcl(1:n_hcl), 0)
         call data_on_model_grid(idf_hcl, lclm, 'mltc')

         if (icl_htop .eq. 1) 
     *   call read_zt (idf_hcl,lclm,npt, 1,1, 'mltc', hclf(1,1,1))

         if (icl_htop .eq. 0) call afill (npt, hclf, h0(1))

         sigk = sigma(3)
         do i = 1, npt
            hclf(i,2,1) = 0.5*hclf(i,1,1) + sigk*(z_begin - 1.5*hclf(i,1,1))
         enddo
         do k = 3, nsig - 1
            sigkp = sigma(k+1)
            do i = 1, npt
               hclf(i,k,1) = (sigk+sigkp) * (z_begin - 1.5*hclf(i,1,1))
            enddo
            sigk = sigkp
         enddo
         k = nsig
         do i = 1, npt
            hclf(i,k,1) = dzin(k+1) + sigma(k)*(z_begin - 1.5*hclf(i,1,1))
         enddo
         call odb_close(idf_hcl)

      elseif (icl_h .eq. 2) then
c........H: sigma structure between Mixed Layer & Thermocline - Dynamic
         call odb_open(idf_hcl, fbhcl(1:n_hcl), 0)
         call data_on_model_grid(idf_hcl, lclm, 'mltc')

         call odb_rddm(idf_hcl, 'T', ntclm)
         call mem_alloc(p_tclm, ntclm, 2, 'tclm')
         call odb_rdgr(idf_hcl, 'T', ntclm, tclm)

         call it_catch (ntclm, tclm, nstart, it1, it2, clm_tscl)
         iclm = it2

         if (icl_htop .eq. 1) then
            call read_zt (idf_hcl, lclm,npt, 1, it1, 'mltc', hclf(1,1,1))
            call read_zt (idf_hcl, lclm,npt, 1, it2, 'mltc', hclf(1,1,2))
         endif

         if (icl_htop .eq. 0) then
            call afill (npt, hclf(1,1,1), h0(1))
            call afill (npt, hclf(1,1,2), h0(1))
         endif

         sigk = sigma(3)
         do i = 1, npt
            hclf(i,2,1) = 0.5*hclf(i,1,1) + sigk*(z_begin - 1.5*hclf(i,1,1))
            hclf(i,2,2) = 0.5*hclf(i,1,2) + sigk*(z_begin - 1.5*hclf(i,1,2))
         enddo
         do k = 3, nsig - 1
            sigkp = sigma(k+1)
            do i = 1, npt
               hclf(i,k,1) = (sigk+sigkp) * (z_begin - 1.5*hclf(i,1,1))
               hclf(i,k,2) = (sigk+sigkp) * (z_begin - 1.5*hclf(i,1,2))
            enddo
            sigk = sigkp
         enddo
         k = nsig
         do i = 1, npt
            hclf(i,k,1) = dzin(k+1) + sigma(k)*(z_begin - 1.5*hclf(i,1,1))
            hclf(i,k,2) = dzin(k+1) + sigma(k)*(z_begin - 1.5*hclf(i,1,2))
         enddo
         call odb_close(idf_hcl)
      endif

      if (icl_psi .eq. 1) then
         call odb_open(idf_psi, fbpsi(1:n_psi), 0)
         call data_on_model_grid(idf_psi, lpsi, 'psi')

         call read_zt (idf_psi,lpsi,npt, 1,1, 'psi', psif(1,1))
         call odb_close(idf_psi)

      elseif (icl_psi .eq. 2) then

         call odb_open(idf_psi, fbpsi(1:n_psi), 0)
         call data_on_model_grid(idf_psi, lpsi, 'psi')

         call odb_rddm(idf_psi, 'T', ntpsi)
         call mem_alloc(p_tpsi, ntpsi, 2, 'psi')
         call odb_rdgr(idf_psi, 'T', ntpsi, tpsi)

         call it_catch (ntpsi, tpsi, nstart, it1, it2, psi_tscl)
         ipsi = it2

         call read_zt (idf_psi, lpsi,npt, 1, it1, 'psi', psif(1,1))
         call read_zt (idf_psi, lpsi,npt, 1, it2, 'psi', psif(1,2))
         call odb_close(idf_psi)
      endif

      if (icl_ts .eq. 0) return

      if (initt.ne.5) then
         call odb_open(idf_t, fbtem(1:n_tem), 0)
         if (icl_h .eq. 0) call data_on_model_grid(idf_t, lclm, 'temp')
         
         call odb_rddm(idf_t, 'Z', nzclm)
         call mem_alloc(p_zclm, nzclm, 2, 'MZ for temp climatology')
         call odb_rdgr (idf_t, 'Z', nzclm, zclm)
         mz = nzclm
         
         
         call odb_open(idf_s, fbsal(1:n_sal), 0)
         call odb_rddm(idf_s, 'Z', mz)
         if (mz .ne. nzclm) 
     *       call perror1('Temp & Salt DATA should be on the same Z grid', 1)
         
         if (icl_ts .eq. 1) then !!.....time independent case:
       call read_linz(idf_t,lclm,npt,mpack,nz,mz,1,hclf,tclf,zclm,tp,'temp') 
       call read_linz(idf_s,lclm,npt,mpack,nz,mz,1,hclf,sclf,zclm,tp,'salt') 
            
            if (iv_bot .eq. 4) then
               do i = 1, npt
                  tclf(i,nz,1) = TEMP_BOT
                  sclf(i,nz,1) = SALT_BOT
               enddo
            endif
            
         elseif (icl_ts .eq. 2) then !!.....time dependent case:
            
            call odb_rddm(idf_t, 'T', i)
            
            if  (icl_h .ne. 2) then !!.....H_clim is time independent
               ntclm = i
               call mem_alloc(p_tclm, ntclm, 2, 'tclm')
               call odb_rdgr(idf_t, 'T', ntclm, tclm)
               call it_catch (ntclm, tclm, nstart, it1, it2, clm_tscl)
               iclm = it2
            elseif (i .ne. ntclm ) then
             call perror1('MxTh & Temp DATA should be on the same T grid', 1)
            endif
            
            call odb_rddm(idf_s, 'T', i)
            if (i .ne. ntclm) 
     *       call perror1('Temp & Salt DATA should be on the same T grid', 1)
            
            call read_linz(idf_t, lclm, npt,mpack,nz,mz, it1,
     *           hclf(1,1,1),tclf(1,1,1),zclm,tp, 'temp') 
            k2_h = 1
            if (icl_h .eq. 2) k2_h = 2
            call read_linz(idf_t, lclm, npt,mpack,nz,mz, it2,
     *           hclf(1,1,k2_h),tclf(1,1,2),zclm,tp, 'temp') 
            
            call read_linz(idf_s, lclm, npt,mpack,nz,mz, it1,
     *           hclf(1,1,1),sclf(1,1,1),zclm,tp, 'salt') 
            call read_linz(idf_s, lclm, npt,mpack,nz,mz, it2,
     *           hclf(1,1,k2_h),sclf(1,1,2),zclm,tp, 'salt') 
            
         endif

         call odb_close(idf_t)
         call odb_close(idf_s)

      endif
      
      
      if (icl_fix.eq.2) then
      n0 = 0
      n1 = npt 
      n2 = 2*npt 
      do ic = 1, min(2,icl_ts)
         do k = nz_fix, nz
            tp(n0+k) = 0.
            tp(n1+k) = 0.
            tp(n2+k) = 0.
         enddo
         do i = 1, npt
            do k = nz_fix, nzi(i)
               tp(n0+k) = tp(n0+k) + tclf(i,k,ic)
               tp(n1+k) = tp(n1+k) + sclf(i,k,ic)
               tp(n2+k) = tp(n2+k) + 1.
            enddo
         enddo
         do k = nz_fix, nz
            tp(n0+k) = tp(n0+k)/tp(n2+k)
            tp(n1+k) = tp(n1+k)/tp(n2+k)
         enddo
         do i = 1, npt
            do k = nz_fix, nzi(i)
               tclf(i,k,ic) = tp(n0+k)
               sclf(i,k,ic) = tp(n1+k)
            enddo
         enddo
      enddo
      endif

      if (icl_rlx .eq. 1) then
         do j = 1, nyp
            if     (ym(j) .gt. clm_no) then
               tp(j) = (ym(j)-clm_no)/(ym(nyp)-clm_no)
            elseif (ym(j) .lt. clm_so) then
               tp(j) = (clm_so-ym(j))/(clm_so-ym(1))
            else
               tp(j) = 0.
            endif
         enddo
         
         do i = 1, npt
            j = (iox(i)-1)/nxp + 1
            sponge(i) = tp(j)
         enddo

      elseif (icl_rlx .eq. 2) then
         do i = 1, npt
            sponge(i) = 1.
         enddo
      elseif (icl_rlx .eq. 3) then
         dsponge = real(ksponge)
         do i = 1, npt
            ixy = iox(i)
            ix = mod (ixy -1 ,nxp) + 1
            iy = (ixy - ix)/nxp + 1
            dmin = nxp+nyp
            do j = 1, nsponge
               kxy = lsponge(j)
               kx = mod (kxy -1 ,nxp) + 1
               ky = (kxy - kx)/nxp + 1
               d = sqrt(float((ix-kx)**2 + (iy-ky)**2))
c               d = abs(ix-kx) + abs(iy-ky)
               dmin = min(d,dmin)
            enddo
c            sponge(i) = max((dsponge-dmin)/dsponge,0.)
            sponge(i) = exp(-dmin/dsponge)
         enddo
      elseif (icl_rlx .eq. 4) then
c.....Read Sponge weights From a File

         call odb_open(idf_sp, fbsp(1:n_sp), 0)
         call data_on_model_grid(idf_sp, lret, 'bath') 
         call read_zt (idf_sp, lret, npt, 1, 1, 'weight', sponge)
         call odb_close(idf_sp)

      endif

      return
      end

c     ----------------------------------------------
      subroutine h_init (npt, nz, nzi, nstart, hmf, hclf)
c     ----------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      dimension hmf(npt,1), hclf(npt,nz,1), nzi(1)

      if     (icl_h .eq. 1) then
         do i = 1, npt
            do k = 1, nzi(i)
               hmf(i,k) = hclf(i,k,1) 
            enddo
         enddo

      elseif (icl_h .eq. 2) then
         do i = 1, npt
            do k = 1, nzi(i)
               hmf(i,k) = hclf(i,k,1) +  clm_tscl*(hclf(i,k,2) - hclf(i,k,1))
            enddo
         enddo

      endif
      
      return
      end

c    -----------------------------------------------------------------
      subroutine temp_init (npt, nzl, nstart, t0, tmf, tclf)
c    -----------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      include 'comm_data.h'
      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc

      dimension t0(1), tmf(npt,1), tclf(npt,nzl,1)
c
      if     (initt .eq. 0) then
c.......constant values for each layer according to T_INIT
         iseed = 10001
         do i = 1, npt
            do k = 1, nzi(i)
               rand = 1.+temp_coef*(2.*ran(iseed)-1.)
               tmf(i,k) = t0(k)*rand
            enddo
         enddo

      elseif (initt .eq. 1) then
c.......linearly decreasing temps with latitude
         y1 = ym(1)
         yn = ym(nyp)
         ys = yn-y1
         do j = 1, nyp
            tp(j) = 1. - (ym(j) - y1)*temp_coef/ys
         enddo

         do i = 1, npt
            j = (iox(i)-1)/nxp + 1
            do k = 1, nzi(i)
               tmf(i,k) = t0(k)*tp(j)
            enddo
         enddo

      elseif (initt .eq. 2) then
c.......front
         y1 = ym(1)
         yn = ym(nyp)
         ys = yn-y1
         do j = 1, nyp/2
            tp(j) = 1. 
         enddo
         do j = nyp/2+1,nyp
            tp(j) = 1. - temp_coef 
         enddo

         do i = 1, npt
            j = (iox(i)-1)/nxp + 1
            do k = 1, nzi(i)
               tmf(i,k) = t0(k)*tp(j)
            enddo
         enddo

      elseif (initt .eq. 3) then
c.....from CLIMATOLOGY Data.
         if     (icl_ts .eq. 2) then
            do i = 1, npt
               do k = 1, nzi(i)
                  tmf(i,k) = tclf(i,k,1)+ clm_tscl*(tclf(i,k,2)- tclf(i,k,1))
               enddo
            enddo
         else
            do i = 1, npt
               do k = 1, nzi(i)
                  tmf(i,k) = tclf(i,k,1)
               enddo
            enddo
         endif
         if (nz_fix.le.nz) then
            n0 = 0
            n1 = npt 
            do k = nz_fix, nz
               tp(n0+k) = 0.
               tp(n1+k) = 0.
            enddo
            do i = 1, npt
               do k = nz_fix, nzi(i)
                  tp(n0+k) = tp(n0+k) + tmf(i,k)
                  tp(n1+k) = tp(n1+k) + 1.
               enddo
            enddo
            do k = nz_fix, nz
               tp(n0+k) = tp(n0+k)/tp(n1+k)
            enddo
            do i = 1, npt
               do k = nz_fix, nzi(i)
                  tmf(i,k) = tp(n0+k)
               enddo
            enddo
         endif
      elseif (initt .eq. 4) then
c.......linearly decreasing temps with longitude
         xw = xm(1)
         xe = xm(nxp)
         xl = xe-xw
         do j = 1, nxp
            tp(j) = 0. + 30.*(xe-xm(j))/xl
         enddo

         do k = 1, npt
            j = (iox(k)-1)/nxp + 1
            i = iox(k) - (j-1)*nxp
            do n = 1, nzi(k)
               tmf(k,n) = tp(i)
            enddo
         enddo
      elseif (initt .eq. 5) then
         do k = 1, npt
            j = (iox(k)-1)/nxp + 1
            i = iox(k) - (j-1)*nxp
c            if (i.gt.nxp-2.and.j.gt.nyp-2) then
            if (j.gt.nyp-2) then
               do n = 1, nzi(k)
                  tmf(k,n) = t0(n)
                  tclf(k,n,1) = t0(n)
               enddo
            else
               do n = 1, nzi(k)
                  tmf(k,n) = t0(1)
                  tclf(k,n,1) = t0(1)
               enddo
            endif
         enddo

      endif

      return
      end

c    -----------------------------------------------------------------
      subroutine salt_init (npt, nz, nstart, s0, smf, sclf)
c    -----------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      include 'comm_data.h'
      dimension s0(1), smf(npt,1), sclf(npt,nz,1)
c
      if     (inits .eq. 0) then
         do i = 1, npt
            do k = 1, nzi(i)
               smf(i,k) = s0(k)
            enddo
         enddo

      elseif (inits .eq. 3) then
c.....from CLIMATOLOGY Data.
         if     (icl_ts .eq. 2) then
            do i = 1, npt
               do k = 1, nzi(i)
                  smf(i,k) = sclf(i,k,1)+ clm_tscl*(sclf(i,k,2)- sclf(i,k,1))
               enddo
            enddo
         else
            do i = 1, npt
               do k = 1, nzi(i)
                  smf(i,k) = sclf(i,k,1)
               enddo
            enddo
         endif
         if (nz_fix.le.nz) then
            n0 = 0
            n1 = npt 
            do k = nz_fix, nz
               tp(n0+k) = 0.
               tp(n1+k) = 0.
            enddo
            do i = 1, npt
               do k = nz_fix, nzi(i)
                  tp(n0+k) = tp(n0+k) + smf(i,k)
                  tp(n1+k) = tp(n1+k) + 1.
               enddo
            enddo
            do k = nz_fix, nz
               tp(n0+k) = tp(n0+k)/tp(n1+k)
            enddo
            do i = 1, npt
               do k = nz_fix, nzi(i)
                  smf(i,k) = tp(n0+k)
               enddo
            enddo
         endif
      endif
      
      return
      end

c----------------------------------------------------
      subroutine tau_init (nstart,npt, t_anom, dtaux, dtauy, 
     *                     dtaux_anom,dtauy_anom,tmonth)
c----------------------------------------------------
c..........initalize the winds according to MTAU
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      include 'comm_data.h'
      include 'comm_nao.h'

      parameter (TAUCON = 10300., API = 3.14159265, TAUINV = 1./TAUCON)

      PARAMETER(JPT=31)
      common/mlevy/ztaux,ztauy,ztot,zsol
      REAL ZTAUX(JPT),ZTAUY(JPT),ZTOT(JPT),ZSOL(JPT)

      common/winds/mtau,matau,tausc,atau,froude
      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch

      dimension dtaux(npt,1),dtauy(npt,1)
      dimension dtaux_anom(npt,2),dtauy_anom(npt,2)

      taub = tausc*TAUINV
      tauc = atau *TAUINV

      do k = 1, npt
         j = (iox(k)-1)/nxp + 1
         i = iox(k) - (j-1)*nxp
         dtaux(k,1) = xm(i)
         dtauy(k,1) = ym(j)
      enddo

      if     (mtau .eq. 1 .or. mtau .eq. 2) then
         call odb_open(idf_tx, fbtau(1:n_tau)//'.x', 0)
         call odb_open(idf_ty, fbtau(1:n_tau)//'.y', 0)

         call odb_rddm(idf_tx, 'T', ntau)
         call mem_alloc(p_ttau, ntau, 2, 'ttau')
         call odb_rdgr(idf_tx, 'T', ntau, ttau)

         call it_catch (ntau, ttau, nstart, it1, it2, tscl)
         itau = it2

         call data_on_model_grid(idf_tx, ltau, 'tau')

         call read_zt (idf_tx, ltau, npt, 1, it1, 'taux', dtaux(1,1)) 
         call read_zt (idf_tx, ltau, npt, 1, it2, 'taux', dtaux(1,2)) 
         call read_zt (idf_ty, ltau, npt, 1, it1, 'tauy', dtauy(1,1)) 
         call read_zt (idf_ty, ltau, npt, 1, it2, 'tauy', dtauy(1,2)) 


         do i = 1, npt
            dtaux(i,1) = taub * dtaux(i,1)
            dtaux(i,2) = taub * dtaux(i,2)
            dtauy(i,1) = tauc * dtauy(i,1)
            dtauy(i,2) = tauc * dtauy(i,2)
         enddo
         if (ipole .eq. 1) then
            call rot_vec (nxp,nyp,npt, iox, xm, ym, dtaux, dtauy)
            call rot_vec (nxp,nyp,npt, iox, xm, ym, dtaux(1,2), dtauy(1,2))
         endif

         if (mtau .eq. 1) then
            do i = 1, npt
               taux(i) = dtaux(i,1) + tscl * (dtaux(i,2) - dtaux(i,1))
               tauy(i) = dtauy(i,1) + tscl * (dtauy(i,2) - dtauy(i,1))
            enddo
         endif

         if (inao.ge.1) then
            call odb_open(idf_txa, fbtau_a(1:n_tnao)//'.x', 0)
            call odb_open(idf_tya, fbtau_a(1:n_tnao)//'.y', 0)

            call data_on_model_grid(idf_txa, ltau, 'tau_anom')

            call read_zt (idf_txa, ltau, npt, 1, 1, 'taux_anom', dtaux_anom) 
            call read_zt (idf_tya, ltau, npt, 1, 1, 'tauy_anom', dtauy_anom) 

            call odb_close(idf_txa)
            call odb_close(idf_tya)

            pi = 3.1415926
            tscl_a = t_anom* max(0,cos(pi*(tmonth-1.5)/6.))
            if (ipole .eq. 1) then
               call rot_vec (nxp,nyp,npt, iox, xm,ym, dtaux_anom,dtauy_anom)
            endif

            do i = 1, npt
               taux(i) = taux(i) + TAUINV*tscl_a * dtaux_anom(i,1)
               tauy(i) = tauy(i) + TAUINV*tscl_a * dtauy_anom(i,1)
            enddo
         endif

         call odb_close(idf_tx)
         call odb_close(idf_ty)

      elseif (mtau .eq. 3) then
c..........3 - annualy averaged climatology
         call odb_open(idf_tx, fbtau(1:n_tau)//'.x', 0)
         call odb_open(idf_ty, fbtau(1:n_tau)//'.y', 0)

         call odb_rddm(idf_tx, 'T', ntau)

         call data_on_model_grid(idf_tx, ltau, 'tau')

         do k = 1, ntau
            call read_zt (idf_tx, ltau, npt, 1, k, 'taux', dtaux(1,2)) 
            call read_zt (idf_ty, ltau, npt, 1, k, 'tauy', dtauy(1,2)) 

            do i = 1, npt
               dtaux(i,1) = dtaux(i,1) + dtaux(i,2)
               dtauy(i,1) = dtauy(i,1) + dtauy(i,2)
            enddo
         enddo
         if (ipole .eq. 1) then
            call rot_vec (nxp,nyp,npt, iox, xm,ym, dtaux,dtauy)
         endif

         coef = 1./real(ntau)
         do i = 1, npt
            taux(i) = taub* coef * dtaux(i,1)
            tauy(i) = tauc* coef * dtauy(i,1)
         enddo
         call odb_close(idf_tx)
         call odb_close(idf_ty)

      elseif (mtau .eq. 4) then

         open(85,file='wind.data',form='unformatted')
         READ (85) ZTAUX, ZTAUY

      elseif (mtau .eq. 5) then
c..........5 - COSINE winds
         if     (itau_cos .eq. 0) then
            do j = 1, nyp
               tp(j) = taub*cos(API*(ym(j))/80.)
            enddo
         elseif (itau_cos .eq. 1 .or. itau_cos .eq. 3) then
            y1 = ym(1)
            y2 = ym(nyp)
            do j = 1, nyp
               tp(j) = taub*cos(2.*API*( (ym(j)-y1)/(y2-y1) - 0.5))
            enddo
         elseif (itau_cos .eq. 2) then
            y1 = ym(1)
            y2 = ym(nyp)
            do j = 1, nyp
c#               tp(j) = taub*cos(API*( (ym(j)-y1)/(y2-y1) - 0.5))
               tp(j) = taub*cos(API*( (ym(j)-y1)/(y2-y1) ))
            enddo
         elseif (itau_cos .eq. 4) then
            y1 = ym(1)
            y2 = ym(nyp)
            do j = 1, nyp
               tp(j) = taub
            enddo
         endif

         do k = 1, npt
            j = (iox(k)-1)/nxp + 1
            tmpx       = tp(j)
            taux(k)    = tmpx
            dtaux(k,1) = tmpx
            tauy(k)    = tauc
            dtauy(k,1) = tauc
         enddo

         if (itau_cos .eq. 3) then
c  just for a beta plane vs. sloped f-plane problem:

            do k = 1, npt
               tmpx       = taux(k)*dept(k)/dept(1)
               taux(k)    = tmpx
               dtaux(k,1) = tmpx
            enddo
         endif
         if (ipole .eq. 1) then
            call rot_vec (nxp,nyp,npt, iox, xm, ym, taux, tauy)
            call rot_vec (nxp,nyp,npt, iox, xm, ym, dtaux, dtauy)
         endif

      endif


      return
      end

c     ------------------------------------------------------------------
      subroutine tau_lin (nstep, inao,t_anom, npt,iox,xm,ym,
     *                    taux,tauy,dtx,dty,taux_anom,tauy_anom,tmonth)
c     ------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      common/winds/mtau,matau,tausc,atau,froude

      parameter (TAUCON = 10300., API = 3.14159265, TAUINV = 1./TAUCON)

      dimension taux(1),tauy(1),dtx(npt,1),dty(npt,1),iox(1),xm(1),ym(1)
      dimension taux_anom(npt,2),tauy_anom(npt,2)

      call it_catch (ntau, ttau, nstep, it1, it2, tscl)

      taub = tausc*TAUINV
      tauc = atau *TAUINV

      if (it2 .ne. itau) then
         itau = it2

         do i = 1, npt
            dtx(i,1) = dtx(i,2)
            dty(i,1) = dty(i,2)
         enddo

         call odb_open(idf_tx, fbtau(1:n_tau)//'.x', 0)
         call odb_open(idf_ty, fbtau(1:n_tau)//'.y', 0)

         call read_zt (idf_tx, ltau, npt, 1, it2, 'taux', dtx(1,2)) 
         call read_zt (idf_ty, ltau, npt, 1, it2, 'tauy', dty(1,2)) 

         call odb_close(idf_tx)
         call odb_close(idf_ty)

         do i = 1, npt
            dtx(i,2) = taub * dtx(i,2)
            dty(i,2) = tauc * dty(i,2)
         enddo
         if (ipole .eq. 1) then
            call rot_vec (nxp,nyp,npt, iox, xm,ym, dtx(1,2),dty(1,2))
         endif
      endif
      
      do i = 1, npt
         taux(i) = dtx(i,1) + tscl * (dtx(i,2) - dtx(i,1))
         tauy(i) = dty(i,1) + tscl * (dty(i,2) - dty(i,1))
      enddo

      if (inao.ge.1) then
c         if (tmonth.gt.10.5.or.tmonth.lt.4.5) then
c            tscl_a = TAUINV* t_anom
c         else
c            tscl_a = 0.
c         endif
         pi = 3.1415926
         tscl_a = t_anom* max(0,cos(pi*(tmonth-1.5)/6.))
         do i = 1, npt
            taux(i) = taux(i) + TAUINV* tscl_a * taux_anom(i,1)
            tauy(i) = tauy(i) + TAUINV* tscl_a * tauy_anom(i,1)
         enddo
      endif
      
      return
      end

c     ------------------------------------------------------------------
      subroutine hflx_init (nstart, npt, nx, ny, temp, sstf, cldf, slrf,
     *                      nrelax, lrelax)
c     ------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      include 'comm_data.h'
      include 'comm_pbl.h'
      PARAMETER(JPT=31)
      common/mlevy/ztaux,ztauy,ztot,zsol
      REAL ZTAUX(JPT),ZTAUY(JPT),ZTOT(JPT),ZSOL(JPT)

      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch

      dimension temp(1), sstf(npt,1), cldf(npt,1), slrf(npt,1), lrelax(1)

      if (initq .eq. 0) then
c.....Haney case  Q = QCOF * (TATM - T) with Heat Transfer Coef.=30 w/m**2/k.  
         do i = 1, npt
            sstf(i,1) = TATM
            sstf(i,3) = TATM
         enddo
         return

      elseif (initq .eq. 1) then
c.....T_atm = initial T(i,1) = const(time)
         do i = 1, npt
            sstf(i,1) = temp(i)
            sstf(i,3) = temp(i)
         enddo

      elseif (initq .eq. 2) then
c.....T_atm = average(SST) = const(time)
         call odb_open(idf_sst, fbsst(1:n_sst), 0)
         call odb_rddm(idf_sst, 'T', nsst)

         call data_on_model_grid(idf_sst, lsst, 'sst') 

         do k = 1, nsst
            call read_zt (idf_sst, lsst, npt, 1, k, 'sst', sstf(1,2)) 

            do i = 1, npt
               sstf(i,1) = sstf(i,1) + sstf(i,2)
            enddo
         enddo

         call odb_close(idf_sst)

         coef = 1./real(nsst)
         do i = 1, npt
            sstf(i,1) = coef * sstf(i,1)
            sstf(i,3) = sstf(i,1)
         enddo

      elseif (initq .eq. 3) then
c.....T_atm = SST(time) - climatology
         call odb_open(idf_sst, fbsst(1:n_sst), 0)

         call odb_rddm(idf_sst, 'T', nsst)

         call mem_alloc(p_tsst, nsst, 2, 'tsst')
         call odb_rdgr(idf_sst, 'T', nsst, tsst)

         call it_catch (nsst, tsst, nstart, it1, it2, tscl)
         isst = it2

         call data_on_model_grid(idf_sst, lsst, 'sst') 

         call read_zt (idf_sst, lsst, npt, 1, it1, 'sst', sstf(1,1)) 
         call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', sstf(1,2)) 

         if (use_wnsp) then
            call odb_open(idf_wsp, fwsp(1:n_wsp), 0)
            call odb_rddm(idf_wsp, 'T', nwsp)
            call read_zt (idf_wsp,lsst,npt,1, it1, 'wndspd', wnsp(1,1)) 
            call read_zt (idf_wsp,lsst,npt,1, it2, 'wndspd', wnsp(1,2)) 
            call odb_close(idf_wsp)
         endif

         call odb_close(idf_sst)

      elseif (initq .eq. 4) then
         open(86,file='hflx.data',form='unformatted')
         READ (86) ZTOT,ZSOL

      elseif (initq .eq. 5) then
c.....*second* Richard-Benno formulation: Q = Q(T, solr, wndsp, clouds)
         call odb_open(idf_sst, fbsst(1:n_sst), 0)
         call odb_open(idf_cld, fbcld(1:n_cld), 0)
         call odb_open(idf_slr, fbslr(1:n_slr), 0)

         call odb_rddm(idf_sst, 'T', nsst)
         call mem_alloc(p_tsst, nsst, 2, 'tsst')
         call odb_rdgr(idf_sst, 'T', nsst, tsst)
         call odb_rddm(idf_cld, 'T', i)
         if (i .ne. nsst) 
     *        call perror1('H.flx & cloud data should be on the same grid',1)
         call odb_rddm(idf_slr, 'T', i)
         if (i .ne. nsst) 
     *        call perror1('H.flx & Sol.Rad. data should be on the same grid',1)

         call it_catch (nsst, tsst, nstart, it1, it2, tscl)
         isst = it2

         call data_on_model_grid(idf_sst, lsst, 'sst')

         call read_zt (idf_sst, lsst, npt, 1, it1, 'sst', sstf(1,1)) 
         call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', sstf(1,2)) 
         call read_zt (idf_cld, lsst, npt, 1, it1, 'cld', cldf(1,1)) 
         call read_zt (idf_cld, lsst, npt, 1, it2, 'cld', cldf(1,2)) 
         call read_zt (idf_slr, lsst, npt, 1, it1, 'solr',slrf(1,1)) 
         call read_zt (idf_slr, lsst, npt, 1, it2, 'solr',slrf(1,2)) 

         call odb_close(idf_sst)
         call odb_close(idf_cld)
         call odb_close(idf_slr)

      elseif (initq .eq. 6) then
         trans_coef_sst = 0.

      elseif (initq .eq. 7) then
c.....SST = Annual Mean SST; Cld = Annual Mean Cloud Cover:
         call odb_open(idf_sst, fbsst(1:n_sst), 0)
         call odb_open(idf_cld, fbcld(1:n_cld), 0)
         call odb_rddm(idf_cld, 'T', nsst)
         call odb_rddm(idf_cld, 'T', ncld)
         if (nsst .ne. ncld) 
     *      call perror1('SST & Cloud Cover DATA should have same grids', 1)

         call data_on_model_grid(idf_sst, lsst, 'sst')

         do k = 1, nsst
            call read_zt (idf_sst, lsst, npt, 1, k, 'sst', sstf(1,2))
            call read_zt (idf_cld, lsst, npt, 1, k, 'cld', cldf(1,2)) 
 
            do i = 1, npt
               sstf(i,1) = sstf(i,1) + sstf(i,2)
               cldf(i,1) = cldf(i,1) + cldf(i,2)
            enddo
         enddo

         coef = 1./real(nsst)
         do i = 1, npt
            sstf(i,1) = coef * sstf(i,1)
            cldf(i,1) = coef * cldf(i,1)
         enddo

         call odb_close(idf_sst)
         call odb_close(idf_cld)

      elseif (initq .eq. 8) then
c.....PBL model Q = Q(T, solr, wndsp, clouds)
         call odb_open(idf_sst, fbsst(1:n_sst), 0)
         call odb_open(idf_cld, fbcld(1:n_cld), 0)
         call odb_open(idf_slr, fbslr(1:n_slr), 0)

         call odb_open(idf_wsp, fwsp(1:n_wsp), 0)
         call odb_open(idf_uwd, fuwd(1:n_uwd), 0)
         call odb_open(idf_vwd, fvwd(1:n_vwd), 0)
         call odb_open(idf_ah,  fah(1:n_ah),   0)
         call odb_open(idf_at,  fat(1:n_at),   0)

         call odb_rddm(idf_sst, 'T', nsst)
         call odb_rddm(idf_cld, 'T', ncld)
         call odb_rddm(idf_slr, 'T', nslr)
         call odb_rddm(idf_wsp, 'T', nwsp)
         call odb_rddm(idf_uwd, 'T', nuwd)
         call odb_rddm(idf_vwd, 'T', nvwd)
         call odb_rddm(idf_ah, 'T', nah)
         call odb_rddm(idf_at, 'T', nat)

         if (initq.eq.8 .and. nslr.ne.nsst) 
     *        call perror1('Solar radiation data is not on PBL grid', 1)
         if (nsst+ncld+nwsp+nuwd+nvwd+nah+nat .ne. 7*nsst)
     *        call perror1('All PBL data should be on the same grid', 1)

         call mem_alloc(p_tsst, nsst, 2, 'tsst')
         call odb_rdgr(idf_sst, 'T', nsst, tsst)

         call it_catch (nsst, tsst, nstart, it1, it2, tscl)
         isst = it2

         call data_on_model_grid(idf_sst, lsst, 'sst')

         call read_zt (idf_sst, lsst, npt, 1, it1, 'sst', sstf(1,1)) 
         call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', sstf(1,2)) 

         call read_zt (idf_cld, lsst, npt, 1, it1, 'cld', cldf(1,1)) 
         call read_zt (idf_cld, lsst, npt, 1, it2, 'cld', cldf(1,2)) 

         call read_zt (idf_slr, lsst, npt, 1, it1, 'solr',slrf(1,1)) 
         call read_zt (idf_slr, lsst, npt, 1, it2, 'solr',slrf(1,2)) 

         call read_zt (idf_wsp, lsst, npt, 1, it1, 'wndspd', wnsp(1,1)) 
         call read_zt (idf_wsp, lsst, npt, 1, it2, 'wndspd', wnsp(1,2)) 

         call read_zt (idf_uwd, lsst, npt, 1, it1, 'uwnd', uwnd(1,1)) 
         call read_zt (idf_uwd, lsst, npt, 1, it2, 'uwnd', uwnd(1,2)) 

         call read_zt (idf_vwd, lsst, npt, 1, it1, 'vwnd', vwnd(1,1)) 
         call read_zt (idf_vwd, lsst, npt, 1, it2, 'vwnd', vwnd(1,2)) 

         if (ipole .eq. 1) then
            call rot_vec (nxp,nyp,npt, iox, xm, ym, uwnd, vwnd)
            call rot_vec (nxp,nyp,npt, iox, xm, ym, uwnd(1,2), vwnd(1,2))
         endif

         NXY = nxp * nyp

         call read_zt (idf_ah,  0,    NXY, 1, it1, 'airhum', ahum(1,1)) 
         call read_zt (idf_ah,  0,    NXY, 1, it2, 'airhum', ahum(1,2)) 
         call read_zt (idf_at,  0,    NXY, 1, it1, 'airtem', atem(1,1)) 
         call read_zt (idf_at,  0,    NXY, 1, it2, 'airtem', atem(1,2)) 

         call odb_close(idf_sst)
         call odb_close(idf_cld)
         call odb_close(idf_slr)

         call odb_close(idf_wsp)
         call odb_close(idf_uwd)
         call odb_close(idf_vwd)
         call odb_close(idf_ah)
         call odb_close(idf_at)

      elseif (initq .eq. 9) then
         call odb_open(idf_sst, fbsst(1:n_sst), 0)
         call odb_rddm(idf_sst, 'T', nsst)
         call mem_alloc(p_tsst, nsst, 2, 'tsst')
         call odb_rdgr(idf_sst, 'T', nsst, tsst)
         call it_catch (nsst, tsst, nstart, it1, it2, tscl)
         isst = it2
         call data_on_model_grid(idf_sst, lsst, 'sst')
         
         call read_zt (idf_sst, lsst, npt, 1, it1, 'sst', sstf(1,1)) 
         call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', sstf(1,2)) 
         
         call odb_close(idf_sst)

         call odb_open(idf_slr, fbslr(1:n_slr), 0)
         call odb_rddm(idf_slr, 'T', nslr)
         call mem_alloc(p_tslr, nslr, 2, 'tslr')
         call odb_rdgr(idf_slr, 'T', nslr, tslr)
         call it_catch (nslr, tslr, nstart, it1, it2, tscl)
         islr = it2
         call data_on_model_grid(idf_slr, lslr, 'slr')
         
         call read_zt (idf_slr, lslr, npt, 1, it1, 'solr',slrf(1,1)) 
         call read_zt (idf_slr, lslr, npt, 1, it2, 'solr',slrf(1,2)) 
         call odb_close(idf_slr)
         
      elseif (initq .eq. 10) then

c                               Uses heat flux parameterization of
c                               Han and Oberhuber
c                                   Q(i) = (dQ/dT)*(sst(i)-t(i)) + Qobs(i)
c                                
c                               Arrays used for storage of forcing fields:
c                                   Qobs  =  slrf(1,1)
c                                   sst   =  sstf(1,1)
c                                   dQdT  =  cldf(1,1)


         call odb_open(idf_slr, fbslr(1:n_slr),0)
         call odb_open(idf_sst, fbsst(1:n_sst),0)
         call odb_open(idf_cld, fbcld(1:n_cld),0)

         call odb_rddm(idf_sst,'T',nsst)
         call mem_alloc(p_tsst,nsst,2,'tsst')
         call odb_rdgr(idf_sst,'T',nsst,tsst)

         call odb_rddm(idf_slr,'T',i)

         if (i .ne. nsst) 
     *       call perror1('H.flx & sst data should be on the same grid',1)

         call odb_rddm(idf_cld,'T',i)

         if (i .ne. nsst)
     *       call perror1('DQDT and sst data should be on same grid',1)
         
         call it_catch(nsst,tsst,nstart,it1,it2,tscl)
         isst = it2

         call data_on_model_grid(idf_sst,lsst,'sst')

         call read_zt(idf_sst,lsst,npt,1,it1,'sst',sstf(1,1))
         call read_zt(idf_sst,lsst,npt,1,it2,'sst',sstf(1,2))
         call read_zt(idf_slr,lsst,npt,1,it1,'heatflux',slrf(1,1))
         call read_zt(idf_slr,lsst,npt,1,it2,'heatflux',slrf(1,2))
         call read_zt(idf_cld,lsst,npt,1,it1,'dqdt',cldf(1,1))
         call read_zt(idf_cld,lsst,npt,1,it2,'dqdt',cldf(1,2))

         if (use_wnsp) then
            call odb_open(idf_wsp,fwsp(1:n_wsp),0)
            call odb_rddm(idf_wsp,'T',nwsp)
            call read_zt(idf_wsp,lsst,npt,1,it1,'wnsp',wnsp(1,1))
            call read_zt(idf_wsp,lsst,npt,1,it2,'wnsp',wnsp(1,2))
            call odb_close(idf_wsp)
         endif

         call odb_close(idf_slr)
         call odb_close(idf_sst)
         call odb_close(idf_cld)

      endif

         drelax = real(krelax)
         do i = 1, npt
            ixy = iox(i)
            ix = mod (ixy -1 ,nxp) + 1
            iy = (ixy - ix)/nxp + 1
            dmin = nxp+nyp
            do j = 1, nrelax
               kxy = lrelax(j)
               kx = mod (kxy -1 ,nxp) + 1
               ky = (kxy - kx)/nxp + 1
               d = sqrt(float((ix-kx)**2 + (iy-ky)**2))
               dmin = min(d,dmin)
            enddo
            relax(i) = max((drelax-dmin)/drelax,0.)
         enddo

      return
      end

c     ---------------------------------------------------
      subroutine ep_init (nstart, npt, salt, sssf, prpf)
c     ---------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      include 'comm_data.h'
      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch

      dimension salt(1), sssf(npt,1), prpf(npt,1)
      
      if (initep .eq. 0) then
c   given by EP = QCOF * (SATM - S) with Salt Transfer Coef.=30 o.e./m**2. 
         do i = 1, npt
            sssf(i,1) = SATM
            sssf(i,3) = SATM
         enddo
         return
         
      elseif (initep .eq. 1) then
c.....S_atm = initial S(i,1) = const(time)
         do i = 1, npt
            sssf(i,1) = salt(i)
            sssf(i,3) = salt(i)
         enddo

      elseif (initep .eq. 2) then
c.....S_atm = average(SSS) = const(time)
         call odb_open(idf_sss, fbsss(1:n_sss), 0)
         call odb_rddm(idf_sss, 'T', nsss)
         if (nsst .ne. nsss) 
     *      call perror1('SST & SSS DATA should be on the same T grid',1)

         do k = 1, nsst
            call read_zt (idf_sss, lsst, npt, 1, k, 'sss', sssf(1,2)) 

            do i = 1, npt
               sssf(i,1) = sssf(i,1) + sssf(i,2)
            enddo
         enddo

         coef = 1./real(nsst)
         do i = 1, npt
            sssf(i,1) = coef * sssf(i,1)
            sssf(i,3) = sssf(i,1)
         enddo
         call odb_close(idf_sss)

      elseif (initep .eq. 3) then
c.....S_atm = SSS(time) - climatology
         call odb_open(idf_sss, fbsss(1:n_sss), 0)
         call odb_rddm(idf_sss, 'T', nsss)
         if (nsst .ne. nsss) 
     *      call perror1('SST & SSS DATA should be on the same T grid',1)

         call it_catch (nsst, tsst, nstart, it1, it2, tscl)
         isss = it2
         call read_zt (idf_sss, lsst, npt, 1, it1, 'sss', sssf(1,1)) 
         call read_zt (idf_sss, lsst, npt, 1, it2, 'sss', sssf(1,2)) 
         call odb_close(idf_sss)

      elseif (initep .eq. 4) then
c.....read EP directly:
         call odb_open(idf_prp, fprp(1:n_prp), 0)

         call odb_rddm(idf_prp, 'T', nprp)
         call mem_alloc(p_tprp, nprp, 2, 'precip')
         call odb_rdgr(idf_prp, 'T', nprp, tprp)

         call it_catch (nprp, tprp, nstart, it1, it2, tscl)
         iprp = it2

         call data_on_model_grid(idf_prp, lprp, 'precip') 

         call read_zt (idf_prp, lprp, npt, 1, it1, 'precip', prpf(1,1)) 
         call read_zt (idf_prp, lprp, npt, 1, it2, 'precip', prpf(1,2)) 

         call odb_close(idf_prp)

      elseif (initep .eq. 6) then
         trans_coef_sss = 0.

      elseif (initep.eq.8) then

         call odb_open(idf_sss, fbsss(1:n_sss), 0)
         call odb_open(idf_prp, fprp(1:n_prp), 0)

         call odb_rddm(idf_sss, 'T', nsss)
         call odb_rddm(idf_prp, 'T', nprp)

         if(nprp.ne.nsst) call perror1('prp should have PBL T grid',1)
         if (nsss .ne. nsst) 
     *      call perror1('SST & SSS DATA should be on the same T grid',1)

         call it_catch (nsst, tsst, nstart, it1, it2, tscl)

         isss = it2
         call read_zt (idf_sss, lsst, npt, 1, it1, 'sss', sssf(1,1)) 
         call read_zt (idf_sss, lsst, npt, 1, it2, 'sss', sssf(1,2)) 

         iprp = it2
         call read_zt (idf_prp, lsst,npt,1,it1,'precip' ,prpf(1,1)) 
         call read_zt (idf_prp, lsst,npt,1,it2,'precip' ,prpf(1,2)) 

         call odb_close(idf_sss)
         call odb_close(idf_prp)
      endif

      return
      end

c--------------------------------------------------------------------
      subroutine qforc(nstep, npt, nx, ny, sstf, cldf, slrf, tpf, qbf) 
c---------------------------------------------------------------------
c     update heat flux using current t(1), and SST(i)
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      include 'comm_data.h'
      include 'comm_pbl.h'
      dimension slrf(npt,1),sstf(npt,1),cldf(npt,1), tpf(npt,1), qbf(npt,1)

      if     (initq .eq. 0) then
         do i = 1, npt
            q(i) = trans_coef_sst * (TATM - t(i))
         enddo

      elseif (initq .eq. 1 .or. initq .eq. 2) then
         do i = 1, npt
            q(i) = trans_coef_sst * (sstf(i,1) - t(i))
         enddo

      elseif (initq .eq. 3) then

         call it_catch (nsst, tsst, nstep, it1, it2, sst_tscl)

         if (it2 .ne. isst) then
            isst = it2
            do i = 1, npt
               sstf(i,1) = sstf(i,2)
            enddo

            call odb_open(idf_sst, fbsst(1:n_sst), 0)
            call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', sstf(1,2)) 
            call odb_close(idf_sst)

            if (use_wnsp) then
               do i = 1, npt
                  wnsp(i,1) = wnsp(i,2)
               enddo
               call odb_open(idf_wsp,fwsp(1:n_wsp),0)
               call read_zt (idf_wsp, lsst, npt, 1, it2, 'wndspd', wnsp(1,2))
               call odb_close(idf_wsp)
            endif
         endif

         do i = 1, npt
            sst_d = sstf(i,1) + sst_tscl * (sstf(i,2) - sstf(i,1))
            q(i) = trans_coef_sst * (sst_d - t(i)) 
            sstf(i,3) = sst_d
         enddo
         if (use_wnsp) then
            do i = 1, npt
               tpf(i,1)  = wnsp(i,1) + sst_tscl * (wnsp(i,2) - wnsp(i,1))
            enddo
         endif

      elseif (initq .eq. 9) then
c.....input heat flux directly from file
         call it_catch (nsst, tsst, nstep, it1, it2, tscl_sst)
         if (it2 .ne. isst) then
            isst = it2
            do i = 1, npt
               sstf(i,1) = sstf(i,2)
            enddo
            call odb_open(idf_sst, fbsst(1:n_sst), 0)
            call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', sstf(1,2)) 
            call odb_close(idf_sst)
         endif

         do i = 1, npt
            sstf(i,3) = sstf(i,1) + tscl_sst * (sstf(i,2) - sstf(i,1))
         enddo

         call it_catch (nslr, tslr, nstep, it1, it2, tscl_slr)
         if (it2 .ne. islr) then
            islr = it2
            do i = 1, npt
               slrf(i,1) = slrf(i,2)
            enddo
            call odb_open(idf_slr, fbslr(1:n_slr), 0)
            call read_zt (idf_slr, lslr, npt, 1, it2, 'solr',slrf(1,2)) 
            call odb_close(idf_slr)
         endif

         do i = 1, npt
            slrf(i,3) = slrf(i,1) + tscl_slr * (slrf(i,2) - slrf(i,1))
         enddo
         
         qcon_inv = 1./QCON
         do i = 1, npt
c...........assumes that slr data in [watts/m^2]:
            qcorr = trans_coef_sst * (sstf(i,3) - t(i))
            rx = relax(i)
            q(i)  = qcon_inv * slrf(i,3) + rx * qcorr
         enddo

      elseif (initq .eq. 5) then
c.....*new* Richard-Benno formulation using Solar Radiation & Clouds:
         call it_catch (nsst, tsst, nstep, it1, it2, tscl)
         if (it2 .ne. isst) then
            isst = it2
            do i = 1, npt
               sstf(i,1) = sstf(i,2)
               cldf(i,1) = cldf(i,2)
               slrf(i,1) = slrf(i,2)
            enddo

            call odb_open(idf_sst, fbsst(1:n_sst), 0)

            call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', sstf(1,2)) 
            call read_zt (idf_cld, lsst, npt, 1, it2, 'cld', cldf(1,2)) 
            call read_zt (idf_slr, lsst, npt, 1, it2, 'solr',slrf(1,2)) 

            call odb_close(idf_sst)
         endif

         
         do i = 1, npt
            sstf(i,3) = sstf(i,1) + tscl * (sstf(i,2) - sstf(i,1))
            cldf(i,3) = cldf(i,1) + tscl * (cldf(i,2) - cldf(i,1))
            slrf(i,3) = slrf(i,1) + tscl * (slrf(i,2) - slrf(i,1))
         enddo
         call hflx_s94(npt,t,taux,tauy,sstf(1,3),cldf(1,3),slrf(1,3),q,qr,qb)
         qcon_inv = 1./QCON
         do i = 1, npt
            qsr(i) = qr(i)/ qcon_gam
         enddo 

      elseif (initq .eq. 6) then
         qcon_gam = solr_gamma/QCON
         qcon_inv = 1./QCON
         do i = 1, npt
            q(i) = - qcon_gam * 150
            qr(i)=   qcon_gam * 150
            qsr(i) = qr(i)/ qcon_gam
         enddo

      elseif (initq .eq. 7) then

         qcon_inv = 1./QCON
         tenso = enso_start + enso_scale * nstep
         call hflx_s89(tenso,npt,iox,t,sstf,cldf,ym,taux,tauy,q,qr,qb,tpf)
         do i = 1, npt
            qsr(i) = qr(i)/ qcon_gam
         enddo 

      elseif (initq .eq. 8) then
c.....PBL model

         call it_catch (nsst, tsst, nstep, it1, it2, tscl)

         NXY = nx*ny

         if (it2 .ne. isst) then
            isst = it2

            do i = 1, npt
               sstf(i,1) = sstf(i,2)
               cldf(i,1) = cldf(i,2)
               slrf(i,1) = slrf(i,2)
               wnsp(i,1) = wnsp(i,2)
               uwnd(i,1) = uwnd(i,2)
               vwnd(i,1) = vwnd(i,2)
            enddo

            do i = 1, NXY
               ahum(i,1) = ahum(i,2)
               atem(i,1) = atem(i,2)
            enddo

            call odb_open(idf_sst, fbsst(1:n_sst), 0)
            call odb_open(idf_cld, fbcld(1:n_cld),0)
            call odb_open(idf_slr, fbslr(1:n_slr), 0)

            call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', sstf(1,2)) 
            call read_zt (idf_cld, lsst, npt, 1, it2, 'cld', cldf(1,2)) 
            call read_zt (idf_slr, lsst, npt, 1, it2, 'solr',slrf(1,2)) 

            call odb_open(idf_wsp,fwsp(1:n_wsp),0)
            call odb_open(idf_uwd, fuwd(1:n_uwd), 0)
            call odb_open(idf_vwd, fvwd(1:n_vwd), 0)

            call read_zt (idf_wsp, lsst, npt, 1, it2, 'wndspd', wnsp(1,2)) 
            call read_zt (idf_uwd, lsst, npt, 1, it2, 'uwnd', uwnd(1,2)) 
            call read_zt (idf_vwd, lsst, npt, 1, it2, 'vwnd', vwnd(1,2)) 
            if (ipole .eq. 1) then
               call rot_vec (nxp,nyp,npt, iox, xm, ym, uwnd(1,2), vwnd(1,2))
            endif

            call odb_open(idf_ah,  fah(1:n_ah),   0)
            call odb_open(idf_at,  fat(1:n_at),   0)

            call read_zt (idf_ah,  0,    NXY, 1, it2, 'airhum', ahum(1,2)) 
            call read_zt (idf_at,  0,    NXY, 1, it2, 'airtem', atem(1,2)) 

            call odb_close(idf_sst)
            call odb_close(idf_cld)
            call odb_close(idf_slr)

            call odb_close(idf_wsp)
            call odb_close(idf_uwd)
            call odb_close(idf_vwd)

            call odb_close(idf_ah)
            call odb_close(idf_at)
         endif

         if (mod(nstep, nstep_pbl).eq.0) then 
            do i = 1, npt
               sstf(i,3) = sstf(i,1) + tscl * (sstf(i,2) - sstf(i,1))
               cldf(i,3) = cldf(i,1) + tscl * (cldf(i,2) - cldf(i,1))
               tpf(i,2)  = uwnd(i,1) + tscl * (uwnd(i,2) - uwnd(i,1))
               tpf(i,3)  = vwnd(i,1) + tscl * (vwnd(i,2) - vwnd(i,1))
               slrf(i,3) = slrf(i,1) + tscl * (slrf(i,2) - slrf(i,1))
            enddo
            
            do i = 1, npt
               wnd_speed = wnsp(i,1) + tscl * (wnsp(i,2) - wnsp(i,1))
               if (wnd_speed .lt. pbl_wmin) then
                  tpf(i,1) = pbl_wmin
               else
                  tpf(i,1) = wnd_speed
               endif
            enddo

            do i = 1, NXY
               ahum(i,3) = ahum(i,1) + tscl * (ahum(i,2) - ahum(i,1))
               atem(i,3) = atem(i,1) + tscl * (atem(i,2) - atem(i,1))
            enddo

            call htflux_pbl (npt, nx, ny, iox, xm, ym, 
     *         t,cldf(1,3), tpf(1,1),tpf(1,2),tpf(1,3), ahum(1,3),atem(1,3), 
     *         qbf(1,2), qbf(1,3), qbf(1,4), amhum, amth)
         endif

         qcon_inv = 1./QCON
         qcon_gam = solr_gamma/QCON
         do i = 1, npt
            qsolr    = slrf(i,3)
            qbf(i,1) = qsolr
            qcorr = trans_coef_sst * (sstf(i,3) - t(i))
            qbf(i,5) = QCON * qcorr

c...........Total Heat Flux at the surface:
c...........Q = (1-gamma)*Q_sol - Q_lh - Q_sh - Q_lw
            qr(i) = qcon_gam * qsolr
            qtot  = qsolr - qbf(i,2) - qbf(i,3) - qbf(i,4)
            rx = relax(i)
            q(i)  = qcon_inv * qtot - qr(i) + rx * qcorr
            qsr(i) = qr(i)/ qcon_gam
         enddo
      elseif (initq .eq. 10) then

c                               Uses heat flux parameterization of
c                               Han and Oberhuber
c                                   Q(i) = (dQ/dT)*(sst(i)-t(i)) + Qobs(i)
c                                
c                               Arrays used for storage of forcing fields:
c                                   Qobs  =  slrf(1,1)
c                                   sst   =  sstf(1,1)
c                                   dQdT  =  cldf(1,1)

         call it_catch(nsst,tsst,nstep,it1,it2,sst_tscl)

         if (it2 .ne. isst) then
            isst = it2
            do i=1,npt
               sstf(i,1) = sstf(i,2)
               cldf(i,1) = cldf(i,2)
               slrf(i,1) = slrf(i,2)
            enddo
            
            call odb_open(idf_sst, fbsst(1:n_sst), 0)
            call odb_open(idf_cld, fbcld(1:n_cld),0)
            call odb_open(idf_slr, fbslr(1:n_slr), 0)

            call read_zt(idf_sst,lsst,npt,1,it2,'sst',sstf(1,2))
            call read_zt(idf_cld,lsst,npt,1,it2,'dqdt',cldf(1,2))
            call read_zt(idf_slr,lsst,npt,1,it2,'heatflux',slrf(1,2))

            if (use_wnsp) then
               do i=1,npt
                  wnsp(i,1) = wnsp(i,2)
               enddo
               call odb_open(idf_wsp,fwsp(1:n_wsp),0)
               call read_zt(idf_wsp,lsst,npt,1,it2,'wndspd',wnsp(1,2))
               call odb_close(idf_wsp)
            endif

            call odb_close(idf_sst)
            call odb_close(idf_cld)
            call odb_close(idf_slr)

         endif

            qcon_inv = 1./QCON
            do i=1,npt
               slrf(i,3) = slrf(i,1) + sst_tscl*(slrf(i,2)-slrf(i,1))
               cldf(i,3) = cldf(i,1) + sst_tscl*(cldf(i,2)-cldf(i,1))
               sstf(i,3) = sstf(i,1) + sst_tscl*(sstf(i,2)-sstf(i,1))

c               q(i) = cldf(i,3)*(sstf(i,3) - t(i)) + slrf(i,3)
               q(i) = cldf(i,3)*(t(i) - sstf(i,3)) + slrf(i,3)
               q(i) = qcon_inv * q(i)
            enddo

      endif

      return
      end

c---------------------------------------------------------
      subroutine epforc(nstep, npt, salt, sssf, prpf, qbf)
c---------------------------------------------------------
c     update EP using current sal(1), and SSS(i)
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_data.h'
      include 'comm_new.h'
      dimension salt(1),sssf(npt,1), qbf(npt,1), prpf(npt,1)
      parameter (R_MMDAY2MSEC = 1./(24. * 3600. * 1000.))
      parameter (CLATHT2EVAP = 1./(2.5e6*1028.))

      if     (initep .eq. 0) then
         do i = 1, npt
            ep(i) = trans_coef_sss * (SATM - salt(i))
         enddo

      elseif (initep .eq. 1 .or. initep .eq. 2) then
         do i = 1, npt
            ep(i) = trans_coef_sss * (sssf(i,1) - salt(i))
         enddo

      elseif (initep .eq. 3) then

         if (isss .ne. isst) then
            isss = isst
            do i = 1, npt
               sssf(i,1) = sssf(i,2)
            enddo

            call odb_open(idf_sss, fbsss(1:n_sss), 0)
            call read_zt (idf_sss, lsst, npt, 1, isss, 'sss', sssf(1,2)) 
            call odb_close(idf_sss)
         endif

         do i = 1, npt
            sss_d = sssf(i,1) + sst_tscl * (sssf(i,2) - sssf(i,1))
            ep(i) = trans_coef_sss * (sss_d - salt(i)) 
            sssf(i,3) = sss_d
         enddo

      elseif (initep .eq. 4) then
         call it_catch (nprp, tprp, nstep, it1, it2, tscl)

         if (it2 .ne. iprp) then
            iprp = it2
            do i = 1, npt
               prpf(i,1) = prpf(i,2)
            enddo

            call odb_open(idf_prp, fprp(1:n_prp), 0)
            call read_zt (idf_prp, lprp, npt, 1, it2, 'precip', prpf(1,2)) 
            call odb_close(idf_prp)
         endif
         
         do i = 1, npt
c...........assumes that prp data in [mm/day]:
            E_P = prpf(i,1) + tscl * (prpf(i,2) - prpf(i,1))
            ep(i) = R_MMDAY2MSEC * salt(i) * E_P 
         enddo

      elseif (initep.eq.8) then
         call it_catch (nsst, tsst, nstep, it1, it2, tscl)

         if (it2 .ne. iprp) then
            iprp = it2
            do i = 1, npt
               prpf(i,1) = prpf(i,2)
               sssf(i,1) = sssf(i,2)
            enddo

            call odb_open(idf_prp, fprp(1:n_prp), 0)
            call odb_open(idf_sss, fbsss(1:n_sss), 0)
            call read_zt (idf_prp, lsst, npt, 1, it2, 'precip',prpf(1,2)) 
            call read_zt (idf_sss, lsst, npt, 1, it2, 'sss', sssf(1,2)) 
            call odb_close(idf_prp)
            call odb_close(idf_sss)

         endif
         
         do i = 1, npt
            precip = R_MMDAY2MSEC *(prpf(i,1) + tscl*(prpf(i,2) - prpf(i,1)))
#ifdef gidon
            evapor = CLATHT2EVAP * qbf(i,2) * 0.7
#else
            evapor = CLATHT2EVAP * qbf(i,2) 
#endif
            sssf(i,3) = sssf(i,1) + tscl * (sssf(i,2) - sssf(i,1))
            sss_d = sssf(i,3)
            ecorr = trans_coef_sss * (sss_d - salt(i))
            rx = relax(i)
            ep(i) = (evapor - precip) * salt(i) + rx * ecorr
         enddo
      endif
      
      return
      end

c---------------------------------------------
      subroutine amlice_flux(nstep, delt, tmonth, t_anom, npt, nx, ny,
     *                    sstf, cldf, slrf, tpf, salt, sssf, prpf, qbf, flor)
c---------------------------------------------
c.....AML(Richard) + ICE(Martin,Bob)
      include 'comm_new.h'
      include 'comm_data.h'
      include 'comm_pbl.h'
      include 'amlice.h'
      include 'comm_nao.h'

      dimension slrf(npt,1),sstf(npt,1),cldf(npt,1), tpf(npt,1), qbf(npt,1)
      dimension salt(1),sssf(npt,1),prpf(npt,1)
      parameter (R_MMDAY2MSEC = 1./(24. * 3600. * 1000.))
      parameter (CLATHT2EVAP = 1./(2.5e6*1028.))
      parameter (D2SEC = 24. * 3600.)

      NXY = nx*ny

      dtpbl   = delt * real(nstep_pbl) * D2SEC
      
      call it_catch (ncld, tcld, nstep, it1, it2, tscl_cld)

      if (it2 .ne. icld) then
         icld = it2
         do i = 1, npt
            cldf(i,1) = cldf(i,2)
         enddo
         call odb_open(idf_cld, fbcld(1:n_cld), 0)
         call read_zt (idf_cld, lcld, npt, 1, it2, 'cld',cldf(1,2)) 
         call odb_close(idf_cld)
      endif

      call it_catch (nslr, tslr, nstep, it1, it2, tscl_slr)
      if (it2 .ne. islr) then
         islr = it2
         do i = 1, npt
            slrf(i,1) = slrf(i,2)
         enddo
         call odb_open(idf_slr, fbslr(1:n_slr), 0)
         call read_zt (idf_slr, lslr, npt, 1, it2, 'solr',slrf(1,2)) 
         call odb_close(idf_slr)
      endif

      call it_catch (nuwd, tuwd, nstep, it1, it2, tscl_uwd)
      if (it2 .ne. iuwd) then
         iuwd = it2
         do i = 1, npt
            uwnd(i,1) = uwnd(i,2)
            vwnd(i,1) = vwnd(i,2)
         enddo
         call odb_open(idf_uwd, fuwd(1:n_uwd), 0)
         call odb_open(idf_vwd, fvwd(1:n_vwd), 0)
         call read_zt (idf_uwd, luwd, npt, 1, it2, 'uwnd', uwnd(1,2)) 
         call read_zt (idf_vwd, luwd, npt, 1, it2, 'vwnd', vwnd(1,2)) 
         call odb_close(idf_uwd)
         call odb_close(idf_vwd)
         if (ipole .eq. 1) then
            call rot_vec (nxp,nyp,npt, iox, xm, ym, uwnd(1,2), vwnd(1,2))
         endif
         if (.not.use_wndspd) then
            do i = 1, npt
               wnsp(i,1) = wnsp(i,2)
               wnsp(i,2) = sqrt(uwnd(i,2)**2 + vwnd(i,2)**2)
            enddo
         endif
      endif
      
      if (use_wndspd) then
         call it_catch (nwsp, twsp, nstep, it1, it2, tscl_wsp)
         if (it2 .ne. iwsp) then
            iwsp = it2
            do i = 1, npt
               wnsp(i,1) = wnsp(i,2)
            enddo
            call odb_open(idf_wsp,fwsp(1:n_wsp),0)
            call read_zt (idf_wsp, lwsp, npt, 1, it2, 'wndspd',wnsp(1,2)) 
            call odb_close(idf_wsp)
         endif
      endif
      
      if (initep.eq.8) then
         call it_catch (nprp, tprp, nstep, it1, it2, tscl_prp)
         if (it2 .ne. iprp) then
            iprp = it2
            do i = 1, npt
               prpf(i,1) = prpf(i,2)
            enddo
            call odb_open(idf_prp, fprp(1:n_prp), 0)
            call read_zt (idf_prp, lprp, npt, 1, it2, 'precip', prpf(1,2)) 
            call odb_close(idf_prp)
         endif
      endif

      nxy = nx*ny

      call it_catch (nah, tah, nstep, it1, it2, tscl_ah)
      if (it2 .ne. iah) then
         iah = it2
         do i = 1, NXY
            ahum(i,1) = ahum(i,2)
         enddo
         call odb_open(idf_ah,  fah(1:n_ah),   0)
         call read_zt (idf_ah,  0, nxy, 1, it2, 'airhum', ahum(1,2)) 
         call odb_close(idf_ah)
      endif

      call it_catch (nat, tat, nstep, it1, it2, tscl_at)
      if (it2 .ne. iat) then
         iat = it2
         do i = 1, NXY
            atem(i,1) = atem(i,2)
         enddo
         call odb_open(idf_at,  fat(1:n_at),   0)
         call read_zt (idf_at,  0, nxy, 1, it2, 'airtem', atem(1,2)) 
         call odb_close(idf_at)
      endif

      call it_catch (nsst, tsst, nstep, it1, it2, tscl_sst)
      if (it2 .ne. isst) then
         isst = it2
         do i = 1, npt
            sstf(i,1) = sstf(i,2)
         enddo
         call odb_open(idf_sst, fbsst(1:n_sst), 0)
         call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', sstf(1,2)) 
         call odb_close(idf_sst)
      endif

      call it_catch (nsss, tsss, nstep, it1, it2, tscl_sss)
      if (it2 .ne. isss) then
         isss = it2
         do i = 1, npt
            sssf(i,1) = sssf(i,2)
         enddo
         call odb_open(idf_sss, fbsss(1:n_sss), 0)
         call read_zt (idf_sss, lsss, npt, 1, it2, 'sss', sssf(1,2)) 
         call odb_close(idf_sss)
      endif

      if (first_step .or. mod(nstep, nstep_pbl).eq.0) then 
         do i = 1, npt
            cldf(i,3) = cldf(i,1) + tscl_cld * (cldf(i,2) - cldf(i,1))
            slrf(i,3) = slrf(i,1) + tscl_slr * (slrf(i,2) - slrf(i,1))
            tpf(i,2)  = uwnd(i,1) + tscl_uwd * (uwnd(i,2) - uwnd(i,1))
            tpf(i,3)  = vwnd(i,1) + tscl_uwd * (vwnd(i,2) - vwnd(i,1))
            prpf(i,3) = prpf(i,1) + tscl_prp * (prpf(i,2) - prpf(i,1))
            prpf(i,3) = R_MMDAY2MSEC * prpf(i,3)
            slrf(i,3) = - slrf(i,3)/(1. - albedoocean)
            sstf(i,3) = sstf(i,1) + tscl_sst * (sstf(i,2) - sstf(i,1))
            sssf(i,3) = sssf(i,1) + tscl_sss * (sssf(i,2) - sssf(i,1))
         enddo
         
         do i = 1, npt
            wnd_speed = wnsp(i,1) + tscl_wsp * (wnsp(i,2) - wnsp(i,1))
            if (wnd_speed .lt. pbl_wmin) then
               tpf(i,1) = pbl_wmin
            else
               tpf(i,1) = wnd_speed
            endif
         enddo
         
         do i = 1, NXY
            ahum(i,3) = ahum(i,1) + tscl_ah * (ahum(i,2) - ahum(i,1))
            atem(i,3) = atem(i,1) + tscl_at * (atem(i,2) - atem(i,1))
         enddo

         if (inao.ge.1) then
            pi = 3.1415926
            tscl_a = t_anom* max(0,cos(pi*(tmonth-1.5)/6.))
            do i = 1, npt
               tpf(i,1)  = tpf(i,1) + tscl_a * wnsp_anom(i)
               tpf(i,2)  = tpf(i,2) + tscl_a * uwnd_anom(i)
               tpf(i,3)  = tpf(i,3) + tscl_a * vwnd_anom(i)
               prpf(i,3)  = prpf(i,3) + R_MMDAY2MSEC * tscl_a * prp_anom(i)
            enddo
            do i = 1, NXY
c**************************
c******WARNING  ahum + ahum_anom  should be positive
c**************************
               ahum(i,3) = max(0,ahum(i,3) + tscl_a * ahum_anom(i))
               atem(i,3) = atem(i,3) + tscl_a * atem_anom(i)
            enddo
         endif
         
         call link2htfluxice (iglob, npt, nx, ny, iox, xm, ym, dtpbl, tmonth,
     *        t,cldf(1,3), tpf(1,1),tpf(1,2),tpf(1,3),ahum(1,3),atem(1,3), 
     *        qbf(1,2), qbf(1,3), qbf(1,4), amhum, amth, rh, 
     *        sal, slrf(1,3), prpf(1,3), qbf, pp, qios, brne,
     *        trs(i_hice), trs(i_cice), trs(i_thice),
     *        tsnw, rlhi, shi, qlwi, qswi)
      endif
      
      qcon_inv = 1./QCON
      qcon_gam = solr_gamma/QCON
      do i = 1, npt
         qcorr = trans_coef_sst * (sstf(i,3) - t(i))

         qbf(i,5) = QCON * qcorr
         
         qsr(i)= -qbf(i,1)
         qr(i) = qcon_gam * qsr(i)

c...........Total Heat Flux at the surface:
c          qbf(*,1:2:3:4) = (*,qsw:rlh:sh:qlw)
         qtot  = qsr(i) - qbf(i,2) - qbf(i,3) - qbf(i,4) - qios(i)

         q(i)  = qcon_inv * qtot + qcorr - qr(i)

      enddo

      if (initep .eq. 8) then
         coef_p = coef_precip
         if (use_budget) coef_p = coef_bud
         do i = 1, npt
            precip = coef_p*pp(i)
            evapor = CLATHT2EVAP * qbf(i,2) 
            ecorr = trans_coef_sss * (sssf(i,3) - salt(i))/salt(i)
            eb = evapor - brne(i)
            ep(i) = (eb - precip + ecorr)* salt_fixed
            gridfac = flor*sxm(i)*sym(i)
            budget_re = budget_re  + ecorr* gridfac
            budget_eb = budget_eb  + eb* gridfac
            budget_pr = budget_pr  - precip* gridfac
         enddo
      elseif (initep .eq. 0) then
         do i = 1, npt
            ecorr = trans_coef_sss * (SATM - salt(i))
            eb = ecorr/salt(i)
            ep(i) = eb * salt_fixed
            gridfac = flor*sxm(i)*sym(i)
            budget_re = budget_re  + eb* gridfac
         enddo

      elseif (initep .eq. 1 .or. initep .eq. 2) then
         do i = 1, npt
            ecorr = trans_coef_sss * (sssf(i,1) - salt(i))
            eb = ecorr/salt(i)
            ep(i) = eb * salt_fixed
            gridfac = flor*sxm(i)*sym(i)
            budget_re = budget_re  + eb* gridfac
         enddo

      elseif (initep .eq. 3) then
         do i = 1, npt
            ecorr = trans_coef_sss * (sssf(i,3) - salt(i))
            eb = ecorr/salt(i)
            ep(i) = eb * salt_fixed

            gridfac = flor*sxm(i)*sym(i)
            budget_re = budget_re  + eb* gridfac
         enddo
      elseif (initep .eq. 4) then
         call it_catch (nprp, tprp, nstep, it1, it2, tscl)

         if (it2 .ne. iprp) then
            iprp = it2
            do i = 1, npt
               prpf(i,1) = prpf(i,2)
            enddo

            call odb_open(idf_prp, fprp(1:n_prp), 0)
            call read_zt (idf_prp, lprp, npt, 1, it2, 'precip', prpf(1,2)) 
            call odb_close(idf_prp)
         endif
         
         do i = 1, npt
            E_P = prpf(i,1) + tscl * (prpf(i,2) - prpf(i,1))
            ep(i) = salt(i)*E_P 
         enddo

      endif

      return
      end
            
c---------------------------------------------
      subroutine amlice_data_init(nstart, npt, nx, ny, 
     *                  temp, sstf, cldf, slrf, salt, sssf, prpf, 
     *                    nrelax, lrelax,istep_pbl)
c---------------------------------------------
      include 'comm_new.h'
      include 'comm_data.h'
      include 'comm_pbl.h'
      include 'comm_nao.h'
      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      
      dimension temp(1), sstf(npt,1), cldf(npt,1), slrf(npt,1), lrelax(1)
      dimension salt(1), sssf(npt,1), prpf(npt,1)

      istep_pbl = nstep_pbl

      call odb_open(idf_cld, fbcld(1:n_cld), 0)
      call odb_rddm(idf_cld, 'T', ncld)
      call mem_alloc(p_tcld, ncld, 2, 'tcld')
      call odb_rdgr(idf_cld, 'T', ncld, tcld)
      call it_catch (ncld, tcld, nstart, it1, it2, tscl)
      icld = it2
      call data_on_model_grid(idf_cld, lcld, 'cld')

      call read_zt (idf_cld, lcld, npt, 1, it1, 'cld', cldf(1,1)) 
      call read_zt (idf_cld, lcld, npt, 1, it2, 'cld', cldf(1,2)) 
      call odb_close(idf_cld)

      call odb_open(idf_slr, fbslr(1:n_slr), 0)
      call odb_rddm(idf_slr, 'T', nslr)
      call mem_alloc(p_tslr, nslr, 2, 'tslr')
      call odb_rdgr(idf_slr, 'T', nslr, tslr)
      call it_catch (nslr, tslr, nstart, it1, it2, tscl)
      islr = it2
      call data_on_model_grid(idf_slr, lslr, 'slr')

      call read_zt (idf_slr, lslr, npt, 1, it1, 'solr',slrf(1,1)) 
      call read_zt (idf_slr, lslr, npt, 1, it2, 'solr',slrf(1,2)) 
      call odb_close(idf_slr)
      
      call odb_open(idf_uwd, fuwd(1:n_uwd), 0)
      call odb_rddm(idf_uwd, 'T', nuwd)
      call odb_open(idf_vwd, fvwd(1:n_vwd), 0)
      call odb_rddm(idf_vwd, 'T', nvwd)
      if (nuwd.ne.nvwd) 
     *     call perror1('uwnd and vwnd should be on same T grid', 1)
      call mem_alloc(p_tuwd, nuwd, 2, 'tuwd')
      call odb_rdgr(idf_uwd, 'T', nuwd, tuwd)
      call it_catch (nuwd, tuwd, nstart, it1, it2, tscl)
      iuwd = it2
      call data_on_model_grid(idf_uwd, luwd, 'uwd')

      call read_zt (idf_uwd, luwd, npt, 1, it1, 'uwnd', uwnd(1,1)) 
      call read_zt (idf_uwd, luwd, npt, 1, it2, 'uwnd', uwnd(1,2)) 
      
      call read_zt (idf_vwd, luwd, npt, 1, it1, 'vwnd', vwnd(1,1)) 
      call read_zt (idf_vwd, luwd, npt, 1, it2, 'vwnd', vwnd(1,2)) 

      call odb_close(idf_uwd)
      call odb_close(idf_vwd)

      if (ipole .eq. 1) then
         call rot_vec (nxp,nyp,npt, iox, xm, ym, uwnd(1,1), vwnd(1,1))
         call rot_vec (nxp,nyp,npt, iox, xm, ym, uwnd(1,2), vwnd(1,2))
      endif

      if (use_wndspd) then
         call odb_open(idf_wsp, fwsp(1:n_wsp), 0)
         call odb_rddm(idf_wsp, 'T', nwsp)
         call mem_alloc(p_twsp, nwsp, 2, 'twsp')
         call odb_rdgr(idf_wsp, 'T', nwsp, twsp)
         call it_catch (nwsp, twsp, nstart, it1, it2, tscl)
         iwsp = it2
         call data_on_model_grid(idf_wsp, lwsp, 'wndspd')
         
         call read_zt (idf_wsp, lwsp, npt, 1, it1, 'wndspd', wnsp(1,1)) 
         call read_zt (idf_wsp, lwsp, npt, 1, it2, 'wndspd', wnsp(1,2)) 
         call odb_close(idf_wsp)
      else
         do i = 1, npt
            wnsp(i,1) = sqrt(uwnd(i,1)**2 + vwnd(i,1)**2)
            wnsp(i,2) = sqrt(uwnd(i,2)**2 + vwnd(i,2)**2)
         enddo
      endif
      
      if (initep.eq.8.or.initep.eq.4) then
         call odb_open(idf_prp, fprp(1:n_prp), 0)
         call odb_rddm(idf_prp, 'T', nprp)
         call mem_alloc(p_tprp, nprp, 2, 'tprp')
         call odb_rdgr(idf_prp, 'T', nprp, tprp)
         call it_catch (nprp, tprp, nstart, it1, it2, tscl)
         iprp = it2
         call data_on_model_grid(idf_prp, lprp, 'prp')

         call read_zt (idf_prp, lprp, npt, 1, it1, 'precip', prpf(1,1)) 
         call read_zt (idf_prp, lprp, npt, 1, it2, 'precip', prpf(1,2)) 
         call odb_close(idf_prp)
      endif

      nxy = nx*ny

      call odb_open(idf_ah,  fah(1:n_ah),   0)
      call odb_rddm(idf_ah, 'T', nah)
      call mem_alloc(p_tah, nah, 2, 'tah')
      call odb_rdgr(idf_ah, 'T', nah, tah)
      call it_catch (nah, tah, nstart, it1, it2, tscl)
      iah = it2

      call read_zt (idf_ah,  0, nxy, 1, it1, 'airhum', ahum(1,1)) 
      call read_zt (idf_ah,  0, nxy, 1, it2, 'airhum', ahum(1,2)) 
      call odb_close(idf_ah)

      call odb_open(idf_at,  fat(1:n_at),   0)
      call odb_rddm(idf_at, 'T', nat)
      call mem_alloc(p_tat, nat, 2, 'tat')
      call odb_rdgr(idf_at, 'T', nat, tat)
      call it_catch (nat, tat, nstart, it1, it2, tscl)
      iat = it2
      
      call read_zt (idf_at,  0, nxy, 1, it1, 'airtem', atem(1,1)) 
      call read_zt (idf_at,  0, nxy, 1, it2, 'airtem', atem(1,2)) 
      call odb_close(idf_at)

      call odb_open(idf_sst, fbsst(1:n_sst), 0)
      call odb_rddm(idf_sst, 'T', nsst)
      call mem_alloc(p_tsst, nsst, 2, 'tsst')
      call odb_rdgr(idf_sst, 'T', nsst, tsst)
      call it_catch (nsst, tsst, nstart, it1, it2, tscl)
      isst = it2
      call data_on_model_grid(idf_sst, lsst, 'sst')

      call read_zt (idf_sst, lsst, npt, 1, it1, 'sst', sstf(1,1)) 
      call read_zt (idf_sst, lsst, npt, 1, it2, 'sst', sstf(1,2)) 

      call odb_close(idf_sst)

      
      call odb_open(idf_sss, fbsss(1:n_sss), 0)
      call odb_rddm(idf_sss, 'T', nsss)
      call mem_alloc(p_tsss, nsss, 2, 'tsss')
      call odb_rdgr(idf_sss, 'T', nsss, tsss)
      call it_catch (nsss, tsss, nstart, it1, it2, tscl)
      isss = it2
      call data_on_model_grid(idf_sss, lsss, 'sss')
      
      call read_zt (idf_sss, lsss, npt, 1, it1, 'sss', sssf(1,1)) 
      call read_zt (idf_sss, lsss, npt, 1, it2, 'sss', sssf(1,2)) 
      
      call odb_close(idf_sss)
      
      if (inao.ge.1) then
         call odb_open(idf_wspa, fwspa(1:n_wspa), 0)
         call odb_open(idf_uwda, fuwda(1:n_uwda), 0)
         call odb_open(idf_vwda, fvwda(1:n_vwda), 0)
         call odb_open(idf_aha,  faha(1:n_aha),   0)
         call odb_open(idf_ata,  fata(1:n_ata),   0)
         call odb_open(idf_prpa,  fprpa(1:n_prpa),   0)
         
         call data_on_model_grid(idf_wspa, lret, 'wsp_anom')
         
         call read_zt (idf_wspa,lret,npt,1,1,'wndspd_anom',wnsp_anom) 
         call read_zt (idf_uwda,lret,npt,1,1,'uwnd_anom', uwnd_anom) 
         call read_zt (idf_vwda,lret,npt,1,1,'vwnd_anom', vwnd_anom) 
         if (ipole .eq. 1) then
            call rot_vec (nxp,nyp,npt, iox, xm, ym, uwnd_anom, vwnd_anom)
         endif
         call read_zt (idf_aha,  0, nxy, 1, 1, 'airhum_anom', ahum_anom) 
         call read_zt (idf_ata,  0, nxy, 1, 1, 'airtem_anom', atem_anom) 

         call read_zt (idf_prpa, lret, npt, 1, 1, 'precip_anom', prp_anom) 

         call odb_close(idf_wspa)
         call odb_close(idf_uwda)
         call odb_close(idf_vwda)
         call odb_close(idf_aha)
         call odb_close(idf_ata)
         call odb_close(idf_prpa)
      endif

      drelax = real(krelax)
      do i = 1, npt
         ixy = iox(i)
         ix = mod (ixy -1 ,nxp) + 1
         iy = (ixy - ix)/nxp + 1
         dmin = nxp+nyp
         do j = 1, nrelax
            kxy = lrelax(j)
            kx = mod (kxy -1 ,nxp) + 1
            ky = (kxy - kx)/nxp + 1
            d = sqrt(float((ix-kx)**2 + (iy-ky)**2))
            dmin = min(d,dmin)
         enddo
c         relax(i) = max((drelax-dmin)/drelax,0.)
         relax(i) = exp(-dmin/drelax)
      enddo

      return
      end
      
c----------------------------------------------------------------------
      subroutine it_catch (NN, tt, nstep, it1, it2, tscl)
c----------------------------------------------------------------------
c.....Returns relative shift & indexes which are bracket nstep. 
      implicit real(a-h,o-z),integer(i-n)
      dimension tt(1)
      include 'comm_new.h'

      denso = enso_start + enso_scale * nstep

      tstep = tt(2) - tt(1)

      if (tt(nn)-tt(1)+tstep .eq. 12.) then
c.....Periodic Climatology Data
c#         denso = denso - 12.*(int(denso)/12)
         denso = mod(denso, 12.)
         if (denso .lt. 0.) denso = denso + 12. 

         do it2 = 1, NN
            if (denso .lt. tt(it2)) goto 100 
         enddo
         denso = denso - 12.
         it2 = 1
  100    if (it2 .eq. 1) then 
            it1  = NN
            tscl = (12. - tt(NN) + denso)/(12. - tt(NN) + tt(1)) 
         else
            it1  = it2 - 1
            tscl = (denso - tt(it1))/(tt(it2) - tt(it1)) 
         endif
      else
c.....Non-periodic Data       
         do it2 = 1, NN
            if (denso .lt. tt(it2)) goto 200 
         enddo
         it2 = NN+1
         tscl = 0.
  200    if    (it2 .eq. 1) then
            it1  = 1
            tscl = 0.
         elseif(it2 .eq. NN+1) then
            it1  = NN
            it2  = NN
            tscl = 0.
         else
            it1  = it2 - 1
            tscl = (denso - tt(it1))/(tt(it2) - tt(it1)) 
         endif
      endif
      return
      end

c--------------------------------------------------------
      subroutine data_on_model_grid (idf, lret, tag)
c--------------------------------------------------------
      character*(*) tag
      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch

      include 'comm_new.h'
      include 'comm_data.h'

      logical grids_equiv

      lret = 1
      if (idatgr .eq. 0) then
c........check if data on the same grid as the model:
         if ( grids_equiv(idf, nxp,nyp,nxyc, nsx,nsy, tp)) lret = 0
         
c........check if data on the same grid as previous data:
      elseif ( grids_equiv(idf, mxp,myp,mpack,msx,msy, tp)) then 
         lret = 1
         
c........check if data on the same grid as the model:
      elseif ( grids_equiv(idf, nxp,nyp,nxyc, nsx,nsy, tp)) then
         lret = 0
         
      else
         write(6, *) tag, ' Only one data GRID allowed! Stop.'
         stop
      endif
      
      return
      end

c-------------------------------------------------------------------
      logical function new_grids_equiv (idf, kxp,kyp,kpack,ksx,ksy, tpp)
c-------------------------------------------------------------------
      dimension tpp(1)
      logical odb_ifatt, odb_ifvar

      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      include 'comm_new.h'
      include 'comm_data.h'

      m_x = 0
      call odb_rdgr(idf, 'X', m_x, tpp)
      xer = 1.e-6*(xm(nxp) - xm(1))  
      if (xm(1) .lt. tpp(1)-xer .or. xm(nxp) .gt. tpp(m_x)+xer) then
         write (6, *) '!!! X grid of DATA must cover the model region'
         stop
      endif

      m_sx = 0
      if ( odb_ifatt(idf, 'X', 'stretched') )
     *     call odb_getiattr(idf, 'X', 'ns', m_sx) 

      m_y = 0
      yer = 1.e-6*(ym(nyp) - ym(1))  
      call odb_rdgr(idf, 'Y', m_y, tpp(m_x+1))
      if (ym(1).lt. tpp(m_x+1)-yer .or. ym(nyp) .gt. tpp(m_x+m_y)+yer) then
         write (6, *) '!!! Y grid of DATA must cover the model region'
         stop
      endif
      m_sy = 0
      if ( odb_ifatt(idf, 'Y', 'stretched') )
     *     call odb_getiattr(idf, 'Y', 'ns', m_sy) 

      m_p   = m_x*m_y
      m_seg = 2
      tpp(1) = 1
      tpp(2) = m_p

      if (odb_ifvar(idf,'MASK')) then
         m_p   = 0
         call odb_rddm(idf, 'NPACK', m_p)
         m_seg = 0
         call odb_rddm (idf, 'NMASK', m_seg)
         call odb_rdvar(idf, 'MASK',  tpp(m_x+m_y+1))
      endif

      if (m_p.ne.kpack .or. 
     *    m_x.ne.kxp   .or. m_y.ne.kyp .or.
     *    m_sx.ne.ksx  .or. m_sy.ne.ksy ) then

         grids_equiv = .FALSE.
         if (idatgr .eq. 0) then
            mpack = m_p
            mxp   = m_x
            msx   = m_sx
            myp   = m_y
            msy   = m_sy
            mseg  = m_seg

            call datagrid_memory(tpp)
            call blin_indx(tpp)
            call blin_coef(tpp(mseg+1))
         endif
      else
         grids_equiv = .TRUE.
      endif

      return
      end

c-----------------------------------------------------------------
      subroutine read_zt (idf, key, npt, iz, it, tag, fdata) 
c-----------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension fdata(1)
      character*(*) tag
      include 'comm_data.h'
      common /odb_loc1/ IDTYPE, idf_curr, idv_curr, iok

      if (key .eq. 0) then      !! data on MODEL grid
         call odb_rd1v3 (idf, iz, it, tag, fdata)
      else                      !! data on a different grid
         call odb_rd1v3 (idf, iz, it, tag, tpdata)
         call blin_intr(npt, ixd, im2d, blcf, tpdata, fdata)
      endif

      return
      end

c-------------------------------------------------------
      integer function nearest(ixy, xd, yd, iseg)
c-------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension xd(1), yd(1), iseg(1)
      include 'comm_new.h'

      iy0 = 1 + (ixy-1)/mxp
      ix0 = ixy - mxp*(iy0-1)
      x0 = xd(ix0)
      y0 = yd(iy0)

      nearest = iseg(1)
      iy = 1       + (nearest-1)/mxp
      ix = nearest - (iy-1)*mxp
      dmin = (xd(ix)-x0)**2 + (yd(iy)-y0)**2
      
      do j = 2, mseg
         i = iseg(j)
         iy = 1 + (i-1)/mxp
         ix = i - mxp*(iy-1)
         d = (xd(ix)-x0)**2 + (yd(iy)-y0)**2
         if (d .lt. dmin) then
            dmin = d
            nearest = i
         endif
      enddo

      return
      end

c--------------------------------------
      subroutine blin_coef(iw)
c--------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension iw(1)
      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      include 'comm_data.h'
      include 'comm_new.h'
      
      call bracket(mxp, xd, nxp, xm, iw)
      call bracket(myp, yd, nyp, ym, iw(nxp+1))
      
      do k = 1, nxyc
         j = 1 + (iox(k)-1)/nxp
         i = iox(k) - (j-1)*nxp
c.....find the i,j location for the four surrounding DATA grid points.
         i1 = iw(i)
         j1 = iw(nxp+j)
c.....x-sort index of four DATA points surrounding MODEL point (i,j).
         im2d(k) = i1 + (j1-1)*mxp

c.....find the interpolation ratios.
         fx = (xm(i)-xd(i1))/(xd(i1+1)-xd(i1))
         fy = (ym(j)-yd(j1))/(yd(j1+1)-yd(j1))

         blcf(k)        = (1.-fx)*(1.-fy)   
         blcf(k+nxyc)   = fx*(1.-fy)        
         blcf(k+2*nxyc) = (1.-fx)*fy        
         blcf(k+3*nxyc) = fx*fy             
      enddo

      return
      end

c     -----------------------------------------------------------------
      subroutine bracket(nx1,x1,nx2,x2,it)
c     -----------------------------------------------------------------
c     find the elements of x1 which bracket each element of x2.
c     returns it(i), for i=1,nx2 such that:
c          x1(it(i)) .le. x2(i)  .and.  x2(i) .le. x1(it(i)+1) 
c     nx1  = (input) length of x1.
c     x1   = (input) must have x1(i+1) .gt. x1(i), i=1,nx1-1.
c     nx2  = (input) length of x2.
c     x2   = (input) must have x2(i+1) .gt. x2(i), i=1,nx2-1.
c     it   = (output) nx2 indices of the lower side of the pair of
c            consecutive elements of x1 which bracket x2(i).
c
c     must input x1(1).le.x2(1) .and. x1(nx1).ge.x2(nx2).
c
      dimension x1(1),x2(1),it(1)
c
      i1 = 1
      do 20 i2=1,nx2
   10 if(x2(i2).ge.x1(i1) .and. x2(i2).le.x1(i1+1)) goto 20
      i1 = i1 + 1
      if(i1.lt.nx1-1) goto 10
   20 it(i2) = i1
      return
c     end of bracket.
      end

c----------------------------------------------------------------
      subroutine clim_updt(npt,nz,nstep, sigma,dzin,hclf,tclf,sclf,dclf)
c----------------------------------------------------------------
      include 'comm_new.h'
      include 'comm_data.h'
      dimension hclf(npt,nz,1),tclf(npt,nz,1),sclf(npt,nz,1)
      dimension dclf(npt,nz),sigma(nz),dzin(nz+1)
      common /coords/ alat,blat,rlat,alon,blon,nsx,nsy,mgrid,z_begin,nsig,nystrch

      if (icl_ts .ne. 2) return

      call it_catch (ntclm, tclm, nstep, it1, it2, clm_tscl)
      
c.....if we are still within bracketing months:      
      if (it2 .eq. iclm) return

      iclm = it2
      
      if (icl_h .eq. 2 ) then   !!..H-clim is time dependent
         do k = 1, nz           !!...save old time
            do i = 1, npt
               hclf(i,k,1) = hclf(i,k,2)
            enddo
         enddo
         

         if (icl_htop .eq. 1) then
            call odb_open(idf_hcl, fbhcl(1:n_hcl), 0)
            call read_zt (idf_hcl, lclm,npt, 1, it2, 'mltc', hclf(1,1,2))
            call odb_close(idf_hcl)
         endif
         
         sigk = sigma(3)
         do i = 1, npt
            hclf(i,2,2) = 0.5*hclf(i,1,2) + sigk*(z_begin - 1.5*hclf(i,1,2))
         enddo
         do k = 3, nsig - 1
            sigkp = sigma(k+1)
            do i = 1, npt
               hclf(i,k,2) = (sigk+sigkp) * (z_begin - 1.5*hclf(i,1,2))
            enddo
            sigk = sigkp
         enddo
         k = nsig
         do i = 1, npt
            hclf(i,k,2) = dzin(k+1) + sigma(k)*(z_begin - 1.5*hclf(i,1,2))
         enddo
      endif
      
      do k = 1, nz              !!...save old time
         do i = 1, npt
            tclf(i,k,1) = tclf(i,k,2)
            sclf(i,k,1) = sclf(i,k,2)
         enddo
      enddo
      
      k2_h = 1
      if (icl_h .eq. 2) k2_h = 2

      call odb_open(idf_t, fbtem(1:n_tem), 0)
      call odb_open(idf_s, fbsal(1:n_sal), 0)

      call read_linz(idf_t, lclm, npt,mpack,nz,nzclm, it2,
     *               hclf(1,1,k2_h),tclf(1,1,2),zclm,tp, 'temp') 
      call read_linz(idf_s, lclm, npt,mpack,nz,nzclm, it2,
     *     hclf(1,1,k2_h),sclf(1,1,2),zclm,tp, 'salt') 

      call odb_close(idf_t)
      call odb_close(idf_s)

      if (icl_fix.eq.2) then
         n0 = 0
         n1 = npt 
         n2 = 2*npt 
         do k = nz_fix, nz
            tp(n0+k) = 0.
            tp(n1+k) = 0.
            tp(n2+k) = 0.
         enddo
         do i = 1, npt
            do k = nz_fix, nzi(i)
               tp(n0+k) = tp(n0+k) + tclf(i,k,2)
               tp(n1+k) = tp(n1+k) + sclf(i,k,2)
               tp(n2+k) = tp(n2+k) + 1.
            enddo
         enddo
         do k = nz_fix, nz
            tp(n0+k) = tp(n0+k)/tp(n2+k)
            tp(n1+k) = tp(n1+k)/tp(n2+k)
         enddo
         do i = 1, npt
            do k = nz_fix, nzi(i)
               tclf(i,k,2) = tp(n0+k)
               sclf(i,k,2) = tp(n1+k)
            enddo
         enddo
      endif
      
      return
      end


c----------------------------------------------------------------
      subroutine psi_updt(npt,nstep,psif)
c----------------------------------------------------------------
      include 'comm_new.h'
      include 'comm_data.h'
      dimension psif(npt,1)

      if (icl_psi .ne. 2) return

      call it_catch (ntpsi, tpsi, nstep, it1, it2, psi_tscl)

      if (it2 .eq. ipsi) return

      ipsi = it2
      
      do i = 1, npt
         psif(i,1) = psif(i,2)
      enddo
         
      call odb_open(idf_psi, fbpsi(1:n_psi), 0)
      call read_zt (idf_psi, lpsi,npt, 1, it2, 'psi', psif(1,2))
      call odb_close(idf_psi)
         
      return
      end

c-----------------------------------------------------------------------
      subroutine clim_relax(npt,nz,ntr,hmf,tmf,smf,umf,vmf,ub,vb,hclf,tclf,sclf,
     *                      trmf,trclf,flor,dnt,t_sp,sal_sp,tr_sp)
c-----------------------------------------------------------------------
      include 'comm_new.h'
      include 'comm_data.h'
      dimension hmf(npt,1),tmf(npt,1),smf(npt,1),
     *          hclf(npt,nz,1),tclf(npt,nz,1),sclf(npt,nz,1),trclf(npt,nz,ntr),
     *          umf(npt,nz),vmf(npt,nz),ub(npt),vb(npt),trmf(npt,nz,ntr),
     *          t_sp(npt,nz),sal_sp(npt,nz),tr_sp(npt,nz+1,ntr)

      if (icl_rlx .eq. 0) return !! NO RELAXATION - EXIT

      if     (icl_h .eq. 1) then !! H-clim time independent
         do i = 1, npt
            coef = clm_coef*sponge(i)
            do k = 1, nzi(i)
               hmf(i,k) = hmf(i,k) - coef*(hmf(i,k)-hclf(i,k,1))
            enddo
         enddo
      elseif (icl_h .eq. 2) then !! h_clim time varying
         do i = 1, npt
            coef = clm_coef*sponge(i)
            if (coef .ne. 0.) then
               do k = 1, nzi(i)
                  hcl = hclf(i,k,1) + clm_tscl*(hclf(i,k,2) - hclf(i,k,1))
                  hmf(i,k) = hmf(i,k) - coef*(hmf(i,k) - hcl)
               enddo
            endif
         enddo
      endif
      
      if     (icl_ts .eq. 1) then 
         do i = 1, npt
            coef = clm_coef*sponge(i)
            coef_v = clm_coef_v*coef
            gridfac = sxm(i)*sym(i)
            if (coef .ne. 0.) then
               do k = 1, nzi(i)
                  dtemp = - coef*(tmf(i,k)-tclf(i,k,1))
                  dsalt = - coef*(smf(i,k)-sclf(i,k,1))
                  tmf(i,k) = tmf(i,k) + dtemp
                  smf(i,k) = smf(i,k) + dsalt
                  c = flor*hmf(i,k)/dnt
                  t_sp(i,k)   = t_sp(i,k)   + c*dtemp
                  sal_sp(i,k) = sal_sp(i,k) + c*dsalt/salt_fixed
                  budget_sp = budget_sp + gridfac*c*dsalt/salt_fixed
                  umf(i,k) = (1.-coef_v)*umf(i,k)+ coef_v*ub(i)
                  vmf(i,k) = (1.-coef_v)*vmf(i,k)+ coef_v*vb(i)
               enddo
            endif
         enddo
         
      elseif (icl_ts .eq. 2) then !!...vary hmix/hthermo
         do i = 1, npt
            coef = clm_coef*sponge(i)
            coef_v = clm_coef_v*coef
            gridfac = sxm(i)*sym(i)
            if (coef .ne. 0.) then
               do k = 1, nzi(i)
                  tcl = tclf(i,k,1) + clm_tscl*(tclf(i,k,2) - tclf(i,k,1))
                  scl = sclf(i,k,1) + clm_tscl*(sclf(i,k,2) - sclf(i,k,1))
                  dtemp = - coef*(tmf(i,k)-tcl)
                  dsalt = - coef*(smf(i,k)-scl)
                  tmf(i,k) = tmf(i,k) + dtemp
                  smf(i,k) = smf(i,k) + dsalt
                  c = flor*hmf(i,k)/dnt
                  t_sp(i,k)   = t_sp(i,k)   + c*dtemp
                  sal_sp(i,k) = sal_sp(i,k) + c*dsalt/salt_fixed
                  budget_sp = budget_sp + gridfac*c*dsalt/salt_fixed
                  umf(i,k) = (1.-coef_v)*umf(i,k)+ coef_v*ub(i)
                  vmf(i,k) = (1.-coef_v)*vmf(i,k)+ coef_v*vb(i)
               enddo
            endif
         enddo
      endif

c    only annual sponges are allowed for tracer sponges

      c = flor/dnt

      do itrac = 1, ntr
         if (inittr(itrac).eq.3) then
            do i = 1, npt
               coef = clm_coef_tr*sponge(i)
               if (coef .ne. 0.) then
                  do k = 1, nzi(i)
                     dtr = - coef*(trmf(i,k,itrac)-trclf(i,k,itrac))
                     trmf(i,k,itrac) = trmf(i,k,itrac) + dtr
                     tr_sp(i,k,itrac) = tr_sp(i,k,itrac)   + c*dtr*hmf(i,k)
                  enddo
               endif
            enddo
         endif
      enddo

      return
      end
      
c-----------------------------------------------------------------------
      subroutine psi_relax(npt,pmf,psif,psib,nbx,lxx,nby,lyy)
c-----------------------------------------------------------------------
      include 'comm_new.h'
      include 'comm_data.h'
      dimension pmf(npt),psif(npt,1),psib(npt),lxx(1),lyy(1)

      if (icl_psi .eq. 0) return

c  keep boundary conditions for psi:

      do i = 1, npt
         psib(i) = pmf(i)
      enddo

      if     (icl_psi .eq. 1) then 
         do i = 1, npt
c#            coef = 1.
            coef = clm_psi*sponge(i)
            pmf(i) = pmf(i) - coef*(pmf(i)-psif(i,1))
         enddo
      elseif (icl_psi .eq. 2) then 
         do i = 1, npt
            coef = clm_psi*sponge(i)
            p = psif(i,1) + psi_tscl*(psif(i,2) - psif(i,1))
            pmf(i) = pmf(i) - coef*(pmf(i) - tcl)
         enddo
      endif

c  restore boundary conditions for psi:
      do i = 1, nbx
         ib = lxx(i)
         pmf(ib) = psib(ib)
      enddo
      do i = 1, nby
         ib = lyy(i)
         pmf(ib) = psib(ib)
      enddo

      
      return
      end
      
c---------------------------------      
      subroutine afill(n, a, v)
c---------------------------------      
      dimension a(1)
      do i = 1, n
         a(i) = v
      enddo
      return
      end


c----------------------------------------------------------------------------- 
      subroutine read_linz(idf,key,NPT,MPT,NZ,MZ,it, hdat,fdat,zvert,fvert,tag)
c----------------------------------------------------------------------------- 
      implicit real(a-h,o-z),integer(i-n)
      dimension hdat(npt,1),fdat(npt,1), zvert(1),fvert(1)
      character*(*) tag
      include 'comm_new.h'
      include 'comm_data.h'

      dimension aa(npt,1), bb(mpt,1)
      pointer   (p_aa, aa), (p_bb, bb)
      save p_aa, p_bb

      if (key .eq. 0) then
         call mem_alloc(p_aa, MZ*npt, 2, 'AA space in read_linz')

         do k = 1, MZ
            call odb_rd1v3(idf, k, it, tag, aa(1,k)) 
         enddo

         do i = 1, npt 
            do k = 1, mz
               fvert(k) = aa(i,k) 
            enddo
            
            call zlin_intrp (i, npt,nz,mz, hdat,fdat,zvert,fvert)
         enddo

         call mem_free(p_aa, MZ*npt, 2)

      else
         call mem_alloc(p_bb, MZ*mpt, 2, 'BB space in read_linz')
      
         do k = 1, MZ
            call odb_rd1v3(idf, k, it, tag, bb(1,k))
         enddo

         call zlin_blin(NPT,MPT,NZ,MZ,ixd,im2d,blcf,bb,hdat,fdat,zvert,fvert)

         call mem_free(p_bb, MZ*mpt, 2)
      endif

      return
      end

c----------------------------------------------------------------------- 
      subroutine zlin_blin(NPT,MPT,NZ,MZ,ixd,im2d,blcf,aa,h,f,zval,fval)
c----------------------------------------------------------------------- 
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      dimension ixd(1),im2d(1),blcf(npt,1)
      dimension aa(mpt,1), h(npt,1),f(npt,1), zval(1),fval(1)
      
      do i = 1, npt 
         i1 = im2d(i)
         n1 = ixd(i1)
         n2 = ixd(i1+1)
         n3 = ixd(i1+mxp)
         n4 = ixd(i1+mxp+1)

         b1 = blcf(i,1)
         b2 = blcf(i,2)
         b3 = blcf(i,3)
         b4 = blcf(i,4)
         do k = 1, mz
            fval(k) = b1*aa(n1,k) + b2*aa(n2,k) + b3*aa(n3,k) + b4*aa(n4,k) 
         enddo

         call zlin_intrp (i, npt,nz,mz, h,f,zval,fval)
      enddo

      return
      end

c------------------------------------------------------------------------
      subroutine blin_intr(npt, ixd, im2d, blcf, fd, f)
c------------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension ixd(1), im2d(1), blcf(npt,1), f(1), fd(1)
      include 'comm_new.h'

      do i = 1, npt
         i1 = im2d(i)
         f(i) = blcf(i,1)*fd(ixd(i1))     + blcf(i,2)*fd(ixd(i1+1))
     *        + blcf(i,3)*fd(ixd(i1+mxp)) + blcf(i,4)*fd(ixd(i1+mxp+1))
      enddo
      return
      end

c-------------------------------------------------------
      subroutine blin_indx (iseg)
c-------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension iseg(2,1)
      include 'comm_data.h'
      include 'comm_new.h'
      
c.....fill in ixd() according with the data compression iseg(2,mseg/2):      
      k = 0
      do i = 1, mseg/2
         do j = iseg(1,i), iseg(2,i)
            k = k + 1
            ixd(j) = k
         enddo
      enddo
c.....continue data to all points
      do k = 1, mxp*myp
         if (ixd(k) .eq. 0) ixd(k) = ixd(nearest(k, xd, yd, iseg))
      enddo

      return
      end
c-------------------------------------------------------------
      subroutine zlin_intrp(i, npt, nz, mz, h, f, zval, fval)
c-------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension h(npt,1), f(npt,1), zval(1), fval(1)
      
      zbot = zval(mz)
      fbot = fval(mz)

      mnext = 2
      scale = (fval(mnext)-fval(mnext-1))/(zval(mnext)-zval(mnext-1))
      shift = fval(mnext-1) - scale*zval(mnext-1)

      dlay = 0.5*h(i,1)
      zlay = dlay
      do k = 1, nz
         if (zlay .gt. zval(mnext)) then
            if (zlay .gt. zbot) then
               shift = fbot
               scale = 0.
            else
               do m = mnext+1, mz
                  if (zlay .le. zval(m)) then
                     scale = (fval(m) - fval(m-1))/(zval(m) - zval(m-1))
                     shift = fval(m-1) - scale*zval(m-1)                  
                     mnext = m
                     goto 100
                  endif 
               enddo
            endif
         endif
  100    f(i,k) = shift + scale * zlay
         dlay = h(i,k) - dlay
         zlay = zlay + 2.*dlay
      enddo

      return
      end



c-------------------------------------------------------------
      subroutine t_limit(npt, nzi, t_min, t_max, t)
c-------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension t(npt,1), nzi(npt)

      do i = 1, npt
         do k = 1, nzi(i)
            t(i,k) = max(t_min,t(i,k))
         enddo
      enddo

      return
      end

c-------------------------------------------------------------
      subroutine ice_limit(hicemax, npt, sst, cice, hice, thice)
c-------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension sst(npt), cice(npt), hice(npt), thice(npt)
      include 'amlice.h'


      hicemax = 1.

      do i = 1, npt
         if (hice(i).lt.1.e-3) then
            hice(i)=0.0
            cice(i)=0.0
            thice(i)=0.0
         else
            cice(i)=min(cice(i),cicemax)
            cice(i)=max(cice(i),0.1)
            if (hice(i).ge.hicemax) then
               dhice = hice(i) - hicemax
               hice(i) = hicemax
               thice(i)=  (1.-dhice/hicemax)*thice(i)
            endif
         endif
      enddo

      return
      end

c-------------------------------------------------------------
      subroutine model_limit(npt, nz, nzi, h,u,v,uc,vc)
c-------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension h(npt,nz),nzi(npt),u(npt,nz),v(npt,nz),uc(npt,nz),vc(npt,nz)

      do i = 1, npt
         do k = 1, nzi(i)
            u(i,k) = min(max(u(i,k),-1.),1.)
            v(i,k) = min(max(v(i,k),-1.),1.)
            uc(i,k) = h(i,k)*u(i,k)
            vc(i,k) = h(i,k)*v(i,k)
         enddo
      enddo

      return
      end


c-------------------------------------------------------------------
      logical function grids_equiv (idf, kxp,kyp,kpack,ksx,ksy, tpp)
c-------------------------------------------------------------------
      dimension tpp(1)
      logical odb_ifatt

      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      include 'comm_new.h'
      include 'comm_data.h'

      m_p = 0
      call odb_rddm(idf, 'NPACK', m_p)
      m_x = 0
      call odb_rdgr(idf, 'X', m_x, tpp)
      xer = 1.e-6*(xm(nxp) - xm(1))  
      if (xm(1) .lt. tpp(1)-xer .or. xm(nxp) .gt. tpp(m_x)+xer) then
         write (6, *) '!!! X grid of DATA must cover the model region'
         stop
      endif

      m_sx = 0
      if ( odb_ifatt(idf, 'X', 'stretched') )
     *     call odb_getiattr(idf, 'X', 'ns', m_sx) 

      m_y = 0
      yer = 1.e-6*(ym(nyp) - ym(1))  
      call odb_rdgr(idf, 'Y', m_y, tpp(m_x+1))
      if (ym(1).lt. tpp(m_x+1)-yer .or. ym(nyp) .gt. tpp(m_x+m_y)+yer) then
         write (6, *) '!!! Y grid of DATA must cover the model region'
         stop
      endif
      m_sy = 0
      if ( odb_ifatt(idf, 'Y', 'stretched') )
     *     call odb_getiattr(idf, 'Y', 'ns', m_sy) 

      m_seg = 0
      call odb_rddm (idf, 'NMASK', m_seg)
      call odb_rdvar(idf, 'MASK',  tpp(m_x+m_y+1))

      if (m_p.ne.kpack .or. 
     *    m_x.ne.kxp   .or. m_y.ne.kyp .or.
     *    m_sx.ne.ksx  .or. m_sy.ne.ksy ) then

         grids_equiv = .FALSE.
         if (idatgr .eq. 0) then
            mpack = m_p
            mxp   = m_x
            msx   = m_sx
            myp   = m_y
            msy   = m_sy
            mseg  = m_seg

            call datagrid_memory(tpp)
            call blin_indx(tpp)
            call blin_coef(tpp(mseg+1))
         endif
      else
         grids_equiv = .TRUE.
      endif

      return
      end


c---------------------------------------------------------
      subroutine add_anom (npt,nx,ny,nz,nzi,xm,ym,iox,h,t,anom_fac)
c---------------------------------------------------------
c INPUT
c    xm           longitude, in degrees
c    ym           latitude , in degrees
c    nx           number of x points 
c    iox          compressed->uncompressed map
c    h            layer depths
c    t            temperature profile
      
c---------------------------------------------------------
c OUTPUT
c    t            temperature profile with added anomaly
c---------------------------------------------------------

      dimension xm(nx),ym(ny),iox(npt),h(npt,nz),nzi(npt),t(npt,nz)
      real inp_flt
      parameter (aq = 8.33333e-6)
      parameter (bq = 7.5e-4)

      ianom_t  = inp_int(c_str(Anom_temp), 0)

      if (ianom_t.eq.1) then
         do i = 1, npt
            ixy = iox(i)
            ix = mod (ixy -1 ,nx) + 1
            iy = (ixy - ix)/nx + 1
            x = (xm(ix) - 200.)/5.
            y = (ym(iy) -  30.)/5.
            hanom = exp(-x**2-y**2)
            dep = 0.
            do k = 1, nzi(i)
               dep = dep + h(i,k)
               zanom  = 0.
               if (dep.lt.400) then
                  zanom = (aq*dep+bq)*(dep-400)
               endif
               t(i,k) = t(i,k) + anom_fac* zanom* hanom
            enddo
         enddo
      endif
      
      return
      end
      
c     ------------------------------------------------
      subroutine Kmap_compute (npt,nz,nxp,nyp)
c     ------------------------------------------------

c   calculate horizontally variable k's from model state

c   see Visbeck et al. JPO 27 381ff

      include 'comm_para.h'
      include 'comm_new.h'
      include 'comm_diff.h'
      include 'comm_data.h'
      include 'diffiso.h'
      include 'comm_topo.h'
      
      const = GRAVTY/(1000.+SITUD_BOT)
      do i = 1, npt
         sum_m2on = 0.
         sum_w    = 0.
         z = 0.
         do k = 1, nzi(i) 
            index = (k-1)*npt + i
            z = z + h(index)
            buoy_x = const*sigx(i,k)
            buoy_y = const*sigy(i,k)
            buoy_z = - const*sigz(i,k)
            buoy_h = sqrt(buoy_x**2 + buoy_y**2)
            buoy_z = max(1e-6,buoy_z)
            bm2on = buoy_h/sqrt(buoy_z)
            ww   = tanh(z/100.)*exp(-z/300.)*h(index)
            sum_m2on = sum_m2on + bm2on*ww
            sum_w    = sum_w    + ww
         enddo
            xl = 250.   !km
            wm2on = 0.015*(xl*1000.)**2*sum_m2on/sum_w
            wm2on = min(1500.,wm2on)
            wm2on = max( 200.,wm2on)
            Kmap(i) = wm2on
         enddo
         
c  smooth Kmap in x and y
c         goto 100
         do i = 0, nxp+1
            do j = 0, nyp+1
               t1(i,j) = 200.
               t2(i,j) = 200.
            enddo
         enddo
         do ipt = 1, npt
            i = inx(ipt)
            j = iny(ipt)
            t1(i,j) = Kmap(ipt)
         enddo
         if (iperx) then
            do j = 1, ny
               t1(0,j)     = t1(nxp,j)
               t1(nxp+1,j) = t1(1,j)
            enddo
         endif
         do j = 1, nyp
            do i = 1, nxp
               ip = i+1
               im = i-1
               t2(i,j)= 0.25* (t1(ip,j) + 2.*t1(i,j) + t1(im,j))
            enddo
         enddo
         do j = 1, nyp
            do i = 1, nxp
               jp = j+1
               jm = j-1
               t1(i,j)= 0.25* (t2(i,jp) + 2.*t2(i,j) + t2(i,jm))
            enddo
         enddo
         do ipt = 1, npt
            i = inx(ipt)
            j = iny(ipt)
            Kmap(ipt) = t1(i,j)
         enddo
  100    continue
         
         return
         end
      
