ice_itd.f90

Go to the documentation of this file.

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

!=======================================================================
!BOP
!
! !MODULE: ice_itd - initialize and redistribute ice in the ITD
!
! !DESCRIPTION:
!
! Routines to initialize the ice thickness distribution and 
! utilities to redistribute ice among categories. These routines 
! are not specific to a particular numerical implementation. \!
! See Bitz, C.M., and W.H. Lipscomb, 1999: 
! An energy-conserving thermodynamic model of sea ice,
! J. Geophys. Res., 104, 15,669--15,677. \!     
! See Bitz, C.M., M.M. Holland, A.J. Weaver, M. Eby, 2001: 
! Simulating the ice-thickness distribution in a climate model,
! J. Geophys. Res., 106, 2441--2464. \!
! !REVISION HISTORY:
!
! author: C. M. Bitz, UW
!         Elizabeth C. Hunke, LANL
!         William H. Lipscomb, LANL
!
! Summer 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb (LANL)
!
! !INTERFACE:
!
      module ice_itd
!
! !USES:
!
      use ice_kinds_mod
      use ice_model_size
      use ice_constants
      use ice_state
      use ice_fileunits
!      use ice_exit

!
!EOP
!
      implicit none
      save

      integer (kind=int_kind) :: &
         kitd              & ! type of itd conversions 
                             !   0 = delta function
                             !   1 = linear remap
      ,  kcatbound         & !   0 = old category boundary formula
                             !   1 = new formula giving round numbers
      ,  ilyr1(ncat)      & ! position of the top layer in each cat
      ,  ilyrn(ncat)        ! position of the bottom layer in each cat

      real (kind=dbl_kind), parameter :: &
         hi_min = p01       ! minimum ice thickness allowed (m)

      real (kind=dbl_kind) :: &
         hin_max(0:ncat)    ! category limits                (m)

!-------------------------------------------------------------------
! a note regarding hi_min and hin_max(0):
! both represent a minimum ice thickness.  hin_max(0) is
! intended to be used for particular numerical implementations
! of category conversions in the ice thickness distribution.
! hi_min is a more general purpose parameter, but is specifically 
! for maintaining stability in the thermodynamics.  Currently,
! hi_min = 0.1 m
! hin_max(0) = 0.1 m for the delta function itd
! hin_max(0) = 0.0 m for linear remapping
!
! similarly, there are two values of minimum snow thickness
! (the other is defined in ice_vthermo.H since it is used only
! for thermo.) 
!
! Also note that the upper limit on the thickest category
! is only used for the linear remapping scheme
! and it is not a true upper limit on the thickness
!-------------------------------------------------------------------

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

      contains

!=======================================================================
!BOP
!
! !IROUTINE: init_itd - initalize area fraction and thickness boundaries for ITD
!
! !INTERFACE:
!
      subroutine init_itd
!
! !DESCRIPTION:
!
! Initialize area fraction and thickness boundaries for the itd model
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
!          Elizabeth C. Hunke LANL
!          C. M. Bitz UW
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
      integer (kind=int_kind) :: &
           ni    ! thickness category index

      real (kind=dbl_kind) :: &
           cc1, cc2, cc3       &        ! parameters for kcatbound = 0
      ,    x1

      real (kind=dbl_kind) :: &
           rn        &! real(n)
      ,    rncat     &! real(ncat)
      ,    d1        &                  ! parameters for kcatbound = 1 (m)
      ,    d2

      if (ncat == 1 .and. kitd == 1) then
         write (nu_diag,*) 'Remapping the ITD is not allowed for ncat=1'
         write (nu_diag,*) 'Use the delta function ITD option instead'
