module SOLARINS More...

Functions/Subroutines
subroutine, public	atmos_solarins_setup (iyear)
	setup solar incidence module More...

subroutine, public	atmos_solarins_orbit (iyear)
	setup solar incidence module More...

subroutine	atmos_solarins_insolation_0d (solins, cosSZA, Re_factor, real_lon, real_lat, now_date, offset_year)
	calc factor of Earths solar insolation More...

subroutine	atmos_solarins_insolation_2d (solins, cosSZA, real_lon, real_lat, now_date, offset_year)
	calc factor of Earths solar insolation More...

Variables
real(rp), public	atmos_solarins_constant = 1360.250117_RP

logical, public	atmos_solarins_fixedlatlon = .false.

logical, public	atmos_solarins_fixeddate = .false.

real(rp), public	atmos_solarins_lon

real(rp), public	atmos_solarins_lat

integer, dimension(6), public	atmos_solarins_date

Detailed Description

module SOLARINS

Description: calculate solar insolation based on berger(1978), berger et al.(1993) The original algorithm is valid to +/- 1,000,000 years from AD1950

Author: Team SCALE

History

2013-01-29 (H.Yashiro) [new]

Function/Subroutine Documentation

◆ atmos_solarins_setup()

subroutine, public scale_atmos_solarins::atmos_solarins_setup ( integer, intent(in) iyear )

setup solar incidence module

Definition at line 547 of file scale_atmos_sub_solarins.F90.

References atmos_solarins_constant, atmos_solarins_date, atmos_solarins_fixeddate, atmos_solarins_fixedlatlon, atmos_solarins_lat, atmos_solarins_lon, atmos_solarins_orbit(), scale_const::const_d2r, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_l, scale_stdio::io_lnml, scale_process::prc_mpistop(), scale_grid_real::real_basepoint_lat, and scale_grid_real::real_basepoint_lon.

