00001
00002
00003
00004
00005 !> @file
00006 !!
00007 !! Determination of the orbital parameters
00008 !! (eccentricity, obliquity, climate parameter CP, anomaly of vernal equinox,
00009 !! mean-annual north- or south-polar insolation).
00010 !!
00011 !! @section Copyright
00012 !!
00013 !! Copyright 2010, 2011 Ralf Greve, Bjoern Grieger, Oliver J. Stenzel
00014 !!
00015 !! @section License
00016 !!
00017 !! This file is part of MAIC-2.
00018 !!
00019 !! MAIC-2 is free software: you can redistribute it and/or modify
00020 !! it under the terms of the GNU General Public License as published by
00021 !! the Free Software Foundation, either version 3 of the License, or
00022 !! (at your option) any later version.
00023 !!
00024 !! MAIC-2 is distributed in the hope that it will be useful,
00025 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
00026 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00027 !! GNU General Public License for more details.
00028 !!
00029 !! You should have received a copy of the GNU General Public License
00030 !! along with MAIC-2. If not, see <http://www.gnu.org/licenses/>.
00031 !<
00032 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00033
00034
00035 !> Determination of the orbital parameters
00036 !! (eccentricity, obliquity, climate parameter CP, anomaly of vernal equinox,
00037 !! mean-annual north- or south-polar insolation).
00038 !<------------------------------------------------------------------------------
00039 subroutine get_orb_par(time, ecc, obl, cp, ave, insol_ma_90NS)
00040
00041 use maic2_types
00042 use maic2_variables
00043
00044 implicit none
00045
00046 real(dp), intent(in) :: time
00047 real(dp), intent(out) :: ecc, obl, cp, ave, insol_ma_90NS
00048
00049 integer(i4b) :: ndata_insol
00050 integer(i4b) :: i_kl, i_gr
00051 real(dp) :: time_kl, time_gr, ave_data_i_kl, ave_data_i_gr
00052
00053 ndata_insol = (insol_time_max-insol_time_min)/insol_time_stp
00054
00055 if (time/YEAR_SEC < real(insol_time_min,dp)) then
00056
00057 ecc = ecc_data(0)
00058 obl = obl_data(0)
00059 cp = cp_data(0)
00060 ave = ave_data(0)
00061 insol_ma_90NS = insol_ma_90(0)
00062
00063 else if (time/YEAR_SEC < real(insol_time_max,dp)) then
00064
00065 i_kl = floor( ((time/YEAR_SEC)-real(insol_time_min,dp))
00066 /real(insol_time_stp,dp) )
00067 i_kl = max(i_kl, 0)
00068
00069 i_gr = ceiling( ((time/YEAR_SEC)-real(insol_time_min,dp))
00070 /real(insol_time_stp,dp) )
00071 i_gr = min(i_gr, ndata_insol)
00072
00073 if (i_kl == i_gr) then
00074
00075 ecc = ecc_data(i_kl)
00076 obl = obl_data(i_kl)
00077 cp = cp_data(i_kl)
00078 ave = ave_data(i_kl)
00079 insol_ma_90NS = insol_ma_90(i_kl)
00080
00081 else
00082
00083 time_kl = (insol_time_min + i_kl*insol_time_stp) *YEAR_SEC
00084 time_gr = (insol_time_min + i_gr*insol_time_stp) *YEAR_SEC
00085
00086 ecc = ecc_data(i_kl) &
00087 +(ecc_data(i_gr)-ecc_data(i_kl)) &
00088 *(time-time_kl)/(time_gr-time_kl)
00089 obl = obl_data(i_kl) &
00090 +(obl_data(i_gr)-obl_data(i_kl)) &
00091 *(time-time_kl)/(time_gr-time_kl)
00092 insol_ma_90NS = insol_ma_90(i_kl) &
00093 +(insol_ma_90(i_gr)-insol_ma_90(i_kl)) &
00094 *(time-time_kl)/(time_gr-time_kl)
00095
00096 if ( abs(ave_data(i_gr)-ave_data(i_kl)) < pi ) then
00097 ave_data_i_kl = ave_data(i_kl)
00098 ave_data_i_gr = ave_data(i_gr)
00099 else
00100 if ( ave_data(i_gr) > ave_data(i_kl) ) then
00101 ave_data_i_kl = ave_data(i_kl) + 2.0_dp*pi
00102 ave_data_i_gr = ave_data(i_gr)
00103 else
00104 ave_data_i_kl = ave_data(i_kl)
00105 ave_data_i_gr = ave_data(i_gr) + 2.0_dp*pi
00106 end if
00107 end if
00108
00109 ave = ave_data_i_kl &
00110 +(ave_data_i_gr-ave_data_i_kl) &
00111 *(time-time_kl)/(time_gr-time_kl)
00112
00113 cp = cp_data(i_kl) &
00114 +(cp_data(i_gr)-cp_data(i_kl)) &
00115 *(time-time_kl)/(time_gr-time_kl)
00116
00117 end if
00118
00119 else
00120
00121 ecc = ecc_data(ndata_insol)
00122 obl = obl_data(ndata_insol)
00123 cp = cp_data(ndata_insol)
00124 ave = ave_data(ndata_insol)
00125 insol_ma_90NS = insol_ma_90(ndata_insol)
00126
00127 end if
00128
00129 end subroutine get_orb_par
00130