!         call abort_ice ('(init_itd) Linear remapping not allowed')
      endif

      rncat = real(ncat, kind=dbl_kind)
      d1 = 3.0_dbl_kind / rncat
      d2 = 0.5_dbl_kind / rncat

      !-----------------------------------------------------------------
      ! Choose category boundaries based on one of two formulas.
      !
      ! The first formula (kcatbound = 0) was used in Lipscomb (2001) 
      !  and in CICE versions 3.0 and 3.1.
      !
      ! The second formula is more user-friendly in the sense that it
      !  is easy to obtain round numbers for category boundaries:
      !
      !    H(n) = n * [d1 + d2*(n-1)] 
      ! 
      ! Default values are d1 = 300/ncat, d2 = 50/ncat.
      ! For ncat = 5, boundaries in cm are 60, 140, 240, 360, which are 
      !  close to the standard values given by the first formula.
      ! For ncat = 10, boundaries in cm are 30, 70, 120, 180, 250, 330,
      !  420, 520, 630.    
      !-----------------------------------------------------------------

      if (kcatbound == 0) then   ! original scheme

         if (kitd == 1) then
            ! linear remapping itd category limits
            cc1 = c3i/rncat
            cc2 = c15*cc1
            cc3 = c3i

            hin_max(0) = c0i     ! minimum ice thickness, m
         else
            ! delta function itd category limits
            cc1 = max(1.1_dbl_kind/rncat,c1i*hi_min)
            cc2 = c25*cc1
            cc3 = 2.25_dbl_kind

            ! hin_max(0) should not be zero
            ! use some caution in making it less than 0.10
            hin_max(0) = hi_min ! minimum ice thickness, m
         endif                  ! kitd

         do ni = 1, ncat
            x1 = real(ni-1,kind=dbl_kind) / rncat
            hin_max(ni) = hin_max(ni-1)   &
                       + cc1 + cc2*(c1i + tanh(cc3*(x1-c1i)))
         enddo

      elseif (kcatbound == 1) then  ! new scheme

         hin_max(0) = c0i
         do ni = 1, ncat
            rn = real(ni, kind=dbl_kind)
            hin_max(ni) = rn * (d1 + (rn-c1i)*d2)
         enddo

      endif

      if (my_task == master_task) then
         write (nu_diag,*) ' '
         write (nu_diag,*) 'The thickness categories are:'
         write (nu_diag,*) ' '
         write (nu_diag,*) 'hin_max(n-1) < Cat n < hin_max(n)'
         do ni = 1, ncat
         write (nu_diag,*) hin_max(ni-1),' < Cat ',ni, ' < ',hin_max(ni)
         enddo
         write (nu_diag,*) ' '
      endif

      !-----------------------------------------------------------------
      ! vectors identifying first and last layer in each category
      !-----------------------------------------------------------------
      ilyr1(1) = 1                       ! if nilyr  = 4
      ilyrn(1) = nilyr                   !   ilyr1 = { 1,5,9 }
      do ni = 2,ncat                      !   ilyrn = { 4,8,12} etc
         ilyr1(ni) = ilyrn(ni-1) + 1
         ilyrn(ni) = ilyrn(ni-1) + nilyr
      enddo

      end subroutine init_itd

!=======================================================================
!BOP
!
! !IROUTINE: aggregate - aggregate ice state variables
!
! !INTERFACE:
!
      subroutine aggregate
!
! !DESCRIPTION:
!
! Aggregate ice state variables over thickness categories.
!
! !REVISION HISTORY:
!
! authors: C. M. Bitz, UW
!          W. H. Lipscomb, LANL
!
! !USES:
!
      use ice_domain
      use ice_flux, only : tf 
      use ice_grid
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
      integer (kind=int_kind) :: i, j, k, ni 

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

      integer (kind=int_kind) :: &
        icells                 & ! number of ocean/ice cells
      , ij                    ! combined i/j horizontal index

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

      !-----------------------------------------------------------------
      ! Initialize
      !-----------------------------------------------------------------
      do j=jlo,jhi
      do i=ilo,ihi
         aice0(i,j) = c1i
         aice(i,j) = c0i
         vice(i,j) = c0i
         vsno(i,j) = c0i
         eice(i,j) = c0i
         esno(i,j) = c0i
         tsfc(i,j) = c0i
      enddo
      enddo

      !-----------------------------------------------------------------
      ! Aggregate
      !-----------------------------------------------------------------

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

      do ni = 1, ncat

!DIR$ CONCURRENT !Cray
!cdir nodep      !NEC
!ocl novrec      !Fujitsu
         do ij = 1, icells
            i = indxi(ij)
            j = indxj(ij)
            aice(i,j) = aice(i,j) + aicen(i,j,ni)
            vice(i,j) = vice(i,j) + vicen(i,j,ni)
            vsno(i,j) = vsno(i,j) + vsnon(i,j,ni)
            esno(i,j) = esno(i,j) + esnon(i,j,ni)
            tsfc(i,j) = tsfc(i,j) + tsfcn(i,j,ni)*aicen(i,j,ni)
         enddo                  ! ij

         do k = 1, nilyr
!DIR$ CONCURRENT !Cray
!cdir nodep      !NEC
!ocl novrec      !Fujitsu
            do ij = 1, icells
               i = indxi(ij)
               j = indxj(ij)
               eice(i,j) = eice(i,j) + eicen(i,j,(ni-1)*nilyr+k)
            enddo
         enddo                  ! k

      enddo                     ! n

      !-----------------------------------------------------------------
      ! Temperature, open water fraction
      !-----------------------------------------------------------------
      do j=jlo,jhi
      do i=ilo,ihi
         if (aice(i,j) > c0i)  then
            aice0(i,j) = max(c1i - aice(i,j), c0i)
!            Tsfc(i,j)  = Tsfc(i,j) / aice(i,j)
!!  ggao change 0616-2008
            tsfc(i,j)  = tsfc(i,j) /(aice(i,j)+epsilon(eps))
         else
            tsfc(i,j)  = tf(i,j)
         endif
      enddo                     ! i
      enddo                     ! j

      end subroutine aggregate

!=======================================================================
!BOP
!
! !IROUTINE: aggregate_area - aggregate ice area
!
! !INTERFACE:
!
      subroutine aggregate_area