Referenced by mod_atmos_driver::atmos_driver_setup().

     use scale_process, only: &
        prc_mpistop
     use scale_const, only: &
        const_d2r
     use scale_grid_real, only: &
        real_basepoint_lon, &
        real_basepoint_lat
     implicit none
 
     integer, intent(in) :: iyear ! year at setup
 
     namelist / param_atmos_solarins / &
        atmos_solarins_constant,    &
        atmos_solarins_fixedlatlon, &
        atmos_solarins_fixeddate,   &
        atmos_solarins_lon,         &
        atmos_solarins_lat,         &
        atmos_solarins_date
 
     real(RP) :: dyear ! delta t [year]
     integer  :: year  ! used year
 
     integer  :: ierr
     !---------------------------------------------------------------------------
 
     if( io_l ) write(io_fid_log,*)
     if( io_l ) write(io_fid_log,*) '++++++ Module[SOLARINS] / Categ[ATMOS SHARE] / Origin[SCALElib]'
 
     atmos_solarins_lon     = real_basepoint_lon
     atmos_solarins_lat     = real_basepoint_lat
     atmos_solarins_date(:) = -1
 
     !--- read namelist
     rewind(io_fid_conf)
     read(io_fid_conf,nml=param_atmos_solarins,iostat=ierr)
     if( ierr < 0 ) then !--- missing
        if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
     elseif( ierr > 0 ) then !--- fatal error
        write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_SOLARINS. Check!'
        call prc_mpistop
     endif
     if( io_lnml ) write(io_fid_log,nml=param_atmos_solarins)
 
     arcsec2d = 1.0_rp / (60.0_rp*60.0_rp)
     arcsec2r = arcsec2d * const_d2r
 
     year = iyear
 
     if ( atmos_solarins_fixeddate ) then
        if( atmos_solarins_date(1) >= 0 ) year = atmos_solarins_date(1)
     endif
 
     call atmos_solarins_orbit( year )
 
     dyear = real( year-year_ref, kind=rp )
 
     !----- report data -----
     if( io_l ) write(io_fid_log,*)
     if( io_l ) write(io_fid_log,'(1x,A)')       '*** insolation parameters'
     if( io_l ) write(io_fid_log,'(1x,A,I7)')    '*** reference year      : ', year_ref
     if( io_l ) write(io_fid_log,'(1x,A,I7)')    '*** current   year      : ', iyear
     if( io_l ) write(io_fid_log,'(1x,A,I7)')    '*** difference from ref.: ', int(dyear)
     if( io_l ) write(io_fid_log,*)
     if( io_l ) write(io_fid_log,'(1x,A,F12.7)') '*** solar constant                 [W/m2]: ', atmos_solarins_constant
     if( io_l ) write(io_fid_log,'(1x,A,F12.7)') '*** obliquity                       [deg]: ', obliquity / const_d2r
     if( io_l ) write(io_fid_log,'(1x,A,F12.7)') '*** eccentricity                         : ', e
     if( io_l ) write(io_fid_log,'(1x,A,F12.7)') '*** longitude of perihelion         [deg]: ', omega     / const_d2r
     if( io_l ) write(io_fid_log,'(1x,A,F12.7)') '*** longitude at the vernal equinox [deg]: ', lambda_m0 / const_d2r
     if( io_l ) write(io_fid_log,*)
 
     if( io_l ) write(io_fid_log,*)              '*** Latitude/Longitude is fixed?         : ', atmos_solarins_fixedlatlon
     if ( atmos_solarins_fixedlatlon ) then
        if( io_l ) write(io_fid_log,*)           '*** Longitude                       [deg]: ', atmos_solarins_lon
        if( io_l ) write(io_fid_log,*)           '*** Latitude                        [deg]: ', atmos_solarins_lat
        atmos_solarins_lon = atmos_solarins_lon * const_d2r
        atmos_solarins_lat = atmos_solarins_lat * const_d2r
     endif
 
     if( io_l ) write(io_fid_log,*)              '*** Date is fixed?                       : ', atmos_solarins_fixeddate
     if ( atmos_solarins_fixeddate ) then
        if( io_l ) write(io_fid_log,*)           '*** Date                                 : ', atmos_solarins_date
     endif
 
     return

Here is the call graph for this function:

Here is the caller graph for this function:

◆ atmos_solarins_orbit()

subroutine, public scale_atmos_solarins::atmos_solarins_orbit ( integer, intent(in) iyear )

setup solar incidence module

Definition at line 638 of file scale_atmos_sub_solarins.F90.

References scale_const::const_d2r, and scale_const::const_pi.

