MAIC-2  Revision 19
 All Classes Files Functions Variables
boundary_maic2.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : b o u n d a r y _ m a i c 2
4 !
5 !> @file
6 !!
7 !! Determination of the surface temperature and of the net mass balance
8 !! (accumulation-ablation rate) for the polar caps of Mars.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2010-2013 Ralf Greve, Bjoern Grieger, Oliver J. Stenzel
13 !!
14 !! @section License
15 !!
16 !! This file is part of MAIC-2.
17 !!
18 !! MAIC-2 is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
22 !!
23 !! MAIC-2 is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
27 !!
28 !! You should have received a copy of the GNU General Public License
29 !! along with MAIC-2. If not, see <http://www.gnu.org/licenses/>.
30 !<
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32 
33 !-------------------------------------------------------------------------------
34 !> Determination of the surface temperature and of the net mass balance
35 !! (accumulation-ablation rate) for the polar caps of Mars.
36 !<------------------------------------------------------------------------------
37 subroutine boundary_maic2(time, ls, psi, dtime)
38 
39 use maic2_types
41 use instemp
42 use evaporation
43 use condensation
44 
45 implicit none
46 
47 real(dp), intent(in) :: time, dtime
48 real(dp), intent(inout) :: ls, psi
49 
50 integer(i4b) :: l, n, n1, n2
51 real(dp) :: temp_co2_mean
52 real(dp) :: ecc, obl, cp, ave, insol_ma_90ns, time_help
53 real(dp) :: evap_coeff, tau_cond
54 real(dp) :: dtime_inv
55 type (ins), save :: temp
56 
57 logical :: first_iteration = .true.
58 real(dp) :: time_of_last_temp_update = 0.0
59 
60 dtime_inv = 1.0_dp/dtime
61 
62 !-------- Surface pressure, CO2 condensation temperature --------
63 
64 p_surf = p_surf
65 
66 temp_co2 = 3182.48_dp/(23.3494_dp-log(0.01_dp*p_surf))
67 temp_co2_mean = sum(temp_co2)/real(size(temp_co2),dp)
68 
69 !-------- Orbital parameters, table of true anomalies
70 ! and seasonal cycle of surface temperatures --------
71 
72 if ( first_iteration &
73  .or. time - time_of_last_temp_update > (1.0e+03_dp*year_sec) ) then
74 
75  print*, 'Surface temperature update ...'
76  first_iteration = .false.
77  time_of_last_temp_update = time
78 
79  call get_orb_par(time, ecc, obl, cp, ave, insol_ma_90ns)
80 
81  call get_psi_tab(ecc, ave)
82 
83  call setinstemp(temp, ecc = ecc, ave = ave*pi_180_inv, &
84  obl = obl*pi_180_inv, &
85  sa = albedo, sac = albedo_co2, op = mars_year, &
86  ct = temp_co2_mean)
87 
88 end if
89 
90 !-------- Solar longitude (orbital position with respect to vernal equinox)
91 ! and true anomaly (orbital position with respect to perihelion) --------
92 
93 time_help = real(ntime,dp)*modulo(time/mars_year, 1.0_dp)
94 
95 n1 = floor(time_help)
96 n2 = n1 + 1
97 
98 if ((n1 < 0).or.(n2 > ntime)) &
99  stop ' Subroutine boundary_maic2: n1 or n2 out of range!'
100 
101 ls = ls_tab(n1) + (time_help-real(n1,dp))*(ls_tab(n2)-ls_tab(n1))
102 ls = modulo(ls, 2.0_dp*pi)
103 
104 psi = psi_tab(n1) + (time_help-real(n1,dp))*(psi_tab(n2)-psi_tab(n1))
105 psi = modulo(psi, 2.0_dp*pi)
106 
107 !-------- Instantaneous surface temperature --------
108 
109 ! ------ Daily mean
110 
111 do l=0, lmax
112  temp_surf(l) = inst(temp, psi*pi_180_inv, phi_node(l)*pi_180_inv) ! in K
113 end do
114 
115 ! ------ Amplitude of daily cycle
116 
117 if (temp_surf_amp_eq < eps) then
118 
119  temp_surf_amp = 0.0_dp ! no daily cycle
120 
121 else
122 
123  do l=0, lmax
124  temp_surf_amp(l) &
125  = temp_surf_amp_eq &
126  * (1.0_dp-(abs(phi_node(l))*2.0_dp*pi_inv)**temp_surf_amp_exp)
127  ! in K
128 ! * max((cos(phi_node(l))),0.0_dp)**TEMP_SURF_AMP_EXP ! in K
129  end do
130 
131 end if
132 
133 !-------- Initializations of evaporation and condensation --------
134 
135 evap = 0.0_dp
136 cond = 0.0_dp
137 
138 !-------- Evaporation (Buoyant-diffusion approach by Ingersoll)--------
139 
140 call setevappar_bd(evap_fact_para=evap_fact, gamma_reg_para=gamma_reg, &
141  g_para=g)
142 call getevap_bd(temp_surf, temp_surf_amp, p_surf, h, evap)
143 
144 !-------- Diffusive transport --------
145 
146 call diff_trans(dtime) ! computed with cond = 0.0_dp;
147  ! condensation will be accounted for later
148 
149 do l=0, lmax
150 
151  if (water_new(l) >= -eps) then
152  water(l) = max(water_new(l),0.0_dp) ! new values -> old values
153  else ! physically meaningless
154  stop ' Negative water content!'
155  end if
156 
157 end do
158 
159 !-------- Condensation --------
160 
161 #if COND==1 /* Removal of water exceeding the surface saturation pressure */
162 
163 call setcondpar(gravity=g)
164 call getcond_1(temp_surf, water, cond, dtime)
165 
166 #elif COND==2 /* Continuous, quadratic dependence on humidity */
167 
168 tau_cond = tau_cond*year_sec ! a -> s
169 call setcondpar(gravity=g, timescale=tau_cond)
170 call getcond_2(temp_surf, water, cond)
171 
172 #endif
173 
174 do l=0, lmax
175 
176  water_new(l) = water(l) - dtime * cond(l)
177 
178  if (water_new(l) >= 0.0_dp) then
179  water(l) = water_new(l) ! new values -> old values
180  else ! physically meaningless, needs to be corrected
181  cond(l) = cond(l) + water_new(l)*dtime_inv
182  water(l) = 0.0_dp
183  end if
184 
185 end do
186 
187 !-------- Net mass balance --------
188 
189 do l=0, lmax
190  a_net(l) = (cond(l)-evap(l))*rho_inv
191 end do
192 
193 end subroutine boundary_maic2
194 !