7 !! Mars Atmosphere-Ice Coupler MAIC-1.5:
8 !! Computation of the surface temperature (must be less than 0 deg C)
9 !! and of the accumulation-ablation rate for the north polar cap of Mars.
10 !! Computation of the geothermal heat flux.
14 !! Copyright 2009-2013 Ralf Greve
18 !! This file is part of SICOPOLIS.
20 !! SICOPOLIS is free software: you can redistribute it and/or modify
21 !! it under the terms of the GNU General Public License as published by
22 !! the Free Software Foundation, either version 3 of the License, or
23 !! (at your option) any later version.
25 !! SICOPOLIS is distributed in the hope that it will be useful,
26 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
27 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28 !! GNU General Public License for more details.
30 !! You should have received a copy of the GNU General Public License
31 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
33 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
36 !> Mars Atmosphere-Ice Coupler MAIC-1.5:
37 !! Computation of the surface temperature (must be less than 0 deg C)
38 !! and of the accumulation-ablation rate for the north polar cap of Mars.
39 !! Computation of the geothermal heat flux.
40 !<------------------------------------------------------------------------------
42 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
53 real(dp),
intent(in) :: time, dtime, dxi, deta
55 real(dp),
intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
56 real(dp),
intent(inout) :: z_sl
65 integer(i4b) :: i_gr, i_kl
66 integer(i4b) :: ndata_insol
68 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
69 real(dp) :: time_gr, time_kl
70 real(dp) :: z_sle_present, z_sle_help
71 real(dp),
dimension(0:JMAX,0:IMAX) :: temp_ma
72 real(dp) :: eld, g_mb, accum_factor
73 real(dp) :: erosion_chasm
74 real(dp),
dimension(0:JMAX,0:IMAX) :: et
75 real(dp) :: t_obliq_main, t_obliq_mod, &
76 obliq_ampl_max, obliq_ampl_min, obliq_ampl, obliq0, obliq, &
78 real(dp) :: insol_ma_90_now, insol_ma_90_present, &
79 obl_now, obl_present, ecc_now, ecc_present, &
80 ave_now, ave_present, cp_now, cp_present, &
81 temp_ma_90, temp_ma_90_present, zs_90
82 real(dp),
dimension(0:JMAX,0:IMAX) :: dist
83 real(dp) :: ave_data_i_kl, ave_data_i_gr
84 real(dp) :: q_geo_chasm
85 logical,
dimension(0:JMAX,0:IMAX) :: check_point
88 type (ins) :: temp_now, temp_present
91 real(dp),
parameter :: &
92 time_present = 0.0_dp, &
93 zs_90_present = -2.0e+03_dp
95 real(dp),
parameter :: &
100 real(dp),
parameter :: &
101 lambda_h2o = 2.86e+06_dp, &
105 real(dp),
parameter :: &
106 temp_s_min = -125.0_dp
126 delta_ts = -sine_amplit &
127 *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
134 t_obliq_main = 1.25e+05_dp *year_sec
135 t_obliq_mod = 1.3e+06_dp *year_sec
137 obliq0 = 25.2_dp *pi_180
138 obliq_ampl_max = 10.0_dp *pi_180
139 obliq_ampl_min = 2.5_dp *pi_180
141 obliq_ampl = 0.5_dp*(obliq_ampl_max+obliq_ampl_min) &
142 -0.5_dp*(obliq_ampl_max-obliq_ampl_min) &
143 *cos(2.0_dp*pi*time/t_obliq_mod)
145 obliq = obliq0 + obliq_ampl*sin(2.0_dp*pi*time/t_obliq_main)
152 insol_ma_90_now = (sol/pi)*sin(obliq)/sqrt(1.0_dp-ecc**2)
154 insol_ma_90_present = (sol/pi)*sin(obliq0)/sqrt(1.0_dp-ecc0**2)
157 #elif ( TSURFACE==5 || TSURFACE==6 )
161 ndata_insol = (insol_time_max-insol_time_min)/insol_time_stp
163 if (time/year_sec.lt.
real(insol_time_min,dp)) then
165 insol_ma_90_now = insol_ma_90(0)
166 obl_now = obl_data(0)
167 ecc_now = ecc_data(0)
168 ave_now = ave_data(0)
171 else if (time/year_sec.lt.
real(insol_time_max,dp)) then
173 i_kl = floor(((time/year_sec) &
174 -
real(insol_time_min,dp))/
real(insol_time_stp,dp))
177 i_gr = ceiling(((time/year_sec) &
178 -
real(insol_time_min,dp))/
real(insol_time_stp,dp))
179 i_gr = min(i_gr, ndata_insol)
181 if (i_kl.eq.i_gr)
then
183 insol_ma_90_now = insol_ma_90(i_kl)
184 obl_now = obl_data(i_kl)
185 ecc_now = ecc_data(i_kl)
186 ave_now = ave_data(i_kl)
187 cp_now = cp_data(i_kl)
191 time_kl = (insol_time_min + i_kl*insol_time_stp) *year_sec
192 time_gr = (insol_time_min + i_gr*insol_time_stp) *year_sec
194 insol_ma_90_now = insol_ma_90(i_kl) &
195 +(insol_ma_90(i_gr)-insol_ma_90(i_kl)) &
196 *(time-time_kl)/(time_gr-time_kl)
197 obl_now = obl_data(i_kl) &
198 +(obl_data(i_gr)-obl_data(i_kl)) &
199 *(time-time_kl)/(time_gr-time_kl)
200 ecc_now = ecc_data(i_kl) &
201 +(ecc_data(i_gr)-ecc_data(i_kl)) &
202 *(time-time_kl)/(time_gr-time_kl)
204 if ( abs(ave_data(i_gr)-ave_data(i_kl)) < pi )
then
205 ave_data_i_kl = ave_data(i_kl)
206 ave_data_i_gr = ave_data(i_gr)
208 if ( ave_data(i_gr) > ave_data(i_kl) )
then
209 ave_data_i_kl = ave_data(i_kl) + 2.0_dp*pi
210 ave_data_i_gr = ave_data(i_gr)
212 ave_data_i_kl = ave_data(i_kl)
213 ave_data_i_gr = ave_data(i_gr) + 2.0_dp*pi
217 ave_now = ave_data_i_kl &
218 +(ave_data_i_gr-ave_data_i_kl) &
219 *(time-time_kl)/(time_gr-time_kl)
221 cp_now = cp_data(i_kl) &
222 +(cp_data(i_gr)-cp_data(i_kl)) &
223 *(time-time_kl)/(time_gr-time_kl)
229 insol_ma_90_now = insol_ma_90(ndata_insol)
230 obl_now = obl_data(ndata_insol)
231 ecc_now = ecc_data(ndata_insol)
232 ave_now = ave_data(ndata_insol)
233 cp_now = cp_data(ndata_insol)
239 if (time_present/year_sec.lt.
real(insol_time_min,dp)) then
241 insol_ma_90_present = insol_ma_90(0)
242 obl_present = obl_data(0)
243 ecc_present = ecc_data(0)
244 ave_present = ave_data(0)
245 cp_present = cp_data(0)
247 else if (time_present/year_sec.lt.
real(insol_time_max,dp)) then
249 i_kl = floor(((time_present/year_sec) &
250 -
real(insol_time_min,dp))/
real(insol_time_stp,dp))
253 i_gr = ceiling(((time_present/year_sec) &
254 -
real(insol_time_min,dp))/
real(insol_time_stp,dp))
255 i_gr = min(i_gr, ndata_insol)
257 if (i_kl.eq.i_gr)
then
259 insol_ma_90_present = insol_ma_90(i_kl)
260 obl_present = obl_data(i_kl)
261 ecc_present = ecc_data(i_kl)
262 ave_present = ave_data(i_kl)
263 cp_present = cp_data(i_kl)
267 time_kl = (insol_time_min + i_kl*insol_time_stp) *year_sec
268 time_gr = (insol_time_min + i_gr*insol_time_stp) *year_sec
270 insol_ma_90_present = insol_ma_90(i_kl) &
271 +(insol_ma_90(i_gr)-insol_ma_90(i_kl)) &
272 *(time_present-time_kl)/(time_gr-time_kl)
273 obl_present = obl_data(i_kl) &
274 +(obl_data(i_gr)-obl_data(i_kl)) &
275 *(time_present-time_kl)/(time_gr-time_kl)
276 ecc_present = ecc_data(i_kl) &
277 +(ecc_data(i_gr)-ecc_data(i_kl)) &
278 *(time_present-time_kl)/(time_gr-time_kl)
280 if ( abs(ave_data(i_gr)-ave_data(i_kl)) < pi )
then
281 ave_data_i_kl = ave_data(i_kl)
282 ave_data_i_gr = ave_data(i_gr)
284 if ( ave_data(i_gr) > ave_data(i_kl) )
then
285 ave_data_i_kl = ave_data(i_kl) + 2.0_dp*pi
286 ave_data_i_gr = ave_data(i_gr)
288 ave_data_i_kl = ave_data(i_kl)
289 ave_data_i_gr = ave_data(i_gr) + 2.0_dp*pi
293 ave_present = ave_data_i_kl &
294 +(ave_data_i_gr-ave_data_i_kl) &
295 *(time_present-time_kl)/(time_gr-time_kl)
297 cp_present = cp_data(i_kl) &
298 +(cp_data(i_gr)-cp_data(i_kl)) &
299 *(time_present-time_kl)/(time_gr-time_kl)
305 insol_ma_90_present = insol_ma_90(ndata_insol)
306 obl_present = obl_data(ndata_insol)
307 ecc_present = ecc_data(ndata_insol)
308 ave_present = ave_data(ndata_insol)
309 cp_present = cp_data(ndata_insol)
315 #if (TSURFACE==4 || TSURFACE==5)
319 temp_ma_90 = sqrt(sqrt(insol_ma_90_now*(1.0_dp-albedo)/(epsilon*sigma))) &
325 = sqrt(sqrt(insol_ma_90_present*(1.0_dp-albedo)/(epsilon*sigma))) &
330 delta_ts = temp_ma_90 - temp_ma_90_present
337 ecc = ecc_now, ave = ave_now*pi_180_inv, &
338 obl = obl_now*pi_180_inv, sa = albedo, ct = 148.7_dp)
340 temp_ma_90 =
instam(temp_now, 90.0_dp) + delta_ts0 &
346 ecc = ecc_present, ave = ave_present*pi_180_inv, &
347 obl = obl_present*pi_180_inv, sa = albedo, ct = 148.7_dp)
349 temp_ma_90_present =
instam(temp_present, 90.0_dp) + delta_ts0 &
354 delta_ts = temp_ma_90 - temp_ma_90_present
366 stop
' boundary: SEA_LEVEL==2 not allowed for nmars application!'
370 stop
' boundary: SEA_LEVEL==3 not allowed for nmars application!'
376 if ( z_sl_old > -999999.9_dp )
then
377 dzsl_dtau = (z_sl-z_sl_old)/dtime
386 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
388 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 )
389 z_mar = fact_z_mar*z_sl
390 #elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 )
391 if (z_sl >= -80.0_dp)
then
394 z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
396 z_mar = fact_z_mar*z_mar
408 check_point(j,i) = .false.
414 if (maske(j,i).ge.2)
then
415 check_point(j ,i ) = .true.
416 check_point(j ,i+1) = .true.
417 check_point(j ,i-1) = .true.
418 check_point(j+1,i ) = .true.
419 check_point(j-1,i ) = .true.
426 if (check_point(j,i))
then
436 if (check_point(j,i))
then
437 maske(j,i) = maske_neu(j,i)
444 call
borehole(zs, 0.0_dp, 0.0_dp, dxi, deta,
'grid', zs_90)
452 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
454 temp_ma_90_present = temp0_ma_90n
456 temp_ma(j,i) = temp_ma_90_present &
457 + (gamma_ma)*(zs(j,i)-zs_90_present) &
458 + (c_ma*pi_180_inv)*((90.0_dp*pi_180)-phi(j,i))
460 #elif (TSURFACE==4 || TSURFACE==5)
462 temp_ma(j,i) = temp_ma_90_present &
463 + (gamma_ma)*(zs(j,i)-zs_90_present) &
464 + (c_ma*pi_180_inv)*((90.0_dp*pi_180)-phi(j,i))
470 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
472 temp_ma(j,i) = temp_ma(j,i) + delta_ts
474 #elif (TSURFACE==4 || TSURFACE==5)
476 temp_ma(j,i) = temp_ma(j,i) + delta_ts &
477 - (gamma_ma)*(zs_90-zs_90_present)
485 temp_ma(j,i) =
instam(temp_now, phi(j,i)*pi_180_inv) + delta_ts0 &
498 accum_factor = 1.0_dp+gamma_s*delta_ts
500 accum_factor = exp(gamma_s*delta_ts)
502 accum_factor = exp(lambda_h2o/(r_h2o*temp0_ref) &
503 -lambda_h2o/(r_h2o*(temp0_ref+delta_ts)))
508 accum(j,i) = accum_present(j,i)*accum_factor
509 if (accum(j,i).lt.0.0_dp) stop
' boundary: Negative accumulation rate!'
515 #if ( ABLSURFACE==1 || ABLSURFACE==2)
517 eld = eld_0 *1.0e+03_dp
520 g_mb = g_0 *(1.0e-06_dp/year_sec)*(rho_w/rho_i) &
521 *(1.0_dp/(1.0_dp-frac_dust))
524 g_mb = g_0 *(1.0e-06_dp/year_sec)
529 eld = eld*(1.0_dp+gamma_eld*delta_ts)
530 g_mb = g_mb*(gamma_g*accum_factor)
532 eld = eld*exp(gamma_eld*delta_ts)
533 g_mb = g_mb*(gamma_g*accum_factor)
536 if (eld.lt.eps) stop
' boundary: Negative equilibrium line distance eld!'
540 dist(j,i) = sqrt( xi(i)**2 + eta(j)**2 )
549 erosion_chasm = erosion_chasm *(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
550 *(1.0_dp/(1.0_dp-frac_dust))
553 erosion_chasm = erosion_chasm *(1.0e-03_dp/year_sec)
562 #if ( ABLSURFACE==1 || ABLSURFACE==2 )
564 as_perp(j,i) = g_mb*(eld-dist(j,i))
565 if (as_perp(j,i).gt.accum(j,i)) as_perp(j,i) = accum(j,i)
569 evap(j,i) = accum(j,i) - as_perp(j,i)
576 if ( (maske_chasm(j,i) == 7) &
577 .and.(time >= time_chasm_init) &
578 .and.(time <= time_chasm_end) )
then
579 evap(j,i) = evap(j,i) + erosion_chasm
580 as_perp(j,i) = accum(j,i) - evap(j,i)
590 temp_s = min(temp_ma, -eps)
592 temp_s = max(temp_s, temp_s_min)
603 q_geo_chasm = q_geo_chasm *1.0e-03_dp
608 if ( (maske_chasm(j,i) == 7) &
609 .and.(time >= time_chasm_init) &
610 .and.(time <= time_chasm_end) )
then
611 q_geo(j,i) = q_geo_chasm
613 q_geo(j,i) = q_geo_normal(j,i)