ice_itd_linear.f90

Go to the documentation of this file.

!/===========================================================================/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

!=======================================================================
!BOP
!
! !MODULE: ice_itd_linear - linear remapping scheme for ITD
!
! !DESCRIPTION:
!
! Linear remapping scheme for the ice thickness distribution
!
! See Lipscomb, W. H.  Remapping the thickness distribution of sea  \!     ice. 2001, J. Geophys. Res., Vol 106, 13989--14000.
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
!          Elizabeth C. Hunke, LANL
!
! Vectorized by Clifford Chen (Fujitsu) and William Lipscomb (LANL)
!
! !INTERFACE:
!
      module ice_itd_linear
!
! !USES:
!
      use ice_model_size
      use ice_kinds_mod
      use ice_domain
      use ice_constants
      use ice_state
      use ice_itd
      use ice_calendar
      use ice_fileunits
!
!EOP
!

      implicit none

!=======================================================================

      contains

!=======================================================================
!BOP
!
! !IROUTINE: linear_itd - ITD scheme that shifts ice among categories
!
! !INTERFACE:
!
      subroutine linear_itd (hicen_old, hicen)
!
! !DESCRIPTION:
!
! Ice thickness distribution scheme that shifts ice among categories. \!
! The default scheme is linear remapping, which works as follows.  See
! Lipscomb (2001) for more details. \!
!   Using the thermodynamic "velocities", interpolate to find the 
!   velocities in thickness space at the category boundaries, and 
!   compute the new locations of the boundaries.  Then for each 
!   category, compute the thickness distribution function,  g(h), 
!   between hL and hR, the left and right boundaries of the category.
!   Assume g(h) is a linear polynomial that satisfies two conditions: \!
!   (1) The ice area implied by g(h) equals aicen(n).
!   (2) The ice volume implied by g(h) equals aicen(n)*hicen(n).
!
! Given g(h), at each boundary compute the ice area and volume lying 
! between the original and new boundary locations.  Transfer area 
! and volume across each boundary in the appropriate direction, thus
! restoring the original boundaries. See Lipscomb (2001) for details. 
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
!          Elizabeth C. Hunke, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat), &
           intent(in) :: &
         hicen_old        ! starting value of hicen (m)

      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat), &
           intent(inout) :: &
         hicen            ! ice thickness for each cat        (m)
!
!EOP
!
      integer (kind=int_kind) :: & 
         i, j           &  ! horizontal indices
      ,  ni, nd         &   ! category indices
      ,  k                ! ice layer index

      real (kind=dbl_kind) :: &
         slope            &! rate of change of dhice with hice
      ,  dh0              &! change in ice thickness at h = 0
      ,  da0              &! area melting from category 1
      ,  damax            &! max allowed reduction in category 1 area
      ,  etamin, etamax   &! left and right limits of integration 
      ,  x1               &! etamax - etamin
      ,  x2               &! (etamax^2 - etamin^2) / 2
      ,  x3               &! (etamax^3 - etamin^3) / 3
      ,  wk1, wk2          ! temporary variables

      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,0:ncat) :: &
         hbnew            ! new boundary locations

      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi) :: &
         hb0            &  ! hin_max(0)
      ,  hb1              ! hin_max(1)

      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat) :: &
         dhicen          & ! thickness change for remapping (m)
      ,  g0              & ! constant coefficient in g(h)
      ,  g1              & ! linear coefficient in g(h) 
      ,  hl              & ! left end of range over which g(h) > 0
      ,  hr               ! right end of range over which g(h) > 0

      real (kind=dbl_kind), dimension(imt_local, jmt_local) :: &
         vice_init, vice_final  &! ice volume summed over categories
      ,  vsno_init, vsno_final  &! snow volume summed over categories
      ,  eice_init, eice_final  &! ice energy summed over categories
      ,  esno_init, esno_final  ! snow energy summed over categories

      ! NOTE: Third index of donor, daice, dvice should be ncat-1,
      !       except that compilers would have trouble when ncat = 1 
      integer (kind=int_kind), dimension(ilo:ihi,jlo:jhi,ncat) :: & 
         donor            ! donor category index

      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat) :: &
         daice          &  ! ice area transferred across boundary
      ,  dvice            ! ice volume transferred across boundary

      logical (kind=log_kind), dimension(ilo:ihi,jlo:jhi) :: &
         remap_flag       ! remap ITD if remap_flag(i,j) is true

      character (len=char_len) :: &
         fieldid           ! field identifier

      logical (kind=log_kind), parameter :: &
         l_conservation_check = .true.   ! if true, check conservation 
                                         ! (useful for debugging)

       integer (kind=int_kind) :: &
         icells         &   ! number of grid cells with ice
      ,  ij                ! combined horizontal index

       integer (kind=int_kind), dimension (1:(ihi-ilo+1)*(jhi-jlo+1)) :: &
        indxi            &  ! compressed i/j indices
      , indxj

      logical (kind=log_kind) :: &
         flag_changed

