stim_lebedev.F90 Source File

!> The Lebedev ice model
!>
!> authors: Adolf Stips and Karsten Bolding

   MODULE stim_lebedev

  !! https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JC012199

   use stim_variables, only: transmissivity, albedo_ice
   use stim_variables, only: rk, Hice, dHis, dHib, Tf, fdd
   use stim_variables, only: rho_ice, L_ice, ocean_ice_flux

   IMPLICIT NONE

   private

   public init_stim_lebedev, do_stim_lebedev, clean_stim_lebedev

   !! Constants provided by Adolf Stips
   real(rk)           :: lebedev_fac=1.33_rk
   real(rk)           :: damp_leb_swr = -1.6_rk
   real(rk)           :: damp_leb_wind = -1.6_rk
   real(rk)           :: damp_leb_shf = -1.6_rk
   real(rk)           :: freeze_fac = -0.0575_rk
   real(rk)           :: lebedev_albedo = 0.545

!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
   SUBROUTINE init_stim_lebedev(ice_cover)

   integer, intent(inout)  :: ice_cover

   integer                   :: rc
!-----------------------------------------------------------------------
   if (Hice .gt. 0._rk) then
      fdd = (100._rk*Hice/lebedev_fac)**(1._rk/0.58_rk)
      ice_cover = 2
   else
      fdd = 0._rk
      ice_cover = 0
   end if
   END SUBROUTINE init_stim_lebedev

!-----------------------------------------------------------------------

   SUBROUTINE do_stim_lebedev(ice_cover,dt,Tw,S,Ta,precip)
     !! Calculate ice thickness according to Lebedev 1938 

   real(rk), intent(in)    :: dt,Ta,S,precip

   integer, intent(inout)  :: ice_cover
   real(rk), intent(inout) :: Tw

   real(rk) :: x
!-----------------------------------------------------------------------
   Tf = freeze_fac*S
   if (ice_cover .eq. 2) then
      ! calculate the cumulative freezing days
      fdd = fdd+(Tf-Ta)*dt/86400._rk
   else ! ice_hi was 0._rk
      if (Ta .lt. Tf .and. Tw .lt. Tf) then
         fdd = fdd+(Tf-Ta)*dt/86400._rk
      end if 
   end if

   ! have at least 1 cm of ice,.. - the melting is not done correct - KB
   if (fdd .gt. 1._rk) then
      x = Hice
      Hice = 0.01_rk*lebedev_fac*fdd**0.58_rk
      dHis = Hice - x
#if 0
      if (dHis .lt. 0._rk) then
         dHib = (-dt*sensible_ice_water)/(rho_ice*L_ice)
      else
         dHib = 0._rk
      end if
#endif
      Tw = Tf
      ice_cover = 2
      albedo_ice = lebedev_albedo
      transmissivity = exp(Hice/damp_leb_swr)
!      transmissivity = exp(Hice*damp_leb_swr)
      if (Hice .lt. 0._rk) then
         albedo_ice = 0._rk
         transmissivity = 1._rk
         Hice = 0._rk
         ice_cover = 0
      end if
   else
      fdd = 0._rk
      Hice = 0._rk
      dHis = 0._rk
      ice_cover = 0
   end if

#if 0
   ! if we have ice change the heatfluxes !
   where (ice_hi .gt. 0._rk)
      ! damp the shortwave radiation
      swr = swr * exp(damp_leb_swr*ice_hi ) 
      albedo = albedo_ice
      ! damp the wind stress
      tausx = tausx * exp(damp_leb_wind*ice_hi ) 
      tausy = tausy * exp(damp_leb_wind*ice_hi )
      ! damp the heatflux ­ this is of course tricky as latent heatflux and
      ! longwave radiation behave different from sensible heatflux
      ! flux routine must be rewritten to treat them all separate
      ! heat = heat * exp(damp_leb_shf * ice_hi)
      ! estimate sensible heat flux between ice and ocean water shf = sens_ice_water * ( freezing_fac * sss ­ sst )
   end where
#endif
   END SUBROUTINE do_stim_lebedev

!-----------------------------------------------------------------------

   SUBROUTINE clean_stim_lebedev()

   IMPLICIT NONE
!-----------------------------------------------------------------------
   LEVEL2 'clean_stim_lebedev'
   END SUBROUTINE clean_stim_lebedev

!-----------------------------------------------------------------------

   END MODULE stim_lebedev

!-----------------------------------------------------------------------
! Copyright by the STIM-team under the GNU Public License - www.gnu.org
!-----------------------------------------------------------------------