Referenced by atmos_solarins_setup().

     use scale_const, only: &
          pi => const_pi, &
          const_d2r
     implicit none
 
     integer, intent(in) :: iyear ! year at setup
 
     real(RP) :: dyear ! delta t [year]
 
     real(RP) :: esinomg     ! e * sin(w) : ecliptic parameter
     real(RP) :: ecosomg     ! e * cos(w) : ecliptic parameter
     real(RP) :: perih_pi    ! pi         : longitude of fixed perihelion [rad]
     real(RP) :: perih_psi   ! psi        : general precession [degree]
     real(RP) :: perih_omega ! omega_bar  : longitude of perihelion [degree]
 
     real(RP) :: temp, beta, ecosomg_mod
     integer  :: i
     !---------------------------------------------------------------------------
 
     ! time from reference year(1950.0 AD)
     dyear = real( iyear-year_ref, kind=rp )
 
     ! obliquity
     temp = 0.0_rp
     do i = 1, nobliq
        temp = temp + obliq_amp(i)*arcsec2d * cos( obliq_rate(i)*arcsec2r*dyear + obliq_phase(i)*const_d2r )
     enddo
     obliquity = ( obliquity_ref + temp ) * const_d2r
 
     ! eccentricity
     esinomg = 0.0_rp
     ecosomg = 0.0_rp
     do i = 1, neclip
        esinomg = esinomg + eclip_amp(i) * sin( eclip_rate(i)*arcsec2r*dyear + eclip_phase(i)*const_d2r )
        ecosomg = ecosomg + eclip_amp(i) * cos( eclip_rate(i)*arcsec2r*dyear + eclip_phase(i)*const_d2r )
     enddo
     e = sqrt( esinomg*esinomg + ecosomg*ecosomg )
 
     ! longitude of fixed perihelion
     if ( ecosomg == 0.0_rp ) then
        ecosomg_mod  = 1.e-30_rp
     else
        ecosomg_mod  = ecosomg
     endif
     perih_pi = atan(esinomg/ecosomg_mod) + pi * ( 0.5_rp - sign(0.5_rp, ecosomg) )
 
 !    if( abs(EcosOMG) < 1.E-8_RP ) then
 !       if    ( EsinOMG == 0.0_RP ) then
 !          Perih_pi = 0.0_RP
 !       elseif( EsinOMG <  0.0_RP ) then
 !          Perih_pi = 1.5_RP * PI
 !       elseif( EsinOMG >  0.0_RP ) then
 !          Perih_pi = 0.5_RP * PI
 !       endif
 !
 !    elseif( EcosOMG < 0.0_RP ) then
 !
 !       Perih_pi = atan(EsinOMG/EcosOMG) + PI
 !
 !    elseif( EcosOMG > 0.0_RP ) then
 !
 !       if ( EsinOMG <  0.0_RP ) then
 !          Perih_pi = atan(EsinOMG/EcosOMG) + 2.0_RP * PI
 !       else
 !          Perih_pi = atan(EsinOMG/EcosOMG)
 !       endif
 !
 !    endif
     perih_pi = perih_pi / const_d2r ! [rad]->[degree]
 
     ! general precession
     temp = 0.0_rp
     do i = 1, nprece
        temp = temp + prece_amp(i)*arcsec2d * sin( prece_rate(i)*arcsec2r*dyear + prece_phase(i)*const_d2r )
     enddo
     perih_psi = psi_bar * arcsec2d * dyear + zeta + temp
 
     ! longitude of perihelion
     perih_omega = perih_pi + perih_psi
 
     do i = 1, 1000
        if    ( perih_omega + 180.0_rp <    0.0_rp ) then
           perih_omega = perih_omega + 360.0_rp
        elseif( perih_omega + 180.0_rp >= 360.0_rp ) then
           perih_omega = perih_omega - 360.0_rp
        else
           exit
        endif
     enddo
 
     ! longitude of perigee (see berger et al.(1993))
     omega = ( perih_omega + 180.0_rp ) * const_d2r
 
     beta = sqrt( 1.0_rp - e*e )
 
     ! longitude at the vernal equinox
     lambda_m0 = 2.0_rp * ( ( 1.0_rp/2.0_rp*e + 1.0_rp/8.0_rp*e*e*e ) * ( 1.0_rp        + beta ) * sin(       omega) &
                          - (                   1.0_rp/4.0_rp*e*e   ) * ( 1.0_rp/2.0_rp + beta ) * sin(2.0_rp*omega) &
                          + (                   1.0_rp/8.0_rp*e*e*e ) * ( 1.0_rp/3.0_rp + beta ) * sin(3.0_rp*omega) )
 
     return

Here is the caller graph for this function:

◆ atmos_solarins_insolation_0d()

subroutine scale_atmos_solarins::atmos_solarins_insolation_0d	(	real(rp), intent(out)	solins,
		real(rp), intent(out)	cosSZA,
		real(rp), intent(out)	Re_factor,
		real(rp), intent(in)	real_lon,
		real(rp), intent(in)	real_lat,
		integer, dimension(6), intent(in)	now_date,
		integer, intent(in)	offset_year
	)

calc factor of Earths solar insolation

Definition at line 751 of file scale_atmos_sub_solarins.F90.