!!  ggao 60162008
      real (kind=dbl_kind) ::eps
!! change end


      !-----------------------------------------------------------------
      ! Compute volume and energy sums that linear remapping should 
      !  conserve.
      !-----------------------------------------------------------------

      if (l_conservation_check) then
         call column_sum (ncat,   vicen, vice_init)
         call column_sum (ncat,   vsnon, vsno_init)
         call column_sum (ntilay, eicen, eice_init)
         call column_sum (ncat,   esnon, esno_init)
      endif

      !-----------------------------------------------------------------
      ! Compute thickness change in each category. 
      !-----------------------------------------------------------------

      dhicen = c0i
      do ni = 1, ncat
         do j = jlo,jhi
         do i = ilo,ihi
            if (aicen(i,j,ni) > puny) then
               dhicen(i,j,ni) = hicen(i,j,ni) - hicen_old(i,j,ni)
            endif               ! aicen > puny
         enddo                  ! i
         enddo                  ! j
      enddo                     ! n

      !-----------------------------------------------------------------
      ! Compute fractional ice area in each grid cell.
      !-----------------------------------------------------------------
      call aggregate_area

      !-----------------------------------------------------------------
      ! Identify grid cells with ice and initialize remapping flag.
      ! Remapping is done wherever remap_flag = .true.
      ! In rare cases the category boundaries may shift too far for the 
      !  remapping algorithm to work, and remap_flag is set to .false.  
      ! In these cases the simpler 'rebin' subroutine will shift ice 
      !  between categories if needed.
      !-----------------------------------------------------------------

      icells = 0
      do j = jlo,jhi
      do i = ilo,ihi
         if (aice(i,j) > puny) then
            remap_flag(i,j) = .true.
            icells = icells + 1
            indxi(icells) = i
            indxj(icells) = j
         else
            remap_flag(i,j) = .false.
         endif
      enddo
      enddo

      !-----------------------------------------------------------------
      ! Compute new category boundaries, Hbnew, based on changes in
      ! ice thickness from vertical thermodynamics.
      !-----------------------------------------------------------------

      hin_max(ncat) = 999.9_dbl_kind ! arbitrary big number
      hbnew = c0i

      do ni = 1, ncat-1 

!DIR$ CONCURRENT !Cray
!cdir nodep      !NEC
!ocl novrec      !Fujitsu
         do ij = 1, icells       ! aice(i,j) > puny
            i = indxi(ij)
            j = indxj(ij)

            if (hicen_old(i,j,ni)   > puny .and.         &
                hicen_old(i,j,ni+1) > puny) then
                 ! interpolate between adjacent category growth rates
               slope = (dhicen(i,j,ni+1)-dhicen(i,j,ni)) /  & 
