31 public :: atmos_solarins_insolation
33 interface atmos_solarins_insolation
36 end interface atmos_solarins_insolation
65 real(
rp),
private :: obliquity
66 real(
rp),
private :: e
67 real(
rp),
private :: omega
68 real(
rp),
private :: lambda_m0
70 integer,
private,
parameter :: year_ref = 1950
71 integer,
private :: ve_date(6)
72 data ve_date / 1950, 3, 21, 0, 0, 0 /
74 real(
rp),
private,
parameter :: obliquity_ref = 23.320556_rp
75 real(
rp),
private,
parameter :: psi_bar = 50.439273_rp
76 real(
rp),
private,
parameter :: zeta = 3.392506_rp
79 integer,
private,
parameter :: nobliq = 47
80 real(
rp),
private :: obliq_amp (nobliq)
81 real(
rp),
private :: obliq_rate (nobliq)
82 real(
rp),
private :: obliq_phase(nobliq)
231 integer,
private,
parameter :: neclip = 19
232 real(
rp),
private :: eclip_amp (neclip)
233 real(
rp),
private :: eclip_rate (neclip)
234 real(
rp),
private :: eclip_phase(neclip)
299 integer,
private,
parameter :: nprece = 78
300 real(
rp),
private :: prece_amp (nprece)
301 real(
rp),
private :: prece_rate (nprece)
302 real(
rp),
private :: prece_phase(nprece)
544 real(
rp),
private :: arcsec2d, arcsec2r
546 logical,
private :: debug = .false.
554 basepoint_lon, basepoint_lat, &
561 integer,
intent(in) :: iyear
562 real(
rp),
intent(in) :: basepoint_lon
563 real(
rp),
intent(in) :: basepoint_lat
565 namelist / param_atmos_solarins / &
585 log_info(
"ATMOS_SOLARINS_setup",*)
'Setup'
593 read(
io_fid_conf,nml=param_atmos_solarins,iostat=ierr)
595 log_info(
"ATMOS_SOLARINS_setup",*)
'Not found namelist. Default used.'
596 elseif( ierr > 0 )
then
597 log_error(
"ATMOS_SOLARINS_setup",*)
'Not appropriate names in namelist PARAM_ATMOS_SOLARINS. Check!'
600 log_nml(param_atmos_solarins)
624 arcsec2d = 1.0_rp / (60.0_rp*60.0_rp)
635 dyear = real( year-year_ref, kind=
rp )
639 log_info(
"ATMOS_SOLARINS_setup",
'(1x,A)')
'Solar insolation parameters '
640 log_info(
"ATMOS_SOLARINS_setup",
'(1x,A,I7)')
'Reference year : ', year_ref
641 log_info(
"ATMOS_SOLARINS_setup",
'(1x,A,I7)')
'Current year : ', year
642 log_info(
"ATMOS_SOLARINS_setup",
'(1x,A,I7)')
'Difference from ref. : ', int(dyear)
645 log_info(
"ATMOS_SOLARINS_setup",
'(1x,A,F12.7)')
'Obliquity [deg] : ', obliquity /
const_d2r
646 log_info(
"ATMOS_SOLARINS_setup",
'(1x,A,F12.7)')
'Eccentricity : ', e
647 log_info(
"ATMOS_SOLARINS_setup",
'(1x,A,F12.7)')
'Longitude of perihelion [deg] : ', omega /
const_d2r
648 log_info(
"ATMOS_SOLARINS_setup",
'(1x,A,F12.7)')
'Longitude at the vernal equinox [deg] : ', lambda_m0 /
const_d2r
677 integer,
intent(in) :: iyear
683 real(
rp) :: perih_psi
684 real(
rp) :: perih_omega
686 real(
rp) :: temp, beta, ecosomg_mod
691 dyear = real( iyear-year_ref, kind=
rp )
696 temp = temp + obliq_amp(i)*arcsec2d * cos( obliq_rate(i)*arcsec2r*dyear + obliq_phase(i)*d2r )
698 obliquity = ( obliquity_ref + temp ) * d2r
704 esinomg = esinomg + eclip_amp(i) * sin( eclip_rate(i)*arcsec2r*dyear + eclip_phase(i)*d2r )
705 ecosomg = ecosomg + eclip_amp(i) * cos( eclip_rate(i)*arcsec2r*dyear + eclip_phase(i)*d2r )
707 e = sqrt( esinomg*esinomg + ecosomg*ecosomg )
715 if ( ecosomg == 0.0_rp )
then
716 ecosomg_mod = 1.e-30_rp
718 ecosomg_mod = ecosomg
720 perih_pi = atan(esinomg/ecosomg_mod) + pi * ( 0.5_rp - sign(0.5_rp, ecosomg) )
744 perih_pi = perih_pi / d2r
749 temp = temp + prece_amp(i)*arcsec2d * sin( prece_rate(i)*arcsec2r*dyear + prece_phase(i)*d2r )
751 perih_psi = psi_bar * arcsec2d * dyear + zeta + temp
754 perih_omega = perih_pi + perih_psi
757 if ( perih_omega + 180.0_rp < 0.0_rp )
then
758 perih_omega = perih_omega + 360.0_rp
759 elseif( perih_omega + 180.0_rp >= 360.0_rp )
then
760 perih_omega = perih_omega - 360.0_rp
767 omega = ( perih_omega + 180.0_rp ) * d2r
769 beta = sqrt( 1.0_rp - e*e )
772 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) &
773 - ( 1.0_rp/4.0_rp*e*e ) * ( 1.0_rp/2.0_rp + beta ) * sin(2.0_rp*omega) &
774 + ( 1.0_rp/8.0_rp*e*e*e ) * ( 1.0_rp/3.0_rp + beta ) * sin(3.0_rp*omega) )
798 real(RP),
intent(out) :: Re_factor
799 real(RP),
intent(out) :: sinDEC
800 real(RP),
intent(out) :: cosDEC
801 real(RP),
intent(out) :: hourangle
802 integer,
intent(in) :: now_date(6)
803 integer,
intent(in) :: offset_year
807 integer :: absday, absday_ve
808 real(DP) :: DayOfYear, abssec
815 date(:) = now_date(:)
848 lambda_m = lambda_m0 + 2.0_rp * pi * real(absday-absday_ve,kind=rp) / dayofyear
850 nu = lambda_m - omega
854 + 2.0_rp*e * cos( nu ) &
855 + 0.5_rp*e*e * cos( 2.0_rp*nu ) &
857 + 4.0_rp*e*e*e * cos( 3.0_rp*nu )
861 + ( 2.0_rp*e - 1.0_rp/ 4.0_rp*e*e*e ) * sin( nu ) &
862 + ( 5.0_rp/ 4.0_rp*e*e ) * sin( 2.0_rp*nu ) &
863 + ( 13.0_rp/12.0_rp*e*e*e ) * sin( 3.0_rp*nu )
866 sindec = sin(lambda) * sin(obliquity)
867 cosdec = sqrt( 1.0_rp - sindec*sindec )
870 hourangle = 2.0_rp * pi * abssec / (24.0_rp*60.0_rp*60.0_rp)
873 log_info(
"ATMOS_SOLARINS_ecliptic_longitude",*) lambda_m, nu, lambda
882 real_lon, real_lat, &
883 now_date, offset_year, &
888 real(RP),
intent(in) :: real_lon
889 real(RP),
intent(in) :: real_lat
890 integer,
intent(in) :: now_date(6)
891 integer,
intent(in) :: offset_year
893 real(RP),
intent(out) :: solins
894 real(RP),
intent(out) :: cosSZA
896 real(RP) :: Re_factor
897 real(RP) :: sinDEC, cosDEC
898 real(RP) :: hourangle
919 cossza = sin(lat)*sindec - cos(lat)*cosdec*cos(lon+hourangle)
928 IA, IS, IE, JA, JS, JE, &
929 real_lon, real_lat, &
930 now_date, offset_year, &
935 integer,
intent(in) :: IA, IS, IE
936 integer,
intent(in) :: JA, JS, JE
938 real(RP),
intent(in) :: real_lon(IA,JA)
939 real(RP),
intent(in) :: real_lat(IA,JA)
940 integer,
intent(in) :: now_date(6)
941 integer,
intent(in) :: offset_year
943 real(RP),
intent(out) :: solins (IA,JA)
944 real(RP),
intent(out) :: cosSZA (IA,JA)
946 real(RP) :: Re_factor
947 real(RP) :: sinDEC, cosDEC
948 real(RP) :: hourangle
950 real(RP) :: lon(IA,JA)
951 real(RP) :: lat(IA,JA)
975 lon(i,j) = real_lon(i,j)
976 lat(i,j) = real_lat(i,j)
984 cossza(i,j) = sin(lat(i,j))*sindec - cos(lat(i,j))*cosdec*cos(lon(i,j)+hourangle)