SICOPOLIS V3.3
sicopolis.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Program : s i c o p o l i s
4 ! (SImulation COde for POLythermal Ice Sheets)
5 !
6 #define MODEL_SICOPOLIS
7 #define VERSION '3.3'
8 #define DATE '2017-08-03'
9 !
10 !> @mainpage
11 !!
12 !! @section Description
13 !!
14 !! SICOPOLIS (SImulation COde for POLythermal Ice Sheets) is a 3-d
15 !! dynamic/thermodynamic model that simulates the evolution of large ice
16 !! sheets and ice caps. It was originally created by Greve (1997a,b) in a
17 !! version for the Greenland ice sheet. Since then, SICOPOLIS has been
18 !! developed continuously and applied to problems of past, present and
19 !! future glaciation of Greenland, Antarctica, the entire northern
20 !! hemisphere, the polar ice caps of the planet Mars and others.
21 !!
22 !! The model is based on the shallow ice approximation for grounded ice, the
23 !! shallow shelf approximation for floating ice (e.g., Greve and Blatter 2009)
24 !! and, optionally, hybrid shallow-ice--shelfy-stream dynamics for ice streams
25 !! (Bernales et al. 2017). It is coded in Fortran and uses finite difference
26 !! discretisation on a staggered (Arakawa C) grid, the velocity
27 !! components being taken between grid points. A variety of different
28 !! thermodynamics solvers are available, namely the polythermal two-layer
29 !! method, two versions of the one-layer enthalpy method, the cold-ice
30 !! method and the isothermal method (Greve and Blatter 2016).
31 !!
32 !! The coding is based on a consequent low-tech philosophy. All
33 !! structures are kept as simple as possible, and advanced coding
34 !! techniques are only employed where it is deemed appropriate. The use
35 !! of external libraries is kept at an absolute minimum. In fact,
36 !! SICOPOLIS can be run without external libraries at all, which makes
37 !! the installation very easy and fast.
38 !!
39 !! Required model forcing:
40 !! @li Surface mass balance (precipitation, evaporation, runoff).
41 !! @li Mean annual air temperature above the ice.
42 !! @li Eustatic sea level.
43 !! @li Geothermal heat flux.
44 !!
45 !! Main output (as functions of position and time):
46 !! @li Extent and thickness of the ice sheet.
47 !! @li Velocity field.
48 !! @li Temperature field.
49 !! @li Water-content field (temperate regions).
50 !! @li Age of the ice.
51 !! @li Isostatic displacement and temperature of the lithosphere.
52 !!
53 !! References:
54 !! @li Bernales, J., I. Rogozhina, R. Greve and M. Thomas. 2017.\n
55 !! Comparison of hybrid schemes for the combination of
56 !! shallow approximations in numerical simulations of the
57 !! Antarctic Ice Sheet.\n
58 !! Cryosphere 11 (1), 247-265.
59 !! @li Greve, R. 1997a.\n
60 !! A continuum-mechanical formulation for shallow polythermal ice sheets.\n
61 !! Phil. Trans. R. Soc. A 355 (1726), 921-974.
62 !! @li Greve, R. 1997b.\n
63 !! Application of a polythermal three-dimensional ice sheet model to the
64 !! Greenland ice sheet: Response to steady-state and transient climate
65 !! scenarios.\n
66 !! J. Climate 10 (5), 901-918.
67 !! @li Greve, R. and H. Blatter. 2009.\n
68 !! Dynamics of Ice Sheets and Glaciers.\n
69 !! Springer, Berlin, Germany etc., 287 pp.
70 !! @li Greve, R. and H. Blatter. 2016.\n
71 !! Comparison of thermodynamics solvers in the polythermal ice sheet model
72 !! SICOPOLIS.\n
73 !! Polar Sci. 10 (1), 11-23.
74 !! @li SICOPOLIS website: <http://www.sicopolis.net/>
75 !! @li Changelog: <http://tinyurl.com/commit-logs-sico-trunk/>
76 !!
77 !! @section Copyright
78 !!
79 !! Copyright 2009-2017 Ralf Greve\n
80 !! (with contributions by Jorge Bernales, Heinz Blatter, Reinhard Calov,
81 !! Thorben Dunse, Ben Galton-Fenzi, Thomas Goelles, Philipp Hancke,
82 !! Nina Kirchner, Sascha Knell, Alex Robinson, Tatsuru Sato, Malte Thoma,
83 !! Roland Warner)
84 !!
85 !! @section License
86 !!
87 !! SICOPOLIS is free software: you can redistribute it and/or modify
88 !! it under the terms of the GNU General Public License as published by
89 !! the Free Software Foundation, either version 3 of the License, or
90 !! (at your option) any later version.
91 !!
92 !! SICOPOLIS is distributed in the hope that it will be useful,
93 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
94 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
95 !! GNU General Public License for more details.
96 !!
97 !! You should have received a copy of the GNU General Public License
98 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
99 !!
100 !! @file
101 !!
102 !! Main program file of SICOPOLIS.
103 !!
104 !! @section Copyright
105 !!
106 !! Copyright 2009-2017 Ralf Greve
107 !!
108 !! @section License
109 !!
110 !! This file is part of SICOPOLIS.
111 !!
112 !! SICOPOLIS is free software: you can redistribute it and/or modify
113 !! it under the terms of the GNU General Public License as published by
114 !! the Free Software Foundation, either version 3 of the License, or
115 !! (at your option) any later version.
116 !!
117 !! SICOPOLIS is distributed in the hope that it will be useful,
118 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
119 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
120 !! GNU General Public License for more details.
121 !!
122 !! You should have received a copy of the GNU General Public License
123 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
124 !<
125 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
126 
127 !-------- Include run specification header --------
128 
129 #include "sico_specs.h"
130 
131 !-------- Include header for the Library of Iterative Solvers Lis
132 ! (only if required) --------
133 
134 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
135 #include "lisf.h"
136 #endif
137 
138 !-------- Include modules --------
139 
140 #include "subroutines/general/sico_types_m.F90"
141 #include "subroutines/general/sico_variables_m.F90"
142 
143 #if (defined(ANT))
144 #include "subroutines/ant/sico_vars_m.F90"
145 #elif (defined(ASF))
146 #include "subroutines/asf/sico_vars_m.F90"
147 #elif (defined(EMTP2SGE))
148 #include "subroutines/emtp2sge/sico_vars_m.F90"
149 #elif (defined(GRL))
150 #include "subroutines/grl/sico_vars_m.F90"
151 #elif (defined(NHEM))
152 #include "subroutines/nhem/sico_vars_m.F90"
153 #elif (defined(NMARS))
154 #include "subroutines/nmars/sico_vars_m.F90"
155 #elif (defined(SCAND))
156 #include "subroutines/scand/sico_vars_m.F90"
157 #elif (defined(SMARS))
158 #include "subroutines/smars/sico_vars_m.F90"
159 #elif (defined(TIBET))
160 #include "subroutines/tibet/sico_vars_m.F90"
161 #elif (defined(XYZ))
162 #include "subroutines/xyz/sico_vars_m.F90"
163 #endif
164 
165 #include "subroutines/general/ice_material_properties_m.F90"
166 #include "subroutines/general/stereo_proj_m.F90"
167 #include "subroutines/general/metric_m.F90"
168 #include "subroutines/general/sico_maths_m.F90"
169 #include "subroutines/general/compare_float_m.F90"
170 
171 #if (NETCDF==2)
172 #include "subroutines/general/nc_check_m.F90"
173 #endif
174 
175 #include "subroutines/general/read_m.F90"
176 #include "subroutines/general/init_temp_water_age_m.F90"
177 
178 #include "subroutines/general/mask_update_sea_level_m.F90"
179 #include "subroutines/general/pdd_m.F90"
180 
181 #if ((MARGIN==2) \
182  && (marine_ice_formation==2) \
183  && (marine_ice_calving==9))
184 #include "subroutines/general/calving_underwater_ice_m.F90"
185 #endif
186 
187 #if (defined(GRL))
188 #if (DISC>0)
189 #include "subroutines/grl/discharge_workers_m.F90"
190 #endif
191 #endif
192 
193 #include "subroutines/general/calc_vxy_m.F90"
194 #include "subroutines/general/calc_vz_m.F90"
195 #include "subroutines/general/calc_dxyz_m.F90"
196 
197 #include "subroutines/general/calc_gia_m.F90"
198 #include "subroutines/general/topograd_m.F90"
199 #include "subroutines/general/calc_thk_m.F90"
200 
201 #include "subroutines/general/enth_temp_omega_m.F90"
202 
203 #if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1)
204 #include "subroutines/general/calc_temp_m.F90"
205 #elif (CALCMOD==2 || CALCMOD==3)
206 #include "subroutines/general/calc_temp_enth_m.F90"
207 #endif
208 
209 #include "subroutines/general/calc_enhance_m.F90"
210 
211 #if (BASAL_HYDROLOGY==1)
212 #include "subroutines/non-foss/hydro_m.F90"
213 #endif
214 
215 #if (defined(NMARS) || defined(SMARS))
216 #include "subroutines/general/mars_instemp_m.f90"
217 #endif
218 
219 #include "subroutines/general/calc_temp_melt_bas_m.F90"
220 #include "subroutines/general/calc_bas_melt_m.F90"
221 #include "subroutines/general/calc_thk_water_bas_m.F90"
222 
223 #include "subroutines/general/output_m.F90"
224 
225 #if (defined(ANT))
226 #include "subroutines/ant/boundary_m.F90"
227 #elif (defined(ASF))
228 #include "subroutines/asf/boundary_m.F90"
229 #elif (defined(EMTP2SGE))
230 #include "subroutines/emtp2sge/boundary_m.F90"
231 #elif (defined(GRL))
232 #include "subroutines/grl/boundary_m.F90"
233 #elif (defined(NHEM))
234 #include "subroutines/nhem/boundary_m.F90"
235 #elif (defined(NMARS))
236 #include "subroutines/nmars/boundary_m.F90"
237 #elif (defined(SCAND))
238 #include "subroutines/scand/boundary_m.F90"
239 #elif (defined(SMARS))
240 #include "subroutines/smars/boundary_m.F90"
241 #elif (defined(TIBET))
242 #include "subroutines/tibet/boundary_m.F90"
243 #elif (defined(XYZ))
244 #include "subroutines/xyz/boundary_m.F90"
245 #endif
246 
247 #if (defined(ANT))
248 #include "subroutines/ant/sico_init_m.F90"
249 #elif (defined(ASF))
250 #include "subroutines/asf/sico_init_m.F90"
251 #elif (defined(EMTP2SGE))
252 #include "subroutines/emtp2sge/sico_init_m.F90"
253 #elif (defined(GRL))
254 #include "subroutines/grl/sico_init_m.F90"
255 #elif (defined(NHEM))
256 #include "subroutines/nhem/sico_init_m.F90"
257 #elif (defined(NMARS))
258 #include "subroutines/nmars/sico_init_m.F90"
259 #elif (defined(SCAND))
260 #include "subroutines/scand/sico_init_m.F90"
261 #elif (defined(SMARS))
262 #include "subroutines/smars/sico_init_m.F90"
263 #elif (defined(TIBET))
264 #include "subroutines/tibet/sico_init_m.F90"
265 #elif (defined(XYZ))
266 #include "subroutines/xyz/sico_init_m.F90"
267 #endif
268 
269 #include "subroutines/general/sico_main_loop_m.F90"
270 #include "subroutines/general/sico_end_m.F90"
271 
272 !-------------------------------------------------------------------------------
273 !> Main program of SICOPOLIS.
274 !<------------------------------------------------------------------------------
275 program sicopolis
276 
277 !-------- Description of variables declared locally --------
278 
279 ! i : Index for coordinate xi
280 ! j : Index for coordinate eta
281 ! kc : Index for coordinate zeta_c
282 ! kt : Index for coordinate zeta_t
283 ! kr : Index for coordinate zeta_r
284 ! (in the bedrock)
285 ! ios : IOSTAT variable for files
286 ! itercount : Counter for the time steps
287 ! ndat2d/3d : Counters for the time-slice files
288 ! n_output : For OUTPUT==2 number of time-slice files
289 ! to be produced
290 ! dH_t_smooth(j,i) : Amount of smoothing of H_t
291 ! delta_ts : Time-dependent surface-temperature variation
292 ! glac_index : Time-dependent glacial index
293 ! forcing_flag : 1 - forcing by delta_ts, 2 - forcing by glac_index
294 ! z_sl : Sea level
295 ! dzsl_dtau : Derivative of z_sl by tau (time)
296 ! precip_mam_present(j,i) : Measured present spring precipitation
297 ! precip_jja_present(j,i) : Measured present summer precipitation
298 ! precip_son_present(j,i) : Measured present autumn precipitation
299 ! precip_djf_present(j,i) : Measured present winter precipitation
300 ! temp_mam_present(j,i) : Present spring surface temperature
301 ! from data
302 ! temp_jja_present(j,i) : Present summer surface temperature
303 ! from data
304 ! temp_son_present(j,i) : Present autumn surface temperature
305 ! from data
306 ! temp_djf_present(j,i) : Present winter surface temperature
307 ! from data
308 ! mean_accum : Mean present accumulation over land
309 ! time : Current time of simulation
310 ! time_init : Initial time of simulation
311 ! time_end : Final time of simulation
312 ! dtime : Time step for computation of velocity and
313 ! topography
314 ! dtime_temp : Time step for computation of temperature,
315 ! water content and age of the ice
316 ! dtime_wss : Time step for computation of isostatic steady-state
317 ! displacement of the lithosphere (ELRA model)
318 ! dtime_out : Time step for output of time-slice files
319 ! dtime_ser : Time step for writing of data in time-series
320 ! file
321 ! time_output(n) : For OUTPUT==2 specified times for output of
322 ! time-slice files
323 ! (.)0 : Quantity (.) before its conversion to MKS units
324 ! dxi : Grid spacing in x-direction
325 ! deta : Grid spacing in y-direction
326 ! dzeta_c : Grid spacing in z-direction in the upper (kc) ice domain
327 ! (in sigma-coordinate zeta_c)
328 ! dzeta_t : Grid spacing in z-direction in the lower (kt) ice domain
329 ! (in sigma-coordinate zeta_t)
330 ! dzeta_r : Grid spacing in z-direction in the bedrock (kr) domain
331 ! (in sigma-coordinate zeta_r)
332 ! runname : Name of simulation
333 ! anfdatname : Name of initial-value file
334 ! (only for ANF_DAT==3)
335 ! datname : Auxiliary string for generation of file names
336 
337 !-------- Declaration of variables --------
338 
339 use sico_types_m
341 use sico_vars_m
342 
343 use sico_init_m
345 use sico_end_m
346 
347 implicit none
348 
349 integer(i4b) :: ndat2d, ndat3d
350 integer(i4b) :: n_output
351 integer(i4b), dimension((IMAX+1)*(JMAX+1)) :: ii, jj
352 integer(i4b), dimension(0:JMAX,0:IMAX) :: nn
353 real(dp) :: delta_ts, glac_index
354 real(dp) :: mean_accum
355 real(dp) :: dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser
356 real(dp) :: time, time_init, time_end
357 real(dp), dimension(100) :: time_output
358 real(dp) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
359 real(dp) :: z_sl, dzsl_dtau, z_mar
360 character(len=100) :: runname
361 
362 !-------- Initialisations --------
363 
364 call sico_init(delta_ts, glac_index, &
365  mean_accum, &
366  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
367  time, time_init, time_end, time_output, &
368  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
369  z_sl, dzsl_dtau, z_mar, &
370  ii, jj, nn, &
371  ndat2d, ndat3d, n_output, &
372  runname)
373 
374 !-------- Main loop --------
375 
376 call sico_main_loop(delta_ts, glac_index, &
377  mean_accum, &
378  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
379  time, time_init, time_end, time_output, &
380  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
381  z_sl, dzsl_dtau, z_mar, &
382  ii, jj, nn, &
383  ndat2d, ndat3d, n_output, &
384  runname)
385 
386 !-------- Endings --------
387 
388 call sico_end()
389 
390 !-------- End of program --------
391 
392 end program sicopolis
393 
394 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
395 ! End of sicopolis.F90
396 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
397 !
subroutine sico_init(delta_ts, glac_index, mean_accum, dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, time, time_init, time_end, time_output, dxi, deta, dzeta_c, dzeta_t, dzeta_r, z_sl, dzsl_dtau, z_mar, ii, jj, nn, ndat2d, ndat3d, n_output, runname)
Main routine of sico_init_m: Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:59
Main loop of SICOPOLIS.
program sicopolis
Main program of SICOPOLIS.
Definition: sicopolis.F90:275
subroutine sico_main_loop(delta_ts, glac_index, mean_accum, dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, time, time_init, time_end, time_output, dxi, deta, dzeta_c, dzeta_t, dzeta_r, z_sl, dzsl_dtau, z_mar, ii, jj, nn, ndat2d, ndat3d, n_output, runname)
Main routine of sico_main_loop_m: Main loop of SICOPOLIS.
subroutine sico_end()
Main routine of sico_end_m: Ending of SICOPOLIS.
Definition: sico_end_m.F90:51
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:35
Declarations of kind types for SICOPOLIS.
Ending of SICOPOLIS.
Definition: sico_end_m.F90:35
Declarations of global variables for SICOPOLIS.