SICOPOLIS V5-dev  Revision 1420
sico_main_loop_iter_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : s i c o _ m a i n _ l o o p _ i t e r _ m
4 !
5 !> @file
6 !!
7 !! A catch-all module for openad-related subroutines.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2017-2018 ...
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
33 
34  implicit none
35 
36  public
37 
38 contains
39 
40 !-------------------------------------------------------------------------------
41 !> The interior of subroutine sico_main_loop. If new subroutines / modules are
42 !! introduced into the main calculation, they must also be inserted here.
43 !<------------------------------------------------------------------------------
44  subroutine sico_main_loop_iter(delta_ts, glac_index, &
45  mean_accum, &
46  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
47  time, time_init, time_end, time_output, &
48  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
49  z_sl, dzsl_dtau, z_mar, &
50  ndat2d, ndat3d, n_output, &
51  runname, &
52  itercount,iter_temp,iter_wss,iter_ser,&
53  iter_out,iter_output)
54 
55  use sico_types_m
57  use sico_vars_m
58  use boundary_m
59 #if (defined(GRL) && DISC>0)
61 #endif
62 #if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1)
63  use calc_temp_m
64 #elif (CALCMOD==2 || CALCMOD==3)
66 #endif
67 
68  use calc_enhance_m
70  use calc_vxy_m
71  use calc_vz_m
72  use calc_dxyz_m
73  use calc_gia_m
74  use calc_thk_m
76  use calc_bas_melt_m
78  use ctrl_m, only: myfloor
79 
80  implicit none
81 
82  integer(i4b), intent(in) :: n_output
83  real(dp), intent(in) :: mean_accum
84  real(dp), intent(in) :: dtime, dtime_temp, dtime_wss, &
85  dtime_out, dtime_ser
86  real(dp), intent(in) :: time_init, time_end, time_output(100)
87  real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
88  character(len=100), intent(in) :: runname
89 
90  integer(i4b), intent(inout) :: ndat2d, ndat3d
91  real(dp), intent(inout) :: delta_ts, glac_index
92  real(dp), intent(inout) :: time
93  real(dp), intent(inout) :: z_sl, dzsl_dtau, z_mar
94 
95  integer(i4b) :: i, j, kc, kt, kr, n
96  integer(i4b), intent(inout) :: itercount, iter_temp, iter_wss
97  integer(i4b), intent(inout) :: iter_ser, iter_out, iter_output(100)
98  real(dp) :: dtime_temp_inv
99  logical :: flag_3d_output
100 
101  real(dp) :: tmp, valmod
102  integer(i4b) :: valmodint
103 
104  write(unit=6, fmt='(2x,i0)') itercount
105 
106  !-------- Update of time --------
107 
108  time = time_init + real(itercount,dp)*dtime
109 
110  !-------- Save old mask --------
111 
112  maske_old = maske
113 
114  !-------- Boundary conditions --------
115 
116  call boundary(time, dtime, dxi, deta, delta_ts, glac_index, &
117  z_sl, dzsl_dtau, z_mar)
118 
119  !-------- Temperature, water content, age, flow enhancement factor --------
120  if ((itercount - (iter_temp*int(real(itercount)/real(iter_temp)))) == 0) then
121 
122  write(unit=6, fmt='(10x,a)') 'Computation of T'
123 
124  ! ------ Temperature, water content, age
125 
126 #if (CALCMOD==1)
127  call calc_temp_poly(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
128  ! polythermal scheme (POLY)
129 
130 #elif (CALCMOD==0)
131  call calc_temp_cold(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
132  ! cold-ice scheme (COLD)
133 
134 #elif (CALCMOD==2 || CALCMOD==3)
135  call calc_temp_enth(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
136  ! enthalpy scheme (ENTC or ENTM)
137 
138 #elif (CALCMOD==-1)
139  call calc_temp_const() ! isothermal scheme (ISOT)
140 
141 #else
142  stop ' >>> sico_main_loop: Parameter CALCMOD must be between -1 and 3!'
143 #endif
144 
145  ! ------ Time derivative of H_t (further time derivatives are
146  ! computed in subroutine calc_thk_xxx)
147 
148  dtime_temp_inv = 1.0_dp/dtime_temp
149 
150  dh_t_dtau = (h_t_neu - h_t)*dtime_temp_inv
151 
152  ! ------ New values --> old values
153 
154  n_cts = n_cts_neu
156  zm = zm_neu
157  h_c = h_c_neu
158  h_t = h_t_neu
160  age_c = age_c_neu
162  age_t = age_t_neu
164 
165 #if (CALCMOD==2 || CALCMOD==3)
169 #endif
170 
171  ! ------ Flow enhancement factor
172 
173 #if (ENHMOD==1)
174  call calc_enhance_1()
175 #elif (ENHMOD==2)
176  call calc_enhance_2()
177 #elif (ENHMOD==3)
178  call calc_enhance_3(time)
179 #elif (ENHMOD==4)
180  call calc_enhance_4()
181 #elif (ENHMOD==5)
182  call calc_enhance_5()
183 #else
184  stop ' >>> sico_main_loop: Parameter ENHMOD must be between 1 and 5!'
185 #endif
186 
187  end if
188  ! End of computation of temperature, water content, age and
189  ! enhancement factor
190 
191  !-------- Velocity --------
192 
193  call flag_update_gf_gl_cf()
194 
195 #if (DYNAMICS==1 || DYNAMICS==2)
196 
197  call calc_vxy_b_sia(time, z_sl)
198  call calc_vxy_sia(dzeta_c, dzeta_t)
199 
200 #if (MARGIN==3 || DYNAMICS==2)
201  call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t)
202 #endif
203 
204  call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
205 
206 #if (MARGIN==3)
207  call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
208 #endif
209 
210 #elif (DYNAMICS==0)
211 
212  call calc_vxy_static()
213  call calc_vz_static()
214 
215 #else
216  stop ' >>> sico_main_loop: DYNAMICS must be either 0, 1 or 2!'
217 #endif
218 
219  call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
220 
221  !-------- Glacial isostatic adjustment and ice topography --------
222 
223  call calc_gia(time, z_sl, dtime, dxi, deta, itercount, iter_wss)
224 
225  call calc_thk_init()
226 
227 #if ((MARGIN==3 || DYNAMICS==2) && (CALCTHK==1 || CALCTHK==2 || CALCTHK==3))
228  write(6, fmt='(a)') ' >>> sico_main_loop:'
229  write(6, fmt='(a)') ' Non-SIA dynamics combined with'
230  write(6, fmt='(a)') ' SIA ice thickness solver!'
231  stop
232 #endif
233 
234 #if (CALCTHK==1)
235  call calc_thk_sia_expl(time, dtime, dxi, deta, z_mar)
236 #elif (CALCTHK==2 || CALCTHK==3)
237  call calc_thk_sia_impl(time, dtime, dxi, deta, z_mar, mean_accum)
238 #elif (CALCTHK==4)
239  call calc_thk_expl(time, dtime, dxi, deta, z_mar)
240 #elif (CALCTHK==5 || CALCTHK==6)
241  call calc_thk_impl(time, dtime, dxi, deta, z_mar, mean_accum)
242 #else
243  stop ' >>> sico_main_loop: Parameter CALCTHK must be between 1 and 6!'
244 #endif
245 
246 #if (MARGIN==3) /* coupled SIA/SSA or SIA/SStA/SSA dynamics */
247  call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 3_i1b)
248 #elif (DYNAMICS==2) /* hybrid SIA/SStA dynamics */
249  call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 2_i1b)
250 #else /* SIA-only dynamics */
251 #if (CALCTHK==1 || CALCTHK==2 || CALCTHK==3)
252  call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 1_i1b)
253 #elif (CALCTHK==4 || CALCTHK==5 || CALCTHK==6)
254  call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 2_i1b)
255 #endif
256 #endif
257 
258  ! ------ New values --> old values
259 
260  call flag_update_gf_gl_cf()
261 
262  zs = zs_neu
263  zm = zm_neu
264  zb = zb_neu
265  zl = zl_neu
266  h_c= h_c_neu
267  h_t= h_t_neu
268 
269  !-------- Melting temperature --------
270 
271  call calc_temp_melt()
272 
273  !-------- Basal temperature --------
274 
275  call calc_temp_bas()
276 
277  !-------- Basal melting rate --------
278 
279  call calc_qbm(time, z_sl, dzeta_c, dzeta_r)
280 
281  !-------- Effective thickness of subglacial water --------
282 
283  call calc_thk_water_bas(z_sl)
284 
285  !-------- Data output --------
286  ! We do none of the original data output in adjoint mode
287 
288  end subroutine sico_main_loop_iter
289 
290 end module sico_main_loop_iter_m
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
subroutine, public calc_thk_water_bas(z_sl)
Main subroutine of calc_thk_water_bas_m: Computation of the thickness of the water column under the i...
Computation of the ice thickness.
Definition: calc_thk_m.F90:35
integer(i1b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_thk_init()
Initialisations for the ice thickness computation.
Definition: calc_thk_m.F90:62
Update of the flags for the land-terminating grounded front, marine-terminating grounded front...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
subroutine, public calc_qbm(time, z_sl, dzeta_c, dzeta_r)
Computation of the basal melting rate Q_bm. Summation of Q_bm and Q_tld (water drainage rate from the...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public flag_update_gf_gl_cf()
Main subroutine of flag_update_gf_gl_cf_m: Update of the flags for the land-terminating grounded fron...
Definition: ctrl_m.F90:9
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t
age_t(kt,j,i): Age in the lower (kt) ice domain
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r
temp_r(kr,j,i): Temperature in the bedrock
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_enhance_2()
Computation of the flow enhancement factor. Case ENHMOD==2: two different values depending on age for...
subroutine, public calc_vz_floating(z_sl, dxi, deta, dzeta_c)
Computation of the vertical velocity vz for floating ice.
Definition: calc_vz_m.F90:259
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t
enth_t(kt,j,i): Enthalpy in the lower (kt) ice domain
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
Computation of the thickness of the water column under the ice base.
integer(i1b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
subroutine, public calc_temp_cold(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
Computation of temperature and age in cold-ice mode.
subroutine, public calc_thk_sia_expl(time, dtime, dxi, deta, z_mar)
Explicit solver for the diffusive SIA ice surface equation.
Definition: calc_thk_m.F90:156
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
Definition: calc_dxyz_m.F90:37
real(dp), dimension(0:jmax, 0:imax) h_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_temp_const()
Isothermal mode: Setting of the temperature and age to constant values.
subroutine, public calc_enhance_5()
Computation of the flow enhancement factor. Case ENHMOD==5: minimal anisotropic enhancement factor fo...
subroutine, public calc_thk_sia_impl(time, dtime, dxi, deta, z_mar, mean_accum)
Over-implicit solver for the diffusive SIA ice surface equation.
Definition: calc_thk_m.F90:234
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
subroutine, public calc_thk_expl(time, dtime, dxi, deta, z_mar)
Explicit solver for the general ice thickness equation.
Definition: calc_thk_m.F90:592
subroutine, public calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
real(dp), dimension(0:jmax, 0:imax) zl_neu
(.)_neu: New value of quantity (.) computed during an integration step
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
Computation of the flow enhancement factor.
Computation of the horizontal velocity vx, vy.
Definition: calc_vxy_m.F90:36
subroutine, public calc_temp_melt()
Computation of the melting temperatures.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz for grounded ice.
Definition: calc_vz_m.F90:52
subroutine, public calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, n_calc_thk_mask_update_aux)
Update of the ice-land-ocean mask etc.
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
subroutine sico_main_loop_iter(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, ndat2d, ndat3d, n_output, runname, itercount, iter_temp, iter_wss, iter_ser, iter_out, iter_output)
The interior of subroutine sico_main_loop. If new subroutines / modules are introduced into the main ...
integer(i4b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
Computation of temperature, water content and age.
Definition: calc_temp_m.F90:35
subroutine, public calc_temp_bas()
Computation of the basal temperatures.
Computation of the melting and basal temperatures.
integer(i1b), dimension(0:jmax, 0:imax) n_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
Computation of the vertical velocity vz.
Definition: calc_vz_m.F90:35
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_gia(time, z_sl, dtime, dxi, deta, itercount, iter_wss)
Main subroutine of calc_gia_m: Computation of the glacial isostatic adjustment of the lithosphere sur...
Definition: calc_gia_m.F90:55
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) zs_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_vz_static()
Computation of the vertical velocity vz for static ice.
Definition: calc_vz_m.F90:374
integer(i4b), dimension(0:jmax, 0:imax) kc_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) h_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_enhance_1()
Computation of the flow enhancement factor. Case ENHMOD==1: constant for grounded ice...
real(dp), dimension(0:jmax, 0:imax) zm
zm(j,i): Coordinate z of the bottom of the upper (kc) ice domain = top of the lower (kt) ice domain (...
integer(i1b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
Computation of temperature, water content and age with the enthalpy method.
subroutine, public calc_vxy_sia(dzeta_c, dzeta_t)
Computation of the shear stresses txz, tyz, the effective shear stress sigma, the depth-averaged flui...
Definition: calc_vxy_m.F90:573
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Main routine of boundary_m: Computation of the surface temperature (must be less than 0 deg C!) and o...
Definition: boundary_m.F90:61
subroutine, public calc_temp_enth(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
Main subroutine of calc_temp_enth_m: Computation of temperature, water content and age with the entha...
subroutine, public calc_enhance_3(time)
Computation of the flow enhancement factor. Case ENHMOD==3: two different values depending on time of...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
subroutine, public calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
Main subroutine of calc_dxyz_m: Computation of all components of the strain-rate tensor, the full effective strain rate and the shear fraction.
Definition: calc_dxyz_m.F90:57
subroutine, public calc_enhance_4()
Computation of the flow enhancement factor. Case ENHMOD==4: minimal anisotropic enhancement factor fo...
Computation of the basal melting rate.
real(dp), dimension(0:jmax, 0:imax) zb_neu
(.)_neu: New value of quantity (.) computed during an integration step
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
Definition: boundary_m.F90:37
real(dp), dimension(0:jmax, 0:imax) zm_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_thk_impl(time, dtime, dxi, deta, z_mar, mean_accum)
Over-implicit solver for the general ice thickness equation.
Definition: calc_thk_m.F90:734
Ice discharge parameterization for the Greenland ice sheet.
integer, parameter dp
Double-precision reals.
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_temp_poly(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
Computation of temperature, water content and age in polythermal mode.
Definition: calc_temp_m.F90:57
Computation of the glacial isostatic adjustment of the lithosphere surface.
Definition: calc_gia_m.F90:36
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
subroutine, public calc_vxy_b_sia(time, z_sl)
Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice approximation.
Definition: calc_vxy_m.F90:55
Declarations of global variables for SICOPOLIS.