ice_dh.f

CCSM Research Tools: Community Atmosphere Model (CAM)

      end subroutine freeboard

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

      subroutine srfsub(  qi,   qs, delti, delts,
     $                   subi, subs,  etop,  enet )

!---!-------------------------------------------------------------------
!---! compute the sea ice and snow thickness changes from  
!---! sublimation/condensation
!---!-------------------------------------------------------------------

      real (kind=dbl_kind), intent(in) :: 
     &   qi (1:plevmx), qs      ! energy of melt of ice/snow per vol  (J/m**3)
     &,  delti(plevmx), delts   ! thickness of ice/snow layer              (m)

      real (kind=dbl_kind), intent(out) :: 
     &   subi, subs  ! subl/cond. amount for ice/snow           (m)

      real (kind=dbl_kind), intent(inout) :: 
     &   etop                 ! energy avail to sub/cond ice/snow   (J/m**2)
     &,  enet                 ! energy needed to melt all ice/snow  (J/m**2)

      real (kind=dbl_kind) :: 
     &   ue                   ! energy of melt + vapor for ice/snow (J/m**3)
     &,  subil                ! subl from this layer

      integer (kind=int_kind) :: layer

      subs = c0
      subi = c0
      if ( delts .gt. 0. ) then
        ! convert etop into snow thickness
        ! positive if cond/negative if subl
        ue = qs - rLvs               
        subs = etop/ue     
        if ( (delts + subs) .ge. 0. ) then
          ! either all condensation becomes snow
          ! or sublimate some snow but no ice
          etop = c0
          enet = enet + subs * qs
          return
        else
          ! sublimate all of the snow and some ice too
          subs = -delts
          etop = etop + delts * ue
          enet = enet + subs  * qs
        endif
      endif
      
      do layer = 1,ni
        ! convert etop into ice thickness
        ! positive if cond/negative if subl
        ue = qi(layer) - rLvi
        subil = etop/ue
        if ( (delti(layer)+subil) .ge. 0. ) then
          ! either all condensation becomes ice 
          ! or sublimate some ice from this layer
          subi = subi + subil
          etop = c0
          enet = enet + subil * qi(layer)
          return
        else
          ! sublimate all of the layer 
          subi = subi - delti(layer)  
          etop = etop + delti(layer) * ue
          enet = enet - delti(layer) * qi(layer)
        endif
      enddo
      
      write(6,*)  'ERROR in srfsub',etop
      errflag = 1

      end subroutine srfsub
      
c=======================================================================

      subroutine srfmelt(   qi,   qs, delti, delts,
     $                    dhit,  dhs,  etop  ) 

!---!-------------------------------------------------------------------
!---! melt ice/snow from the top srf
!---!-------------------------------------------------------------------

      real (kind=dbl_kind), intent(in) :: 
     &   qi (1:plevmx), qs      ! energy of melt of ice/snow per vol (J/m**3)
     &,  delti(plevmx), delts   ! thickness of ice/snow layer             (m)

      real (kind=dbl_kind), intent(out) :: dhit, dhs ! ice/snow thickness change               (m)

      real (kind=dbl_kind), intent(inout) :: etop    ! energy avail to melt ice and snow  (J/m**2)

      real (kind=dbl_kind) :: dhitl  ! melt from this layer          (m)
      integer (kind=int_kind) :: layer

      dhit = c0
      dhs  = c0

      if ( delts .gt. 0. ) then
        ! convert etop into snow thickness
        dhs  =  etop/qs
        if ( (delts + dhs) .ge. 0. ) then
          ! melt only some of the snow
          etop = c0
          return
        else
          ! melt all of the snow and some ice too
          dhs  = -delts
          etop = etop+qs*delts
        endif
      endif

      do layer = 1,ni
        ! convert etop into ice thickness
        dhitl = etop/qi(layer)
        if ( (delti(layer)+dhitl) .ge. 0. ) then
          ! melt some ice from this layer
          dhit = dhit + dhitl
          etop = c0
          return
        else
          ! melt all of the ice in this layer
          dhit = dhit - delti(layer)   
          etop = etop + delti(layer) * qi(layer)
        endif
      enddo
      
      write(6,*)  'ERROR in srfmelt',etop
      errflag = 1

      end subroutine srfmelt
      
c=======================================================================

      subroutine botmelt(   qi,   qs, delti, delts,
     $                    dhib,  dhs,  ebot  )

!---!-------------------------------------------------------------------
!---! melt from bottom
!---!-------------------------------------------------------------------

      real (kind=dbl_kind), intent(in) :: 
     &   qi (1:plevmx), qs      ! energy of melt of ice/snow per vol (J/m**3)
     &,  delti(plevmx), delts   ! thickness of ice/snow layer             (m)

      real (kind=dbl_kind), intent(out) :: dhib  ! ice thickness change (m)

      real (kind=dbl_kind), intent(inout) :: 
     &   ebot                 ! energy avail to melt ice and snow  (J/m**2)
     &,  dhs                  ! snow thickness change                   (m)

      real (kind=dbl_kind) :: 
     &   dhibl                ! melt from this ice layer                (m)
     &,  dhsl                 ! melt from this snow layer               (m)
      integer (kind=int_kind) :: layer

      dhib = c0

      do layer = ni,1,-1
         dhibl = ebot/qi(layer)
         if ( (delti(layer)+dhibl) .ge. 0. ) then
            dhib = dhib + dhibl
            ebot = c0
            return
         else
            dhib = dhib - delti(layer)   
            ebot = ebot + delti(layer) * qi(layer)
         endif
      enddo

      ! finally melt snow if necessary
      dhsl = ebot/qs
      if ( (delts + dhsl) .ge. 0. ) then
         dhs = dhs + dhsl
         ebot = c0
         return
      endif

      write(6,*)  'ERROR in botmelt',ebot
      errflag = 1

      end subroutine botmelt

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

      subroutine adjust(hi0,dhib,dhit,dhif,dhsf,qiflood,qigrow,qi_tw)

