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