MAIC-2  Revision 19
 All Classes Files Functions Variables
get_orb_par.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : g e t _ o r b _ p a r
4 !
5 !> @file
6 !!
7 !! Determination of the orbital parameters
8 !! (eccentricity, obliquity, climate parameter CP, anomaly of vernal equinox,
9 !! mean-annual north- or south-polar insolation).
10 !!
11 !! @section Copyright
12 !!
13 !! Copyright 2010-2013 Ralf Greve, Bjoern Grieger, Oliver J. Stenzel
14 !!
15 !! @section License
16 !!
17 !! This file is part of MAIC-2.
18 !!
19 !! MAIC-2 is free software: you can redistribute it and/or modify
20 !! it under the terms of the GNU General Public License as published by
21 !! the Free Software Foundation, either version 3 of the License, or
22 !! (at your option) any later version.
23 !!
24 !! MAIC-2 is distributed in the hope that it will be useful,
25 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
26 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27 !! GNU General Public License for more details.
28 !!
29 !! You should have received a copy of the GNU General Public License
30 !! along with MAIC-2. If not, see <http://www.gnu.org/licenses/>.
31 !<
32 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33 
34 !-------------------------------------------------------------------------------
35 !> Determination of the orbital parameters
36 !! (eccentricity, obliquity, climate parameter CP, anomaly of vernal equinox,
37 !! mean-annual north- or south-polar insolation).
38 !<------------------------------------------------------------------------------
39 subroutine get_orb_par(time, ecc, obl, cp, ave, insol_ma_90NS)
40 
41 use maic2_types
43 
44 implicit none
45 
46 real(dp), intent(in) :: time
47 real(dp), intent(out) :: ecc, obl, cp, ave, insol_ma_90ns
48 
49 integer(i4b) :: ndata_insol
50 integer(i4b) :: i_kl, i_gr
51 real(dp) :: time_kl, time_gr, ave_data_i_kl, ave_data_i_gr
52 
53 ndata_insol = (insol_time_max-insol_time_min)/insol_time_stp
54 
55 if (time/year_sec < real(insol_time_min,dp)) then
56 
57  ecc = ecc_data(0)
58  obl = obl_data(0)
59  cp = cp_data(0)
60  ave = ave_data(0)
61  insol_ma_90ns = insol_ma_90(0)
62 
63 else if (time/year_sec < real(insol_time_max,dp)) then
64 
65  i_kl = floor( ((time/year_sec)-real(insol_time_min,dp)) &
66  /real(insol_time_stp,dp) )
67  i_kl = max(i_kl, 0)
68 
69  i_gr = ceiling( ((time/year_sec)-real(insol_time_min,dp)) &
70  /real(insol_time_stp,dp) )
71  i_gr = min(i_gr, ndata_insol)
72 
73  if (i_kl == i_gr) then
74 
75  ecc = ecc_data(i_kl)
76  obl = obl_data(i_kl)
77  cp = cp_data(i_kl)
78  ave = ave_data(i_kl)
79  insol_ma_90ns = insol_ma_90(i_kl)
80 
81  else
82 
83  time_kl = (insol_time_min + i_kl*insol_time_stp) *year_sec
84  time_gr = (insol_time_min + i_gr*insol_time_stp) *year_sec
85 
86  ecc = ecc_data(i_kl) &
87  +(ecc_data(i_gr)-ecc_data(i_kl)) &
88  *(time-time_kl)/(time_gr-time_kl)
89  obl = obl_data(i_kl) &
90  +(obl_data(i_gr)-obl_data(i_kl)) &
91  *(time-time_kl)/(time_gr-time_kl)
92  insol_ma_90ns = insol_ma_90(i_kl) &
93  +(insol_ma_90(i_gr)-insol_ma_90(i_kl)) &
94  *(time-time_kl)/(time_gr-time_kl)
95 
96  if ( abs(ave_data(i_gr)-ave_data(i_kl)) < pi ) then ! regular case
97  ave_data_i_kl = ave_data(i_kl)
98  ave_data_i_gr = ave_data(i_gr)
99  else
100  if ( ave_data(i_gr) > ave_data(i_kl) ) then
101  ave_data_i_kl = ave_data(i_kl) + 2.0_dp*pi
102  ave_data_i_gr = ave_data(i_gr)
103  else
104  ave_data_i_kl = ave_data(i_kl)
105  ave_data_i_gr = ave_data(i_gr) + 2.0_dp*pi
106  end if
107  end if
108 
109  ave = ave_data_i_kl &
110  +(ave_data_i_gr-ave_data_i_kl) &
111  *(time-time_kl)/(time_gr-time_kl)
112 
113  cp = cp_data(i_kl) &
114  +(cp_data(i_gr)-cp_data(i_kl)) &
115  *(time-time_kl)/(time_gr-time_kl)
116  ! linear interpolation of the data
117  end if
118 
119 else
120 
121  ecc = ecc_data(ndata_insol)
122  obl = obl_data(ndata_insol)
123  cp = cp_data(ndata_insol)
124  ave = ave_data(ndata_insol)
125  insol_ma_90ns = insol_ma_90(ndata_insol)
126 
127 end if
128 
129 end subroutine get_orb_par
130 !