!---!-------------------------------------------------------------------
!---! Adjusts temperature profile to account for changing
!---! the layer spacing due to growth/melt (incl. subl/cond, flooding)
!---! At start the energy of melting was computed
!---! after updating tiz from the heat equation
!---! hi is the thickness prior to changes from dhib and dhit
!---! hi_tw is the thickness after making these changes
!---! dhib<0 if there is melt at the bottom
!---! dhit<0 if there is melt at the top
!---! generally _tw is a suffix to label the adjusted variables
!---!-------------------------------------------------------------------

      real (kind=dbl_kind), intent(in) :: 
     &   hi0                  ! initial ice thickness                    (m)
     &,  dhib                 ! ice bot, dhib<0 if melt                  (m)
     &,  dhit                 ! ice top, dhit<=0                         (m)
     &,  dhif                 ! ice top from flooding, dhif>0            (m)
     &,  dhsf                 ! snow from flooding, dhsf<0               (m)

     &,  qiflood              ! qi for flooded ice                  (J/m**3)
     &,  qigrow               ! qi for ice growing on bot           (J/m**3)

      real (kind=dbl_kind), intent(inout) :: 
     &   qi_tw(plevmx)          ! energy of melt of ice per unit vol. (J/m**3)


      ! the following pairs are before & after adjusting
      real (kind=dbl_kind) ::     hi, hi_tw          ! ice thickness   (m)
     &,      z(plevmx+2)        ! vertical layer position        (m)
     &,      z_tw(plevmx+1)     ! vertical layer position        (m)
     &,      D, D_tw          ! layer thickness                (m)
      integer (kind=int_kind) ::  k, k_tw          ! index of layer

      real (kind=dbl_kind) :: qi(plevmx)  ! qi of initial layers  (J/m**3)
     &,  fract(plevmx,plevmx)     ! fract of layer k that makes up layer k_tw
      integer (kind=int_kind) :: layer             

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

      ! there is ice left so allow snow-ice conversion from flooding
      hi_tw = hi_tw+dhif

      if ( (abs(dhib) .gt. puny) .or. (abs(dhit) .gt. puny)
     $     .or. (dhif.gt.0)) then
         do layer = 1,ni
            qi(layer) = qi_tw(layer)
         enddo
      
         ! layer thickness      
         D = hi/ni 
         D_tw = hi_tw/ni 

         ! z is positive down and zero is relative to the top 
         ! of the ice from the old time step
         z(1) = -dhit          ! necessary
         z_tw(1) = -dhit-dhif
         do layer = 2,ni
           z(layer) = D*(layer-1)
           z_tw(layer) = z_tw(1)+D_tw*(layer-1) 
         enddo
         z(ni+1) = hi + min(dhib,c0)
         z_tw(ni+1) = z_tw(1)+hi_tw

         do k_tw = 1,ni
           qi_tw(k_tw) = c0
           do k = 1,ni
             fract(k,k_tw) = ( min( z_tw(k_tw+1), z(k+1)) -
     $                         max( z_tw(k_tw)  , z(k)  )   )
             if (fract(k,k_tw).gt.0.)
     $            qi_tw(k_tw) = qi_tw(k_tw)+fract(k,k_tw)*qi(k)
           enddo
         enddo

         if (dhif.gt.D_tw) then
           do k = 1,ni
             qi_tw(k) = qi_tw(k) + qiflood*max(c0,
     $                   min(D_tw,dhif-D_tw*(k-1)))
           enddo
         else
           qi_tw(1) = qi_tw(1) + dhif*qiflood
         endif

         if (dhib.gt.D_tw) then
           z(ni+2) = hi + dhib
           do k = 1,ni
             qi_tw(k) = qi_tw(k) +  qigrow*max(c0,
     $                   (min(z(ni+2),Z_tw(k+1))-max(hi,Z_tw(k))))
           enddo
         else
           qi_tw(ni) = qi_tw(ni) + max(dhib,c0)*qigrow
         endif

         do k_tw = 1,ni
           qi_tw(k_tw) = qi_tw(k_tw)/D_tw
         enddo

      endif

      end subroutine adjust

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

      real function energ(Tmp ,sal)

!---!-------------------------------------------------------------------
!---! compute the energy of melting per unit volume  (J/m**3)
!---! relative to melting (negative quantity)
!---!-------------------------------------------------------------------

      real (kind=dbl_kind), intent(in) ::
     &    Tmp          ! midpt temperature of ice layer       (C)
     &,   sal          ! midpt salinity of ice layer        (ppt)  

      energ = -rLfi - rcpi*(-depressT*sal-Tmp)-rLfidepressT*sal/Tmp

      end function energ

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

      end module ice_dh

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