!
! !DESCRIPTION:
!
! Aggregate ice area (but not other state variables) over thickness categories
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
!          modified Jan 2004 by Clifford Chen, Fujitsu
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!

      use ice_grid, only : tlat,tlon
!     ggao
      integer (kind=int_kind) :: i, j, ni

      logical (kind=log_kind) :: &
         outbound          ! = .true. if aggregate ice area out of bounds


      !-----------------------------------------------------------------
      ! Initialize
      !-----------------------------------------------------------------
      do j=jlo,jhi
      do i=ilo,ihi
         aice(i,j) = c0i
      enddo
      enddo

      !-----------------------------------------------------------------
      ! Aggregate
      !-----------------------------------------------------------------
      do ni = 1, ncat
         do j=jlo,jhi
         do i=ilo,ihi
            aice(i,j) = aice(i,j) + aicen(i,j,ni)
         enddo                  ! i
         enddo                  ! j
      enddo                     ! n

      outbound = .false.

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

      !-----------------------------------------------------------------
      ! Bug check
      !-----------------------------------------------------------------
         if (aice(i,j) > c1i+puny .or.  &
             aice(i,j) < -puny) then
            outbound = .true.
         endif

      !-----------------------------------------------------------------
      ! open water fraction
      !-----------------------------------------------------------------
         aice0(i,j) = max(c1i - aice(i,j), c0i)

      enddo                     ! i
      enddo                     ! j

      if ( outbound ) then      ! area out of bounds
        do j = jlo, jhi
        do i = ilo, ihi

         if (aice(i,j) > c1i+puny .or.  &
             aice(i,j) < -puny) then
!            write(nu_diag,*) ' '
!            write(nu_diag,*) 'aggregate ice area out of bounds'
!            write(nu_diag,*) 'my_task, i, j, aice:',     &
!                              my_task, i, j, aice(i,j)

!!            write(nu_diag,*) 'aicen =', aicen(i,j,:) 
!            call abort_ice('ice_itd: aggregate_area')
         endif

        enddo                   ! i
        enddo                   ! j
      endif                     ! outbound

      end subroutine aggregate_area

!=======================================================================
!BOP
!
! !IROUTINE: bound_aggregate - bound calls for aggregate ice state
!
! !INTERFACE:
!
      subroutine bound_aggregate
!
! !DESCRIPTION:
!
! Get ghost cell values for aggregate ice state variables
! NOTE: This subroutine is called only at initialization.  It could be
!       eliminated if the aggregate variables were defined only in
!       physical grid cells.
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
      use ice_grid

!      call bound (aice0)
!      call bound (aice)
!      call bound (vice)
!      call bound (vsno)
!      call bound (eice)
!      call bound (esno)
!      call bound (Tsfc)

!    need inter processors exchange 
!   ggao


      end subroutine bound_aggregate

!=======================================================================
!BOP
!
! !IROUTINE: bound_state - bound calls for ice state variables
!
! !INTERFACE:
!
      subroutine bound_state
!
! !DESCRIPTION:
!
! Get ghost cell values for ice state variables in each thickness category
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
      use ice_grid

!      call bound_narr (ncat,aicen)
!      call bound_narr (ncat,vicen)
!      call bound_narr (ncat,vsnon)
!      call bound_narr (ntilay,eicen)
!      call bound_narr (ncat,esnon)
!      call bound_narr (ncat,Tsfcn)
!    need inter processors exchange 
!   ggao


      end subroutine bound_state

!=======================================================================
!BOP
!
! !IROUTINE: check_state - require certain fields to be monotone
!
! !INTERFACE:
!
      subroutine check_state
!
! !DESCRIPTION:
!
!  Insist that certain fields are monotone.
!  Should not be necessary if all is well, 
!  but best to keep going. Model will not conserve
!  energy and water if fields are zeroed here.
!
! !REVISION HISTORY:
!
! author: C. M. Bitz, UW
!
! !USES:
!
      use ice_flux
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
      integer (kind=int_kind) :: &
         i, j         &   ! horizontal indices
      ,  ni           &    ! thickness category index
      ,  k               ! ice layer index

      integer (kind=int_kind), dimension (ilo:ihi,jlo:jhi) :: &
         zerout    ! 0=false, 1=true

      do ni = 1, ncat

         do j = jlo, jhi
         do i = ilo, ihi
            if (aicen(i,j,ni) < puny .or. vicen(i,j,ni) < puny) then
               zerout(i,j) = 1
            else
               zerout(i,j) = 0
            endif
         enddo
         enddo

         do k = 1, nilyr
            do j = jlo, jhi
            do i = ilo, ihi
               if (eicen(i,j,ilyr1(ni)+k-1) > -puny) zerout(i,j) = 1
            enddo               ! i
            enddo               ! j
         enddo                  ! k

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

            if (zerout(i,j)==1) then
               aice0(i,j) = aice0(i,j) + aicen(i,j,ni) 
               aicen(i,j,ni) = c0i
               vicen(i,j,ni) = c0i
               vsnon(i,j,ni) = c0i
               esnon(i,j,ni) = c0i
               tsfcn(i,j,ni) = tf(i,j)
            elseif (vsnon(i,j,ni) <= puny) then
               vsnon(i,j,ni) = c0i
               esnon(i,j,ni) = c0i
            endif

            if (vsnon(i,j,ni) > puny) then
               if (-esnon(i,j,ni)/vsnon(i,j,ni)-lfresh*rhos < eps04)  &
                    esnon(i,j,ni) = -vsnon(i,j,ni)*(lfresh*rhos + eps04)
            endif

         enddo                  ! i
         enddo                  ! j

         do k = 1, nilyr
            do j = jlo, jhi
            do i = ilo, ihi
               if (zerout(i,j)==1) eicen(i,j,ilyr1(ni)+k-1) = c0i
            enddo
            enddo
         enddo                  ! k

      enddo                     ! n

      do j=jlo,jhi
      do i=ilo,ihi
         if (aice0(i,j) < puny) aice0(i,j) = c0i
      enddo
      enddo

      end subroutine check_state

