ice_dh.f

CCSM Research Tools: Community Atmosphere Model (CAM)

c=======================================================================
!---! Energy-conserving sea ice model
!---! Routines to grow/melt ice and adjust temperature profile
!---!
!---! author C. M. Bitz
!---!
!---! 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. 
!---!     
!---! The author grants permission to the public to copy and use this
!---! software without charge, provided that this Notice and any statement
!---! of authorship are reproduced on all copies and any publications that
!---! result from the use of this software must (1) refer to the publications 
!---! listed above and (2) acknowledge the origin and author of the model.
!---! This software is without warranty, expressed or implied, and the
!---! author assumes no liability or responsibility for its use. 
c=======================================================================

      module ice_dh

      use ice_kinds_mod
      use ice_constants
!      use ice_itd

      implicit none

      integer (kind=int_kind) :: errflag

! ADD FLAG FOR ATMOSPHERE MODEL VERSION OF THIS ROUTINE
! Skips the last half of this routine
! Maybe someday this model will be used with a mixed layer ocean
! and then this could be set to false
      logical (kind=log_kind), parameter :: fixice = .true.
!
! Add flag if prognostic snow on ice, and flag to reset csim iceprops.
!
      logical (kind=log_kind) :: prognostic_icesnow,reset_csim_iceprops

c=======================================================================

      contains

c=======================================================================

      logical function icemodel_is( name )
!
! Input arguments
!
      character(len=*) :: name
      
      if (name == 'csim' .or. name == 'CSIM' ) then
         icemodel_is = .true.
      else
         icemodel_is = .false.
      end if
      
      return

      end function icemodel_is

c=======================================================================

      subroutine dh(  dtsub,    sal1d,    tiz
     $              ,  tbot,       hi,       hs,   fbot
     $              ,  fnet,    condb,      flh
     $              ,  dhib,     dhit,      dhs,   subi
     $              ,  subs,     dhif,     dhsf,     qi
     $              ,  focn, i,j )

!---!-------------------------------------------------------------------
!---! Computes the thickness changes at the top and bottom
!---! and adjusts layer energy of melt
!---! does not allow h<0  
!---! Focn= actual flux of heat from the ocean layer under sea ice 
!---! (equal to fbot unless all the ice melts away)
!---! compensates for rare case of melting entire slab through 
!---!-------------------------------------------------------------------

!      use ice_state
!      use ice_diagnostics

      real (kind=dbl_kind), intent(in) :: 
     &   dtsub                ! timestep
     &,  sal1d   (plevmx+1)     ! ice salinity                           (ppt)
     &,  tiz   (0:plevmx)       ! snow/ice internal temp                   (C)
     &,  Tbot                 ! ice bottom in                            (C)
     &,  hi                   ! initial ice thickness                    (m)
     &,  hs                   ! initial snow thickness                   (m)
     &,  fbot                 ! flx from ocean, potent.             (W/m**2)
     &,  fnet                 ! net flx at top srf incl. cond. flx  (W/m**2)
     &,  condb                ! cond. flx at bot.                   (W/m**2)
     &,  flh                  ! latent heat flx                     (w/m**2)

      integer (kind=int_kind), intent(in) :: 
     &   i,j ! grid location for debugging

      ! thickness changes from grow/melt (default) or sublimate/flooding      
      real (kind=dbl_kind), intent(out) :: 
     &   dhib                 ! ice bot, dhib<0 if melt                  (m)
     &,  dhit                 ! ice top, dhit<=0                         (m)
     &,  dhs                  ! snow top, dhit<=0                        (m)
     &,  subi                 ! ice top, subi<0 if sublimating           (m)
     &,  subs                 ! snow, subs<0 if sublimating              (m)
     &,  dhif                 ! ice top from flooding, dhif>0            (m)
     &,  dhsf                 ! snow from flooding, dhsf<0               (m)

     &,  qi(plevmx)             ! energy of melt of ice per unit vol. (J/m**3)
     &,  focn                 ! actual flx of heat used from ocn    (w/m**2)

      real (kind=dbl_kind) :: 
     &   delti(plevmx)          ! evolving ice layer thickness 
     &,  delts                ! evolving snow thickness
     &,  sumr                 ! dummy for adjusting delti
     &,  rmvi                 ! another dummy for adjusting delti
     &,  dtop                 ! dhit+subi for adjust
      
     &,  qigrow               ! energy of melt of ice that grows 
     &,  ebot                 ! heat available to grow/melt at bottom
     &,  etop                 ! heat available to melt at top

     &,  qs                   ! energy of melt of snow per unit vol. (J/m**3)
     &,  enet                 ! sum of energy of melt of ice and snow (J/m**2)
     &,  evnet                ! sum of energy of melt and vapor. of ice and snow (J/m**2)

     &,  qiflood              ! energy of melt of flooded ice      (W/m**2)

      real (kind=dbl_kind) ::     hi_tw          ! ice thickness   (m)

      integer (kind=int_kind) :: layer

      logical (kind=log_kind) :: verbos
c      verbos = .true.
      verbos = .false.

      errflag = 0

      dhib  = c0
      dhit  = c0
      dhs   = c0
      dhif  = c0
      dhsf  = c0
      subi  = c0
      subs  = c0
      qiflood = c0

      ! thickness of snow and each ice layer 
      delti(1) = hi/ni
      do layer = 2,ni
        delti(layer) = delti(1)
      enddo
      delts = hs

      ! energy of melt per unit vol for snow and ice for each layer and sum
      qs = -rLfs
      enet = c0
      do layer = 1,ni
        qi(layer) = energ(tiz(layer),sal1d(layer))
        enet = enet + qi(layer)
      enddo
      enet = enet*delti(1) + qs*hs

      !-----------------------------------------------------------------
      ! sublimate/condense 
      !-----------------------------------------------------------------
      etop  = -flh*dtsub           !  etop>0
      evnet = -rLvs*hs-rLvi*hi+enet+etop
      if ( evnet .ge. 0. ) then    !  should never happen
        subi = -hi
        subs = -hs
        focn = condb + evnet/dtsub   
        write(6,*)  flh,dtsub,-rLvs*hs,-rLvi*hi,enet,etop
        write(6,*) 'sublimate away all sea ice'
        write(6,*) 'something is probably seriously wrong'
