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