!=======================================================================
!BOP
!
! !IROUTINE: rebin - rebins thicknesses into defined categories
!
! !INTERFACE:
!
      subroutine rebin
!
! !DESCRIPTION:
!
! Rebins thicknesses into defined categories
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
!          Elizabeth C. Hunke LANL 
!
! !USES:
!
      use ice_grid
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
      integer (kind=int_kind) :: &
         i,j            &    ! horizontal indices
      ,  ni                  ! category index

      logical (kind=log_kind) :: &
         shiftflag          ! = .true. if ice must be shifted

      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat) :: &
         hicen              ! ice thickness for each cat        (m)

      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
      ,  dvice              ! ice volume transferred

      !-----------------------------------------------------------------
      ! Compute ice thickness.
      !-----------------------------------------------------------------
      do ni = 1, ncat
         do j = jlo, jhi
         do i = ilo, ihi
            if (aicen(i,j,ni) > puny) then
               hicen(i,j,ni) = vicen(i,j,ni) / aicen(i,j,ni)
            else
               hicen(i,j,ni) = c0i
            endif
         enddo                  ! i
         enddo                  ! j
      enddo                     ! n

      !-----------------------------------------------------------------
      ! make sure thickness of cat 1 is at least hin_max(0)
      !-----------------------------------------------------------------
      do j = jlo, jhi
      do i = ilo, ihi
         if (aicen(i,j,1) > puny) then
            if (hicen(i,j,1) <= hin_max(0) .and. hin_max(0) > c0i ) then
               aicen(i,j,1) = vicen(i,j,1) / hin_max(0)
               hicen(i,j,1) = hin_max(0)
            endif
         endif
      enddo                     ! i
      enddo                     ! j

      !-----------------------------------------------------------------
      ! If a category thickness is not in bounds, shift the
      ! entire area, volume, and energy to the neighboring category
      !-----------------------------------------------------------------

      !-----------------------------------------------------------------
      ! Initialize shift arrays
      !-----------------------------------------------------------------

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

      !-----------------------------------------------------------------
      ! Move thin categories up
      !-----------------------------------------------------------------

      do ni = 1, ncat-1          ! loop over category boundaries

      !-----------------------------------------------------------------
      ! identify thicknesses that are too big
      !-----------------------------------------------------------------
         shiftflag = .false.
         do j = jlo, jhi
         do i = ilo, ihi
            if (aicen(i,j,ni) > puny .and.  & 
                hicen(i,j,ni) > hin_max(ni)) then
               shiftflag = .true.
               donor(i,j,ni) = ni
               daice(i,j,ni) = aicen(i,j,ni)
               dvice(i,j,ni) = vicen(i,j,ni)
            endif
         enddo                  ! i
         enddo                  ! j

         if (shiftflag) then

      !-----------------------------------------------------------------
      ! shift ice between categories
      !-----------------------------------------------------------------
            call shift_ice (donor, daice, dvice, hicen)
             
      !-----------------------------------------------------------------
      ! reset shift parameters
      !-----------------------------------------------------------------
            do j = jlo, jhi
            do i = ilo, ihi
               donor(i,j,ni) = 0
               daice(i,j,ni) = c0i
               dvice(i,j,ni) = c0i
            enddo
            enddo

         endif                  ! shiftflag

      enddo                     ! n

      !-----------------------------------------------------------------
      ! Move thick categories down
      !-----------------------------------------------------------------

      do ni = ncat-1, 1, -1      ! loop over category boundaries

      !-----------------------------------------------------------------
      ! identify thicknesses that are too small
      !-----------------------------------------------------------------
         shiftflag = .false.
         do j = jlo, jhi
         do i = ilo, ihi
            if (aicen(i,j,ni+1) > puny .and.  &
                hicen(i,j,ni+1) <= hin_max(ni)) then
               shiftflag = .true.
               donor(i,j,ni) = ni+1
               daice(i,j,ni) = aicen(i,j,ni+1)
               dvice(i,j,ni) = vicen(i,j,ni+1)
            endif
         enddo                  ! i
         enddo                  ! j

         if (shiftflag) then

      !-----------------------------------------------------------------
      ! shift ice between categories
      !-----------------------------------------------------------------
            call shift_ice (donor, daice, dvice, hicen)

      !-----------------------------------------------------------------
      ! reset shift parameters
      !-----------------------------------------------------------------
            do j = jlo, jhi
            do i = ilo, ihi
               donor(i,j,ni) = 0
               daice(i,j,ni) = c0i
               dvice(i,j,ni) = c0i
            enddo
            enddo

         endif                  ! shiftflag

      enddo                     ! n


      end subroutine rebin