References atmos_solarins_constant, atmos_solarins_date, atmos_solarins_fixeddate, atmos_solarins_fixedlatlon, atmos_solarins_lat, atmos_solarins_lon, scale_calendar::calendar_getdayofyear(), scale_calendar::calendar_hms2abssec(), scale_calendar::calendar_ymd2absday(), scale_const::const_eps, scale_const::const_pi, scale_calendar::i_day, scale_calendar::i_hour, scale_calendar::i_min, scale_calendar::i_month, scale_calendar::i_sec, and scale_calendar::i_year.

     use scale_const, only: &
        pi  => const_pi,  &
        eps => const_eps
     use scale_calendar, only: &
        calendar_getdayofyear,  &
        calendar_ymd2absday,    &
        calendar_hms2abssec,    &
        i_year, i_month, i_day, &
        i_hour, i_min, i_sec
     implicit none
 
     real(RP), intent(out) :: solins      ! solar insolation
     real(RP), intent(out) :: cossza      ! cos(Solar Zenith Angle)
     real(RP), intent(out) :: re_factor   ! factor of the distance of Earth from the sun (1/rho2)
     real(RP), intent(in)  :: real_lon    ! longitude
     real(RP), intent(in)  :: real_lat    ! latitude
     integer,  intent(in)  :: now_date(6) ! date(yyyy,mm,dd,hh,mm,ss)
     integer,  intent(in)  :: offset_year ! year offset
 
     real(RP) :: lon
     real(RP) :: lat
     integer  :: date(6)
     integer  :: oyear
 
     real(RP) :: lambda_m       ! mean longitude from vernal equinox
     real(RP) :: lambda         !
     real(RP) :: sindec, cosdec ! sin/cos(solar declination)
     real(RP) :: hourangle      ! hour angle: relative longitude of subsolar point
 
     integer  :: absday, absday_ve
     real(DP) :: dayofyear, abssec
     real(RP) :: nu
     !---------------------------------------------------------------------------
 
     lon = real_lon
     lat = real_lat
     if ( atmos_solarins_fixedlatlon ) then
        lon = atmos_solarins_lon
        lat = atmos_solarins_lat
     endif
 
     date(:) = now_date(:)
     oyear   = offset_year
     if ( atmos_solarins_fixeddate ) then
        if( atmos_solarins_date(1) >= 0 ) oyear   = 0
        if( atmos_solarins_date(1) >= 0 ) date(1) = atmos_solarins_date(1)
        if( atmos_solarins_date(2) >= 1 ) date(2) = atmos_solarins_date(2)
        if( atmos_solarins_date(3) >= 1 ) date(3) = atmos_solarins_date(3)
        if( atmos_solarins_date(4) >= 0 ) date(4) = atmos_solarins_date(4)
        if( atmos_solarins_date(5) >= 0 ) date(5) = atmos_solarins_date(5)
        if( atmos_solarins_date(6) >= 0 ) date(6) = atmos_solarins_date(6)
     endif
 
     call calendar_getdayofyear( dayofyear, date(i_year) )
 
     call calendar_ymd2absday( absday,        & ! [OUT]
                               date(i_year),  & ! [IN]
                               date(i_month), & ! [IN]
                               date(i_day),   & ! [IN]
                               oyear          ) ! [IN]
 
     call calendar_ymd2absday( absday_ve,        & ! [OUT]
                               date(i_year),  & ! [IN]
                               ve_date(i_month), & ! [IN]
                               ve_date(i_day),   & ! [IN]
                               oyear             ) ! [IN]
 
     call calendar_hms2abssec( abssec,        & ! [OUT]
                               date(i_hour),  & ! [IN]
                               date(i_min),   & ! [IN]
                               date(i_sec),   & ! [IN]
                               0.0_dp         ) ! [IN]
 
     lambda_m = lambda_m0 + 2.0_rp * pi * real(absday-absday_ve,kind=RP) / dayofyear
 
     nu = lambda_m - omega
 
     ! actual longitude from vernal equinox
     lambda = lambda_m &
            + ( 2.0_rp*e -  1.0_rp/ 4.0_rp*e*e*e ) * sin(        nu ) &
            + (             5.0_rp/ 4.0_rp*e*e   ) * sin( 2.0_rp*nu ) &
            + (            13.0_rp/12.0_rp*e*e*e ) * sin( 3.0_rp*nu )
 
     ! solar declination
     sindec = sin(lambda) * sin(obliquity)
     cosdec = sqrt( 1.0_rp - sindec*sindec )
 
     ! hour angle
     hourangle = lon + 2.0_rp * pi * abssec / (24.0_rp*60.0_rp*60.0_rp)
 
     ! cos(Solar Zenith Angle)
     cossza = sin(lat)*sindec - cos(lat)*cosdec*cos(hourangle)
 
     ! 1 / (rho*rho)
     re_factor = 1.0_rp                          &
               + 2.0_rp*e     * cos(        nu ) &
               + 0.5_rp*e*e   * cos( 2.0_rp*nu ) &
               + 2.5_rp*e*e                      &
               + 4.0_rp*e*e*e * cos( 3.0_rp*nu )
 
 !    solins = ATMOS_SOLARINS_constant * Re_factor * max(cosSZA,0.0_RP)
     solins = atmos_solarins_constant * re_factor * ( 0.5_rp + sign(0.5_rp,cossza-eps) )
 
     return

