SICOPOLIS V5-dev  Revision 1368
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 '5-dev'
8 #define DATE '2018-10-31'
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-2018 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-2018 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 !@ begin openad_extract @
127 
128 !-------- Include run specification header --------
129 
130 #include RUN_SPECS_HEADER
131 
132 !-------- Include header for the Library of Iterative Solvers Lis
133 ! (only if required) --------
134 
135 #ifndef ALLOW_OPENAD /* Normal */
136 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
137 #include "lisf.h"
138 #endif
139 #endif /* Normal */
140 
141 !-------- Include modules --------
142 
143 #include "subroutines/general/sico_types_m.F90"
144 #include "subroutines/general/sico_variables_m.F90"
145 
146 #if (defined(ANT))
147 #include "subroutines/ant/sico_vars_m.F90"
148 #elif (defined(ASF))
149 #include "subroutines/asf/sico_vars_m.F90"
150 #elif (defined(EMTP2SGE))
151 #include "subroutines/emtp2sge/sico_vars_m.F90"
152 #elif (defined(GRL))
153 #include "subroutines/grl/sico_vars_m.F90"
154 #elif (defined(NHEM))
155 #include "subroutines/nhem/sico_vars_m.F90"
156 #elif (defined(NMARS))
157 #include "subroutines/nmars/sico_vars_m.F90"
158 #elif (defined(SCAND))
159 #include "subroutines/scand/sico_vars_m.F90"
160 #elif (defined(SMARS))
161 #include "subroutines/smars/sico_vars_m.F90"
162 #elif (defined(TIBET))
163 #include "subroutines/tibet/sico_vars_m.F90"
164 #elif (defined(XYZ))
165 #include "subroutines/xyz/sico_vars_m.F90"
166 #endif
167 
168 #include "subroutines/general/error_m.F90"
169 
170 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD)) /* OpenAD */
171 #include "subroutines/openad/ctrl_m.F90"
172 #endif /* OpenAD */
173 
174 #include "subroutines/general/ice_material_properties_m.F90"
175 #include "subroutines/general/stereo_proj_m.F90"
176 #include "subroutines/general/metric_m.F90"
177 
178 #if (!defined(ALLOW_OPENAD) || defined(ALLOW_OPENAD_DIFFERENTIATE)) /* Normal */
179 #include "subroutines/general/sico_maths_m.F90"
180 #endif /* Normal */
181 
182 #if (!defined(ALLOW_OPENAD)) /* Normal */
183 #include "subroutines/general/compare_float_m.F90"
184 #endif /* Normal */
185 
186 #if (NETCDF==2)
187 #include "subroutines/general/nc_check_m.F90"
188 #endif
189 
190 #include "subroutines/general/read_m.F90"
191 
192 #include "subroutines/general/mask_update_sea_level_m.F90"
193 #include "subroutines/general/flag_update_gf_gl_cf_m.F90"
194 #include "subroutines/general/pdd_m.F90"
195 
196 #if ((MARGIN==2) \
197  && (marine_ice_formation==2) \
198  && (marine_ice_calving==9))
199 #include "subroutines/general/calving_underwater_ice_m.F90"
200 #endif
201 
202 #if (defined(GRL))
203 #if (DISC>0)
204 #include "subroutines/grl/discharge_workers_m.F90"
205 #endif
206 #endif
207 
208 #include "subroutines/general/calc_vxy_m.F90"
209 #include "subroutines/general/calc_vz_m.F90"
210 #include "subroutines/general/calc_dxyz_m.F90"
211 
212 #include "subroutines/general/calc_gia_m.F90"
213 #include "subroutines/general/topograd_m.F90"
214 #include "subroutines/general/calc_thk_m.F90"
215 
216 #include "subroutines/general/enth_temp_omega_m.F90"
217 
218 #if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1)
219 #include "subroutines/general/calc_temp_m.F90"
220 #elif (CALCMOD==2 || CALCMOD==3)
221 #include "subroutines/general/calc_temp_enth_m.F90"
222 #endif
223 
224 #include "subroutines/general/calc_enhance_m.F90"
225 
226 #if (BASAL_HYDROLOGY==1)
227 #include "subroutines/non-foss/hydro_m.F90"
228 #endif
229 
230 #if (defined(NMARS) || defined(SMARS))
231 #include "subroutines/general/mars_instemp_m.f90"
232 #endif
233 
234 #include "subroutines/general/calc_temp_melt_bas_m.F90"
235 #include "subroutines/general/calc_bas_melt_m.F90"
236 #include "subroutines/general/calc_thk_water_bas_m.F90"
237 
238 #ifndef ALLOW_OPENAD /* Normal */
239 #include "subroutines/general/output_m.F90"
240 #endif /* Normal */
241 
242 #if (defined(ANT))
243 #include "subroutines/ant/boundary_m.F90"
244 #elif (defined(ASF))
245 #include "subroutines/asf/boundary_m.F90"
246 #elif (defined(EMTP2SGE))
247 #include "subroutines/emtp2sge/boundary_m.F90"
248 #elif (defined(GRL))
249 #include "subroutines/grl/boundary_m.F90"
250 #elif (defined(NHEM))
251 #include "subroutines/nhem/boundary_m.F90"
252 #elif (defined(NMARS))
253 #include "subroutines/nmars/boundary_m.F90"
254 #elif (defined(SCAND))
255 #include "subroutines/scand/boundary_m.F90"
256 #elif (defined(SMARS))
257 #include "subroutines/smars/boundary_m.F90"
258 #elif (defined(TIBET))
259 #include "subroutines/tibet/boundary_m.F90"
260 #elif (defined(XYZ))
261 #include "subroutines/xyz/boundary_m.F90"
262 #endif
263 
264 !@ end openad_extract @
265 
266 #include "subroutines/general/init_temp_water_age_m.F90"
267 #if (defined(ANT))
268 #include "subroutines/ant/sico_init_m.F90"
269 #elif (defined(ASF))
270 #include "subroutines/asf/sico_init_m.F90"
271 #elif (defined(EMTP2SGE))
272 #include "subroutines/emtp2sge/sico_init_m.F90"
273 #elif (defined(GRL))
274 #include "subroutines/grl/sico_init_m.F90"
275 #elif (defined(NHEM))
276 #include "subroutines/nhem/sico_init_m.F90"
277 #elif (defined(NMARS))
278 #include "subroutines/nmars/sico_init_m.F90"
279 #elif (defined(SCAND))
280 #include "subroutines/scand/sico_init_m.F90"
281 #elif (defined(SMARS))
282 #include "subroutines/smars/sico_init_m.F90"
283 #elif (defined(TIBET))
284 #include "subroutines/tibet/sico_init_m.F90"
285 #elif (defined(XYZ))
286 #include "subroutines/xyz/sico_init_m.F90"
287 #endif
288 
289 !@ begin openad_extract @
290 
291 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD)) /* OpenAD */
292 #include "subroutines/openad/sico_main_loop_iter_m.F90"
293 #include "subroutines/openad/sico_main_loop_wrapper_m.F90"
294 #endif /* OpenAD */
295 
296 #include "subroutines/general/sico_main_loop_m.F90"
297 #include "subroutines/general/sico_end_m.F90"
298 
299 !@ end openad_extract @
300 
301 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD)) /* OpenAD */
302 #include "subroutines/openad/openad_m.F90"
303 #endif /* OpenAD */
304 !-------------------------------------------------------------------------------
305 !> Main program of SICOPOLIS.
306 !<------------------------------------------------------------------------------
307 program sicopolis
308 
309 !-------- Description of variables declared locally --------
310 
311 ! i : Index for coordinate xi
312 ! j : Index for coordinate eta
313 ! kc : Index for coordinate zeta_c
314 ! kt : Index for coordinate zeta_t
315 ! kr : Index for coordinate zeta_r
316 ! (in the bedrock)
317 ! ios : IOSTAT variable for files
318 ! itercount : Counter for the time steps
319 ! ndat2d/3d : Counters for the time-slice files
320 ! n_output : For OUTPUT==2 number of time-slice files
321 ! to be produced
322 ! dH_t_smooth(j,i) : Amount of smoothing of H_t
323 ! delta_ts : Time-dependent surface-temperature variation
324 ! glac_index : Time-dependent glacial index
325 ! forcing_flag : 1 - forcing by delta_ts, 2 - forcing by glac_index
326 ! z_sl : Sea level
327 ! dzsl_dtau : Derivative of z_sl by tau (time)
328 ! precip_mam_present(j,i) : Measured present spring precipitation
329 ! precip_jja_present(j,i) : Measured present summer precipitation
330 ! precip_son_present(j,i) : Measured present autumn precipitation
331 ! precip_djf_present(j,i) : Measured present winter precipitation
332 ! temp_mam_present(j,i) : Present spring surface temperature
333 ! from data
334 ! temp_jja_present(j,i) : Present summer surface temperature
335 ! from data
336 ! temp_son_present(j,i) : Present autumn surface temperature
337 ! from data
338 ! temp_djf_present(j,i) : Present winter surface temperature
339 ! from data
340 ! mean_accum : Mean present accumulation over land
341 ! time : Current time of simulation
342 ! time_init : Initial time of simulation
343 ! time_end : Final time of simulation
344 ! dtime : Time step for computation of velocity and
345 ! topography
346 ! dtime_temp : Time step for computation of temperature,
347 ! water content and age of the ice
348 ! dtime_wss : Time step for computation of isostatic steady-state
349 ! displacement of the lithosphere (ELRA model)
350 ! dtime_out : Time step for output of time-slice files
351 ! dtime_ser : Time step for writing of data in time-series
352 ! file
353 ! time_output(n) : For OUTPUT==2 specified times for output of
354 ! time-slice files
355 ! (.)0 : Quantity (.) before its conversion to MKS units
356 ! dxi : Grid spacing in x-direction
357 ! deta : Grid spacing in y-direction
358 ! dzeta_c : Grid spacing in z-direction in the upper (kc) ice domain
359 ! (in sigma-coordinate zeta_c)
360 ! dzeta_t : Grid spacing in z-direction in the lower (kt) ice domain
361 ! (in sigma-coordinate zeta_t)
362 ! dzeta_r : Grid spacing in z-direction in the bedrock (kr) domain
363 ! (in sigma-coordinate zeta_r)
364 ! runname : Name of simulation
365 ! anfdatname : Name of initial-value file
366 ! (only for ANF_DAT==3)
367 ! datname : Auxiliary string for generation of file names
368 
369 !-------- Declaration of variables --------
370 
371 !@ begin openad_extract @
372 !openad begin subroutine sicopolis_openad
373 use sico_types_m
375 use sico_vars_m
376 !@ end openad_extract @
377 
378 use sico_init_m
379 !@ begin openad_extract @
381 use sico_end_m
382 !@ end openad_extract @
383 
384 #ifdef ALLOW_GRDCHK /* OpenAD */
385 use openad_m, only: adjoint_master, grdchk_main
386 #endif /* OpenAD */
387 
388 #ifdef ALLOW_OPENAD /* OpenAD */
389 use oad_tape
390 use oad_rev
391 use openad_m, only: adjoint_master
392 #endif /* OpenAD */
393 
394 implicit none
395 
396 !@ begin openad_extract @
397 integer(i4b) :: ndat2d, ndat3d
398 integer(i4b) :: n_output
399 integer(i4b), dimension((IMAX+1)*(JMAX+1)) :: ii, jj
400 integer(i4b), dimension(0:JMAX,0:IMAX) :: nn
401 real(dp) :: delta_ts, glac_index
402 real(dp) :: mean_accum
403 real(dp) :: dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser
404 real(dp) :: time, time_init, time_end
405 real(dp), dimension(100) :: time_output
406 real(dp) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
407 real(dp) :: z_sl, dzsl_dtau, z_mar
408 character(len=100) :: runname
409 !openad sicopolis_independents_cost
410 !@ end openad_extract @
411 
412 #ifdef ALLOW_OPENAD /* OpenAD */
413 integer(i4b) :: mode
414 character(len=32) :: arg
415 logical :: ISPLAIN, ISTAPE, ISADJOINT
416 #endif /* OpenAD */
417 
418 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
419 !-------- Initialisations --------
420 
421 call sico_init(delta_ts, glac_index, &
422  mean_accum, &
423  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
424  time, time_init, time_end, time_output, &
425  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
426  z_sl, dzsl_dtau, z_mar, &
427  ii, jj, nn, &
428  ndat2d, ndat3d, n_output, &
429  runname)
430 
431 !@ begin openad_extract @
432 !-------- Main loop --------
433 
434 call sico_main_loop(delta_ts, glac_index, &
435  mean_accum, &
436  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
437  time, time_init, time_end, time_output, &
438  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
439  z_sl, dzsl_dtau, z_mar, &
440  ii, jj, nn, &
441  ndat2d, ndat3d, n_output, &
442  runname)
443 
444 !openad end subroutine sicopolis_openad
445 
446 !@ end openad_extract @
447 !-------- Endings --------
448 
449 call sico_end()
450 
451 !-------- End of program --------
452 
453 #else /* OpenAD */
454 
455 call adjoint_master
456 
457 #endif /* Normal vs. OpenAD */
458 
459 end program sicopolis
460 
461 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
462 ! End of sicopolis.F90
463 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
464 !
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:60
Main loop of SICOPOLIS.
program sicopolis
Main program of SICOPOLIS.
Definition: sicopolis.F90:307
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
Module for all openad-related subroutines.
Definition: openad_m.F90:37
Declarations of kind types for SICOPOLIS.
Ending of SICOPOLIS.
Definition: sico_end_m.F90:35
Declarations of global variables for SICOPOLIS.
subroutine, public adjoint_master
Adjoint master is the main tool by which sicopolis.F90 invokes the adjoint code. Its job is to figure...
Definition: openad_m.F90:63