!=======================================================================
!BOP
!
! !IROUTINE: reduce_area - reduce area when ice melts for special case ncat=1
!
! !INTERFACE:
!
      subroutine reduce_area(hice1_old, hice1)
!
! !DESCRIPTION:
!
! Reduce area when ice melts for special case of ncat=1
!
! Use CSM 1.0-like method of reducing ice area
! when melting occurs: assume only half the ice volume
! change goes to thickness decrease, the other half
! to reduction in ice fraction
!
! !REVISION HISTORY:
!
! authors: C. M. Bitz, UW 
! modified by: Elizabeth C. Hunke, LANL
!
! !USES:
!
      use ice_grid
!
! !INPUT/OUTPUT PARAMETERS:
!
      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), &
           intent(in) :: &
         hice1_old   ! old ice thickness for category 1 (m)

      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi), &
           intent(inout) :: &
         hice1       ! new ice thickness for category 1 (m)
!
!EOP
!
      integer (kind=int_kind) :: &
         i, j        ! horizontal indices

      real (kind=dbl_kind) :: &
         hi0       & ! current hi for ice fraction adjustment
      ,  dai0      & ! change in aice for ice fraction adjustment
      ,  dhi         ! hice1 - hice1_old

      do j=jlo,jhi
      do i=ilo,ihi
         if (tmask(i,j)) then

      !-----------------------------------------------------------------
      ! make sure thickness of cat 1 is at least hin_max(0)
      !-----------------------------------------------------------------

            if (hice1(i,j) <= hin_max(0) .and. hin_max(0) > c0i ) then
               aicen(i,j,1) = vicen(i,j,1) / hin_max(0)
               hice1(i,j) = hin_max(0)
            endif

            if (aicen(i,j,1) > c0i) then
               dhi = hice1(i,j) - hice1_old(i,j)
               if (dhi < c0i) then  
                  hi0  = vicen(i,j,1) / aicen(i,j,1)
                  dai0 = vicen(i,j,1) / (hi0-p5*dhi) &
                       - aicen(i,j,1)
                  aicen(i,j,1) = aicen(i,j,1) + dai0 
               endif
            endif

         endif                  ! tmask
      enddo                     ! i
      enddo                     ! j

      end subroutine reduce_area

!=======================================================================
!BOP
!
! !IROUTINE: shift_ice - shift ice across category boundaries 
!
! !INTERFACE:
!
      subroutine shift_ice (donor, daice, dvice, hicen)
!
! !DESCRIPTION:
!
! Shift ice across category boundaries, conserving area, volume, and
! energy.
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
!          Elizabeth C. Hunke, LANL
!
! !USES:
!
      use ice_flux
      use ice_work, only: worka
