MAIC-2  Revision 19
 All Classes Files Functions Variables
maic2.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Program : m a i c 2
4 !
5 !> @mainpage
6 !!
7 !! @section Description
8 !!
9 !! The Mars Atmosphere-Ice Coupler MAIC-2 is a simple, latitudinal model,
10 !! which consists of a set of parameterisations for the surface temperature,
11 !! the atmospheric water transport and the surface mass balance (condensation
12 !! minus evaporation) of water ice. It is driven directly by orbital
13 !! parameters. A detailed description of the model is given in the publication
14 !! by Greve et al. (2010).
15 !!
16 !! The model equations of MAIC-2 are discretised by a
17 !! finite-difference/finite-volume scheme.
18 !! Coding is done in the programming language Fortran 90.
19 !!
20 !! Required model forcing (as functions of time):
21 !! @li Obliquity (axial tilt).
22 !! @li Orbital eccentricity.
23 !! @li Solar longitude (Ls) of perihelion.
24 !!
25 !! Output (as functions of latitude and time):
26 !! @li Surface temperature.
27 !! @li Evaporation rate of water ice.
28 !! @li Condensation rate of water ice.
29 !! @li Atmospheric water content.
30 !! @li Surface mass balance of water ice.
31 !! @li Ice thickness.\n
32 !!
33 !! References:
34 !! @li Greve, R., B. Grieger and O. J. Stenzel. 2010.\n
35 !! MAIC-2, a latitudinal model for the Martian surface temperature,
36 !! atmospheric water transport and surface glaciation.\n
37 !! Planetary and Space Science 58 (6), 931-940.
38 !! @li MAIC-2 website: http://maic2.greveweb.net/
39 !!
40 !! @section Copyright
41 !!
42 !! Copyright 2010-2013 Ralf Greve, Bjoern Grieger, Oliver J. Stenzel
43 !!
44 !! @section License
45 !!
46 !! MAIC-2 is free software: you can redistribute it and/or modify
47 !! it under the terms of the GNU General Public License as published by
48 !! the Free Software Foundation, either version 3 of the License, or
49 !! (at your option) any later version.
50 !!
51 !! MAIC-2 is distributed in the hope that it will be useful,
52 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
53 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
54 !! GNU General Public License for more details.
55 !!
56 !! You should have received a copy of the GNU General Public License
57 !! along with MAIC-2. If not, see <http://www.gnu.org/licenses/>.
58 !!
59 !! @file
60 !!
61 !! Main program file of MAIC-2.
62 !!
63 !! @section Copyright
64 !!
65 !! Copyright 2010-2013 Ralf Greve, Bjoern Grieger, Oliver J. Stenzel
66 !!
67 !! @section License
68 !!
69 !! This file is part of MAIC-2.
70 !!
71 !! MAIC-2 is free software: you can redistribute it and/or modify
72 !! it under the terms of the GNU General Public License as published by
73 !! the Free Software Foundation, either version 3 of the License, or
74 !! (at your option) any later version.
75 !!
76 !! MAIC-2 is distributed in the hope that it will be useful,
77 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
78 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
79 !! GNU General Public License for more details.
80 !!
81 !! You should have received a copy of the GNU General Public License
82 !! along with MAIC-2. If not, see <http://www.gnu.org/licenses/>.
83 !<
84 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
85 
86 !-------- Inclusion of specification header --------
87 
88 #include "maic2_specs.h"
89 
90 !-------- Inclusion of kind-type and global-variable modules --------
91 
92 #include "subroutines/maic2_types.F90"
93 #include "subroutines/maic2_variables.F90"
94 
95 !-------- Inclusion of further modules --------
96 
97 #include "subroutines/instemp.f90"
98 #include "subroutines/evaporation.f90"
99 #include "subroutines/condensation.f90"
100 
101 !-------------------------------------------------------------------------------
102 !> Main program of MAIC-2.
103 !<------------------------------------------------------------------------------
104 program maic2
105 
106 !-------- Modules, variables --------
107 
108 use maic2_types
109 use maic2_variables
110 
111 implicit none
112 integer(i4b) :: l, n
113 integer(i4b) :: ios
114 integer(i4b) :: itercount_max
115 integer(i4b) :: ndata_insol
116 real(dp) :: time, time_init, time_end, dtime
117 real(dp) :: ls, psi
118 real(dp) :: dphi_equi
119 real(dp) :: d_dummy
120 character (len=100) :: runname
121 character :: ch_dummy
122 logical :: output_flag
123 
124 !-------- Initializations --------
125 
126 ! ------ Setting of physical parameters
127 
128  rho_i = 9.1e+02_dp
129 ! Density of ice = 910 kg/m3
130 !
131  rho_w = 1.0e+03_dp
132 ! Density of pure water = 1000 kg/m3
133 !
134  g = 3.72_dp
135 ! Gravity acceleration = 3.72 m/s2
136 !
137  r = 3.396e+06_dp
138 ! Radius of Mars = 3396 km
139 
140 ! ------ Density of ice-dust mixture
141 
142 rho = rho_i ! so far no dust considered
143 
144 rho_inv = 1.0_dp/rho
145 
146 ! ------ Specification of current simulation
147 
148 runname = runname
149 
150 ! ------ Conversion of time quantities
151 
152 time_init = time_init0*year_sec ! a --> s
153 time_end = time_end0*year_sec ! a --> s
154 dtime = dtime0*year_sec ! a --> s
155 
156 #if OUTPUT==1
157 dtime_out = dtime_out0 * year_sec ! a --> s
158 #elif OUTPUT==2
159 n_output = n_output
160 time_output( 1) = time_out0_01 * year_sec ! a --> s
161 time_output( 2) = time_out0_02 * year_sec ! a --> s
162 time_output( 3) = time_out0_03 * year_sec ! a --> s
163 time_output( 4) = time_out0_04 * year_sec ! a --> s
164 time_output( 5) = time_out0_05 * year_sec ! a --> s
165 time_output( 6) = time_out0_06 * year_sec ! a --> s
166 time_output( 7) = time_out0_07 * year_sec ! a --> s
167 time_output( 8) = time_out0_08 * year_sec ! a --> s
168 time_output( 9) = time_out0_09 * year_sec ! a --> s
169 time_output(10) = time_out0_10 * year_sec ! a --> s
170 time_output(11) = time_out0_11 * year_sec ! a --> s
171 time_output(12) = time_out0_12 * year_sec ! a --> s
172 time_output(13) = time_out0_13 * year_sec ! a --> s
173 time_output(14) = time_out0_14 * year_sec ! a --> s
174 time_output(15) = time_out0_15 * year_sec ! a --> s
175 time_output(16) = time_out0_16 * year_sec ! a --> s
176 time_output(17) = time_out0_17 * year_sec ! a --> s
177 time_output(18) = time_out0_18 * year_sec ! a --> s
178 time_output(19) = time_out0_19 * year_sec ! a --> s
179 time_output(20) = time_out0_20 * year_sec ! a --> s
180 #endif
181 
182 ! ------ Reading of data for orbital parameters
183 
184 insol_ma_90 = 0.0_dp ! Assignment of dummy values
185 obl_data = 0.0_dp
186 ecc_data = 0.0_dp
187 ave_data = 0.0_dp
188 cp_data = 0.0_dp
189 
190 open(21, iostat=ios, &
191  file=inpath//'/'//insol_ma_90n_file, &
192  status='old')
193 if (ios /= 0) stop ' Error when opening the data file for orbital parameters!'
194 
195 read(21,*) ch_dummy, insol_time_min, insol_time_stp, insol_time_max
196 
197 if (ch_dummy /= '#') then
198  write(6,*) 'insol_time_min, insol_time_stp, insol_time_max not defined in'
199  write(6,*) 'data file!'
200 end if
201 
202 ndata_insol = (insol_time_max-insol_time_min)/insol_time_stp
203 
204 if (ndata_insol > 100000) &
205  stop 'Too many data in orbital-parameter-data file!'
206 
207 do n=0, ndata_insol
208  read(21,*) d_dummy, ecc_data(n), obl_data(n), cp_data(n), ave_data(n), insol_ma_90(n)
209  obl_data(n) = obl_data(n) *pi_180 ! deg -> rad
210  ave_data(n) = ave_data(n) *pi_180 ! deg -> rad
211 end do
212 
213 close(21, status='keep')
214 
215 ! ------ Numerical grid
216 
217 ! ---- Nodes (grid points)
218 
219 dphi_equi = 180.0_dp/lmax *pi_180 ! deg -> rad
220 
221 do l=0, lmax
222  phi_node(l) = -90.0_dp*pi_180 + dble(l)*dphi_equi
223  ! Grid points laid out between -90 deg (90S, south pole)
224  ! and +90 deg (90N, north pole).
225  ! So far only equidistant grid points implemented.
226 end do
227 
228 ! ---- Spacing
229 
230 do l=1, lmax
231  dphi(l) = phi_node(l) - phi_node(l-1)
232  dphi_inv(l) = 1.0_dp/dphi(l)
233 end do
234 
235 ! ---- Lower cell boundaries
236 
237 phi_cb1(0) = phi_node(0)
238 do l=1, lmax
239  phi_cb1(l) = 0.5_dp*(phi_node(l-1)+phi_node(l))
240 end do
241 
242 cos_phi_cb1 = cos(phi_cb1)
243 sin_phi_cb1 = sin(phi_cb1)
244 
245 ! ---- Upper cell boundaries
246 
247 do l=0, lmax-1
248  phi_cb2(l) = 0.5_dp*(phi_node(l)+phi_node(l+1))
249 end do
250 phi_cb2(lmax) = phi_node(lmax)
251 
252 cos_phi_cb2 = cos(phi_cb2)
253 sin_phi_cb2 = sin(phi_cb2)
254 
255 ! ---- Auxiliary quantity needed for the diffusional transport
256 
257 do l=0, lmax
258  diff_aux(l) = diff_water_maic/(r**2*(sin_phi_cb2(l)-sin_phi_cb1(l)))
259 end do
260 
261 ! ------ Definition of initial values
262 
263 #if H_INIT==1
264 
265 h = thick_init
266 
267 #elif H_INIT==2
268 
269 do l=0, lmax
270  h(l) = 3000.0_dp &
271  * ( 1.0_dp-(90.0_dp*pi_180-abs(phi_node(l)))**2/(10.0_dp*pi_180)**2 )
272  h(l) = max(h(l), -0.1_dp)
273 end do
274 
275 #endif
276 
277 water = water_init
278 
279 ! ------ Output file
280 
281 #if OUTPUT==1
282 iter_out = nint(dtime_out/dtime)
283 #elif OUTPUT==2
284 do n=1, n_output
285  iter_output(n) = nint((time_output(n)-time_init)/dtime)
286 end do
287 #endif
288 
289 open(12, iostat=ios, &
290  file=outpath//'/'//trim(runname)//'.out', &
291  status='replace')
292 if (ios /= 0) stop ' Error when opening the output file runname.out!'
293 
294 !-------- Main loop --------
295 
296 ! ------ Begin of main loop
297 
298 write(6,*) ' '
299 
300 itercount=1; write(6,'(i10)') itercount
301 
302 itercount_max = nint((time_end-time_init)/dtime)
303 
304 main_loop : do itercount=1, itercount_max
305 
306  if ( mod(itercount, 1000) == 0 ) &
307  write(6,'(i10)') itercount
308 
309 ! ------ Update of time
310 
311  time = time_init + real(itercount,dp)*dtime
312 
313 ! ------ Boundary conditions
314 
315  call boundary_maic2(time, ls, psi, dtime)
316 
317 ! ------ Topography
318 
319  call calc_top_maic2(time, dtime)
320 
321  h = h_new ! new values -> old values
322 
323 ! ------ Data output
324 
325 output_flag = .false.
326 
327 #if OUTPUT == 1
328 
329 if ( mod(itercount, iter_out) == 0 ) output_flag = .true.
330 
331 #elif OUTPUT == 2
332 
333 do n=1, n_output
334  if (itercount == iter_output(n)) output_flag = .true.
335 end do
336 
337 #endif
338 
339 if ( output_flag ) call output(time, ls)
340 
341 end do main_loop ! End of main loop
342 
343 !-------- End of main program --------
344 
345 close(12, status='keep')
346 
347 write(6,'(a)') ' '
348 write(6,'(a)') ' '
349 write(6,'(a)') &
350 ' * * * maic2.F90 r e a d y * * *'
351 write(6,'(a)') ' '
352 write(6,'(a)') ' '
353 
354 end program maic2
355 
356 !-------- Inclusion of subroutines --------
357 
358 #include "subroutines/boundary_maic2.F90"
359 #include "subroutines/get_orb_par.F90"
360 #include "subroutines/get_psi_tab.F90"
361 #include "subroutines/p_sat.f90"
362 #include "subroutines/diff_trans.F90"
363 #include "subroutines/tri_sle.F90"
364 #include "subroutines/calc_top_maic2.F90"
365 #include "subroutines/output.F90"
366 
367 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
368 ! End of maic2.F90
369 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
370 !