scale_atm_phy_rd_solarins_simple.F90

Go to the documentation of this file.

!> module FElib / Atmosphere / Physics radiation / Solar insolation / Simple gray-radiation scheme
!!
!! @par Description
!!      A module for idealized solar insolation scheme
!!
!! @author Yuta Kawai, Team SCALE
!!
!! @par Reference
!!  - Frierson, D.M.W., Held, I.M., and Zurita-Gotor, P. 2006:
!!    A gray-radiation aquaplanet moist GCM. Part I: Static stability and eddy scale.
!!    J. Atmos. Sci., 63, 2548–2566.
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_phy_rd_solarins_simple
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
 
  use scale_const, only: &
    undef => const_undef, &
    pi => const_pi
 
  !-----------------------------------------------------------------------------
  implicit none
  private  
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
  public :: atm_phy_rd_solarins_simple_setup
  public :: atm_phy_rd_solarins_simple_update
  public :: atm_phy_rd_solarins_simple_get
 
  !-----------------------------------------------------------------------------
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !++ Private procedure
  !
  private :: annual_mean_insol
  private :: daily_mean_insol
  private :: days_in_month
  private :: is_leap_year
  private :: seconds_to_hms
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer :: SOLARINS_TYPE_ID !< Solar insolation type ID
 
  integer, parameter :: SOLARINS_SIMPLE_TYPE_ID_CONST        = 1 !< Type ID for constant solar insolation
  integer, parameter :: SOLARINS_SIMPLE_TYPE_ID_ANNUAL_MEAN  = 2 !< Type ID for annual mean solar insolation
  integer, parameter :: SOLARINS_SIMPLE_TYPE_ID_FRIERSON2006 = 3 !< Type ID for an idealized solar insolation based on Frierson et al. (2006)
 
 
  real(RP) :: CONST_FLUX            = 340.0_rp !< Constant solar flux value [W/m2]
  
  integer :: ORBIT_REFERENCE_YEAR   = 2000   !< Reference year for orbital calculations
  integer :: ANNUAL_MEAN_YEAR       = 2000   !< Year for computing annual mean insolation
  integer :: ANNUAL_SAMPLES_PER_DAY = 1      !< Number of samples per day for annual mean calculation
  logical :: CACHE_ANNUAL_MEAN      = .true. !< Flag to cache annual mean insolation
 
  integer, parameter :: DATE_YEAR   = 1
  integer, parameter :: DATE_MONTH  = 2
  integer, parameter :: DATE_DAY    = 3
  integer, parameter :: DATE_HOUR   = 4
  integer, parameter :: DATE_MINUTE = 5
  integer, parameter :: DATE_SECOND = 6
  integer, parameter :: DATE_SIZE = 6
 
  !-
  real(RP) :: FRIERSON2006_DELTA_s = 1.4e0_rp !< Parameter for  latitudinal variation (Delta_s) in Frierson et al. (2006) insolation scheme
 
contains
  !> Setup the simplified solar insolation module. 
  subroutine atm_phy_rd_solarins_simple_setup()
    use scale_atmos_solarins, only: &
        atmos_solarins_setup    
    implicit none
    character(len=H_SHORT) :: solarins_type = 'ANNUAL_MEAN'
    namelist / param_atmos_solarins_simple / &
        solarins_type,                    &
        const_flux,                       &
        orbit_reference_year,             &
        annual_mean_year,                 &
        annual_samples_per_day,           &
        cache_annual_mean,                &
        frierson2006_delta_s
 
    integer :: ierr
    !------------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_RD_SOLARINS_SIMPLE_setup",*) 'Setup'
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_solarins_simple,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_PHY_RD_SOLARINS_SIMPLE_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_PHY_RD_SOLARINS_SIMPLE_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_SOLARINS_SIMPLE. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_solarins_simple)
 
    select case (solarins_type)
    case ('CONST')
       solarins_type_id = solarins_simple_type_id_const
    case ('ANNUAL_MEAN')
       solarins_type_id = solarins_simple_type_id_annual_mean
    case ('FRIERSON2006')
       solarins_type_id = solarins_simple_type_id_frierson2006
    case default
       log_error("ATMOS_PHY_RD_SOLARINS_SIMPLE_setup",*) 'Not appropriate SOLARINS_TYPE. Check!'
       call prc_abort
    end select
 
    !- Initializes orbital parameters, the vernal-equinox reference date, and related internal state.
    !  Here, longitude and latitude are irrelevant because this module uses only the orbital state returned by ecliptic_longitude.
    call atmos_solarins_setup( 0.0_rp, 0.0_rp, orbit_reference_year )
 
    return
  end subroutine atm_phy_rd_solarins_simple_setup
 