!
! !INPUT/OUTPUT PARAMETERS:
! 
      ! 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), &
         intent(in) :: &
         donor             ! donor category index

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

      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              & ! thickness category index
      ,  n2              & ! receiver category
      ,  n1              & ! donor category
      ,  k                 ! ice layer index

      real (kind=dbl_kind), dimension(ilo:ihi,jlo:jhi,ncat) :: &
         atsfn

      real (kind=dbl_kind) :: &
         dvsnow          & ! snow volume transferred
      ,  desnow          & ! snow energy transferred
      ,  deice           & ! ice energy transferred
      ,  datsf             ! aicen*Tsfcn transferred

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

      integer (kind=int_kind) :: &
        icells        & ! number of cells with ice to transfer
      , ij              ! combined i/j horizontal index

      logical (kind=log_kind) :: &
        daice_negative      & ! true if daice < -puny
      , dvice_negative      & ! true if dvice < -puny
      , daice_greater_aicen & ! true if daice > aicen
      , dvice_greater_vicen  ! true if dvice > vicen

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


      !-----------------------------------------------------------------
      ! Define a variable equal to aicen*Tsfcn
      !-----------------------------------------------------------------
      do ni = 1, ncat
         do j = jlo,jhi
         do i = ilo,ihi
            atsfn(i,j,ni) = aicen(i,j,ni)*tsfcn(i,j,ni)
         enddo                  ! i
         enddo                  ! j
      enddo                     ! n

      !-----------------------------------------------------------------
      ! Check for daice or dvice out of range, allowing for roundoff error
      !-----------------------------------------------------------------

      do ni = 1, ncat-1

         daice_negative = .false.
         dvice_negative = .false.
         daice_greater_aicen = .false.
         dvice_greater_vicen = .false.

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

            if (donor(i,j,ni) > 0) then 
               n1 = donor(i,j,ni)

               if (daice(i,j,ni) < c0i) then
                  if (daice(i,j,ni) > -puny*aicen(i,j,n1)) then   
                     daice(i,j,ni) = c0i ! shift no ice
                     dvice(i,j,ni) = c0i
                  else
                     daice_negative = .true.
                  endif
               endif
         
               if (dvice(i,j,ni) < c0i) then
                  if (dvice(i,j,ni) > -puny*vicen(i,j,n1)) then   
                     daice(i,j,ni) = c0i ! shift no ice
                     dvice(i,j,ni) = c0i
                  else
                     dvice_negative = .true.
                  endif
               endif

               if (daice(i,j,ni) > aicen(i,j,n1)*(c1i-puny)) then
                  if (daice(i,j,ni) < aicen(i,j,n1)*(c1i+puny)) then
                     daice(i,j,ni) = aicen(i,j,n1)
                     dvice(i,j,ni) = vicen(i,j,n1)
                  else
                     daice_greater_aicen = .true.
                  endif
               endif    

               if (dvice(i,j,ni) > vicen(i,j,n1)*(c1i-puny)) then
                  if (dvice(i,j,ni) < vicen(i,j,n1)*(c1i+puny)) then
                     daice(i,j,ni) = aicen(i,j,n1)
                     dvice(i,j,ni) = vicen(i,j,n1)
                  else
                     dvice_greater_vicen = .true.
                  endif
               endif
               
            endif               ! donor > 0 
         enddo                  ! i
         enddo                  ! j

      !-----------------------------------------------------------------
      ! error messages
      !-----------------------------------------------------------------

         if (daice_negative) then
            do j = jlo,jhi
            do i = ilo,ihi
               if (donor(i,j,ni) > 0 .and.                            &
                   daice(i,j,ni) <= -puny*aicen(i,j,n1)) then          
                  write(nu_diag,*) my_task,':',i,j,                   &
                       'ITD Neg daice =',daice(i,j,ni),' boundary',ni   
!                  call abort_ice ('ice: ITD Neg daice')                
               endif                                                   
            enddo                                                      
            enddo                                                      
         endif                                                         
                                                                       
         if (dvice_negative) then                                      
            do j = jlo,jhi                                             
            do i = ilo,ihi                                             
               if (donor(i,j,ni) > 0 .and.                            &
                   dvice(i,j,ni) <= -puny*vicen(i,j,n1)) then          
                  write(nu_diag,*) my_task,':',i,j,                   &
                       'ITD Neg dvice =',dvice(i,j,ni),' boundary',ni    
!                  call abort_ice ('ice: ITD Neg dvice')                
               endif                                                   
            enddo                                                      
            enddo                                                      
         endif                                                         
                                                                       
         if (daice_greater_aicen) then                                 
            do j = jlo,jhi                                             
            do i = ilo,ihi                                             
               if (donor(i,j,ni) > 0) then                              
                  n1 = donor(i,j,ni)                                    
                  if (daice(i,j,ni) >= aicen(i,j,n1)*(c1i+puny)) then    
                     write(nu_diag,*) my_task,':',i,j,                &
                          'ITD daice > aicen, cat',n1                  
                     write(nu_diag,*) my_task,':',i,j,                &
                          'daice =', daice(i,j,ni),                    &
                          'aicen =', aicen(i,j,n1)
!                     call abort_ice ('ice: ITD daice > aicen')
                  endif
               endif
            enddo
            enddo
         endif

         if (dvice_greater_vicen) then
            do j = jlo,jhi
            do i = ilo,ihi
               if (donor(i,j,ni) > 0) then
                  n1 = donor(i,j,ni)
                  if (dvice(i,j,ni) >= vicen(i,j,n1)*(c1i+puny)) then
                     write(nu_diag,*) my_task,':',i,j,          &
                          'ITD dvice > vicen, cat',n1            
                     write(nu_diag,*) my_task,':',i,j,          &
                          'dvice =', dvice(i,j,ni),              &
                          'vicen =', vicen(i,j,n1)
!                     call abort_ice ('ice: ITD dvice > vicen')
                  endif
               endif
            enddo
            enddo
         endif

      !-----------------------------------------------------------------
      ! transfer volume and energy between categories
      !-----------------------------------------------------------------

         icells = 0
         do j = jlo, jhi
         do i = ilo, ihi
           if (daice(i,j,ni) > c0i) then ! daice(n) can be < puny
             icells = icells + 1
             indxi(icells) = i
             indxj(icells) = j
           endif   ! tmask
         enddo
         enddo

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

            n1 = donor(i,j,ni)
           !worka(i,j) = dvice(i,j,ni) / vicen(i,j,n1)
!! ggao change
           worka(i,j) = dvice(i,j,ni) / (vicen(i,j,n1)+epsilon(eps))