!                       (hicen_old(i,j,ni+1)-hicen_old(i,j,ni))
!!  ggao change 0616-2008
                       (hicen_old(i,j,ni+1)-hicen_old(i,j,ni)+epsilon(eps))
               hbnew(i,j,ni) = hin_max(ni) + dhicen(i,j,ni)  &
                            + slope * (hin_max(ni) - hicen_old(i,j,ni))
            elseif (hicen_old(i,j,ni) > puny) then ! hicen_old(ni+1)=0
               hbnew(i,j,ni) = hin_max(ni) + dhicen(i,j,ni)
            elseif (hicen_old(i,j,ni+1) > puny) then ! hicen_old(ni)=0
               hbnew(i,j,ni) = hin_max(ni) + dhicen(i,j,ni+1)
            else
               hbnew(i,j,ni) = hin_max(ni)
            endif
         enddo                  ! ij

      !-----------------------------------------------------------------
      ! Check that each boundary lies between adjacent values of hicen.
      ! If not, set remap_flag = .false.
      !-----------------------------------------------------------------

         flag_changed = .false.
         do ij = 1, icells       ! aice(i,j) > puny
            i = indxi(ij)
            j = indxj(ij)

            if (aicen(i,j,ni) > puny .and.            &
                hicen(i,j,ni) >= hbnew(i,j,ni)) then
               remap_flag(i,j) = .false.
               flag_changed = .true.
            elseif (aicen(i,j,ni+1) > puny .and.      &
                    hicen(i,j,ni+1) <= hbnew(i,j,ni)) then
               remap_flag(i,j) = .false.
               flag_changed = .true.
            endif

      !-----------------------------------------------------------------
      ! Check that Hbnew(n) lies between hin_max(n-1) and hin_max(n+1).
      ! If not, set remap_flag = .false.
      ! (In principle we could allow this, but it would make the code 
      ! more complicated.)
      !-----------------------------------------------------------------

            if (hbnew(i,j,ni) > hin_max(ni+1)) then
               remap_flag(i,j) = .false.
               flag_changed = .true.
            endif

            if (hbnew(i,j,ni) < hin_max(ni-1)) then
               remap_flag(i,j) = .false.
               flag_changed = .true.
            endif

         enddo                  ! ij

      !-----------------------------------------------------------------
      ! Write diagnosis outputs if remap_flag was changed to false
      !-----------------------------------------------------------------

         if (flag_changed) then
            do ij = 1, icells   ! aice(i,j) > puny

               i = indxi(ij)
               j = indxj(ij)

               if (aicen(i,j,ni) > puny .and.                  &
                   hicen(i,j,ni) >= hbnew(i,j,ni)) then         
                  write(nu_diag,*) 'istep1 = ',istep1           
                  write(nu_diag,*) my_task,':',i,j,            &
                       'ITD: hicen(ni) > Hbnew(ni)'             
                  write(nu_diag,*) 'cat ',ni                    
                  write(nu_diag,*) my_task,':',i,j,            &
                       'hicen(ni) =', hicen(i,j,ni)             
                  write(nu_diag,*) my_task,':',i,j,            &
                       'Hbnew(ni) =', hbnew(i,j,ni)             
               elseif (aicen(i,j,ni+1) > puny .and.            &
                       hicen(i,j,ni+1) <= hbnew(i,j,ni)) then   
                  write(nu_diag,*) 'istep1 = ',istep1           
                  write(nu_diag,*) my_task,':',i,j,            &
                       'ITD: hicen(ni+1) < Hbnew(ni)'           
                  write(nu_diag,*) 'cat ',ni                     
                  write(nu_diag,*) my_task,':',i,j,            &
                       'hicen(ni+1) =', hicen(i,j,ni+1)         
                  write(nu_diag,*) my_task,':',i,j,            &
                       'Hbnew(ni) =', hbnew(i,j,ni)             
               endif                                            
                                                                
               if (hbnew(i,j,ni) > hin_max(ni+1)) then            
                  write(nu_diag,*) 'istep1 = ',istep1           
                  write(nu_diag,*) my_task,':',i,j,            &
                       'ITD Hbnew(ni) > hin_max(ni+1)'
                  write(nu_diag,*) 'cat ',ni
                  write(nu_diag,*) my_task,':',i,j,           &
                       'Hbnew(ni) =', hbnew(i,j,ni)              
                  write(nu_diag,*) my_task,':',i,j,           &
                       'hin_max(ni+1) =', hin_max(ni+1)          
               endif                                           
                                                               
               if (hbnew(i,j,ni) < hin_max(ni-1)) then           
                  write(nu_diag,*) 'istep1 = ',istep1          
                  write(nu_diag,*) my_task,':',i,j,           &
                       'ITD: Hbnew(ni) < hin_max(ni-1)'          
                  write(nu_diag,*) 'cat ',ni                    
                  write(nu_diag,*) my_task,':',i,j,           &
                       'Hbnew(ni) =', hbnew(i,j,ni)              
                  write(nu_diag,*) my_task,':',i,j,           &
                       'hin_max(ni-1) =', hin_max(ni-1)
               endif

            enddo               ! ij
         endif                  ! flag_changed

      enddo                     ! boundaries, 1 to ncat-1

      !-----------------------------------------------------------------
      ! Identify cells where the ITD is to be remapped
      !-----------------------------------------------------------------

      icells = 0
      do j = jlo,jhi
      do i = ilo,ihi
         if (remap_flag(i,j)) then
            icells = icells + 1
            indxi(icells) = i
            indxj(icells) = j
         endif
      enddo
      enddo

      !-----------------------------------------------------------------
      ! Fill arrays with initial boundaries of category 1
      ! Prescribe Hbnew(0) and Hbnew(ncat)
      !-----------------------------------------------------------------

      do j = jlo,jhi
      do i = ilo,ihi
         hb0(i,j) = hin_max(0)
         hb1(i,j) = hin_max(1)

         hbnew(i,j,0) = c0i

         if (aicen(i,j,ncat) > puny) then
           hbnew(i,j,ncat) = c3i*hicen(i,j,ncat) - c2i*hbnew(i,j,ncat-1)
         else
            hbnew(i,j,ncat) = hin_max(ncat)
         endif

         if (hbnew(i,j,ncat) < hin_max(ncat-1))   &
              hbnew(i,j,ncat) = hin_max(ncat-1)
      enddo
      enddo

      !-----------------------------------------------------------------
      ! Compute g(h) for category 1 at start of time step 
      ! (hicen = hicen_old)
      !-----------------------------------------------------------------

      call fit_line(1,         hb0,       hb1,       hicen_old(:,:,1), &
                    g0(:,:,1), g1(:,:,1), hl(:,:,1), hr(:,:,1),        &
                    remap_flag)
 
      !-----------------------------------------------------------------
      ! Find area lost due to melting of thin (category 1) ice
      !-----------------------------------------------------------------

