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_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              = ''
    logical                :: ocean_phy_ice_nudging_basename_add_num      = .false.
    integer                :: ocean_phy_ice_nudging_number_of_files       = 1
    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_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_basename_add_num,      &
!        OCEAN_PHY_ICE_nudging_number_of_files,       &
!        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_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 == '' ) 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_phy_ice_nudging_basename_add_num,      & ! [IN]
                                        ocean_phy_ice_nudging_number_of_files,       & ! [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_defval,                & ! [IN]
                                        check_coordinates = ocean_phy_ice_nudging_check_coordinates, & ! [IN]
                                        step_limit        = 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
    !---------------------------------------------------------------------------
    !$acc data copyin(ICE_MASS) copyout(ICE_FRAC)
 
    !$acc kernels
    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
    !$acc end kernels
 
    !$acc end data
    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,      &
       MASS_FLUX,     &
       ENGI_FLUX,     &
       MASS_SUPL,     &
       ENGI_SUPL      )
    use scale_const, only: &
       dwatr => const_dwatr, &
       eps   => const_eps
    use scale_atmos_hydrometeor, only: &
       cv_water, &
       cv_ice,   &
       lhf
    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), intent(out)   :: mass_flux(oia,oja)
    real(rp), intent(out)   :: engi_flux(oia,oja)
    real(rp), intent(out)   :: mass_supl (oia,oja)
    real(rp), intent(out)   :: engi_supl (oia,oja)
 
    real(rp) :: c_w
    real(rp) :: ice_mass_frz
    real(rp) :: ice_mass_prev
 
    integer  :: i, j
    !---------------------------------------------------------------------------
    !$acc data copy(OCEAN_TEMP,ICE_TEMP,ICE_MASS) &
    !$acc      copyin(calc_flag) &
    !$acc      copyout(MASS_FLUX,ENGI_FLUX,MASS_SUPL,ENGI_SUPL)
 
    c_w = cv_water * dwatr * ocean_depth
 
    !$acc kernels
    !$omp parallel do &
    !$omp private(ICE_MASS_frz,ICE_MASS_prev)
    do j = ojs, oje
    do i = ois, oie
       if (       calc_flag(i,j) &
            .and. ocean_temp(i,j) < ocean_phy_ice_freezetemp .and. ice_mass(i,j) < ocean_phy_ice_mass_limit ) then
          ice_mass_frz = c_w * ( ocean_phy_ice_freezetemp - ocean_temp(i,j) ) &
                       / ( ( cv_water - cv_ice ) * ocean_phy_ice_freezetemp + lhf )
          ice_mass_frz = min( ice_mass_frz, dwatr * ocean_depth )
 
          ! update ice mass
          ice_mass_prev = ice_mass(i,j)
          ice_mass(i,j) = ice_mass(i,j) + ice_mass_frz
          ice_mass(i,j) = min( ice_mass(i,j), ocean_phy_ice_mass_limit ) ! apply limiter
          ice_mass_frz = ice_mass(i,j) - ice_mass_prev
 
          ! update ice temperature
          if ( ice_mass(i,j) > eps ) then
             ice_temp(i,j) = ice_temp(i,j) &
                  + ( ocean_phy_ice_freezetemp - ice_temp(i,j) ) * ice_mass_frz / ice_mass(i,j)
          else ! all ice melt
             ice_temp(i,j) = ocean_phy_ice_freezetemp ! dummy
          end if
 
          ! update ocean temperature
          if ( c_w - cv_water * ice_mass_frz > eps ) then
             ocean_temp(i,j) = ocean_temp(i,j) &
                  + ( cv_water * ocean_temp(i,j) - cv_ice * ocean_phy_ice_freezetemp + lhf ) * ice_mass_frz &
                  / ( c_w - cv_water * ice_mass_frz )
          else ! all water freeze
             ocean_temp(i,j) = ocean_phy_ice_freezetemp
          end if
 
          mass_flux(i,j) = ice_mass_frz
          engi_flux(i,j) = ( cv_ice * ocean_phy_ice_freezetemp - lhf ) * ice_mass_frz
          mass_supl(i,j) = ice_mass_frz
          engi_supl(i,j) = ice_mass_frz * cv_water * ocean_temp(i,j)
       else
          mass_flux(i,j) = 0.0_rp
          engi_flux(i,j) = 0.0_rp
          mass_supl(i,j) = 0.0_rp
          engi_supl(i,j) = 0.0_rp
       endif
    enddo
    enddo
    !$acc end kernels
 
    !$acc end data
    return
  end subroutine ocean_phy_ice_adjustment
 
  !-----------------------------------------------------------------------------
  subroutine ocean_phy_ice_simple( &
       OIA, OIS, OIE, &
       OJA, OJS, OJE, &
       iflx_water,    &
       iflx_hbalance, &
       subsfc_temp,   &
       TC_dz,         &
       ICE_TEMP,      &
       ICE_MASS,      &
       ICE_FRAC,      &
       calc_flag,     &
       dt,            &
       ICE_TEMP_t,    &
       ICE_MASS_t,    &
       sflx_G,        &
       sflx_water,    &
       sflx_RHOE      )
    use scale_prc, only: &
       prc_abort
    use scale_atmos_hydrometeor, only: &
       cv_water, &
       cv_ice,   &
       lhf
    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)  :: iflx_water   (oia,oja) ! input mass flux [kg/m2/s] (downward)
    real(rp), intent(in)  :: iflx_hbalance(oia,oja) ! input heat flux [J/m2/s]  (downward)
    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]
    real(rp), intent(in)  :: ice_frac     (oia,oja) ! sea ice fraction    [0-1]
    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) ! mass flux from sea ice to subsurface
    real(rp), intent(out) :: sflx_rhoe    (oia,oja) ! internal energy flux from sea ice to subsurface
 
    real(rp) :: ice_mass_new    ! [kg/m2]
    real(rp) :: ice_temp_new    ! [K]
    real(rp) :: mass_budget     ! [kg/m2/s]
    real(rp) :: heat_budget     ! [J/m2/s]
    real(rp) :: g               ! [J/m2/s]
    real(rp) :: dm              ! [kg/m2]
    real(rp) :: de              ! [J/m2]
    real(rp) :: m_mlt           ! [kg/m2]
    real(rp) :: dt_rp
 
    integer  :: i, j
    !---------------------------------------------------------------------------
    !$acc data copyin (iflx_water,iflx_hbalance,subsfc_temp, &
    !$acc              TC_dz,ICE_TEMP,ICE_MASS,ICE_FRAC,calc_flag) &
    !$acc      copyout(ICE_TEMP_t,ICE_MASS_t,SFLX_G,SFLX_water,SFLX_RHOE)
 
    log_progress(*) 'ocean / physics / seaice'
 
    dt_rp = real(dt,kind=rp)
 
    !$acc kernels
    !$omp parallel do &
    !$omp private(mass_budget,heat_budget,dM,dE,G,M_mlt, &
    !$omp         ICE_TEMP_new,ICE_MASS_new)
    do j = ojs, oje
    do i = ois, oie
       if ( calc_flag(i,j) ) then
 
          ! mass change
          dm  = iflx_water(i,j) * ice_frac(i,j) * dt_rp
          ice_mass_new = ice_mass(i,j) + dm
 
          if ( ice_mass_new > 0.0_rp ) then
             ! internal energy change
             g = ( subsfc_temp(i,j) - ice_temp(i,j) ) * tc_dz(i,j) ! heat flux from ocean
             de = ( iflx_hbalance(i,j) + g ) * ice_frac(i,j) * dt_rp
             ice_temp_new = ice_temp(i,j) &
                          + ( de - ( cv_ice * ice_temp(i,j) - lhf ) * dm ) / ( cv_ice * ice_mass_new )
 
             ! melting ice
             m_mlt = cv_ice * ( ice_temp_new - ocean_phy_ice_freezetemp ) * ice_mass_new &
                   / ( ( cv_water - cv_ice ) * ocean_phy_ice_freezetemp + lhf )
             m_mlt = min( max( m_mlt, 0.0_rp ), ice_mass_new )
 
             ice_mass_new = ice_mass_new - m_mlt
             ice_temp_new = ice_temp_new &
                          + ( cv_ice * ice_temp_new - lhf - cv_water * ocean_phy_ice_freezetemp ) * m_mlt &
                          / ( cv_ice * ice_mass_new )
 
             ! ice to ocean flux
             mass_budget = m_mlt / dt_rp
             sflx_rhoe(i,j) = cv_water * ocean_phy_ice_freezetemp * mass_budget
             sflx_g(i,j) = - g * ice_frac(i,j)
             sflx_water(i,j) = mass_budget
 
          else
 
             ice_mass_new = 0.0_rp
             ice_temp_new = ocean_phy_ice_freezetemp ! dummy
 
             sflx_rhoe(i,j) = cv_water * subsfc_temp(i,j) * ice_mass_new / dt_rp
             sflx_g(i,j) = ( cv_ice * ice_temp(i,j) - lhf ) * ice_mass(i,j) &
                             - sflx_rhoe(i,j) + iflx_hbalance(i,j) * ice_frac(i,j)
             sflx_water(i,j) = ice_mass_new ! (negative)
 
          endif
 
          ice_mass_t(i,j) = ( ice_mass_new - 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_rhoe(i,j) = 0.0_rp
          ice_mass_t(i,j) = 0.0_rp
          ice_temp_t(i,j) = 0.0_rp
       endif
    enddo
    enddo
    !$acc end kernels
 
    !$acc end data
    return
  end subroutine ocean_phy_ice_simple
 
end module scale_ocean_phy_ice_simple