00001
00002
00003
00004
00005 !> @file
00006 !!
00007 !! Computation of the daily mean surface temperature of Mars
00008 !! based on obliquity, eccentricity and the anomaly of vernal
00009 !! equinox (local insolation temperature = LIS scheme).
00010 !!
00011 !! @section Copyright
00012 !!
00013 !! Copyright 2010, 2011 Ralf Greve, Bjoern Grieger, Oliver J. Stenzel
00014 !!
00015 !! @section License
00016 !!
00017 !! This file is part of MAIC-2.
00018 !!
00019 !! MAIC-2 is free software: you can redistribute it and/or modify
00020 !! it under the terms of the GNU General Public License as published by
00021 !! the Free Software Foundation, either version 3 of the License, or
00022 !! (at your option) any later version.
00023 !!
00024 !! MAIC-2 is distributed in the hope that it will be useful,
00025 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
00026 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00027 !! GNU General Public License for more details.
00028 !!
00029 !! You should have received a copy of the GNU General Public License
00030 !! along with MAIC-2. If not, see <http://www.gnu.org/licenses/>.
00031 !<
00032 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00033
00034
00035 !> Computation of the daily mean surface temperature of Mars
00036 !! based on obliquity, eccentricity and the anomaly of vernal
00037 !! equinox (local insolation temperature = LIS scheme).
00038 !<------------------------------------------------------------------------------
00039 module instemp
00040
00041 use maic2_types
00042
00043 implicit none
00044
00045 !> Surface temperatures
00046 type ins
00047 real(dp) :: t(0:360,-90:90)
00048 real(dp) :: tam(-90:90)
00049 real(dp) :: tmax(-90:90)
00050 end type ins
00051
00052 contains
00053
00054
00055 !> Main subroutine of module instemp.
00056 !<------------------------------------------------------------------------------
00057 subroutine setinstemp ( o, ecc, ave, obl, sma, sa, sac, op, ct )
00058
00059 implicit none
00060
00061 type(ins) :: o
00062 real(dp), optional :: ecc
00063 real(dp), optional :: ave
00064 real(dp), optional :: obl
00065 real(dp), optional :: sma
00066 real(dp), optional :: sa
00067 real(dp), optional :: sac
00068 real(dp), optional :: op
00069 real(dp), optional :: ct
00070
00071 logical :: co2layer(-90:90)
00072 integer(i4b) :: iphi, ipsi, year
00073 integer(i4b), parameter :: yearmax = 2
00074 real(dp), parameter :: pi = 3.141592653589793_dp, pi_180 = pi/180.0_dp
00075 real(dp), parameter :: SB = 5.67051e-8_dp
00076
00077 real(dp) :: tptd, seps, sdelta, cdelta
00078 real(dp) :: albact, albact_co2, delta, du, e, eps, f,
00079 lambda, phi, psi, psi0, r
00080 real(dp) :: tau0, teq, ufac, usum, w, wg
00081 real(dp), dimension(-90:90) :: co2, t
00082
00083 real(dp) :: j0
00084 real(dp) :: a
00085 real(dp) :: alb, alb_co2
00086 real(dp) :: u
00087 real(dp) :: tco2
00088
00089 j0 = 1367.6_dp
00090
00091 if ( present(ecc) ) then
00092 e = ecc
00093 else
00094 e = 0.0935_dp
00095 end if
00096 if ( present(ave) ) then
00097 psi0 = ave
00098 else
00099 psi0 = 109.13_dp
00100 end if
00101 if ( present(obl) ) then
00102 eps = obl
00103 else
00104 eps = 25.19_dp
00105 end if
00106 if ( present(sma) ) then
00107 a = sma
00108 else
00109 a = 1.524_dp
00110 end if
00111 if ( present(sa) ) then
00112 alb = sa
00113 else
00114 alb = 0.3_dp
00115 end if
00116 if ( present(sac) ) then
00117 alb_co2 = sac
00118 else
00119 alb_co2 = 0.7_dp
00120 end if
00121 if ( present(op) ) then
00122 u = op
00123 else
00124 u = 686.95_dp*24._dp*3600._dp
00125 end if
00126 if ( present(ct) ) then
00127 tco2 = ct
00128 else
00129 tco2 = 146._dp
00130 end if
00131
00132 j0 = j0 / a**2
00133 f = a * e
00134 psi0 = psi0 * pi_180
00135 eps = eps * pi_180
00136
00137 seps = sin(eps)
00138
00139 usum = 0.0_dp
00140 do ipsi = 0, 360, 1
00141 psi = ipsi * pi_180
00142 r = ( a**2 - f**2 ) / ( a + f*cos(psi) )
00143 du = r**2
00144 if (ipsi<360) usum = usum + du
00145 end do
00146 ufac = u / usum
00147
00148 t = 273.15_dp
00149 co2 = 0.0_dp
00150 co2layer = .false.
00151 do year = 1, yearmax
00152 usum = 0.0_dp
00153 do ipsi = 0, 360, 1
00154 psi = ipsi * pi_180
00155 r = ( a**2 - f**2 ) / ( a + f*cos(psi) )
00156 du = ufac * r**2
00157 if (ipsi<360) usum = usum + du
00158 lambda = psi - psi0
00159 delta = asin( seps * sin(lambda) )
00160 sdelta = sin(delta)
00161 cdelta = cos(delta)
00162 do iphi = -89, 89, 1
00163 phi = iphi * pi_180
00164 tptd = -tan(phi) * tan(delta)
00165 if ( tptd .le. -1.0_dp ) then
00166 tau0 = pi
00167 else if ( tptd .ge. 1.0_dp ) then
00168 tau0 = 0.0_dp
00169 else
00170 tau0 = acos( tptd )
00171 end if
00172 w = ( tau0*sin(phi)*sdelta + &
00173 cos(phi)*cdelta*sin(tau0) ) &
00174 * j0 * a**2 / ( pi * r**2 )
00175 wg = 0.0_dp
00176 if ( t(iphi) < -999._dp ) then
00177 albact = 0.95_dp
00178 albact_co2 = 0.95_dp
00179 else
00180 albact = alb
00181 albact_co2 = alb_co2
00182 end if
00183 teq = ( ( wg + (1._dp-albact) * w ) / SB )**.25_dp
00184 if ( teq < tco2 ) then
00185 co2layer(iphi) = .true.
00186 co2(iphi) = co2(iphi) + SB * tco2**4 * du &
00187 - (1._dp-albact_co2) * w * du
00188 t(iphi) = tco2
00189 else
00190 if ( co2layer(iphi) ) then
00191 co2(iphi) = co2(iphi) + SB * tco2**4 * du &
00192 - (1._dp-albact_co2) * w * du
00193 if ( co2(iphi) .le. 0.0_dp ) then
00194 co2(iphi) = 0.0_dp
00195 co2layer(iphi) = .false.
00196 t(iphi) = tco2
00197 end if
00198 else
00199 t(iphi) = teq
00200 end if
00201 end if
00202 if ( year .eq. yearmax ) o%t(ipsi,iphi) = t(iphi)
00203 end do
00204 if ( year .eq. yearmax .and. ipsi .eq. 0 ) then
00205 o%tam = 0.0_dp
00206 o%tmax = 0.0_dp
00207 else
00208 o%tam = o%tam + t
00209 do iphi = -89, 89
00210 o%tmax(iphi) = max( o%tmax(iphi), t(iphi) )
00211 end do
00212 end if
00213 end do
00214 o%tam = o%tam / 360._dp
00215 end do
00216 o%t(:,-90) = o%t(:,-89) + ( o%t(:,-89) - o%t(:,-88) ) / 2._dp
00217 o%t(:, 90) = o%t(:, 89) + ( o%t(:, 89) - o%t(:, 88) ) / 2._dp
00218 o%tam(-90) = o%tam(-89) + ( o%tam(-89) - o%tam(-88) ) / 2._dp
00219 o%tam( 90) = o%tam( 89) + ( o%tam( 89) - o%tam( 88) ) / 2._dp
00220 o%tmax(-90) = o%tmax(-89) + ( o%tmax(-89) - o%tmax(-88) ) / 2._dp
00221 o%tmax( 90) = o%tmax( 89) + ( o%tmax( 89) - o%tmax( 88) ) / 2._dp
00222
00223 end subroutine
00224
00225
00226 !> Annual mean temperature at latitude phi.
00227 !<------------------------------------------------------------------------------
00228 real(dp) function instam ( o, phi )
00229
00230 implicit none
00231 type(ins) :: o
00232 real(dp) :: phi
00233 integer(i4b) :: iphi1, iphi2
00234
00235 iphi1 = nint( phi - 0.5_dp )
00236 if ( iphi1 < -90 ) then
00237 iphi1 = -90
00238 else if ( iphi1 > 89 ) then
00239 iphi1 = 89
00240 end if
00241 iphi2 = iphi1 + 1
00242 instam = o%tam(iphi1) + ( o%tam(iphi2) - o%tam(iphi1) ) &
00243 * ( phi - iphi1 )
00244
00245 end function instam
00246
00247
00248 !> Annual maximum temperature at latitude phi.
00249 !<------------------------------------------------------------------------------
00250 real(dp) function instmax ( o, phi )
00251
00252 implicit none
00253 type(ins) :: o
00254 real(dp) :: phi
00255 integer(i4b) :: iphi1, iphi2
00256
00257 iphi1 = nint( phi - 0.5_dp )
00258 if ( iphi1 < -90 ) then
00259 iphi1 = -90
00260 else if ( iphi1 > 89 ) then
00261 iphi1 = 89
00262 end if
00263 iphi2 = iphi1 + 1
00264 instmax = o%tmax(iphi1) + ( o%tmax(iphi2) - o%tmax(iphi1) ) &
00265 * ( phi - iphi1 )
00266
00267 end function instmax
00268
00269
00270 !> Temperature at orbit position psi and latitude phi.
00271 !<------------------------------------------------------------------------------
00272 real(dp) function inst ( o, psi, phi )
00273
00274 implicit none
00275 type(ins) :: o
00276 real(dp) :: psi, phi
00277 integer(i4b) :: ipsi1, ipsi2
00278
00279 ipsi1 = nint( psi - 0.5_dp )
00280 if ( ipsi1 < 0 ) then
00281 ipsi1 = 0
00282 else if ( ipsi1 > 359 ) then
00283 ipsi1 = 359
00284 end if
00285 ipsi2 = ipsi1 + 1
00286 inst = inst1(o,ipsi1,phi) + ( inst1(o,ipsi2,phi) - inst1(o,ipsi1,phi) ) &
00287 * ( psi - ipsi1 )
00288
00289 end function inst
00290
00291
00292 !> Temperature at orbit position ipsi (integer value) and latitude phi.
00293 !<------------------------------------------------------------------------------
00294 real(dp) function inst1 ( o, ipsi, phi )
00295
00296 implicit none
00297 type(ins) :: o
00298 integer(i4b) :: ipsi
00299 real(dp) :: phi
00300 integer(i4b) :: iphi1, iphi2
00301
00302 iphi1 = nint( phi - 0.5_dp )
00303 if ( iphi1 < -90 ) then
00304 iphi1 = -90
00305 else if ( iphi1 > 89 ) then
00306 iphi1 = 89
00307 end if
00308 iphi2 = iphi1 + 1
00309 inst1 = o%t(ipsi,iphi1) + ( o%t(ipsi,iphi2) - o%t(ipsi,iphi1) ) &
00310 * ( phi - iphi1 )
00311
00312 end function inst1
00313
00314 end module instemp
00315