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).
13 !! Copyright 2009-2013 Ralf Greve, Bjoern Grieger, Oliver J. Stenzel
17 !! This file is part of SICOPOLIS.
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.
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.
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/>.
32 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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 !<------------------------------------------------------------------------------
45 !> Martian surface temperatures
47 !> Temperature as a function of time and latitude
48 real(dp) :: t(0:360,-90:90)
50 real(dp) :: tam(-90:90)
52 real(dp) :: tmax(-90:90)
58 !> Main subroutine of module instemp
59 !<------------------------------------------------------------------------------
60 subroutine setinstemp ( o, ecc, ave, obl, sma, sa, sac, op, ct )
65 real(dp),
optional :: ecc
66 real(dp),
optional :: ave
67 real(dp),
optional :: obl
68 real(dp),
optional :: sma
69 real(dp),
optional :: sa
70 real(dp),
optional :: sac
71 real(dp),
optional :: op
72 real(dp),
optional :: ct
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
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
88 real(dp) :: alb, alb_co2
94 if ( present(ecc) )
then
99 if ( present(ave) )
then
104 if ( present(obl) )
then
109 if ( present(sma) )
then
114 if ( present(sa) )
then
119 if ( present(sac) )
then
124 if ( present(op) )
then
127 u = 686.95_dp*24._dp*3600._dp
129 if ( present(ct) )
then
145 r = ( a**2 - f**2 ) / ( a + f*cos(psi) )
147 if (ipsi<360) usum = usum + du
158 r = ( a**2 - f**2 ) / ( a + f*cos(psi) )
160 if (ipsi<360) usum = usum + du
162 delta = asin( seps * sin(lambda) )
167 tptd = -tan(phi) * tan(delta)
168 if ( tptd .le. -1.0_dp )
then
170 else if ( tptd .ge. 1.0_dp )
then
175 w = ( tau0*sin(phi)*sdelta + &
176 cos(phi)*cdelta*sin(tau0) ) &
177 * j0 * a**2 / ( pi * r**2 )
179 if ( t(iphi) < -999._dp )
then
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
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
198 co2layer(iphi) = .false.
205 if ( year .eq. yearmax ) o%t(ipsi,iphi) = t(iphi)
207 if ( year .eq. yearmax .and. ipsi .eq. 0 )
then
213 o%tmax(iphi) = max( o%tmax(iphi), t(iphi) )
217 o%tam = o%tam / 360._dp
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
229 !> Annual mean temperature at latitude phi
230 !<------------------------------------------------------------------------------
238 integer(i4b) :: iphi1, iphi2
240 iphi1 = nint( phi - 0.5_dp )
241 if ( iphi1 < -90 )
then
243 else if ( iphi1 > 89 )
then
247 instam = o%tam(iphi1) + ( o%tam(iphi2) - o%tam(iphi1) ) &
253 !> Annual maximum temperature at latitude phi
254 !<------------------------------------------------------------------------------
262 integer(i4b) :: iphi1, iphi2
264 iphi1 = nint( phi - 0.5_dp )
265 if ( iphi1 < -90 )
then
267 else if ( iphi1 > 89 )
then
271 instmax = o%tmax(iphi1) + ( o%tmax(iphi2) - o%tmax(iphi1) ) &
277 !> Temperature at orbit position psi and latitude phi
278 !<------------------------------------------------------------------------------
286 integer(i4b) :: ipsi1, ipsi2
288 ipsi1 = nint( psi - 0.5_dp )
289 if ( ipsi1 < 0 )
then
291 else if ( ipsi1 > 359 )
then
301 !> Temperature at orbit position ipsi (integer) and latitude phi
302 !<------------------------------------------------------------------------------
311 integer(i4b) :: iphi1, iphi2
313 iphi1 = nint( phi - 0.5_dp )
314 if ( iphi1 < -90 )
then
316 else if ( iphi1 > 89 )
then
320 inst1 = o%t(ipsi,iphi1) + ( o%t(ipsi,iphi2) - o%t(ipsi,iphi1) ) &