!DIR$ CONCURRENT !Cray
!cdir nodep      !NEC
!ocl novrec      !Fujitsu
      do ij = 1, icells    ! remap_flag = .true.
         i = indxi(ij)
         j = indxj(ij)

         if (aicen(i,j,1) > puny) then

            dh0 = dhicen(i,j,1)         

            if (dh0 < c0i) then   ! remove area from category 1   
               dh0 = min(-dh0,hin_max(1))   ! dh0 --> |dh0| 
    
      !-----------------------------------------------------------------
      ! Integrate g(1) from 0 to dh0 to estimate area melted
      !-----------------------------------------------------------------

               ! right integration limit (left limit = 0)
               etamax = min(dh0,hr(i,j,1)) - hl(i,j,1) 
                                              
               if (etamax > c0i) then
                  x1 = etamax
                  x2 = p5 * etamax*etamax
                  da0 = g1(i,j,1)*x2 + g0(i,j,1)*x1 ! ice area removed

               ! constrain new thickness <= hicen_old
                  damax = aicen(i,j,1)                       &
                        * (c1i-hicen(i,j,1)/hicen_old(i,j,1)) ! damax > 0
                  da0 = min(da0, damax)
      
               ! remove area, conserving volume
                  hicen(i,j,1) = hicen(i,j,1)                    &
!                               * aicen(i,j,1) / (aicen(i,j,1)-da0)
!!  gao change 0616-2008
                               * aicen(i,j,1) / (aicen(i,j,1)-da0+epsilon(eps))
                  aicen(i,j,1) = aicen(i,j,1) - da0
               endif            ! etamax > 0

            else                ! dh0 >= 0
               hbnew(i,j,0) = min(dh0,hin_max(1))  ! shift Hbnew(0) to right
            endif

         endif                  ! aicen(i,j,1) > puny
      enddo                     ! ij

      !-----------------------------------------------------------------
      ! Compute g(h) for each ice thickness category.
      !-----------------------------------------------------------------

      do ni = 1, ncat
        call fit_line(ni, hbnew(:,:,ni-1), hbnew(:,:,ni), hicen(:,:,ni),& 
                      g0(:,:,ni), g1(:,:,ni), hl(:,:,ni),    hr(:,:,ni),& 
                      remap_flag)
      enddo

      !-----------------------------------------------------------------
      ! Compute area and volume to be shifted across each boundary.
      !-----------------------------------------------------------------

      donor(:,:,:) = 0
      daice(:,:,:) = c0i
      dvice(:,:,:) = c0i

      do ni = 1, ncat-1