!OCL SERIAL
  subroutine atm_phy_rd_solarins_simple_update()
  end subroutine atm_phy_rd_solarins_simple_update
 
  !> Get solar insolation and cosine of the solar zenith angle for the given latitude array.
!OCL SERIAL
  subroutine atm_phy_rd_solarins_simple_get( solins, cosSZA, &
    lat, Np )
    use scale_atmos_solarins, only: &
        atmos_solarins_constant
    implicit none
    integer, intent(in) :: np           !< Number of points
    real(rp), intent(out) :: solins(np) !< Solar insolation [W/m2] for each latitude point
    real(rp), intent(out) :: cossza(np) !< Cosine of the solar zenith angle for each latitude point
    real(rp), intent(in) :: lat(np)     !< Array storing latitude values [radians] for each point
    !------------------------------------------------------------------------------
 
    select case (solarins_type_id)
    case (solarins_simple_type_id_const)
      solins(:) = const_flux
      cossza(:) = undef
    case (solarins_simple_type_id_annual_mean)
      call annual_mean_insol( solins, & ! (out)
        lat, np )                       ! (in)
      cossza(:) = undef
    case (solarins_simple_type_id_frierson2006)
      solins(:) = 0.25_rp * atmos_solarins_constant &
                * ( 1.0_rp + 0.25_rp * frierson2006_delta_s * ( 1.0_rp - 3.0_rp * sin(lat(:))**2 ) )
      cossza(:) = undef
    end select
    return
  end subroutine atm_phy_rd_solarins_simple_get
 
!- Private subroutines ------------------------------------
  subroutine annual_mean_insol( solins, &
    lat, Np )
    use scale_atmos_solarins, only: &
        atmos_solarins_ecliptic_longitude
    implicit none
    integer, intent(in) :: np
    real(rp), intent(out) :: solins(np)
    real(rp), intent(in) :: lat(np)
 
    integer :: date(date_size)    
    integer :: month, day
    integer :: ndays_month
    integer :: nsample
 
    real(rp) :: re_factor
    real(rp) :: sin_decl
    real(rp) :: cos_decl
    real(rp) :: hour_angle
    integer, parameter :: offset_year = 0
 
    integer :: i
 
    real(rp) :: sample_second
 
    real(rp) :: daily_solins(np)
    real(rp) :: sum_flux(np)
    !------------------------------------------------------------------------------
 
    date(date_year)   = annual_mean_year
    nsample = 1
    sum_flux(:) = 0.0_rp
 
    do month=1, 12
      ndays_month = days_in_month( annual_mean_year, month )
      date(date_month) = month
      do day=1, ndays_month
        date(date_day) = day
        do i=1, annual_samples_per_day
          ! Midpoint sampling within each day.
          sample_second = 86400.0_rp                                       &
               * ( real(i,rp) - 0.5_rp ) / real(annual_samples_per_day,rp)
        
          call seconds_to_hms( sample_second,                        & ! (in)
               date(date_hour), date(date_minute), date(date_second) ) ! (out)
        
          call atmos_solarins_ecliptic_longitude( &
            re_factor, sin_decl, cos_decl, hour_angle, & ! (out)
            date, offset_year                          ) ! (in)
        
          call daily_mean_insol( daily_solins,     & ! (out)
            lat, re_factor, sin_decl, cos_decl, np ) ! (in)
    
          sum_flux(:) = sum_flux(:) + daily_solins(:)
          nsample = nsample + 1
        end do
      end do
    end do
 
    solins(:) = sum_flux(:) / real(nsample, rp)
    return
  end subroutine annual_mean_insol
 