Here is the call graph for this function:

◆ atmos_solarins_insolation_2d()

subroutine scale_atmos_solarins::atmos_solarins_insolation_2d	(	real(rp), dimension (ia,ja), intent(out)	solins,
		real(rp), dimension (ia,ja), intent(out)	cosSZA,
		real(rp), dimension(ia,ja), intent(in)	real_lon,
		real(rp), dimension(ia,ja), intent(in)	real_lat,
		integer, dimension(6), intent(in)	now_date,
		integer, intent(in)	offset_year
	)

calc factor of Earths solar insolation

Definition at line 866 of file scale_atmos_sub_solarins.F90.

References atmos_solarins_constant, atmos_solarins_date, atmos_solarins_fixeddate, atmos_solarins_fixedlatlon, atmos_solarins_lat, atmos_solarins_lon, scale_calendar::calendar_getdayofyear(), scale_calendar::calendar_hms2abssec(), scale_calendar::calendar_ymd2absday(), scale_const::const_eps, scale_const::const_pi, scale_calendar::i_day, scale_calendar::i_hour, scale_calendar::i_min, scale_calendar::i_month, scale_calendar::i_sec, scale_calendar::i_year, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, and scale_grid_index::js.

     use scale_grid_index
     use scale_const, only: &
        pi  => const_pi,  &
        eps => const_eps
     use scale_calendar, only: &
        calendar_getdayofyear,  &
        calendar_ymd2absday,    &
        calendar_hms2abssec,    &
        i_year, i_month, i_day, &
        i_hour, i_min, i_sec
     implicit none
 
     real(RP), intent(out) :: solins  (ia,ja) ! solar insolation
     real(RP), intent(out) :: cossza  (ia,ja) ! cos(Solar Zenith Angle)
     real(RP), intent(in)  :: real_lon(ia,ja) ! longitude [rad]
     real(RP), intent(in)  :: real_lat(ia,ja) ! latitude  [rad]
     integer,  intent(in)  :: now_date(6)     ! date(yyyy,mm,dd,hh,mm,ss)
     integer,  intent(in)  :: offset_year     ! year offset
 
     real(RP) :: lon(ia,ja)
     real(RP) :: lat(ia,ja)
     integer  :: date(6)
     integer  :: oyear
 
     real(RP) :: lambda_m       ! mean longitude from vernal equinox
     real(RP) :: lambda         !
     real(RP) :: sindec, cosdec ! sin/cos(solar declination)
     real(RP) :: hourangle(ia,ja) ! hour angle: relative longitude of subsolar point
 
     integer  :: absday, absday_ve
     real(DP) :: dayofyear, abssec
     real(RP) :: nu
     real(RP) :: re_factor   ! factor of the distance of Earth from the sun (1/rho2)
 
     integer  :: i, j
     !---------------------------------------------------------------------------
 
     lon(:,:) = real_lon(:,:)
     lat(:,:) = real_lat(:,:)
     if ( atmos_solarins_fixedlatlon ) then
        lon(:,:) = atmos_solarins_lon
        lat(:,:) = atmos_solarins_lat
     endif
 
     date(:) = now_date(:)
     oyear   = offset_year
     if ( atmos_solarins_fixeddate ) then
        if( atmos_solarins_date(1) >= 0 ) oyear   = 0
        if( atmos_solarins_date(1) >= 0 ) date(1) = atmos_solarins_date(1)
        if( atmos_solarins_date(2) >= 1 ) date(2) = atmos_solarins_date(2)
        if( atmos_solarins_date(3) >= 1 ) date(3) = atmos_solarins_date(3)
        if( atmos_solarins_date(4) >= 0 ) date(4) = atmos_solarins_date(4)
        if( atmos_solarins_date(5) >= 0 ) date(5) = atmos_solarins_date(5)
        if( atmos_solarins_date(6) >= 0 ) date(6) = atmos_solarins_date(6)
     endif
 
     call calendar_getdayofyear( dayofyear, date(i_year) )
 
     call calendar_ymd2absday( absday,        & ! [OUT]
                               date(i_year),  & ! [IN]
                               date(i_month), & ! [IN]
                               date(i_day),   & ! [IN]
                               oyear          ) ! [IN]
 
     call calendar_ymd2absday( absday_ve,        & ! [OUT]
                               date(i_year),  & ! [IN]
                               ve_date(i_month), & ! [IN]
                               ve_date(i_day),   & ! [IN]
                               oyear             ) ! [IN]
 
     call calendar_hms2abssec( abssec,       & ! [OUT]
                               date(i_hour), & ! [IN]
                               date(i_min),  & ! [IN]
                               date(i_sec),  & ! [IN]
                               0.0_dp        ) ! [IN]
 
     lambda_m = lambda_m0 + 2.0_rp * pi * real(absday-absday_ve,kind=RP) / dayofyear
 
     nu = lambda_m - omega
 
     ! 1 / (rho*rho)
     re_factor = 1.0_rp                          &
               + 2.0_rp*e     * cos(        nu ) &
               + 0.5_rp*e*e   * cos( 2.0_rp*nu ) &
               + 2.5_rp*e*e                      &
               + 4.0_rp*e*e*e * cos( 3.0_rp*nu )
 
     ! actual longitude from vernal equinox
     lambda = lambda_m &
            + ( 2.0_rp*e -  1.0_rp/ 4.0_rp*e*e*e ) * sin(        nu ) &
            + (             5.0_rp/ 4.0_rp*e*e   ) * sin( 2.0_rp*nu ) &
            + (            13.0_rp/12.0_rp*e*e*e ) * sin( 3.0_rp*nu )
 
     ! solar declination
     sindec = sin(lambda) * sin(obliquity)
     cosdec = sqrt( 1.0_rp - sindec*sindec )
 
     do j = js, je
     do i = is, ie
        ! hour angle
        hourangle(i,j) = lon(i,j) + 2.0_rp * pi * abssec / (24.0_rp*60.0_rp*60.0_rp)
     enddo
     enddo
 
     do j = js, je
     do i = is, ie
        ! cos(Solar Zenith Angle)
        cossza(i,j) = sin(lat(i,j))*sindec - cos(lat(i,j))*cosdec*cos(hourangle(i,j))
     enddo
     enddo
 
     do j = js, je
     do i = is, ie
 !       solins(i,j) = ATMOS_SOLARINS_constant * Re_factor * max(cosSZA(i,j),0.0_RP)
        solins(i,j) = atmos_solarins_constant * re_factor * ( 0.5_rp + sign(0.5_rp,cossza(i,j)-eps) )
     enddo
     enddo
 
     return