scale_ocean_phy_ice_simple.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_ocean_phy_ice_simple
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: ocean_phy_ice_setup
  public :: ocean_phy_ice_fraction
  public :: ocean_phy_ice_adjustment
  public :: ocean_phy_ice_simple

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), public :: ocean_phy_ice_freezetemp      =  271.35_rp ! freezing temperature of sea ice [K]
  real(RP), public :: ocean_phy_ice_density         =  1000.0_rp ! density of sea ice              [kg/m3]

  real(RP), private :: ocean_phy_ice_mass_critical  =  1600.0_rp ! ice amount for fraction = 1     [kg/m2]
  real(RP), private :: ocean_phy_ice_mass_limit     = 50000.0_rp ! maximum ice amount              [kg/m2]
  real(RP), private :: ocean_phy_ice_fraction_limit =     1.0_rp ! maximum ice fraction            [1]
  real(RP), private :: ocean_phy_ice_dt_max         =   5.e-2_rp ! maximum delta ice temperature   [K/s]

  logical,  private :: ocean_phy_ice_nudging        = .false.    
  real(DP), private :: ocean_phy_ice_nudging_tausec
  logical,  private :: ocean_phy_ice_offline_mode   = .false.    

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine ocean_phy_ice_setup
    use scale_prc, only: &
       prc_abort
    use scale_const, only: &
       undef => const_undef
    use scale_calendar, only: &
       calendar_unit2sec
    use scale_file_external_input, only: &
       file_external_input_file_limit, &
       file_external_input_regist
    implicit none

    real(DP)               :: OCEAN_PHY_ICE_nudging_tau                                      = 0.0_dp  ! Relaxation time
    character(len=H_SHORT) :: OCEAN_PHY_ICE_nudging_tau_unit                                 = "SEC"
    character(len=H_LONG)  :: OCEAN_PHY_ICE_nudging_basename(file_external_input_file_limit) = ''
    logical                :: OCEAN_PHY_ICE_nudging_enable_periodic_year                     = .false.
    logical                :: OCEAN_PHY_ICE_nudging_enable_periodic_month                    = .false.
    logical                :: OCEAN_PHY_ICE_nudging_enable_periodic_day                      = .false.
    integer                :: OCEAN_PHY_ICE_nudging_step_fixed                               = 0
    real(RP)               :: OCEAN_PHY_ICE_nudging_offset                                   = 0.0_rp
    real(RP)               :: OCEAN_PHY_ICE_nudging_defval                                  != UNDEF
    logical                :: OCEAN_PHY_ICE_nudging_check_coordinates                        = .true.
    integer                :: OCEAN_PHY_ICE_nudging_step_limit                               = 0

    namelist / param_ocean_phy_ice / &
       ocean_phy_ice_density,                       &
       ocean_phy_ice_mass_critical,                 &
       ocean_phy_ice_mass_limit,                    &
       ocean_phy_ice_fraction_limit,                &
       ocean_phy_ice_dt_max!,                        &
