scale_ocean_sfc_slab.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
module scale_ocean_sfc_slab
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_grid_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: ocean_sfc_slab_setup
  public :: ocean_sfc_slab
  public :: ocean_sfc_slab_simplealbedo

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
  logical, private :: sst_update

  logical, allocatable, private :: is_ocn(:,:)

contains
  !-----------------------------------------------------------------------------
  subroutine ocean_sfc_slab_setup( OCEAN_TYPE )
    use scale_process, only: &
       prc_mpistop
    use scale_landuse, only: &
       landuse_fact_ocean
    implicit none

    character(len=*), intent(in) :: OCEAN_TYPE

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

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[SLAB] / Categ[OCEAN SFC] / Origin[SCALElib]'

    select case( ocean_type )
    case ( 'CONST' )
       sst_update = .false.
    case ( 'SLAB', 'FILE' )
       sst_update = .true.
    case default
       write(*,*) 'xxx wrong OCEAN_TYPE. Check!'
       call prc_mpistop
    end select

    ! judge to run slab ocean model
    allocate( is_ocn(ia,ja) )

    do j = js, je
    do i = is, ie
      if( landuse_fact_ocean(i,j) > 0.0_rp ) then
        is_ocn(i,j) = .true.
      else
        is_ocn(i,j) = .false.
      end if
    end do
    end do

    return
  end subroutine ocean_sfc_slab_setup

  !-----------------------------------------------------------------------------
  subroutine ocean_sfc_slab( &
        SST_t,  & ! [OUT]
        ZMFLX,  & ! [OUT]
        XMFLX,  & ! [OUT]
        YMFLX,  & ! [OUT]
        SHFLX,  & ! [OUT]
        LHFLX,  & ! [OUT]
        WHFLX,  & ! [OUT]
        U10,    & ! [OUT]
        V10,    & ! [OUT]
        T2,     & ! [OUT]
        Q2,     & ! [OUT]
        TMPA,   & ! [IN]
        PRSA,   & ! [IN]
        WA,     & ! [IN]
        UA,     & ! [IN]
        VA,     & ! [IN]
        RHOA,   & ! [IN]
        QVA,    & ! [IN]
        Z1,     & ! [IN]
        PBL,    & ! [IN]
        PRSS,   & ! [IN]
        LWD,    & ! [IN]
        SWD,    & ! [IN]
        TW,     & ! [IN]
        SST,    & ! [IN]
        ALB_LW, & ! [IN]
        ALB_SW, & ! [IN]
        Z0M,    & ! [IN]
        Z0H,    & ! [IN]
        Z0E,    & ! [IN]
        dt      ) ! [IN]
    use scale_grid_index
    use scale_process, only: &
       prc_mpistop
    use scale_const, only: &
      cpdry => const_cpdry, &
      stb   => const_stb
    use scale_bulkflux, only: &
      bulkflux
    use scale_atmos_thermodyn, only: &
       atmos_thermodyn_templhv
    use scale_atmos_saturation, only: &
      qsat => atmos_saturation_pres2qsat_all
    implicit none

    ! arguments
    real(RP), intent(out) :: SST_t(ia,ja) ! tendency of sea surface temperature
    real(RP), intent(out) :: ZMFLX(ia,ja) ! z-momentum flux at the surface [kg/m2/s]
    real(RP), intent(out) :: XMFLX(ia,ja) ! x-momentum flux at the surface [kg/m2/s]
    real(RP), intent(out) :: YMFLX(ia,ja) ! y-momentum flux at the surface [kg/m2/s]
    real(RP), intent(out) :: SHFLX(ia,ja) ! sensible heat flux at the surface [W/m2]
    real(RP), intent(out) :: LHFLX(ia,ja) ! latent heat flux at the surface [W/m2]
    real(RP), intent(out) :: WHFLX(ia,ja) ! water heat flux at the surface [W/m2]
    real(RP), intent(out) :: U10  (ia,ja) ! velocity u at 10m [m/s]
    real(RP), intent(out) :: V10  (ia,ja) ! velocity v at 10m [m/s]
    real(RP), intent(out) :: T2   (ia,ja) ! temperature at 2m [K]
    real(RP), intent(out) :: Q2   (ia,ja) ! water vapor at 2m [kg/kg]

    real(RP), intent(in) :: TMPA(ia,ja) ! temperature at the lowest atmospheric layer [K]
    real(RP), intent(in) :: PRSA(ia,ja) ! pressure at the lowest atmospheric layer [Pa]
    real(RP), intent(in) :: WA  (ia,ja) ! velocity w at the lowest atmospheric layer [m/s]
    real(RP), intent(in) :: UA  (ia,ja) ! velocity u at the lowest atmospheric layer [m/s]
    real(RP), intent(in) :: VA  (ia,ja) ! velocity v at the lowest atmospheric layer [m/s]
    real(RP), intent(in) :: RHOA(ia,ja) ! density at the lowest atmospheric layer [kg/m3]
    real(RP), intent(in) :: QVA (ia,ja) ! ratio of water vapor mass to total mass at the lowest atmospheric layer [kg/kg]
    real(RP), intent(in) :: Z1  (ia,ja) ! cell center height at the lowest atmospheric layer [m]
    real(RP), intent(in) :: PBL (ia,ja) ! the top of atmospheric mixing layer [m]
    real(RP), intent(in) :: PRSS(ia,ja) ! pressure at the surface [Pa]
    real(RP), intent(in) :: LWD (ia,ja) ! downward long-wave radiation flux at the surface [W/m2]
    real(RP), intent(in) :: SWD (ia,ja) ! downward short-wave radiation flux at the surface [W/m2]

    real(RP), intent(in) :: TW    (ia,ja) ! water temperature [K]
    real(RP), intent(in) :: SST   (ia,ja) ! sea surface temperature [K]
    real(RP), intent(in) :: ALB_LW(ia,ja) ! surface albedo for LW [0-1]
    real(RP), intent(in) :: ALB_SW(ia,ja) ! surface albedo for SW [0-1]
    real(RP), intent(in) :: Z0M   (ia,ja) ! roughness length for momentum [m]
    real(RP), intent(in) :: Z0H   (ia,ja) ! roughness length for heat [m]
    real(RP), intent(in) :: Z0E   (ia,ja) ! roughness length for vapor [m]
    real(DP), intent(in) :: dt            ! delta time

    ! works
    real(RP) :: SST1(ia,ja)

    real(RP) :: Ustar        ! friction velocity [m]
    real(RP) :: Tstar        ! friction temperature [K]
    real(RP) :: Qstar        ! friction mixing rate [kg/kg]
    real(RP) :: Uabs         ! modified absolute velocity [m/s]

    real(RP) :: SQV          ! saturation water vapor mixing ratio at surface [kg/kg]
    real(RP) :: LHV(ia,ja)   ! latent heat for vaporization depending on temperature [J/kg]

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

    ! update surface temperature
    do j = js, je
    do i = is, ie
      sst1(i,j) = tw(i,j) ! assumed well-mixed condition
    end do
    end do

    call atmos_thermodyn_templhv( lhv, tmpa )

    do j = js, je
    do i = is, ie

      if( is_ocn(i,j) ) then

        ! saturation at the surface
        call qsat( sqv,       & ! [OUT]
                   sst1(i,j), & ! [IN]
                   prss(i,j)  ) ! [IN]

        call bulkflux( &
            ustar,     & ! [OUT]
            tstar,     & ! [OUT]
            qstar,     & ! [OUT]
            uabs,      & ! [OUT]
            tmpa(i,j), & ! [IN]
            sst1(i,j), & ! [IN]
            prsa(i,j), & ! [IN]
            prss(i,j), & ! [IN]
            qva(i,j), & ! [IN]
            sqv,       & ! [IN]
            ua(i,j), & ! [IN]
            va(i,j), & ! [IN]
            z1(i,j), & ! [IN]
            pbl(i,j), & ! [IN]
            z0m(i,j), & ! [IN]
            z0h(i,j), & ! [IN]
            z0e(i,j)  ) ! [IN]

        zmflx(i,j) = -rhoa(i,j) * ustar**2 / uabs * wa(i,j)
        xmflx(i,j) = -rhoa(i,j) * ustar**2 / uabs * ua(i,j)
        ymflx(i,j) = -rhoa(i,j) * ustar**2 / uabs * va(i,j)

        shflx(i,j) = -cpdry    * rhoa(i,j) * ustar * tstar
        lhflx(i,j) = -lhv(i,j) * rhoa(i,j) * ustar * qstar

        ! calculation for residual
        whflx(i,j) = ( 1.0_rp - alb_sw(i,j) ) * swd(i,j) * (-1.0_rp) &
                   - ( 1.0_rp - alb_lw(i,j) ) * ( lwd(i,j) - stb * sst1(i,j)**4 ) &
                   + shflx(i,j) + lhflx(i,j)

        ! diagnositc variables
        u10(i,j) = ua(i,j) * log( 10.0_rp / z0m(i,j) ) / log( z1(i,j) / z0m(i,j) )
        v10(i,j) = va(i,j) * log( 10.0_rp / z0m(i,j) ) / log( z1(i,j) / z0m(i,j) )
        t2(i,j) = sst1(i,j) + ( tmpa(i,j) - sst1(i,j) ) * ( log(  2.0_rp / z0m(i,j) ) * log(  2.0_rp / z0h(i,j) ) ) &
                                                         / ( log( z1(i,j) / z0m(i,j) ) * log( z1(i,j) / z0h(i,j) ) )
        q2(i,j) = sqv       + (  qva(i,j) - sqv       ) * ( log(  2.0_rp / z0m(i,j) ) * log(  2.0_rp / z0e(i,j) ) ) &
                                                         / ( log( z1(i,j) / z0m(i,j) ) * log( z1(i,j) / z0e(i,j) ) )

      else
        ! not calculate surface flux
        zmflx(i,j) = 0.0_rp
        xmflx(i,j) = 0.0_rp
        ymflx(i,j) = 0.0_rp
        shflx(i,j) = 0.0_rp
        lhflx(i,j) = 0.0_rp
        whflx(i,j) = 0.0_rp
        u10(i,j) = 0.0_rp
        v10(i,j) = 0.0_rp
        t2(i,j) = 0.0_rp
        q2(i,j) = 0.0_rp
      end if

    enddo
    enddo

    ! calculate tendency
    if( sst_update ) then
      do j = js, je
      do i = is, ie
        sst_t(i,j) = ( sst1(i,j) - sst(i,j) ) / dt
      enddo
      enddo
    else
      do j = js, je
      do i = is, ie
        sst_t(i,j) = 0.0_rp
      enddo
      enddo
    end if

    return
  end subroutine ocean_sfc_slab

  !-----------------------------------------------------------------------------
  subroutine ocean_sfc_slab_simplealbedo( &
       SFC_albedo_t, &
       SFC_albedo,   &
       cosSZA,       &
       dt            )
    use scale_grid_index
    use scale_const, only: &
       i_sw  => const_i_sw, &
       i_lw  => const_i_lw
    implicit none

    ! arguments
    real(RP), intent(out) :: SFC_albedo_t(ia,ja,2) ! tendency of sea surface albedo [0-1]
    real(RP), intent(in)  :: SFC_albedo  (ia,ja,2) ! sea surface                    [0-1]
    real(RP), intent(in)  :: cosSZA      (ia,ja)   ! cos(solar zenith angle)        [0-1]
    real(DP), intent(in)  :: dt                    ! delta time

    ! works
    real(RP) :: SFC_albedo1(ia,ja,2)
    real(RP) :: am1

    real(RP), parameter :: c_ocean_albedo(3) = (/ -0.747900_rp, &
                                                  -4.677039_rp, &
                                                   1.583171_rp  /)

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

    do j = js, je
    do i = is, ie
       am1 = max( min( cossza(i,j), 0.961_rp ), 0.0349_rp )

       ! SFC_albedo1(i,j,I_LW) = 0.5_RP do nothing
       sfc_albedo1(i,j,i_sw) = exp( ( c_ocean_albedo(3)*am1 + c_ocean_albedo(2) )*am1 + c_ocean_albedo(1) )
    enddo
    enddo

    ! calculate tendency
    do j = js, je
    do i = is, ie
       sfc_albedo_t(i,j,i_lw) = 0.0_rp
       sfc_albedo_t(i,j,i_sw) = ( sfc_albedo1(i,j,i_sw) - sfc_albedo(i,j,i_sw) ) / dt
    enddo
    enddo

    return
  end subroutine ocean_sfc_slab_simplealbedo

end module scale_ocean_sfc_slab