SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
instemp.f90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : i n s t e m p
4 !
5 !> @file
6 !!
7 !! Computation of the daily mean surface temperature of Mars based on
8 !! obliquity, eccentricity and the anomaly of vernal equinox (local
9 !! insolation temperature scheme = LIS scheme).
10 !!
11 !! @section Copyright
12 !!
13 !! Copyright 2009-2013 Ralf Greve, Bjoern Grieger, Oliver J. Stenzel
14 !!
15 !! @section License
16 !!
17 !! This file is part of SICOPOLIS.
18 !!
19 !! SICOPOLIS 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 !! SICOPOLIS 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 SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
31 !<
32 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33 
34 !-------------------------------------------------------------------------------
35 !> Computation of the daily mean surface temperature of Mars based on
36 !! obliquity, eccentricity and the anomaly of vernal equinox (local
37 !! insolation temperature scheme = LIS scheme).
38 !<------------------------------------------------------------------------------
39 module instemp
40 
41  use sico_types
42 
43  implicit none
44 
45  !> Martian surface temperatures
46  type ins
47  !> Temperature as a function of time and latitude
48  real(dp) :: t(0:360,-90:90)
49  !> Annual mean temperature as a function of latitude
50  real(dp) :: tam(-90:90)
51  !> Annual maximum temperature as a function of latitude
52  real(dp) :: tmax(-90:90)
53  end type ins
54 
55 contains
56 
57 !-------------------------------------------------------------------------------
58 !> Main subroutine of module instemp
59 !<------------------------------------------------------------------------------
60  subroutine setinstemp ( o, ecc, ave, obl, sma, sa, sac, op, ct )
61 
62  implicit none
63 
64  type(ins) :: o
65  real(dp), optional :: ecc ! eccentricity
66  real(dp), optional :: ave ! anom. of vernal equinox in degree
67  real(dp), optional :: obl ! obliquity in degree
68  real(dp), optional :: sma ! semi-major axis in AU
69  real(dp), optional :: sa ! surface albedo (planetary mean)
70  real(dp), optional :: sac ! surface albedo (seasonal CO2 ice caps)
71  real(dp), optional :: op ! orbital period in seconds
72  real(dp), optional :: ct ! CO2 condensation temperature in K
73 
74  logical :: co2layer(-90:90)
75  integer(i4b) :: iphi, ipsi, year
76  integer(i4b), parameter :: yearmax = 2
77  real(dp), parameter :: pi = 3.141592653589793_dp, pi_180 = pi/180.0_dp
78  real(dp), parameter :: sb = 5.67051e-8_dp
79  ! Wikipedia gives 5.670400E-8 +- 0.000040E-8
80  real(dp) :: tptd, seps, sdelta, cdelta
81  real(dp) :: albact, albact_co2, delta, du, e, eps, f, &
82  lambda, phi, psi, psi0, r
83  real(dp) :: tau0, teq, ufac, usum, w, wg
84  real(dp), dimension(-90:90) :: co2, t
85 
86  real(dp) :: j0
87  real(dp) :: a
88  real(dp) :: alb, alb_co2
89  real(dp) :: u
90  real(dp) :: tco2
91 
92  j0 = 1367.6_dp ! solar constant for Earth in W/m**2
93 
94  if ( present(ecc) ) then
95  e = ecc
96  else
97  e = 0.0935_dp
98  end if
99  if ( present(ave) ) then
100  psi0 = ave
101  else
102  psi0 = 109.13_dp
103  end if
104  if ( present(obl) ) then
105  eps = obl
106  else
107  eps = 25.19_dp
108  end if
109  if ( present(sma) ) then
110  a = sma
111  else
112  a = 1.524_dp
113  end if
114  if ( present(sa) ) then
115  alb = sa
116  else
117  alb = 0.3_dp
118  end if
119  if ( present(sac) ) then
120  alb_co2 = sac
121  else
122  alb_co2 = alb
123  end if
124  if ( present(op) ) then
125  u = op
126  else
127  u = 686.95_dp*24._dp*3600._dp
128  end if
129  if ( present(ct) ) then
130  tco2 = ct
131  else
132  tco2 = 146._dp
133  end if
134 
135  j0 = j0 / a**2
136  f = a * e
137  psi0 = psi0 * pi_180
138  eps = eps * pi_180
139 
140  seps = sin(eps)
141 
142  usum = 0.0_dp
143  do ipsi = 0, 360, 1
144  psi = ipsi * pi_180
145  r = ( a**2 - f**2 ) / ( a + f*cos(psi) )
146  du = r**2
147  if (ipsi<360) usum = usum + du
148  end do
149  ufac = u / usum
150 
151  t = 273.15_dp
152  co2 = 0.0_dp
153  co2layer = .false.
154  do year = 1, yearmax
155  usum = 0.0_dp
156  do ipsi = 0, 360, 1
157  psi = ipsi * pi_180
158  r = ( a**2 - f**2 ) / ( a + f*cos(psi) )
159  du = ufac * r**2
160  if (ipsi<360) usum = usum + du
161  lambda = psi - psi0
162  delta = asin( seps * sin(lambda) )
163  sdelta = sin(delta)
164  cdelta = cos(delta)
165  do iphi = -89, 89, 1
166  phi = iphi * pi_180
167  tptd = -tan(phi) * tan(delta)
168  if ( tptd .le. -1.0_dp ) then
169  tau0 = pi
170  else if ( tptd .ge. 1.0_dp ) then
171  tau0 = 0.0_dp
172  else
173  tau0 = acos( tptd )
174  end if
175  w = ( tau0*sin(phi)*sdelta + &
176  cos(phi)*cdelta*sin(tau0) ) &
177  * j0 * a**2 / ( pi * r**2 )
178  wg = 0.0_dp ! 0.045 * j0
179  if ( t(iphi) < -999._dp ) then
180  albact = 0.95_dp
181  albact_co2 = 0.95_dp
182  else
183  albact = alb
184  albact_co2 = alb_co2
185  end if
186  teq = ( ( wg + (1._dp-albact) * w ) / sb )**.25_dp
187  if ( teq < tco2 ) then
188  co2layer(iphi) = .true.
189  co2(iphi) = co2(iphi) + sb * tco2**4 * du &
190  - (1._dp-albact_co2) * w * du
191  t(iphi) = tco2
192  else
193  if ( co2layer(iphi) ) then
194  co2(iphi) = co2(iphi) + sb * tco2**4 * du &
195  - (1._dp-albact_co2) * w * du
196  if ( co2(iphi) .le. 0.0_dp ) then
197  co2(iphi) = 0.0_dp
198  co2layer(iphi) = .false.
199  t(iphi) = tco2
200  end if
201  else
202  t(iphi) = teq
203  end if
204  end if
205  if ( year .eq. yearmax ) o%t(ipsi,iphi) = t(iphi)
206  end do
207  if ( year .eq. yearmax .and. ipsi .eq. 0 ) then
208  o%tam = 0.0_dp
209  o%tmax = 0.0_dp
210  else
211  o%tam = o%tam + t
212  do iphi = -89, 89
213  o%tmax(iphi) = max( o%tmax(iphi), t(iphi) )
214  end do
215  end if
216  end do
217  o%tam = o%tam / 360._dp
218  end do
219  o%t(:,-90) = o%t(:,-89) + ( o%t(:,-89) - o%t(:,-88) ) / 2._dp
220  o%t(:, 90) = o%t(:, 89) + ( o%t(:, 89) - o%t(:, 88) ) / 2._dp
221  o%tam(-90) = o%tam(-89) + ( o%tam(-89) - o%tam(-88) ) / 2._dp
222  o%tam( 90) = o%tam( 89) + ( o%tam( 89) - o%tam( 88) ) / 2._dp
223  o%tmax(-90) = o%tmax(-89) + ( o%tmax(-89) - o%tmax(-88) ) / 2._dp
224  o%tmax( 90) = o%tmax( 89) + ( o%tmax( 89) - o%tmax( 88) ) / 2._dp
225 
226  end subroutine
227 
228 !-------------------------------------------------------------------------------
229 !> Annual mean temperature at latitude phi
230 !<------------------------------------------------------------------------------
231  function instam ( o, phi )
232 
233  implicit none
234  real(dp) :: instam
235  type(ins) :: o
236  real(dp) :: phi
237 
238  integer(i4b) :: iphi1, iphi2
239 
240  iphi1 = nint( phi - 0.5_dp )
241  if ( iphi1 < -90 ) then
242  iphi1 = -90
243  else if ( iphi1 > 89 ) then
244  iphi1 = 89
245  end if
246  iphi2 = iphi1 + 1
247  instam = o%tam(iphi1) + ( o%tam(iphi2) - o%tam(iphi1) ) &
248  * ( phi - iphi1 )
249 
250  end function instam
251 
252 !-------------------------------------------------------------------------------
253 !> Annual maximum temperature at latitude phi
254 !<------------------------------------------------------------------------------
255  function instmax ( o, phi )
256 
257  implicit none
258  real(dp) :: instmax
259  type(ins) :: o
260  real(dp) :: phi
261 
262  integer(i4b) :: iphi1, iphi2
263 
264  iphi1 = nint( phi - 0.5_dp )
265  if ( iphi1 < -90 ) then
266  iphi1 = -90
267  else if ( iphi1 > 89 ) then
268  iphi1 = 89
269  end if
270  iphi2 = iphi1 + 1
271  instmax = o%tmax(iphi1) + ( o%tmax(iphi2) - o%tmax(iphi1) ) &
272  * ( phi - iphi1 )
273 
274  end function instmax
275 
276 !-------------------------------------------------------------------------------
277 !> Temperature at orbit position psi and latitude phi
278 !<------------------------------------------------------------------------------
279  function inst ( o, psi, phi )
280 
281  implicit none
282  real(dp) :: inst
283  type(ins) :: o
284  real(dp) :: psi, phi
285 
286  integer(i4b) :: ipsi1, ipsi2
287 
288  ipsi1 = nint( psi - 0.5_dp )
289  if ( ipsi1 < 0 ) then
290  ipsi1 = 0
291  else if ( ipsi1 > 359 ) then
292  ipsi1 = 359
293  end if
294  ipsi2 = ipsi1 + 1
295  inst = inst1(o,ipsi1,phi) + ( inst1(o,ipsi2,phi) - inst1(o,ipsi1,phi) ) &
296  * ( psi - ipsi1 )
297 
298  end function inst
299 
300 !-------------------------------------------------------------------------------
301 !> Temperature at orbit position ipsi (integer) and latitude phi
302 !<------------------------------------------------------------------------------
303  function inst1 ( o, ipsi, phi )
304 
305  implicit none
306  real(dp) :: inst1
307  type(ins) :: o
308  integer(i4b) :: ipsi
309  real(dp) :: phi
310 
311  integer(i4b) :: iphi1, iphi2
312 
313  iphi1 = nint( phi - 0.5_dp )
314  if ( iphi1 < -90 ) then
315  iphi1 = -90
316  else if ( iphi1 > 89 ) then
317  iphi1 = 89
318  end if
319  iphi2 = iphi1 + 1
320  inst1 = o%t(ipsi,iphi1) + ( o%t(ipsi,iphi2) - o%t(ipsi,iphi1) ) &
321  * ( phi - iphi1 )
322 
323  end function inst1
324 
325 end module instemp
326 !
Generated on Thu Sep 12 2013 09:45:20 for SICOPOLIS V3.1 by   doxygen 1.8.2