!OCL SERIAL
  subroutine daily_mean_insol( solins, &
    lat, re_factor, sin_decl, cos_decl, Np )
    use scale_atmos_solarins, only: &
        atmos_solarins_constant
    implicit none
    integer, intent(in) :: np
    real(rp), intent(out) :: solins(np)
    real(rp), intent(in) :: lat(np)
    real(rp), intent(in) :: re_factor
    real(rp), intent(in) :: sin_decl
    real(rp), intent(in) :: cos_decl
 
    real(rp) :: sin_lat
    real(rp) :: cos_lat
    real(rp) :: cos_h0
    real(rp) :: h0
 
    real(rp), parameter :: pole_eps = 100.0_rp * epsilon(1.0_rp)
 
    integer :: i
    !------------------------------------------------------------------------------
 
    !$omp parallel do private(sin_lat, cos_lat, cos_h0, h0)
    do i=1, np
      sin_lat = sin(lat(i))
      cos_lat = cos(lat(i))
      if ( abs(cos_lat) <= pole_eps ) then
        if ( sin_lat * sin_decl > 0.0_rp ) then
          h0 = pi
        else
          h0 = 0.0_rp
        end if
      else if ( abs(cos_decl) <= pole_eps ) then
        if ( sin_lat * sin_decl > 0.0_rp ) then
          h0 = pi
        else
          h0 = 0.0_rp
        end if
      else
        cos_h0 = - ( sin_lat * sin_decl ) / ( cos_lat * cos_decl )
        if ( cos_h0 >= 1.0_rp ) then ! Polar night
          h0 = 0.0_rp
        else if ( cos_h0 <= -1.0_rp ) then ! Polar day
          h0 = pi
        else
          h0 = acos(cos_h0)
        end if
      end if
 
      solins(i) = atmos_solarins_constant * re_factor / pi             &
            * ( h0 * sin_lat * sin_decl + sin(h0) * cos_lat * cos_decl )
      
      solins(i) = max( solins(i), 0.0_rp )
    end do
    return
  end subroutine daily_mean_insol
 
  !> Number of days in a month
  pure function days_in_month( YEAR, MONTH ) result(NDAYS)
    implicit none
    integer, intent(in) :: year
    integer, intent(in) :: month
 
    integer :: ndays
    integer, parameter :: ndays_normal(12) = [ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 ]
    !---------------------------------------------------------------------------
 
    ndays = ndays_normal(month)
    if ( month == 2 .and. is_leap_year(year) ) ndays = 29
    return
  end function days_in_month
 
  !> Gregorian leap-year test
!OCL SERIAL
  pure function is_leap_year( YEAR ) result(IS_LEAP)
    implicit none
    integer, intent(in) :: year
    logical :: is_leap
    !---------------------------------------------------------------------------
 
    is_leap = mod(year,4) == 0
 
    if ( mod(year,100) == 0 ) is_leap = .false.
    if ( mod(year,400) == 0 ) is_leap = .true.
    return
  end function is_leap_year  
 
  !> Convert seconds from start of day to hour/minute/second
!OCL SERIAL
  pure subroutine seconds_to_hms( DAY_SECOND, &
    HOUR, MINUTE, SECOND                      )
    implicit none
    real(rp), intent(in)  :: day_second
    integer, intent(out) :: hour
    integer, intent(out) :: minute
    integer, intent(out) :: second
 
    integer :: total_second
    !---------------------------------------------------------------------------
 
    total_second = int(day_second)
    total_second = max(0, min(total_second, 86399))
 
    hour = total_second / 3600
    total_second = total_second - hour * 3600
 
    minute = total_second / 60
    second = total_second - minute * 60
    return
  end subroutine seconds_to_hms  
 
end module scale_atm_phy_rd_solarins_simple