00001
00002
00003
00004
00005 !> @file
00006 !!
00007 !! Determination of the surface temperature and of the net mass balance
00008 !! (accumulation-ablation rate) for the polar caps of Mars.
00009 !!
00010 !! @section Copyright
00011 !!
00012 !! Copyright 2010, 2011 Ralf Greve, Bjoern Grieger, Oliver J. Stenzel
00013 !!
00014 !! @section License
00015 !!
00016 !! This file is part of MAIC-2.
00017 !!
00018 !! MAIC-2 is free software: you can redistribute it and/or modify
00019 !! it under the terms of the GNU General Public License as published by
00020 !! the Free Software Foundation, either version 3 of the License, or
00021 !! (at your option) any later version.
00022 !!
00023 !! MAIC-2 is distributed in the hope that it will be useful,
00024 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
00025 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00026 !! GNU General Public License for more details.
00027 !!
00028 !! You should have received a copy of the GNU General Public License
00029 !! along with MAIC-2. If not, see <http://www.gnu.org/licenses/>.
00030 !<
00031 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00032
00033
00034 !> Determination of the surface temperature and of the net mass balance
00035 !! (accumulation-ablation rate) for the polar caps of Mars.
00036 !<------------------------------------------------------------------------------
00037 subroutine boundary_maic2(time, dtime)
00038
00039 use maic2_types
00040 use maic2_variables
00041 use instemp
00042 use evaporation
00043 use condensation
00044
00045 implicit none
00046
00047 real(dp), intent(in) :: time, dtime
00048
00049 integer(i4b) :: l, n, n1, n2
00050 real(dp) :: temp_co2_mean
00051 real(dp) :: ecc, obl, cp, ave, insol_ma_90NS, time_help, psi
00052 real(dp) :: evap_coeff, tau_cond
00053 real(dp) :: dtime_inv
00054 type (ins), save :: temp
00055
00056 logical :: first_iteration = .true.
00057 real(dp) :: time_of_last_temp_update = 0.0
00058
00059 dtime_inv = 1.0_dp/dtime
00060
00061
00062
00063 p_surf = P_SURF
00064
00065 temp_co2 = 3182.48_dp/(23.3494_dp-log(0.01_dp*p_surf))
00066 temp_co2_mean = sum(temp_co2)/real(size(temp_co2),dp)
00067
00068
00069
00070
00071 if ( first_iteration &
00072 .or. time - time_of_last_temp_update > (1.0e+03_dp*YEAR_SEC) ) then
00073
00074 print*, 'Surface temperature update ...'
00075 first_iteration = .false.
00076 time_of_last_temp_update = time
00077
00078 call get_orb_par(time, ecc, obl, cp, ave, insol_ma_90NS)
00079
00080 call get_psi_tab(ecc, ave)
00081
00082 call setinstemp(temp, ecc = ecc, ave = ave*pi_180_inv, &
00083 obl = obl*pi_180_inv, &
00084 sa = ALBEDO, sac = ALBEDO_CO2, op = MARS_YEAR, &
00085 ct = temp_co2_mean)
00086
00087 end if
00088
00089
00090
00091 time_help = real(NTIME,dp)*modulo(time/MARS_YEAR, 1.0_dp)
00092
00093 n1 = floor(time_help)
00094 n2 = n1 + 1
00095
00096 if ((n1 < 0).or.(n2 > NTIME)) &
00097 stop ' Subroutine boundary_maic2: n1 or n2 out of range!'
00098
00099 psi = psi_tab(n1) + (time_help-real(n1,dp))*(psi_tab(n2)-psi_tab(n1))
00100 psi = modulo(psi, 2.0_dp*pi)
00101
00102
00103
00104
00105
00106 do l=0, LMAX
00107 temp_surf(l) = inst(temp, psi*pi_180_inv, phi_node(l)*pi_180_inv)
00108 end do
00109
00110
00111
00112 if (TEMP_SURF_AMP_EQ < eps) then
00113
00114 temp_surf_amp = 0.0_dp
00115
00116 else
00117
00118 do l=0, LMAX
00119 temp_surf_amp(l) &
00120 = TEMP_SURF_AMP_EQ &
00121 * (1.0_dp-(abs(phi_node(l))*2.0_dp*pi_inv)**TEMP_SURF_AMP_EXP)
00122
00123
00124 end do
00125
00126 end if
00127
00128
00129
00130 evap = 0.0_dp
00131 cond = 0.0_dp
00132
00133
00134
00135 call setevappar_bd(evap_fact_para=EVAP_FACT, gamma_reg_para=GAMMA_REG, &
00136 g_para=G)
00137 call getevap_bd(temp_surf, temp_surf_amp, p_surf, H, evap)
00138
00139
00140
00141 call diff_trans(dtime)
00142
00143
00144 do l=0, LMAX
00145
00146 if (water_new(l) >= -eps) then
00147 water(l) = max(water_new(l),0.0_dp)
00148 else
00149 stop ' Negative water content!'
00150 end if
00151
00152 end do
00153
00154
00155
00156 #if COND==1 /* Removal of water exceeding the surface saturation pressure */
00157
00158 call setcondpar(gravity=G)
00159 call getcond_1(temp_surf, water, cond, dtime)
00160
00161 #elif COND==2 /* Continuous, quadratic dependence on humidity */
00162
00163 tau_cond = TAU_COND*YEAR_SEC
00164 call setcondpar(gravity=G, timescale=tau_cond)
00165 call getcond_2(temp_surf, water, cond)
00166
00167 #endif
00168
00169 do l=0, LMAX
00170
00171 water_new(l) = water(l) - dtime * cond(l)
00172
00173 if (water_new(l) >= 0.0_dp) then
00174 water(l) = water_new(l)
00175 else
00176 cond(l) = cond(l) + water_new(l)*dtime_inv
00177 water(l) = 0.0_dp
00178 end if
00179
00180 end do
00181
00182
00183
00184 do l=0, LMAX
00185 a_net(l) = (cond(l)-evap(l))*rho_inv
00186 end do
00187
00188 end subroutine boundary_maic2
00189