!        OCEAN_PHY_ICE_nudging,                       &
!        OCEAN_PHY_ICE_nudging_tau,                   &
!        OCEAN_PHY_ICE_nudging_tau_unit,              &
!        OCEAN_PHY_ICE_nudging_basename,              &
!        OCEAN_PHY_ICE_nudging_enable_periodic_year,  &
!        OCEAN_PHY_ICE_nudging_enable_periodic_month, &
!        OCEAN_PHY_ICE_nudging_enable_periodic_day,   &
!        OCEAN_PHY_ICE_nudging_step_fixed,            &
!        OCEAN_PHY_ICE_nudging_offset,                &
!        OCEAN_PHY_ICE_nudging_defval,                &
!        OCEAN_PHY_ICE_nudging_check_coordinates,     &
!        OCEAN_PHY_ICE_nudging_step_limit

    integer :: ierr
    !---------------------------------------------------------------------------

    log_newline
    log_info("OCEAN_PHY_ICE_setup",*) 'Setup'

    ocean_phy_ice_nudging_defval = undef

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_ocean_phy_ice,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("OCEAN_PHY_ICE_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("OCEAN_PHY_ICE_setup",*) 'Not appropriate names in namelist PARAM_OCEAN_PHY_ICE. Check!'
       call prc_abort
    endif
    log_nml(param_ocean_phy_ice)

    log_newline
    log_info("OCEAN_PHY_ICE_setup",*) 'Ice amount for frac. = 1 [kg/m2] : ', ocean_phy_ice_mass_critical
    if ( ocean_phy_ice_nudging ) then
       call calendar_unit2sec( ocean_phy_ice_nudging_tausec, ocean_phy_ice_nudging_tau, ocean_phy_ice_nudging_tau_unit )

       log_info("OCEAN_PHY_ICE_setup",*) 'Use nudging for sea ice fraction : ON'
       log_info("OCEAN_PHY_ICE_setup",*) 'Relaxation time Tau [sec]        : ', ocean_phy_ice_nudging_tausec

       if ( ocean_phy_ice_nudging_tausec == 0.0_rp ) then
          ocean_phy_ice_offline_mode = .true.
          log_info("OCEAN_PHY_ICE_setup",*) 'Tau=0 means that sea ice is completely replaced by the external data.'
       endif

       if ( ocean_phy_ice_nudging_basename(1) == '' ) then
          log_error("OCEAN_PHY_ICE_setup",*) 'OCEAN_PHY_ICE_nudging_basename is necessary. STOP'
          call prc_abort
       endif
    else
       log_info("OCEAN_PHY_ICE_setup",*) 'Use nudging for sea ice fraction : OFF'
    endif

    if ( ocean_phy_ice_nudging ) then
       call file_external_input_regist( ocean_phy_ice_nudging_basename(:),           & ! [IN]
                                        'OCEAN_ICE_FRAC',                            & ! [IN]
                                        'XY',                                        & ! [IN]
                                        ocean_phy_ice_nudging_enable_periodic_year,  & ! [IN]
                                        ocean_phy_ice_nudging_enable_periodic_month, & ! [IN]
                                        ocean_phy_ice_nudging_enable_periodic_day,   & ! [IN]
                                        ocean_phy_ice_nudging_step_fixed,            & ! [IN]
                                        ocean_phy_ice_nudging_offset,                & ! [IN]
                                        ocean_phy_ice_nudging_defval,                & ! [IN]
                                        ocean_phy_ice_nudging_check_coordinates,     & ! [IN]
                                        ocean_phy_ice_nudging_step_limit             ) ! [IN]
    endif

    return
  end subroutine ocean_phy_ice_setup

  !-----------------------------------------------------------------------------
  subroutine ocean_phy_ice_fraction( &
       OIA, OIS, OIE, &
       OJA, OJS, OJE, &
       ICE_MASS,      &
       ICE_FRAC       )
    implicit none

    integer,  intent(in)  :: OIA, OIS, OIE
    integer,  intent(in)  :: OJA, OJS, OJE
    real(RP), intent(in)  :: ICE_MASS(oia,oja) ! sea ice amount        [kg/m2]
    real(RP), intent(out) :: ICE_FRAC(oia,oja) ! sea ice area fraction [1]

    integer :: i, j
    !---------------------------------------------------------------------------

    do j = ojs, oje
    do i = ois, oie
       ice_frac(i,j) = ice_mass(i,j) / ocean_phy_ice_mass_critical

       ice_frac(i,j) = min( sqrt( max( ice_frac(i,j), 0.0_rp ) ), ocean_phy_ice_fraction_limit )
    enddo
    enddo

    return
  end subroutine ocean_phy_ice_fraction

  !-----------------------------------------------------------------------------
  subroutine ocean_phy_ice_adjustment( &
       OIA, OIS, OIE, &
       OJA, OJS, OJE, &
       calc_flag,     &
       OCEAN_DEPTH,   &
       OCEAN_TEMP,    &
       ICE_TEMP,      &
       ICE_MASS       )
    use scale_const, only: &
       cl    => const_cl,    &
       ci    => const_ci,    &
       emelt => const_emelt, &
       dwatr => const_dwatr
    implicit none

    integer,  intent(in)    :: OIA, OIS, OIE
    integer,  intent(in)    :: OJA, OJS, OJE
    logical,  intent(in)    :: calc_flag (oia,oja) ! to decide calculate or not
    real(RP), intent(in)    :: OCEAN_DEPTH         ! depth of the first layer of the ocean
    real(RP), intent(inout) :: OCEAN_TEMP(oia,oja) ! ocean temperature   [K]
    real(RP), intent(inout) :: ICE_TEMP  (oia,oja) ! sea ice temperature [K]
    real(RP), intent(inout) :: ICE_MASS  (oia,oja) ! sea ice amount      [kg/m2]

    real(RP) :: ICE_MASS_diff
    real(RP) :: ICE_MASS_prev
    real(RP) :: factor
    real(RP) :: sw

    integer  :: i, j
    !---------------------------------------------------------------------------

    !$omp parallel do &
    !$omp private(ICE_MASS_diff,ICE_MASS_prev,factor,sw)
    do j = ojs, oje
    do i = ois, oie
       sw = 0.5_rp + sign(0.5_rp, ocean_temp(i,j) - ocean_phy_ice_freezetemp) ! 1: melt, 0: freeze
       factor = cl * dwatr * ocean_depth &
              / ( emelt + sw * ci * ( ocean_phy_ice_freezetemp - ice_temp(i,j) ) )

       if ( calc_flag(i,j) ) then
          ! update ice mass
          ice_mass_diff = ( ocean_phy_ice_freezetemp - ocean_temp(i,j) ) * factor ! [kg/m2/s], positive is freezing
          ice_mass_prev = ice_mass(i,j)
          ice_mass(i,j) = min( max( ice_mass(i,j) + ice_mass_diff, 0.0_rp ), ocean_phy_ice_mass_limit ) ! update mass w/ limiter
          ice_mass_diff = ice_mass(i,j) - ice_mass_prev

          ! update ocean temperature
          ocean_temp(i,j) = ocean_temp(i,j) + ice_mass_diff / factor

          ! update ice temperature
          if ( ice_mass(i,j) > 0.0_rp ) then
             if ( ice_mass_diff > 0.0_rp ) then ! when ice increases, new ice have the temperature at freezing point (when the limiter does not work)
                ice_temp(i,j) = ( ice_temp(i,j) * ice_mass_prev &
                                + ocean_phy_ice_freezetemp * ice_mass_diff &
                                + ( ocean_phy_ice_freezetemp - ocean_temp(i,j) ) * ice_mass_diff * cl / ci & ! due to the limiter
                                ) / ice_mass(i,j)
             endif
          else ! fill value
             ice_temp(i,j) = ocean_phy_ice_freezetemp
          endif
       endif
    enddo
    enddo

    return
  end subroutine ocean_phy_ice_adjustment

  !-----------------------------------------------------------------------------
  subroutine ocean_phy_ice_simple( &
       OIA, OIS, OIE, &
       OJA, OJS, OJE, &
       sflx_QV,       &
       sflx_rain,     &
       sflx_snow,     &
       sflx_hbalance, &
       subsfc_temp,   &
       TC_dz,         &
       ICE_TEMP,      &
       ICE_MASS,      &
       calc_flag,     &
       dt,            &
       ICE_TEMP_t,    &
       ICE_MASS_t,    &
       sflx_G,        &
       sflx_water,    &
       sflx_ice       )
    use scale_prc, only: &
       prc_abort
    use scale_const, only: &
       const_ci,    &
       const_emelt
    use scale_file_external_input, only: &
       file_external_input_update
    implicit none

    integer,  intent(in)  :: OIA,   OIS, OIE
    integer,  intent(in)  :: OJA,   OJS, OJE
    real(RP), intent(in)  :: sflx_QV      (oia,oja) ! water vapor flux [kg/m2/s]
    real(RP), intent(in)  :: sflx_rain    (oia,oja) ! rain        flux [kg/m2/s]
    real(RP), intent(in)  :: sflx_snow    (oia,oja) ! snow        flux [kg/m2/s]
    real(RP), intent(in)  :: sflx_hbalance(oia,oja) ! surface heat flux balance [J/m2/s] (negative)
    real(RP), intent(in)  :: subsfc_temp  (oia,oja) ! subsurface temperature [K]
    real(RP), intent(in)  :: TC_dz        (oia,oja) ! Thermal conductance [K/m]
    real(RP), intent(in)  :: ICE_TEMP     (oia,oja) ! sea ice temperature [K]
    real(RP), intent(in)  :: ICE_MASS     (oia,oja) ! sea ice amount      [kg/m2]
    logical,  intent(in)  :: calc_flag    (oia,oja) ! to decide calculate or not
    real(DP), intent(in)  :: dt
    real(RP), intent(out) :: ICE_TEMP_t   (oia,oja) ! tendency of sea ice temperature [K/s]
    real(RP), intent(out) :: ICE_MASS_t   (oia,oja) ! tendency of sea ice amount      [kg/m2/s]
    real(RP), intent(out) :: SFLX_G       (oia,oja) ! heat         flux from sea ice to subsurface
    real(RP), intent(out) :: SFLX_water   (oia,oja) ! liquid water flux from sea ice to subsurface
    real(RP), intent(out) :: SFLX_ice     (oia,oja) ! ice    water flux from sea ice to subsurface

    real(RP) :: mass_budget     ! [kg/m2/s]
    real(RP) :: ICE_MASS_new    ! [kg/m2]
    real(RP) :: ICE_MASS_new2   ! [kg/m2]
    real(RP) :: heat_budget     ! [J/m2/s]
    real(RP) :: heat_budget_new ! [J/m2/s]
    real(RP) :: gflx            ! [J/m2/s]
    real(RP) :: ICE_TEMP_new    ! [K]
    real(RP) :: heatcapacity    ! [J/m2/K]
    real(RP) :: heating         ! [J/m2/s]
    real(RP) :: cooling         ! [J/m2/s]
    real(RP) :: sflx_melt       ! [kg/m2/s]
    real(RP) :: dt_RP

    integer  :: i, j
    !---------------------------------------------------------------------------

    log_progress(*) 'ocean / physics / seaice'

    dt_rp = real(dt,kind=rp)

    !$omp parallel do &
    !$omp private(mass_budget,ICE_MASS_new,ICE_MASS_new2, &
    !$omp         heat_budget,heat_budget_new,gflx,ICE_TEMP_new,heatcapacity,heating,cooling, &
    !$omp         sflx_melt)
    do j = ojs, oje
    do i = ois, oie
       if ( calc_flag(i,j) ) then
          ! mass balance
          mass_budget     = -sflx_qv(i,j) + sflx_rain(i,j) + sflx_snow(i,j)
          ice_mass_new    = max( ice_mass(i,j) + mass_budget * dt_rp, 0.0_rp )
          mass_budget     = mass_budget - ( ice_mass_new - ice_mass(i,j) ) / dt_rp ! residual

          ! heat balance
          gflx = ( subsfc_temp(i,j) - ice_temp(i,j) ) * tc_dz(i,j)
          heat_budget     = sflx_hbalance(i,j)           & ! SWD - SWU + LWD - LWU - SH - LH
                          + gflx                         & ! heat flux from ocean
                          + sflx_rain(i,j) * const_emelt   ! rain freezing

          ice_temp_new    = ice_temp(i,j)
          heatcapacity    = const_ci * ice_mass_new

          ! heating ice
          if ( heat_budget > 0.0_rp .AND. ice_mass_new > 0.0_rp ) then
             heating         = max( ocean_phy_ice_freezetemp - ice_temp(i,j), 0.0_rp ) / dt_rp * heatcapacity
             heat_budget_new = max( heat_budget - heating, 0.0_rp )
             ice_temp_new    = ice_temp(i,j) + ( heat_budget - heat_budget_new ) / heatcapacity * dt_rp
             heat_budget     = heat_budget_new
          endif

          ! cooling ice
          if ( heat_budget < 0.0_rp .AND. ice_mass_new > 0.0_rp ) then
             !! cooling limiter: tiny ice does not cooled faster than thermal conductance of ice, residual heat flux goes ocean
             !ICE_TEMP_new    = max( ICE_TEMP(i,j) + heat_budget / heatcapacity * dt_RP, 0.0_RP )
             !cooling         = max( heat_budget, ( ICE_TEMP_new - ICE_TEMP(i,j) ) * TC_dZ(i,j) )

             ! cooling limiter: tiny ice does not cooled faster than the limit (-0.05K/s), residual heat flux goes ocean
             cooling         = max( heat_budget, -ocean_phy_ice_dt_max * heatcapacity )

             ice_temp_new    = ice_temp(i,j) + cooling / heatcapacity * dt_rp
             heat_budget     = heat_budget - cooling
          endif

          ! melting ice
          sflx_melt       = max( heat_budget / const_emelt, 0.0_rp ) ! only for positive heat flux
          ice_mass_new2   = max( ice_mass_new - sflx_melt * dt_rp, 0.0_rp )
          sflx_melt       = ( ice_mass_new - ice_mass_new2 ) / dt_rp
          heat_budget     = heat_budget - sflx_melt * const_emelt
          mass_budget     = mass_budget + sflx_melt

          ! ice to ocean flux
          sflx_g(i,j) = heat_budget - gflx
          sflx_water(i,j) = mass_budget
          sflx_ice(i,j) = 0.0_rp

          ice_mass_t(i,j) = ( ice_mass_new2 - ice_mass(i,j) ) / dt_rp
          ice_temp_t(i,j) = ( ice_temp_new  - ice_temp(i,j) ) / dt_rp
       else
          sflx_g(i,j) = 0.0_rp
          sflx_water(i,j) = 0.0_rp
          sflx_ice(i,j) = 0.0_rp
          ice_mass_t(i,j) = 0.0_rp
          ice_temp_t(i,j) = 0.0_rp
       endif
    enddo
    enddo

    return
  end subroutine ocean_phy_ice_simple

end module scale_ocean_phy_ice_simple