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