!DIR$ CONCURRENT !Cray
!cdir nodep      !NEC
!ocl novrec      !Fujitsu
         do ij = 1, icells   ! remap_flag = .true.
            i = indxi(ij)
            j = indxj(ij)

            if (hbnew(i,j,ni) > hin_max(ni)) then ! transfer from n to n+1

               ! left and right integration limits in eta space
               etamin = max(hin_max(ni), hl(i,j,ni)) - hl(i,j,ni)
               etamax = min(hbnew(i,j,ni), hr(i,j,ni)) - hl(i,j,ni) 
               donor(i,j,ni) = ni

            else             ! Hbnew(n) <= hin_max(n); transfer from n+1 to n

               ! left and right integration limits in eta space
               etamin = c0i                                 
               etamax = min(hin_max(ni), hr(i,j,ni+1)) - hl(i,j,ni+1)
               donor(i,j,ni) = ni+1

            endif            ! Hbnew(n) > hin_max(n)

            if (etamax > etamin) then
               x1  = etamax - etamin
               wk1 = etamin*etamin
               wk2 = etamax*etamax
               x2  = p5 * (wk2 - wk1)
               wk1 = wk1*etamin
               wk2 = wk2*etamax
               x3  = p333 * (wk2 - wk1)
               nd  = donor(i,j,ni)
               daice(i,j,ni) = g1(i,j,nd)*x2 + g0(i,j,nd)*x1        
               if (daice(i,j,ni) > c0i) then
                  dvice(i,j,ni) = g1(i,j,nd)*x3 + g0(i,j,nd)*x2 &
                               + daice(i,j,ni)*hl(i,j,nd)
               else
                  daice(i,j,ni) = c0i
                  donor(i,j,ni) = 0
               endif
            endif


            ! If daice or dvice is very small, shift no ice.

            nd = donor(i,j,ni)

!!!===============================
            IF(nd>0) THEN
!!!===============================

            if (daice(i,j,ni) < aicen(i,j,nd)*puny) then
               daice(i,j,ni) = c0i
               dvice(i,j,ni) = c0i
               donor(i,j,ni) = 0
            endif 

            if (dvice(i,j,ni) < vicen(i,j,nd)*puny) then
               daice(i,j,ni) = c0i
               dvice(i,j,ni) = c0i
               donor(i,j,ni) = 0
            endif

            ! If daice is close to aicen or dvice is close to vicen,
            ! shift entire category

            if (daice(i,j,ni) > aicen(i,j,nd)*(c1i-puny)) then
               daice(i,j,ni) = aicen(i,j,nd)
               dvice(i,j,ni) = vicen(i,j,nd)
            endif

            if (dvice(i,j,ni) > vicen(i,j,nd)*(c1i-puny)) then
               daice(i,j,ni) = aicen(i,j,nd)
               dvice(i,j,ni) = vicen(i,j,nd)
            endif
!!!   ggao fix the bug 02102008
            ELSE

               daice(i,j,ni) = c0i
               dvice(i,j,ni) = c0i
               donor(i,j,ni) = 0
            ENDIF
!!! change end

         enddo                  ! ij
      enddo                     ! boundaries, 1 to ncat-1

      !-----------------------------------------------------------------
      ! Shift ice between categories as necessary  
      !-----------------------------------------------------------------

      call shift_ice (donor, daice, dvice, hicen)

      !-----------------------------------------------------------------
      ! Make sure hice(i,j,1) >= minimum ice thickness hi_min.
      !-----------------------------------------------------------------

!DIR$ CONCURRENT !Cray
!cdir nodep      !NEC
!ocl novrec      !Fujitsu
      do ij = 1, icells          ! remap_flag = .true.
         i = indxi(ij)
         j = indxj(ij)
         if (hi_min > c0i .and.       &
              aicen(i,j,1) > puny .and. hicen(i,j,1) < hi_min) then
            aicen(i,j,1) = aicen(i,j,1) * hicen(i,j,1)/hi_min
            hicen(i,j,1) = hi_min
         endif
      enddo                     ! ij

      !-----------------------------------------------------------------
      ! Update ice and open water area.
      !-----------------------------------------------------------------
      call aggregate_area

      !-----------------------------------------------------------------
      ! Check volume and energy conservation.
      !-----------------------------------------------------------------

      if (l_conservation_check) then

         call column_sum (ncat,   vicen, vice_final)
         fieldid = 'vice, ITD remap'
         call column_conservation_check (vice_init, vice_final,      &
                                         puny,      fieldid)          
                                                                      
         call column_sum (ncat,   vsnon, vsno_final)                  
         fieldid = 'vsno, ITD remap'                                  
         call column_conservation_check (vsno_init, vsno_final,      &
                                         puny,      fieldid)          
                                                                      
         call column_sum (ntilay, eicen, eice_final)                  
         fieldid = 'eice, ITD remap'                                  
         call column_conservation_check (eice_init, eice_final,      &
                                         puny*lfresh*rhoi, fieldid)   
                                                                      
         call column_sum (ncat,   esnon, esno_final)                  
         fieldid = 'esno, ITD remap'                                  
         call column_conservation_check (esno_init, esno_final,      &
                                         puny*lfresh*rhos, fieldid)

      endif                     ! conservation check

      end subroutine linear_itd

