mod_ocean_driver.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module mod_ocean_driver
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_ocean_grid_cartesc_index
  use scale_tracer
  use scale_cpl_sfc_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: ocean_driver_setup
  public :: ocean_driver_calc_tendency
  public :: ocean_driver_update
  public :: ocean_surface_get
  public :: ocean_surface_set
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), private, allocatable :: WSTR   (:,:)
  real(RP), private, allocatable :: QVEF   (:,:)
  real(RP), private, allocatable :: SR     (:,:)
  real(RP), private, allocatable :: ATMOS_W(:,:)
 
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine ocean_driver_setup
    use scale_prc, only: &
       prc_abort
    use scale_const, only: &
       huge => const_huge
    use mod_ocean_admin, only: &
       ocean_do,       &
       ocean_dyn_type, &
       ocean_sfc_type, &
       ocean_ice_type, &
       ocean_alb_type, &
       ocean_rgn_type
    use scale_cpl_phy_sfc_fixed_temp, only: &
       cpl_phy_sfc_fixed_temp_setup
    use scale_atmos_hydrometeor, only: &
       atmos_hydrometeor_dry
    use scale_ocean_dyn_slab, only: &
       ocean_dyn_slab_setup
    use scale_ocean_dyn_offline, only: &
       ocean_dyn_offline_setup
    use scale_ocean_phy_ice_simple, only: &
       ocean_phy_ice_setup
    use scale_ocean_phy_albedo, only: &
       ocean_phy_albedo_const_setup, &
       ocean_phy_albedo_seaice_setup
    use scale_ocean_phy_albedo_nakajima00, only: &
       ocean_phy_albedo_nakajima00_setup
    use scale_ocean_phy_roughness, only: &
       ocean_phy_roughness_const_setup, &
       ocean_phy_roughness_seaice_setup
    use scale_ocean_phy_roughness_miller92, only: &
       ocean_phy_roughness_miller92_setup
    use scale_ocean_phy_roughness_moon07, only: &
       ocean_phy_roughness_moon07_setup
    use scale_ocean_phy_tc, only: &
       ocean_phy_tc_seaice_setup
    use scale_ocean_grid_cartesc, only: &
       cdz => ocean_grid_cartesc_cdz
    implicit none
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("OCEAN_driver_setup",*) 'Setup'
 
    if ( ocean_do ) then
 
       select case ( ocean_dyn_type )
       case ( 'SLAB' )
          call ocean_dyn_slab_setup( cdz(oks) )
       case ( 'OFFLINE' )
          call ocean_dyn_offline_setup
       case ( 'INIT' )
          ! do nothing
       case default
          log_error("OCEAN_driver_setup",*) 'OCEAN_DYN_TYPE is invalid: ', trim(ocean_dyn_type)
          call prc_abort
       end select
 
       select case ( ocean_sfc_type )
       case ( 'FIXED-TEMP' )
          call cpl_phy_sfc_fixed_temp_setup
       case default
          log_error("OCEAN_driver_setup",*) 'OCEAN_SFC_TYPE is invalid: ', trim(ocean_sfc_type)
          call prc_abort
       end select
 
       select case ( ocean_ice_type )
       case ( 'NONE' )
       case ( 'SIMPLE' )
          call ocean_phy_ice_setup
       case ( 'INIT' )
       case default
          log_error("OCEAN_driver_setup",*) 'OCEAN_ICE_TYPE is invalid: ', trim(ocean_ice_type)
          call prc_abort
       end select
 
       ! surface albedo
       select case ( ocean_alb_type )
       case ( 'NAKAJIMA00' )
          call ocean_phy_albedo_const_setup
          call ocean_phy_albedo_nakajima00_setup
          call ocean_phy_albedo_seaice_setup
       case ( 'CONST' )
          call ocean_phy_albedo_const_setup
       case ( 'INIT' )
          ! do nothing
       case default
          log_error("OCEAN_driver_setup",*) 'OCEAN_ALB_TYPE is invalid: ', trim(ocean_alb_type)
          call prc_abort
       end select
 
       ! surface roughness length
       select case ( ocean_rgn_type )
       case ( 'MILLER92' )
          call ocean_phy_roughness_miller92_setup
          call ocean_phy_roughness_seaice_setup
       case ( 'MOON07' )
          call ocean_phy_roughness_moon07_setup
          call ocean_phy_roughness_seaice_setup
       case ( 'CONST' )
          call ocean_phy_roughness_const_setup
       case ( 'INIT' )
          ! do nothing
       case default
          log_error("OCEAN_driver_setup",*) 'OCEAN_RGN_TYPE is invalid: ', trim(ocean_rgn_type)
          call prc_abort
       end select
 
       ! thermal conductivity
       call ocean_phy_tc_seaice_setup
 
       allocate( wstr(oia,oja) )
       allocate( qvef(oia,oja) )
       allocate( sr(oia,oja) )
       allocate( atmos_w(oia,oja) )
       wstr(:,:) = huge
       if ( atmos_hydrometeor_dry ) then
          qvef(:,:) = 0.0_rp
       else
          qvef(:,:) = 1.0_rp
       end if
       sr(:,:) = 0.0_rp
       atmos_w(:,:) = 0.0_rp ! slope of the sea surface is zero
 
    endif
 
    return
  end subroutine ocean_driver_setup
 
  !-----------------------------------------------------------------------------
  subroutine ocean_driver_calc_tendency( force )
    use scale_const, only: &
       eps => const_eps
    use scale_time, only: &
       dt => time_dtsec_ocean
    use scale_statistics, only: &
       statistics_checktotal, &
       statistics_total
    use scale_landuse, only: &
       exists_ocean => landuse_exists_ocean
    use scale_atmos_grid_cartesc_real, only: &
       real_z1 => atmos_grid_cartesc_real_z1
    use scale_atmos_hydrometeor, only: &
       hydrometeor_lhv => atmos_hydrometeor_lhv, &
       hydrometeor_lhs => atmos_hydrometeor_lhs, &
       atmos_hydrometeor_dry, &
       cv_water,              &
       cv_ice,                &
       lhf,                   &
       i_qv
    use scale_ocean_grid_cartesc_real, only: &
       ocean_grid_cartesc_real_vol,    &
       ocean_grid_cartesc_real_totvol, &
       ocean_grid_cartesc_real_area,   &
       ocean_grid_cartesc_real_totarea
    use scale_ocean_grid_cartesc_real, only: &
       ocean_grid_cartesc_real_area, &
       ocean_grid_cartesc_real_totarea
    use scale_cpl_phy_sfc_fixed_temp, only: &
       cpl_phy_sfc_fixed_temp
    use scale_ocean_phy_albedo, only: &
       ocean_phy_albedo_const, &
       ocean_phy_albedo_seaice
    use scale_ocean_phy_albedo_nakajima00, only: &
       ocean_phy_albedo_nakajima00
    use scale_ocean_phy_roughness, only: &
       ocean_phy_roughness_const, &
       ocean_phy_roughness_seaice
    use scale_ocean_phy_roughness_miller92, only: &
       ocean_phy_roughness_miller92
    use scale_ocean_phy_roughness_moon07, only: &
       ocean_phy_roughness_moon07
    use scale_ocean_phy_tc, only: &
       ocean_phy_tc_seaice
    use scale_ocean_phy_ice_simple, only: &
       ocean_phy_ice_simple
    use scale_bulkflux, only: &
       bulkflux_diagnose_scales
    use mod_atmos_admin, only: &
       atmos_sw_phy_ch
    use mod_atmos_phy_ch_driver, only: &
       atmos_phy_ch_driver_ocean_flux
    use mod_ocean_admin, only: &
       ocean_sfc_type, &
       ocean_ice_type, &
       ocean_alb_type, &
       ocean_rgn_type
    use mod_ocean_vars, only: &
       ice_flag,          &
       ocean_temp,        &
       ocean_ocn_z0m,     &
       ocean_ice_temp,    &
       ocean_ice_mass,    &
       ocean_sfc_temp,    &
       ocean_sfc_albedo,  &
       ocean_sfc_z0m,     &
       ocean_sfc_z0h,     &
       ocean_sfc_z0e,     &
       ocean_temp_t,      &
       ocean_salt_t,      &
       ocean_uvel_t,      &
       ocean_vvel_t,      &
       ocean_ice_temp_t,  &
       ocean_ice_mass_t,  &
       atmos_temp,        &
       atmos_pres,        &
       atmos_u,           &
       atmos_v,           &
       atmos_dens,        &
       atmos_qv,          &
       atmos_pbl,         &
       atmos_cossza,      &
       atmos_sfc_dens,    &
       atmos_sfc_pres,    &
       atmos_sflx_rad_dn, &
       atmos_sflx_water,  &
       atmos_sflx_engi,   &
       ocean_sflx_mw,     &
       ocean_sflx_mu,     &
       ocean_sflx_mv,     &
       ocean_sflx_sh,     &
       ocean_sflx_lh,     &
       ocean_sflx_qtrc,   &
       ocean_u10,         &
       ocean_v10,         &
       ocean_t2,          &
       ocean_q2,          &
       ocean_ustar,       &
       ocean_tstar,       &
       ocean_qstar,       &
       ocean_wstar,       &
       ocean_rlmo,        &
       ocean_ocn_ustar,   &
       ocean_ocn_tstar,   &
       ocean_ocn_qstar,   &
       ocean_ocn_wstar,   &
       ocean_ocn_rlmo,    &
       ocean_ice_ustar,   &
       ocean_ice_tstar,   &
       ocean_ice_qstar,   &
       ocean_ice_wstar,   &
       ocean_ice_rlmo,    &
       ocean_sflx_gh,     &
       ocean_sflx_water,  &
       ocean_sflx_engi,   &
       ocean_oflx_gh,     &
       ocean_oflx_water,  &
       ocean_oflx_engi,   &
       ocean_ice_frac
    use scale_file_history, only: &
       file_history_in
    implicit none
 
    logical, intent(in) :: force
 
    real(rp) :: lhv          (oia,oja)
    real(rp) :: lhs          (oia,oja)
    real(rp) :: atmos_uabs   (oia,oja)
    real(rp) :: sfc_temp     (oia,oja)
    real(rp) :: sfc_albedo   (oia,oja,n_rad_dir,n_rad_rgn)
    real(rp) :: sfc_z0m      (oia,oja)
    real(rp) :: sfc_z0h      (oia,oja)
    real(rp) :: sfc_z0e      (oia,oja)
    real(rp) :: subsfc_temp  (oia,oja)
    real(rp) :: tc_dz        (oia,oja)
    real(rp) :: sflx_mw      (oia,oja)
    real(rp) :: sflx_mu      (oia,oja)
    real(rp) :: sflx_mv      (oia,oja)
    real(rp) :: sflx_sh      (oia,oja)
    real(rp) :: sflx_lh      (oia,oja)
    real(rp) :: sflx_qv      (oia,oja)
    real(rp) :: ocean_sflx_qv(oia,oja)
    real(rp) :: u10          (oia,oja)
    real(rp) :: v10          (oia,oja)
    real(rp) :: t2           (oia,oja)
    real(rp) :: q2           (oia,oja)
    real(rp) :: sflx_hbalance(oia,oja)
    real(rp) :: sflx_gh      (oia,oja)
    real(rp) :: sflx_water   (oia,oja)
    real(rp) :: sflx_engi    (oia,oja)
    real(rp) :: ice_mass     (oia,oja)
    logical  :: exists_ice   (oia,oja)
    real(rp) :: sw
 
    real(rp) :: sfc_frac
 
    integer  :: k, i, j, iq, idir, irgn
    !---------------------------------------------------------------------------
 
    call prof_rapstart('OCN_CalcTend', 1)
 
    !########## Get Surface Boundary from coupler ##########
    call ocean_surface_get
 
    !$omp parallel do
    do j = ojs, oje
    do i = ois, oie
       atmos_uabs(i,j) = sqrt( atmos_u(i,j)**2 + atmos_v(i,j)**2 )
    enddo
    enddo
 
    !$omp parallel do
    do j = ojs, oje
    do i = ois, oie
       exists_ice(i,j) = .false.
       if( exists_ocean(i,j) .AND. ocean_ice_frac(i,j) > 0.0_rp ) exists_ice(i,j) = .true.
    enddo
    enddo
 
    !########## reset tendencies ##########
 
    !$omp parallel do
    do j = ojs, oje
    do i = ois, oie
    do k = oks, oke
       ocean_temp_t(k,i,j) = 0.0_rp
       ocean_salt_t(k,i,j) = 0.0_rp
       ocean_uvel_t(k,i,j) = 0.0_rp
       ocean_vvel_t(k,i,j) = 0.0_rp
    enddo
    enddo
    enddo
 
    do iq = 1, qa
    !$omp parallel do
    do j  = ojs, oje
    do i  = ois, oie
       ocean_sflx_qtrc(i,j,iq) = 0.0_rp
    enddo
    enddo
    enddo
 
    if ( ice_flag ) then
       !$omp parallel do
       do j = ojs, oje
       do i = ois, oie
          ocean_ice_temp_t(i,j) = 0.0_rp
          ocean_ice_mass_t(i,j) = 0.0_rp
       enddo
       enddo
    end if
 
 
 
    !########## surface process (ice-free ocean) ##########
 
    !$omp parallel do
    do j = ojs, oje
    do i = ois, oie
       ocean_sfc_temp(i,j) = ocean_temp(oks,i,j)
    enddo
    enddo
 
    ! albedo
    select case ( ocean_alb_type )
    case ( 'NAKAJIMA00' )
       ! for Near-IR, IR
       call ocean_phy_albedo_const     ( oia, ois, oie,        & ! [IN]
                                         oja, ojs, oje,        & ! [IN]
                                         sfc_albedo(:,:,:,:) ) ! [OUT]
       ! for VIS (overwrite)
       call ocean_phy_albedo_nakajima00( oia, ois, oie,              & ! [IN]
                                         oja, ojs, oje,              & ! [IN]
                                         atmos_cossza(:,:),          & ! [IN]
                                         sfc_albedo(:,:,:,i_r_vis) ) ! [OUT]
    case ( 'CONST' )
       call ocean_phy_albedo_const     ( oia, ois, oie,        & ! [IN]
                                         oja, ojs, oje,        & ! [IN]
                                         sfc_albedo(:,:,:,:) ) ! [OUT]
    case ( 'INIT' )
       ! Never update OCEAN_SFC_albedo from initial condition
       do irgn = i_r_ir, i_r_vis
       do idir = i_r_direct, i_r_diffuse
       do j    = ojs, oje
       do i    = ois, oie
       if ( exists_ocean(i,j) ) then
          sfc_albedo(i,j,idir,irgn) = ocean_sfc_albedo(i,j,idir,irgn)
       end if
       enddo
       enddo
       enddo
       enddo
    end select
 
    ! roughness length
    select case ( ocean_rgn_type )
    case ( 'MILLER92' )
       call ocean_phy_roughness_miller92( oia, ois, oie, & ! [IN]
                                          oja, ojs, oje, & ! [IN]
                                          atmos_uabs(:,:),    & ! [IN]
                                          ocean_ocn_z0m(:,:), & ! [OUT]
                                          ocean_sfc_z0h(:,:), & ! [OUT]
                                          ocean_sfc_z0e(:,:)  ) ! [OUT]
    case ( 'MOON07' )
       call ocean_phy_roughness_moon07  ( oia, ois, oie, & ! [IN]
                                          oja, ojs, oje, & ! [IN]
                                          atmos_uabs(:,:), & ! [IN]
                                          real_z1(:,:), & ! [IN]
                                          exists_ocean(:,:), & ! [IN]
                                          ocean_ocn_z0m(:,:), & ! [INOUT]
                                          ocean_sfc_z0h(:,:), & ! [OUT]
                                          ocean_sfc_z0e(:,:)  ) ! [OUT]
    case ( 'CONST' )
       call ocean_phy_roughness_const   ( oia, ois, oie, & ! [IN]
                                          oja, ojs, oje, & ! [IN]
                                          ocean_ocn_z0m(:,:), & ! [OUT]
                                          ocean_sfc_z0h(:,:), & ! [OUT]
                                          ocean_sfc_z0e(:,:)  ) ! [OUT]
    case ( 'INIT' )
       ! Never update from initial condition
    end select
 
    ! tendency
    select case ( ocean_sfc_type )
    case ( 'FIXED-TEMP' )
 
       call hydrometeor_lhv( oia, ois, oie, oja, ojs, oje, & ! [IN]
                             ocean_sfc_temp(:,:), & ! [IN]
                             lhv(:,:)             ) ! [OUT]
 
       call cpl_phy_sfc_fixed_temp( oia, ois, oie,              & ! [IN]
                                    oja, ojs, oje,              & ! [IN]
                                    atmos_temp(:,:),     & ! [IN]
                                    atmos_pres(:,:),     & ! [IN]
                                    atmos_w(:,:),     & ! [IN]
                                    atmos_u(:,:),     & ! [IN]
                                    atmos_v(:,:),     & ! [IN]
                                    atmos_dens(:,:),     & ! [IN]
                                    atmos_qv(:,:),     & ! [IN]
                                    lhv(:,:),     & ! [IN]
                                    real_z1(:,:),     & ! [IN]
                                    atmos_pbl(:,:),     & ! [IN]
                                    atmos_sfc_dens(:,:),     & ! [IN]
                                    atmos_sfc_pres(:,:),     & ! [IN]
                                    atmos_sflx_rad_dn(:,:,:,:), & ! [IN]
                                    ocean_sfc_temp(:,:),     & ! [IN]
                                    wstr(:,:),     & ! [IN]
                                    qvef(:,:),     & ! [IN]
                                    sfc_albedo(:,:,:,:), & ! [IN]
                                    sr(:,:),     & ! [IN]
                                    ocean_ocn_z0m(:,:),     & ! [IN]
                                    ocean_sfc_z0h(:,:),     & ! [IN]
                                    ocean_sfc_z0e(:,:),     & ! [IN]
                                    exists_ocean(:,:),     & ! [IN]
                                    dt,                         & ! [IN]
                                    ocean_sflx_mw(:,:),     & ! [OUT]
                                    ocean_sflx_mu(:,:),     & ! [OUT]
                                    ocean_sflx_mv(:,:),     & ! [OUT]
                                    ocean_sflx_sh(:,:),     & ! [OUT]
                                    ocean_sflx_lh(:,:),     & ! [OUT]
                                    ocean_sflx_qv(:,:),     & ! [OUT]
                                    ocean_sflx_gh(:,:),     & ! [OUT]
                                    ocean_ocn_ustar(:,:),     & ! [OUT]
                                    ocean_ocn_tstar(:,:),     & ! [OUT]
                                    ocean_ocn_qstar(:,:),     & ! [OUT]
                                    ocean_ocn_wstar(:,:),     & ! [OUT]
                                    ocean_ocn_rlmo(:,:),     & ! [OUT]
                                    ocean_u10(:,:),     & ! [OUT]
                                    ocean_v10(:,:),     & ! [OUT]
                                    ocean_t2(:,:),     & ! [OUT]
                                    ocean_q2(:,:)      ) ! [OUT]
    end select
 
    !$omp parallel do
    do j = ojs, oje
    do i = ois, oie
    if ( exists_ocean(i,j) ) then
       ocean_sflx_water(i,j) = atmos_sflx_water(i,j) - ocean_sflx_qv(i,j)
       ocean_sflx_engi(i,j) = atmos_sflx_engi(i,j)                              & ! internal energy of precipitation
                             - ocean_sflx_qv(i,j) * cv_water * ocean_sfc_temp(i,j) ! internal energy of evaporation water
    end if
    enddo
    enddo
 
    ! weighted average
 
    if ( ice_flag ) then
 
       ! history (open ocean)
       call file_history_in( ocean_u10(:,:), 'OCEAN_OCN_U10', 'Wind velocity u at 10 m on open ocean surface', 'm/s'   )
       call file_history_in( ocean_v10(:,:), 'OCEAN_OCN_V10', 'Wind velocity v at 10 m on open ocean surface', 'm/s'   )
       call file_history_in( ocean_t2(:,:), 'OCEAN_OCN_T2',  'Air temperature at 2m on open ocean surface',   'K'     )
       call file_history_in( ocean_q2(:,:), 'OCEAN_OCN_Q2',  'Specific humidity at 2m on open ocean surface', 'kg/kg' )
 
       !$omp parallel do &
       !$omp private(sfc_frac)
       do j = ojs, oje
       do i = ois, oie
       if ( exists_ocean(i,j) ) then
          sfc_frac = 1.0_rp - ocean_ice_frac(i,j)
 
          ocean_sfc_temp(i,j) = ocean_sfc_temp(i,j) * sfc_frac
          ocean_sfc_z0m(i,j) = ocean_ocn_z0m(i,j) * sfc_frac
          ocean_sfc_z0h(i,j) = ocean_sfc_z0h(i,j) * sfc_frac
          ocean_sfc_z0e(i,j) = ocean_sfc_z0e(i,j) * sfc_frac
          ocean_sflx_mw(i,j) = ocean_sflx_mw(i,j) * sfc_frac
          ocean_sflx_mu(i,j) = ocean_sflx_mu(i,j) * sfc_frac
          ocean_sflx_mv(i,j) = ocean_sflx_mv(i,j) * sfc_frac
          ocean_sflx_sh(i,j) = ocean_sflx_sh(i,j) * sfc_frac
          ocean_sflx_lh(i,j) = ocean_sflx_lh(i,j) * sfc_frac
          ocean_sflx_qv(i,j) = ocean_sflx_qv(i,j) * sfc_frac
          ocean_u10(i,j) = ocean_u10(i,j) * sfc_frac
          ocean_v10(i,j) = ocean_v10(i,j) * sfc_frac
          ocean_t2(i,j) = ocean_t2(i,j) * sfc_frac
          ocean_q2(i,j) = ocean_q2(i,j) * sfc_frac
 
          ocean_sflx_gh(i,j) = ocean_sflx_gh(i,j) * sfc_frac
          ocean_sflx_water(i,j) = ocean_sflx_water(i,j) * sfc_frac
          ocean_sflx_engi(i,j) = ocean_sflx_engi(i,j) * sfc_frac
 
          do irgn = i_r_ir, i_r_vis
          do idir = i_r_direct, i_r_diffuse
             ocean_sfc_albedo(i,j,idir,irgn) = sfc_albedo(i,j,idir,irgn) * sfc_frac
          enddo
          enddo
       end if
       enddo
       enddo
 
    end if
 
 
    !########## surface process (ice) ##########
 
    if ( ice_flag ) then
 
       !$omp parallel do
       do j = ojs, oje
       do i = ois, oie
          subsfc_temp(i,j) = ocean_temp(oks,i,j)
       enddo
       enddo
 
       ! albedo
       select case ( ocean_alb_type )
       case ( 'NAKAJIMA00' )
          call ocean_phy_albedo_seaice( oia, ois, oie,      & ! [IN]
                                        oja, ojs, oje,      & ! [IN]
                                        sfc_albedo(:,:,:,:) ) ! [OUT]
       case ( 'CONST' )
          call ocean_phy_albedo_const ( oia, ois, oie,      & ! [IN]
                                        oja, ojs, oje,      & ! [IN]
                                        sfc_albedo(:,:,:,:) ) ! [OUT]
       case ( 'INIT' )
          ! Never update OCEAN_SFC_albedo from initial condition
          do irgn = i_r_ir, i_r_vis
          do idir = i_r_direct, i_r_diffuse
          do j    = ojs, oje
          do i    = ois, oie
             sfc_albedo(i,j,idir,irgn) = ocean_sfc_albedo(i,j,idir,irgn)
          enddo
          enddo
          enddo
          enddo
       end select
 
       ! roughness length
       select case ( ocean_rgn_type )
       case ( 'MILLER92', 'MOON07' )
          call ocean_phy_roughness_seaice( oia, ois, oie, & ! [IN]
                                           oja, ojs, oje, & ! [IN]
                                           sfc_z0m(:,:),  & ! [OUT]
                                           sfc_z0h(:,:),  & ! [OUT]
                                           sfc_z0e(:,:)   ) ! [OUT]
       case ( 'CONST' )
          call ocean_phy_roughness_const ( oia, ois, oie, & ! [IN]
                                           oja, ojs, oje, & ! [IN]
                                           sfc_z0m(:,:),  & ! [OUT]
                                           sfc_z0h(:,:),  & ! [OUT]
                                           sfc_z0e(:,:)   ) ! [OUT]
       case ( 'INIT' )
          ! Never update OCEAN_SFC_Z0M/H/E from initial condition
          !$omp parallel do
          do j = ojs, oje
          do i = ois, oie
             sfc_z0m(i,j) = ocean_sfc_z0m(i,j)
             sfc_z0h(i,j) = ocean_sfc_z0h(i,j)
             sfc_z0e(i,j) = ocean_sfc_z0e(i,j)
          enddo
          enddo
       end select
 
       ! thermal conductivity / depth
       call ocean_phy_tc_seaice( oia, ois, oie,       & ! [IN]
                                 oja, ojs, oje,       & ! [IN]
                                 ocean_ice_mass(:,:), & ! [IN]
                                 ocean_ice_frac(:,:), & ! [IN]
                                 exists_ice(:,:), & ! [IN]
                                 tc_dz(:,:)  ) ! [OUT]
 
       ! tendency
       select case ( ocean_sfc_type )
       case ( 'FIXED-TEMP' )
          !$omp parallel do private(sw)
          do j = ojs, oje
          do i = ois, oie
          if ( exists_ocean(i,j) ) then
             sw = 0.5_rp + sign(0.5_rp, ocean_ice_frac(i,j)-eps)
             ice_mass(i,j) = ocean_ice_mass(i,j) * sw / ( ocean_ice_frac(i,j) + 1.0_rp - sw )
          end if
          end do
          end do
 
          call hydrometeor_lhs( oia, ois, oie, oja, ojs, oje, & ! [IN]
                                ocean_ice_temp(:,:), & ! [IN]
                                lhs(:,:)             ) ! [OUT]
 
          call cpl_phy_sfc_fixed_temp( oia, ois, oie,              & ! [IN]
                                       oja, ojs, oje,              & ! [IN]
                                       atmos_temp(:,:),     & ! [IN]
                                       atmos_pres(:,:),     & ! [IN]
                                       atmos_w(:,:),     & ! [IN]
                                       atmos_u(:,:),     & ! [IN]
                                       atmos_v(:,:),     & ! [IN]
                                       atmos_dens(:,:),     & ! [IN]
                                       atmos_qv(:,:),     & ! [IN]
                                       lhs(:,:),     & ! [IN]
                                       real_z1(:,:),     & ! [IN]
                                       atmos_pbl(:,:),     & ! [IN]
                                       atmos_sfc_dens(:,:),     & ! [IN]
                                       atmos_sfc_pres(:,:),     & ! [IN]
                                       atmos_sflx_rad_dn(:,:,:,:), & ! [IN]
                                       ocean_ice_temp(:,:),     & ! [IN]
                                       ice_mass(:,:),     & ! [IN]
                                       qvef(:,:),     & ! [IN]
                                       sfc_albedo(:,:,:,:), & ! [IN]
                                       sr(:,:),     & ! [IN]
                                       sfc_z0m(:,:),     & ! [IN]
                                       sfc_z0h(:,:),     & ! [IN]
                                       sfc_z0e(:,:),     & ! [IN]
                                       exists_ice(:,:),     & ! [IN]
                                       dt,                         & ! [IN]
                                       sflx_mw(:,:),     & ! [OUT]
                                       sflx_mu(:,:),     & ! [OUT]
                                       sflx_mv(:,:),     & ! [OUT]
                                       sflx_sh(:,:),     & ! [OUT]
                                       sflx_lh(:,:),     & ! [OUT]
                                       sflx_qv(:,:),     & ! [OUT]
                                       sflx_gh(:,:),     & ! [OUT]
                                       ocean_ice_ustar(:,:),     & ! [OUT]
                                       ocean_ice_tstar(:,:),     & ! [OUT]
                                       ocean_ice_qstar(:,:),     & ! [OUT]
                                       ocean_ice_wstar(:,:),     & ! [OUT]
                                       ocean_ice_rlmo(:,:),     & ! [OUT]
                                       u10(:,:),     & ! [OUT]
                                       v10(:,:),     & ! [OUT]
                                       t2(:,:),     & ! [OUT]
                                       q2(:,:)      ) ! [OUT]
       end select
 
       !$omp parallel do
       do j = ojs, oje
       do i = ois, oie
       if ( exists_ocean(i,j) ) then
          sflx_water(i,j) = atmos_sflx_water(i,j) - sflx_qv(i,j)
          sflx_engi(i,j) = atmos_sflx_engi(i,j)                           & ! internal energy of precipitation
                          - sflx_qv(i,j) * ( cv_ice * ocean_ice_temp(i,j) - lhf ) ! internal energy of evaporation water
       end if
       enddo
       enddo
 
       ! weighted average
       !$omp parallel do
       do j = ojs, oje
       do i = ois, oie
       if ( exists_ocean(i,j) ) then
          ocean_sfc_temp(i,j) = ocean_sfc_temp(i,j) + ocean_ice_temp(i,j) * ocean_ice_frac(i,j)
          ocean_sfc_z0m(i,j) = ocean_sfc_z0m(i,j) + sfc_z0m(i,j) * ocean_ice_frac(i,j)
          ocean_sfc_z0h(i,j) = ocean_sfc_z0h(i,j) + sfc_z0h(i,j) * ocean_ice_frac(i,j)
          ocean_sfc_z0e(i,j) = ocean_sfc_z0e(i,j) + sfc_z0e(i,j) * ocean_ice_frac(i,j)
          ocean_sflx_mw(i,j) = ocean_sflx_mw(i,j) + sflx_mw(i,j) * ocean_ice_frac(i,j)
          ocean_sflx_mu(i,j) = ocean_sflx_mu(i,j) + sflx_mu(i,j) * ocean_ice_frac(i,j)
          ocean_sflx_mv(i,j) = ocean_sflx_mv(i,j) + sflx_mv(i,j) * ocean_ice_frac(i,j)
          ocean_sflx_sh(i,j) = ocean_sflx_sh(i,j) + sflx_sh(i,j) * ocean_ice_frac(i,j)
          ocean_sflx_qv(i,j) = ocean_sflx_qv(i,j) + sflx_qv(i,j) * ocean_ice_frac(i,j)
          ocean_sflx_lh(i,j) = ocean_sflx_lh(i,j) + sflx_lh(i,j) * ocean_ice_frac(i,j)
          ocean_u10(i,j) = ocean_u10(i,j) + u10(i,j) * ocean_ice_frac(i,j)
          ocean_v10(i,j) = ocean_v10(i,j) + v10(i,j) * ocean_ice_frac(i,j)
          ocean_t2(i,j) = ocean_t2(i,j) + t2(i,j) * ocean_ice_frac(i,j)
          ocean_q2(i,j) = ocean_q2(i,j) + q2(i,j) * ocean_ice_frac(i,j)
 
          ocean_oflx_gh(i,j) = ocean_sflx_gh(i,j)
          ocean_oflx_water(i,j) = ocean_sflx_water(i,j)
          ocean_oflx_engi(i,j) = ocean_sflx_engi(i,j)
 
          ocean_sflx_gh(i,j) = ocean_sflx_gh(i,j) + sflx_gh(i,j) * ocean_ice_frac(i,j)
          ocean_sflx_water(i,j) = ocean_sflx_water(i,j) + sflx_water(i,j) * ocean_ice_frac(i,j)
          ocean_sflx_engi(i,j) = ocean_sflx_engi(i,j) + sflx_engi(i,j) * ocean_ice_frac(i,j)
       end if
       end do
       end do
 
       !$omp parallel do
       do irgn = i_r_ir, i_r_vis
       do idir = i_r_direct, i_r_diffuse
       do j    = ojs, oje
       do i    = ois, oie
       if ( exists_ocean(i,j) ) then
          ocean_sfc_albedo(i,j,idir,irgn) = ocean_sfc_albedo(i,j,idir,irgn) + sfc_albedo(i,j,idir,irgn) * ocean_ice_frac(i,j)
       end if
       enddo
       enddo
       enddo
       enddo
 
 
       ! seaice
       select case ( ocean_ice_type )
       case ( 'SIMPLE' )
 
          !$omp parallel do
          do j = ojs, oje
          do i = ois, oie
          if ( exists_ocean(i,j) ) then
             sflx_hbalance(i,j) = sflx_gh(i,j) + sflx_engi(i,j)
          end if
          enddo
          enddo
 
          call ocean_phy_ice_simple( oia, ois, oie,         & ! [IN]
                                     oja, ojs, oje,         & ! [IN]
                                     sflx_water(:,:), & ! [IN]
                                     sflx_hbalance(:,:), & ! [IN]
                                     subsfc_temp(:,:), & ! [IN]
                                     tc_dz(:,:), & ! [IN]
                                     ocean_ice_temp(:,:), & ! [IN]
                                     ocean_ice_mass(:,:), & ! [IN]
                                     ocean_ice_frac(:,:), & ! [IN]
                                     exists_ice(:,:), & ! [IN]
                                     dt,                    & ! [IN]
                                     ocean_ice_temp_t(:,:), & ! [OUT]
                                     ocean_ice_mass_t(:,:), & ! [OUT]
                                     sflx_gh(:,:), & ! [OUT]
                                     sflx_water(:,:), & ! [OUT]
                                     sflx_engi(:,:)  ) ! [OUT]
       case ( 'INIT' )
          !$omp parallel do
          do j = ojs, oje
          do i = ois, oie
          if ( exists_ocean(i,j) ) then
             sflx_gh(i,j) = sflx_gh(i,j) * ocean_ice_frac(i,j)
             sflx_water(i,j) = 0.0_rp ! no flux from seaice to ocean
             sflx_engi(i,j) = 0.0_rp ! no flux from seaice to ocean
          end if
          enddo
          enddo
       end select
 
       ! history (sea ice)
       call file_history_in( u10(:,:), 'OCEAN_ICE_U10', 'Wind velocity u at 10 m on sea ice surface', 'm/s'   )
       call file_history_in( v10(:,:), 'OCEAN_ICE_V10', 'Wind velocity v at 10 m on sea ice surface', 'm/s'   )
       call file_history_in( t2(:,:), 'OCEAN_ICE_T2',  'Air temperature at 2m on sea ice surface',   'K'     )
       call file_history_in( q2(:,:), 'OCEAN_ICE_Q2',  'Specific humidity at 2m on sea ice surface', 'kg/kg' )
 
       ! weighted average
       !$omp parallel do
       do j = ojs, oje
       do i = ois, oie
       if ( exists_ocean(i,j) ) then
          ocean_oflx_gh(i,j) = ocean_oflx_gh(i,j) + sflx_gh(i,j)
          ocean_oflx_water(i,j) = ocean_oflx_water(i,j) + sflx_water(i,j)
          ocean_oflx_engi(i,j) = ocean_oflx_engi(i,j) + sflx_engi(i,j)
       end if
       enddo
       enddo
 
       call bulkflux_diagnose_scales( oia, ois, oie, oja, ojs, oje, &
                                      ocean_sflx_mw(:,:), ocean_sflx_mu(:,:), ocean_sflx_mv(:,:), & ! [IN]
                                      ocean_sflx_sh(:,:), ocean_sflx_qv(:,:),                     & ! [IN]
                                      atmos_sfc_dens(:,:), ocean_sfc_temp(:,:), atmos_pbl(:,:),   & ! [IN]
                                      ocean_ustar(:,:), ocean_tstar(:,:), ocean_qstar(:,:),       & ! [OUT]
                                      ocean_wstar(:,:), ocean_rlmo(:,:),                          & ! [OUT]
                                      mask = exists_ocean(:,:)                                    ) ! [IN]
 
 
 
 
    endif ! ICE process?
 
    if ( .NOT. atmos_hydrometeor_dry ) then
       !$omp parallel do
       do j = ojs, oje
       do i = ois, oie
       if ( exists_ocean(i,j) ) then
          ocean_sflx_qtrc(i,j,i_qv) = ocean_sflx_qv(i,j)
       end if
       enddo
       enddo
    endif
 
 
    ! Surface flux for chemical tracers
    if ( atmos_sw_phy_ch ) then
       call atmos_phy_ch_driver_ocean_flux( ocean_sflx_qtrc(:,:,:) ) ! [INOUT]
    endif
 
    if ( statistics_checktotal ) then
       call statistics_total( oka, oks, oke, oia, ois, oie, oja, ojs, oje, &
                              ocean_temp_t(:,:,:), 'OCEAN_TEMP_t',         &
                              ocean_grid_cartesc_real_vol(:,:,:),          &
                              ocean_grid_cartesc_real_totvol               )
       if ( ice_flag ) then
          call statistics_total( oia, ois, oie, oja, ojs, oje,                &
                                 ocean_ice_temp_t(:,:), 'OCEAN_ICE_TEMP_t',   &
                                 ocean_grid_cartesc_real_area(:,:),           &
                                 ocean_grid_cartesc_real_totarea              )
          call statistics_total( oia, ois, oie, oja, ojs, oje,                &
                                 ocean_ice_mass_t(:,:), 'OCEAN_ICE_MASS_t',   &
                                 ocean_grid_cartesc_real_area(:,:),           &
                                 ocean_grid_cartesc_real_totarea              )
       end if
    endif
 
    call prof_rapend  ('OCN_CalcTend', 1)
 
    !########## Set Surface Boundary to coupler ##########
    call ocean_surface_set( countup=.true. )
 
    return
  end subroutine ocean_driver_calc_tendency
 
  !-----------------------------------------------------------------------------
  subroutine ocean_driver_update
    use scale_time, only: &
       nowdaysec => time_nowdaysec,  &
       dt        => time_dtsec_ocean
    use scale_landuse, only: &
       exists_ocean => landuse_exists_ocean
    use mod_ocean_admin, only: &
       ocean_dyn_type, &
       ocean_ice_type
    use mod_ocean_vars, only: &
       ocean_temp,       &
       ocean_salt,       &
       ocean_uvel,       &
       ocean_vvel,       &
       ocean_ice_temp,   &
       ocean_ice_mass,   &
       ocean_ice_frac,   &
       ocean_temp_t,     &
       ocean_salt_t,     &
       ocean_uvel_t,     &
       ocean_vvel_t,     &
       ocean_ice_temp_t, &
       ocean_ice_mass_t, &
       ocean_oflx_gh,    &
       ocean_oflx_water, &
       ocean_oflx_engi,  &
       ocean_mass_supl,  &
       ocean_engi_supl,  &
       ocean_vars_check
    use scale_ocean_dyn_slab, only: &
       ocean_dyn_slab
    use scale_ocean_dyn_offline, only: &
       ocean_dyn_offline
    use scale_ocean_phy_ice_simple, only: &
       ocean_phy_ice_adjustment, &
       ocean_phy_ice_fraction
    use scale_ocean_grid_cartesc, only: &
       cdz => ocean_grid_cartesc_cdz
    implicit none
 
    real(rp) :: mass_flux(oia,oja)
    real(rp) :: engi_flux(oia,oja)
    real(rp) :: mass_supl(oia,oja)
    real(rp) :: engi_supl(oia,oja)
 
    real(rp) :: sflx_gh(oia,oja)
 
    integer :: i, j
    !---------------------------------------------------------------------------
 
    call prof_rapstart('OCN_Update', 2)
 
    !########## Get Surface Boundary from coupler ##########
    call ocean_surface_get
 
    !########## Dynamics / Update ##########
    !$omp parallel do
    do j = ojs, oje
    do i = ois, oie
       ocean_mass_supl(i,j) = 0.0_rp
       ocean_engi_supl(i,j) = 0.0_rp
    end do
    end do
 
    select case ( ocean_dyn_type )
    case ( 'SLAB' )
 
       !$omp parallel do
       do j = ojs, oje
       do i = ois, oie
       if ( exists_ocean(i,j) ) then
          sflx_gh(i,j) = ocean_oflx_gh(i,j) + ocean_oflx_engi(i,j)
       end if
       end do
       end do
 
       call ocean_dyn_slab( okmax, oks, oke,         & ! [IN]
                            oia,   ois, oie,         & ! [IN]
                            oja,   ojs, oje,         & ! [IN]
                            ocean_temp_t(:,:,:), & ! [IN]
                            sflx_gh(:,:),   & ! [IN]
                            ocean_oflx_water(:,:),   & ! [IN]
                            exists_ocean(:,:),   & ! [IN]
                            dt, nowdaysec,           & ! [IN]
                            ocean_temp(:,:,:), & ! [INOUT]
                            mass_supl(:,:),   & ! [OUT]
                            engi_supl(:,:)    ) ! [OUT]
 
       !$omp parallel do
       do j = ojs, oje
       do i = ois, oie
       if ( exists_ocean(i,j) ) then
          ocean_mass_supl(i,j) = ocean_mass_supl(i,j) + mass_supl(i,j)
          ocean_engi_supl(i,j) = ocean_engi_supl(i,j) + engi_supl(i,j)
       end if
       end do
       end do
 
    case ( 'OFFLINE' )
 
       call ocean_dyn_offline( okmax, oks, oke,         & ! [IN]
                               oia,   ois, oie,         & ! [IN]
                               oja,   ojs, oje,         & ! [IN]
                               exists_ocean(:,:),   & ! [IN]
                               dt, nowdaysec,           & ! [IN]
                               ocean_temp(:,:,:)  ) ! [INOUT]
 
    case ( 'INIT' )
       ! Never update OCEAN_TEMP from initial condition
    end select
 
    !########## Ice / Update ##########
    select case ( ocean_ice_type )
    case ( 'SIMPLE' )
 
       !$omp parallel do
       do j = ojs, oje
       do i = ois, oie
       if ( exists_ocean(i,j) ) then
          ocean_ice_temp(i,j) = ocean_ice_temp(i,j) + ocean_ice_temp_t(i,j) * dt
          ocean_ice_mass(i,j) = ocean_ice_mass(i,j) + ocean_ice_mass_t(i,j) * dt
       end if
       enddo
       enddo
 
       ! ice adjustment
       call ocean_phy_ice_adjustment( oia, ois, oie,           & ! [IN]
                                      oja, ojs, oje,           & ! [IN]
                                      exists_ocean(:,:),     & ! [IN]
                                      cdz(oks),                & ! [IN]
                                      ocean_temp(oks,:,:), & ! [INOUT]
                                      ocean_ice_temp(:,:),     & ! [INOUT]
                                      ocean_ice_mass(:,:),     & ! [INOUT]
                                      mass_flux(:,:),     & ! [OUT]
                                      engi_flux(:,:),     & ! [OUT]
                                      mass_supl(:,:),     & ! [OUT]
                                      engi_supl(:,:)      ) ! [OUT]
 
       !$omp parallel do
       do j = ojs, oje
       do i = ois, oie
       if ( exists_ocean(i,j) ) then
          ocean_oflx_water(i,j) = ocean_oflx_water(i,j) - mass_flux(i,j) / dt
          ocean_oflx_engi(i,j) = ocean_oflx_engi(i,j) - engi_flux(i,j) / dt
          ocean_mass_supl(i,j) = ocean_mass_supl(i,j) + mass_supl(i,j) / dt
          ocean_engi_supl(i,j) = ocean_engi_supl(i,j) + engi_supl(i,j) / dt
       end if
       end do
       end do
 
 
       ! update ice fraction
       call ocean_phy_ice_fraction( oia, ois, oie,       & ! [IN]
                                    oja, ojs, oje,       & ! [IN]
                                    ocean_ice_mass(:,:), & ! [IN]
                                    ocean_ice_frac(:,:)  ) ! [OUT]
 
    case ( 'INIT' )
       ! Never update OCEAN_ICE_TEMP, OCEAN_ICE_MASS, OCEAN_ICE_FRAC from initial condition
    end select
 
    call ocean_vars_check
 
    call prof_rapend  ('OCN_Update', 2)
 
    return
  end subroutine ocean_driver_update
 
  !-----------------------------------------------------------------------------
  subroutine ocean_surface_get
    use mod_ocean_admin, only: &
       ocean_do
    use mod_ocean_vars, only: &
       atmos_temp,        &
       atmos_pres,        &
       atmos_w,           &
       atmos_u,           &
       atmos_v,           &
       atmos_dens,        &
       atmos_qv,          &
       atmos_pbl,         &
       atmos_sfc_dens,    &
       atmos_sfc_pres,    &
       atmos_sflx_rad_dn, &
       atmos_cossza,      &
       atmos_sflx_water,  &
       atmos_sflx_engi
    use mod_cpl_vars, only: &
       cpl_getatm_ocn
    implicit none
    !---------------------------------------------------------------------------
 
    call prof_rapstart('OCN_SfcExch', 2)
 
    if ( ocean_do ) then
       call cpl_getatm_ocn( atmos_temp(:,:),     & ! [OUT]
                            atmos_pres(:,:),     & ! [OUT]
                            atmos_w(:,:),     & ! [OUT]
                            atmos_u(:,:),     & ! [OUT]
                            atmos_v(:,:),     & ! [OUT]
                            atmos_dens(:,:),     & ! [OUT]
                            atmos_qv(:,:),     & ! [OUT]
                            atmos_pbl(:,:),     & ! [OUT]
                            atmos_sfc_dens(:,:),     & ! [OUT]
                            atmos_sfc_pres(:,:),     & ! [OUT]
                            atmos_sflx_rad_dn(:,:,:,:), & ! [OUT]
                            atmos_cossza(:,:),     & ! [OUT]
                            atmos_sflx_water(:,:),     & ! [OUT]
                            atmos_sflx_engi(:,:)      ) ! [OUT]
    endif
 
    call prof_rapend  ('OCN_SfcExch', 2)
 
    return
  end subroutine ocean_surface_get
 
  !-----------------------------------------------------------------------------
  subroutine ocean_surface_set( countup )
    use mod_ocean_admin, only: &
       ocean_do
    use mod_ocean_vars, only: &
       ocean_sfc_temp,   &
       ocean_sfc_albedo, &
       ocean_sfc_z0m,    &
       ocean_sfc_z0h,    &
       ocean_sfc_z0e,    &
       ocean_sflx_mw,    &
       ocean_sflx_mu,    &
       ocean_sflx_mv,    &
       ocean_sflx_sh,    &
       ocean_sflx_lh,    &
       ocean_sflx_qtrc,  &
       ocean_u10,        &
       ocean_v10,        &
       ocean_t2,         &
       ocean_q2,         &
       ocean_sflx_gh
    use mod_cpl_vars, only: &
       cpl_putocn
    use scale_landuse, only: &
       exists_ocean => landuse_exists_ocean
    implicit none
 
    logical, intent(in) :: countup
    !---------------------------------------------------------------------------
 
    call prof_rapstart('OCN_SfcExch', 2)
 
    if ( ocean_do ) then
       call cpl_putocn( ocean_sfc_temp(:,:),     & ! [IN]
                        ocean_sfc_albedo(:,:,:,:), & ! [IN]
                        ocean_sfc_z0m(:,:),     & ! [IN]
                        ocean_sfc_z0h(:,:),     & ! [IN]
                        ocean_sfc_z0e(:,:),     & ! [IN]
                        ocean_sflx_mw(:,:),     & ! [IN]
                        ocean_sflx_mu(:,:),     & ! [IN]
                        ocean_sflx_mv(:,:),     & ! [IN]
                        ocean_sflx_sh(:,:),     & ! [IN]
                        ocean_sflx_lh(:,:),     & ! [IN]
                        ocean_sflx_gh(:,:),     & ! [IN]
                        ocean_sflx_qtrc(:,:,:),   & ! [IN]
                        ocean_u10(:,:),     & ! [IN]
                        ocean_v10(:,:),     & ! [IN]
                        ocean_t2(:,:),     & ! [IN]
                        ocean_q2(:,:),     & ! [IN]
                        exists_ocean(:,:),     & ! [IN]
                        countup                    ) ! [IN]
    endif
 
    call prof_rapend  ('OCN_SfcExch', 2)
 
    return
  end subroutine ocean_surface_set
 
end module mod_ocean_driver