!! ggao change end 0616-2008

            if (n1  ==  ni) then
               n2 = n1+1
            else                ! n1 = n+1
               n2 = ni
            endif
            
            aicen(i,j,n1) = aicen(i,j,n1) - daice(i,j,ni)
            aicen(i,j,n2) = aicen(i,j,n2) + daice(i,j,ni)
            vicen(i,j,n1) = vicen(i,j,n1) - dvice(i,j,ni)
            vicen(i,j,n2) = vicen(i,j,n2) + dvice(i,j,ni)
            
            dvsnow = vsnon(i,j,n1) * worka(i,j)
            vsnon(i,j,n1) = vsnon(i,j,n1) - dvsnow
            vsnon(i,j,n2) = vsnon(i,j,n2) + dvsnow
            
            datsf = daice(i,j,ni)*tsfcn(i,j,n1)
            atsfn(i,j,n1) = atsfn(i,j,n1) - datsf
            atsfn(i,j,n2) = atsfn(i,j,n2) + datsf 

            desnow = esnon(i,j,n1) * worka(i,j)
            esnon(i,j,n1) = esnon(i,j,n1) - desnow
            esnon(i,j,n2) = esnon(i,j,n2) + desnow

         enddo                  ! ij

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

               n1 = donor(i,j,ni)
               if (n1  ==  ni) then
                  n2 = n1+1
               else             ! n1 = n+1
                  n2 = ni
               endif

               deice = eicen(i,j,ilyr1(n1)+k-1) * worka(i,j)
               eicen(i,j,ilyr1(n1)+k-1) =           &
                    eicen(i,j,ilyr1(n1)+k-1) - deice
               eicen(i,j,ilyr1(n2)+k-1) =           &
                    eicen(i,j,ilyr1(n2)+k-1) + deice
            enddo               ! ij
         enddo                  ! k

      enddo                     ! boundaries, 1 to ncat-1

      !-----------------------------------------------------------------
      ! Update ice thickness and temperature
      !-----------------------------------------------------------------

      do ni = 1, ncat
         do j = jlo,jhi
         do i = ilo,ihi
            if (aicen(i,j,ni) > puny) then
!               hicen(i,j,ni) = vicen(i,j,ni) / aicen(i,j,ni)
!               Tsfcn(i,j,ni) = aTsfn(i,j,ni) / aicen(i,j,ni)
!!  ggao change
               hicen(i,j,ni) = vicen(i,j,ni) /( aicen(i,j,ni)+epsilon(eps))
               tsfcn(i,j,ni) = atsfn(i,j,ni) /( aicen(i,j,ni)+epsilon(eps))
!! ggao change end 06162008

            else
               hicen(i,j,ni) = c0i
               tsfcn(i,j,ni) = tf(i,j)
            endif
         enddo                  ! i
         enddo                  ! j
      enddo                     ! n

      end subroutine shift_ice

!=======================================================================
!BOP
!
! !IROUTINE: column_sum - sum field over all ice categories
!
! !INTERFACE:
!
      subroutine column_sum (nsum, xin, xout)
!
! !DESCRIPTION:
!
! For each grid cell, sum field over all ice categories.
!
! !REVISION HISTORY:
!
! author: William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
      integer (kind=int_kind), intent(in) :: &
           nsum                              ! number of categories/layers

      real (kind=dbl_kind), intent(in) :: &
           xin(imt_local,jmt_local,nsum)   ! input field


      real (kind=dbl_kind), intent(out) :: &
           xout(imt_local,jmt_local)        ! output field
!
!EOP
!
      integer (kind=int_kind) :: &
           i, j               &  ! horizontal indices
      ,    ni                    ! category/layer index

      do j = 1, jmt_local
      do i = 1, imt_local
         xout(i,j) = c0i
      enddo
      enddo

      do ni = 1, nsum
         do j = 1, jmt_local
         do i = 1, imt_local
            xout(i,j) = xout(i,j) + xin(i,j,ni)
         enddo                  ! i
         enddo                  ! j
      enddo                     ! n

      end subroutine column_sum

!=======================================================================
!BOP
!
! !IROUTINE: column_conservation_check
!
! !INTERFACE:
!
      subroutine column_conservation_check (x1, x2, max_err, fieldid)
!
! !DESCRIPTION:
!
! For each physical grid cell, check that initial and final values
! of a conserved field are equal to within a small value.
!
! !REVISION HISTORY:
!
! author: William H. Lipscomb, LANL
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
      real (kind=dbl_kind), intent(in) :: &
           x1(imt_local,jmt_local)    ! initial field

      real (kind=dbl_kind), intent(in) :: &
           x2(imt_local,jmt_local)    ! final field

      real (kind=dbl_kind), intent(in) :: & 
          max_err                     ! max allowed error

      character (len=char_len), intent(in) :: &
           fieldid                     ! field identifier