!        call print_state('ice state at dh stop',i,j)
        stop
      endif
!      call t_startf ('srfsub')
      call srfsub(    qi,   qs, delti, delts,
     $              subi, subs,  etop,  enet )
!      call t_stopf ('srfsub')

      ! adjust the layer thickness to reflect subl/cond changes
      delts = hs + subs
      rmvi = subi
      do layer = 1,ni
        sumr = max( -delti(layer), rmvi )
        rmvi = rmvi-sumr
        delti(layer) = delti(layer) + sumr
      enddo

      !-----------------------------------------------------------------
      ! melt at top srf melt      
      ! may melt when Tsf < melting, but alway a neglible amount
      !-----------------------------------------------------------------
      if ( fnet .gt. 0. ) then
        etop = fnet * dtsub
        enet = enet + etop
        if ( enet .ge. 0. ) then          !    remotely possible
          dhit = -( hi + subi )
          dhs  = -( hs + subs )
          focn = condb + enet/dtsub
          return
        endif
!        call t_startf ('srfmelt')
        call srfmelt(   qi,   qs, delti, delts,
     $                dhit,  dhs,  etop  )
!        call t_stopf ('srfmelt')
        ! adjust the layer thickness to reflect melt changes
        delts = delts + dhs
        rmvi = dhit
        do layer = 1,ni
          sumr = max( -delti(layer), rmvi )
          rmvi = rmvi - sumr
          delti(layer) = delti(layer) + sumr
        enddo
      endif

      if (fixice) return

      dtop = dhit + subi

      !-----------------------------------------------------------------
      ! melt/grow at bot srf
      !-----------------------------------------------------------------
      focn = fbot
      ebot = dtsub * ( condb - fbot )
      if (ebot .le. 0. ) then            !    grow at bottom
        qigrow = energ( tbot, salnew ) 
        dhib = ebot/qigrow
      else                               !    melt at bottom
        qigrow = c0                      !    on purpose
        if ( (enet+ebot) .ge. 0. ) then  !    remotely possible
          dhib = -( hi + dtop )
          dhs  = -( hs + subs )
          focn = condb + enet / dtsub
          return
        endif
!        call t_startf ('botmelt')
        call botmelt(   qi,   qs, delti, delts, 
     $                dhib,  dhs,  ebot  )
!        call t_stopf ('botmelt')
      endif

      !-----------------------------------------------------------------
      ! stop if error occurred in srfsub, srfmelt or botmelt
      !-----------------------------------------------------------------
      if (errflag.ne.0) then
!         call print_state('state at dh error   ',i,j)
         stop
      endif

      !-----------------------------------------------------------------
      ! check to see if there is any ice left after top/bottom 
      ! melt/growth
      !-----------------------------------------------------------------
      hi_tw = hi+dhib+dtop

      if ( hi_tw .le. 0. ) then
         dhib = -(hi+dtop)
         ! convert any snow that might be left into sea ice
         dhsf = min(-(hs+dhs),c0)
         dhif = -rhos / rhoi * dhsf
         do layer = 1,ni
           qi(layer) = qs * rhoi / rhos
         enddo
      else

      !-----------------------------------------------------------------
      ! flooding
      !-----------------------------------------------------------------
!        call t_startf ('freeboard')
        call freeboard(hs,hi,dhs,qs,dhsf,dhif,qiflood)
!        call t_stopf ('freeboard')
        if (verbos) write(6,*) 'fld',dhib,dtop,dhif,dhsf

      !-----------------------------------------------------------------
      ! adjust layers
      !-----------------------------------------------------------------
!        call t_startf ('adjust')
        call adjust(hi,dhib,dtop,dhif,dhsf,qiflood,qigrow,qi) 
!        call t_stopf ('adjust')

      endif

      end subroutine dh

c=======================================================================

      subroutine freeboard(hs,hi,dhs,qs,dhsf,dhif,qiflood)

!---!-------------------------------------------------------------------
!---! freeboard adjustment due to flooding ... snow-ice formation
!---!-------------------------------------------------------------------

      real (kind=dbl_kind), intent(in) ::
     &    hi  ! initial ice thickness                    (m)
     &,   hs  ! initial snow thickness                   (m)
     &,   dhs ! snow top, dhit<=0                        (m)
     &,   qs  ! energy of melt of snow per unit vol. (J/m**3)

      real (kind=dbl_kind), intent(out) ::
     &    dhif ! ice top from flooding, dhif>0            (m)
     &,   dhsf ! snow from flooding, dhsf<0               (m)
     &,   qiflood ! energy of melt of flooded ice    (W/m**2)

      real (kind=dbl_kind) :: zintfc  ! height of snow/ice interf. wrt ocn (m)

      zintfc  =  hi - (rhos*hs + rhoi*hi)/rhow
      if (( zintfc .lt. 0. ).and.(hs+dhs.gt.0.)) then
        dhsf    =  rhoi / rhos * zintfc
        dhsf    =  max( dhsf, -(hs+dhs) )
        dhif    =  -rhos / rhoi * dhsf
        qiflood =  qs * rhoi / rhos
      endif