!=======================================================================
!BOP
!
! !IROUTINE: fit_line - fit g(h) with a line using area, volume constraints
!
! !INTERFACE:
!
      subroutine fit_line (ni,  HbL, HbR, hice,  &
                          g0, g1,  hL,  hR, remap_flag)
!
! !DESCRIPTION:
!
! Fit g(h) with a line, satisfying area and volume constraints.
! To reduce roundoff errors caused by large values of g0 and g1,
! we actually compute g(eta), where eta = h - hL, and hL is the
! left boundary.
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
!          Elizabeth C. Hunke, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
      integer (kind=int_kind), intent(in) :: ni      ! category index

      real (kind=dbl_kind), dimension (ilo:ihi,jlo:jhi), &
           intent(in) :: & 
         hbl, hbr        ! left and right category boundaries

      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(in) :: & 
         hice            ! ice thickness

      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), intent(out) :: &
         g0, g1       &   ! coefficients in linear equation for g(eta)
      ,  hl           &   ! min value of range over which g(h) > 0
      ,  hr              ! max value of range over which g(h) > 0

      logical (kind=log_kind), dimension (ilo:ihi,jlo:jhi),      &
         intent(in) :: &
         remap_flag
!
!EOP
!
      integer (kind=int_kind) :: &
         i,j             ! horizontal indices

      real  (kind=dbl_kind) :: &
         h13           &  ! HbL + 1/3 * (HbR - HbL)
      ,  h23           &  ! HbL + 2/3 * (HbR - HbL)
      ,  dhr           &  ! 1 / (hR - hL)
      ,  wk1, wk2        ! temporary variables

!!  ggao 60162008
      real (kind=dbl_kind) ::eps
!! change end


      do j = jlo,jhi
      do i = ilo,ihi

         if (remap_flag(i,j) .and. aicen(i,j,ni) > puny    &
                      .and. hbr(i,j) - hbl(i,j) > puny) then

         ! Initialize hL and hR

            hl(i,j) = hbl(i,j)
            hr(i,j) = hbr(i,j)

         ! Change hL or hR if hicen(n) falls outside central third of range

            h13 = p333 * (c2i*hl(i,j) + hr(i,j))
            h23 = p333 * (hl(i,j) + c2i*hr(i,j))
            if (hice(i,j) < h13) then
               hr(i,j) = c3i*hice(i,j) - c2i*hl(i,j)
            elseif (hice(i,j) > h23) then
               hl(i,j) = c3i*hice(i,j) - c2i*hr(i,j)
            endif

         ! Compute coefficients of g(eta) = g0 + g1*eta

!            dhr = c1i / (hR(i,j) - hL(i,j))
!!  ggao change 0616-2008
            dhr = c1i / (hr(i,j) - hl(i,j)+epsilon(eps))
!! change end
            wk1 = c6i * aicen(i,j,ni) * dhr
            wk2 = (hice(i,j) - hl(i,j)) * dhr
            g0(i,j) = wk1 * (p666 - wk2)
            g1(i,j) = c2i*dhr * wk1 * (wk2 - p5)
               
         else                   ! remap_flag = .false. or aicen < puny
                                ! or hbR <= hbL
            hl(i,j) = c0i
            hr(i,j) = c0i
            g0(i,j) = c0i
            g1(i,j) = c0i
            
         endif                  ! aicen > puny

      enddo                     ! i
      enddo                     ! j

      end subroutine fit_line

!=======================================================================

      end module ice_itd_linear

!=======================================================================