!
!EOP
!
      integer (kind=int_kind) :: &
           i, j                        ! horizontal indices      

      logical (kind=log_kind) :: &
         conserv_err          ! = .true. if conservation check failed

      conserv_err = .false.

      do j = jlo, jhi
      do i = ilo, ihi
         if (abs(x2(i,j) - x1(i,j)) > max_err) then
            conserv_err = .true.
         endif
      enddo
      enddo

      if ( conserv_err ) then
        do j = jlo, jhi
        do i = ilo, ihi
         if (abs(x2(i,j) - x1(i,j)) > max_err) then
            write (nu_diag,*) ' '
            write (nu_diag,*) 'Conservation error: ', fieldid
            write (nu_diag,*)  my_task, ':', i, j
            write (nu_diag,*) 'Initial value =', x1(i,j)
            write (nu_diag,*) 'Final value =',   x2(i,j)
            write (nu_diag,*) 'Difference =', x2(i,j) - x1(i,j)
!            call abort_ice ('ice: Conservation error')
         endif
        enddo
        enddo
      endif

      end subroutine column_conservation_check

!=======================================================================
!BOP
!
! !IROUTINE: zap_small_areas - eliminate very small ice areas
!
! !INTERFACE:
!
      subroutine zap_small_areas
!
! !DESCRIPTION:
!
! For each ice category in each grid cell, remove ice if the fractional
! area is less than puny.
!
! !REVISION HISTORY:
!
! author: William H. Lipscomb, LANL
! Nov 2003:  Modified by Julie Schramm to conserve volume and energy
! Sept 2004: Modified by William Lipscomb; replaced normalize_state with
!            additions to local freshwater, salt, and heat fluxes
!            
!
! !USES:
!
      use ice_flux
!      use ice_calendar, only: dt
      use ice_calendar, only: dtice

!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
      integer (kind=int_kind) :: &
           i,j      & ! horizontal indices
      ,    ni       &  ! ice category index
      ,    k        & ! ice layer index
      ,    icells   & ! number of cells with ice to zap
      ,    ij        ! combined i/j horizontal index

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

      real (kind=dbl_kind) :: &
           xtmp      ! temporary variable

      do ni = 1, ncat

      !-----------------------------------------------------------------
      ! Count categories to be zapped.
      ! Abort model in case of negative area.
      !-----------------------------------------------------------------

         icells = 0
         do j = jlo, jhi
         do i = ilo, ihi
            if (aicen(i,j,ni) < -puny) then
               write (nu_diag,*) 'Negative ice area: i, j, n:', i, j, ni
               write (nu_diag,*) 'aicen =', aicen(i,j,ni)
!               call abort_ice ('zap: negative ice area')
            elseif ((aicen(i,j,ni) >= -puny .and. aicen(i,j,ni) < c0i) &
                                         .or.                         &
                    (aicen(i,j,ni) > c0i .and. aicen(i,j,ni) <= puny)) then
               icells = icells + 1
               indxi(icells) = i
               indxj(icells) = j
            endif
         enddo
         enddo

      !-----------------------------------------------------------------
      ! Zap ice energy and use ocean heat to melt ice
      !-----------------------------------------------------------------

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

!               xtmp = eicen(i,j,ilyr1(ni)+k-1) / dt ! < 0
               xtmp = eicen(i,j,ilyr1(ni)+k-1) / dtice ! < 0
               fhnet(i,j)      = fhnet(i,j)      + xtmp
               fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp
               eicen(i,j,ilyr1(ni)+k-1) = c0i

            enddo               ! ij
         enddo                  ! k

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

      !-----------------------------------------------------------------
      ! Zap snow energy and use ocean heat to melt snow
      !-----------------------------------------------------------------

!            xtmp = esnon(i,j,ni) / dt ! < 0
            xtmp = esnon(i,j,ni) / dtice ! < 0
            fhnet(i,j)      = fhnet(i,j)      + xtmp
            fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp
            esnon(i,j,ni) = c0i

      !-----------------------------------------------------------------
      ! zap ice and snow volume, add water and salt to ocean
      !-----------------------------------------------------------------

!            xtmp = (rhoi*vicen(i,j,ni) + rhos*vsnon(i,j,ni)) / dt
            xtmp = (rhoi*vicen(i,j,ni) + rhos*vsnon(i,j,ni)) / dtice
            fresh(i,j)      = fresh(i,j)      + xtmp
            fresh_hist(i,j) = fresh_hist(i,j) + xtmp

!            xtmp = rhoi*vicen(i,j,n)*ice_ref_salinity*p001/ dt
            xtmp = rhoi*vicen(i,j,ni)*ice_ref_salinity*p001/ dtice
            fsalt(i,j)      = fsalt(i,j)      + xtmp
            fsalt_hist(i,j) = fsalt_hist(i,j) + xtmp

            aice0(i,j) = aice0(i,j) + aicen(i,j,ni)
            aicen(i,j,ni) = c0i
            vicen(i,j,ni) = c0i
            vsnon(i,j,ni) = c0i
            tsfcn(i,j,ni) = tf(i,j)

         enddo                  ! ij
      enddo                     ! n

      end subroutine zap_small_areas

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

      end module ice_itd

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