SICOPOLIS V5-dev  Revision 1264
sico_main_loop_iter_m.F90
Go to the documentation of this file.
2  implicit none
3 
4  public
5 
6 contains
7 
8  subroutine sico_main_loop_iter(delta_ts, glac_index, &
9  mean_accum, &
10  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
11  time, time_init, time_end, time_output, &
12  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
13  z_sl, dzsl_dtau, z_mar, &
14  ii, jj, nn, &
15  ndat2d, ndat3d, n_output, &
16  runname, &
17  itercount,iter_temp,iter_wss,iter_ser,&
18  iter_out,iter_output)
19 
20  use sico_types_m
22  use sico_vars_m
23  use boundary_m
24 #if (defined(GRL) && DISC>0)
26 #endif
27 #if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1)
28  use calc_temp_m
29 #elif (CALCMOD==2 || CALCMOD==3)
31 #endif
32 
33  use calc_enhance_m
34  use calc_vxy_m
35  use calc_vz_m
36  use calc_dxyz_m
37  use calc_gia_m
38  use calc_thk_m
40  use calc_bas_melt_m
42 #ifdef ALLOW_GRDCHK
43  use output_m
44 #endif
45  use ctrl_m, only: myfloor
46 
47  implicit none
48 
49  integer(i4b), intent(inout) :: ii((imax+1)*(jmax+1)), &
50  jj((imax+1)*(jmax+1)), &
51  nn(0:jmax,0:imax)
52  integer(i4b), intent(in) :: n_output
53  real(dp), intent(in) :: mean_accum
54  real(dp), intent(in) :: dtime, dtime_temp, dtime_wss, &
55  dtime_out, dtime_ser
56  real(dp), intent(in) :: time_init, time_end, time_output(100)
57  real(dp), intent(in) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
58  character(len=100), intent(in) :: runname
59 
60  integer(i4b), intent(inout) :: ndat2d, ndat3d
61  real(dp), intent(inout) :: delta_ts, glac_index
62  real(dp), intent(inout) :: time
63  real(dp), intent(inout) :: z_sl, dzsl_dtau, z_mar
64 
65  integer(i4b) :: i, j, kc, kt, kr, n
66  integer(i4b), intent(inout) :: itercount, iter_temp, iter_wss
67  integer(i4b), intent(inout) :: iter_ser, iter_out, iter_output(100)
68  real(dp) :: dtime_temp_inv
69  logical :: flag_3d_output
70 
71  real(dp) :: tmp, valmod
72  integer(i4b) :: valmodint
73 
74  write(unit=6, fmt='(2x,i0)') itercount
75 
76  !-------- Update of time --------
77 
78  time = time_init + real(itercount,dp)*dtime
79 
80  !-------- Save old mask --------
81 
83 
84  !-------- Boundary conditions --------
85 
86  call boundary(time, dtime, dxi, deta, delta_ts, glac_index, &
87  z_sl, dzsl_dtau, z_mar)
88 
89  !-------- Temperature, water content, age, flow enhancement factor --------
90 
91  if ( mod(itercount, iter_temp) == 0 ) then
92 
93  write(unit=6, fmt='(10x,a)') 'Computation of T'
94 
95  ! ------ Temperature, water content, age
96 
97 #if (CALCMOD==1)
98  call calc_temp_poly(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
99  ! polythermal scheme (POLY)
100 
101 #elif (CALCMOD==0)
102  call calc_temp_cold(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
103  ! cold-ice scheme (COLD)
104 
105 #elif (CALCMOD==2 || CALCMOD==3)
106  call calc_temp_enth(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
107  ! enthalpy scheme (ENTC or ENTM)
108 
109 #elif (CALCMOD==-1)
110  call calc_temp_const() ! isothermal scheme (ISOT)
111 
112 #else
113  stop ' >>> sico_main_loop: Parameter CALCMOD must be between -1 and 3!'
114 #endif
115 
116  ! ------ Time derivative of H_t (further time derivatives are
117  ! computed in subroutine calc_thk_xxx)
118 
119  dtime_temp_inv = 1.0_dp/dtime_temp
120 
121  dh_t_dtau = (h_t_neu - h_t)*dtime_temp_inv
122 
123  ! ------ New values --> old values
124 
125  n_cts = n_cts_neu
127  zm = zm_neu
128  h_c = h_c_neu
129  h_t = h_t_neu
131  age_c = age_c_neu
133  age_t = age_t_neu
135 
136 #if (CALCMOD==2 || CALCMOD==3)
140 #endif
141 
142  ! ------ Flow enhancement factor
143 
144 #if (ENHMOD==1)
145  call calc_enhance_1()
146 #elif (ENHMOD==2)
147  call calc_enhance_2()
148 #elif (ENHMOD==3)
149  call calc_enhance_3(time)
150 #elif (ENHMOD==4)
151  call calc_enhance_4()
152 #elif (ENHMOD==5)
153  call calc_enhance_5()
154 #else
155  stop ' >>> sico_main_loop: Parameter ENHMOD must be between 1 and 5!'
156 #endif
157 
158  end if
159  ! End of computation of temperature, water content, age and
160  ! enhancement factor
161 
162  !-------- Velocity --------
163 
164 #if (DYNAMICS==1 || DYNAMICS==2)
165 
166  call calc_vxy_b_sia(time, z_sl)
167  call calc_vxy_sia(dzeta_c, dzeta_t)
168 
169 #if (MARGIN==3 || DYNAMICS==2)
170  call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
171 #endif
172 
173  call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
174 
175 #if (MARGIN==3)
176  call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
177 #endif
178 
179 #elif (DYNAMICS==0)
180 
181  call calc_vxy_static()
182  call calc_vz_static()
183 
184 #else
185  stop ' >>> sico_main_loop: DYNAMICS must be either 0, 1 or 2!'
186 #endif
187 
188  call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
189 
190  !-------- Glacial isostatic adjustment and ice topography --------
191 
192  call calc_gia(time, z_sl, dtime, dxi, deta, itercount, iter_wss)
193 
194  call calc_thk_init()
195 
196 #if ((MARGIN==3 || DYNAMICS==2) && (CALCTHK==1 || CALCTHK==2 || CALCTHK==3))
197  write(6, fmt='(a)') ' >>> sico_main_loop:'
198  write(6, fmt='(a)') ' Non-SIA dynamics combined with'
199  write(6, fmt='(a)') ' SIA ice thickness solver!'
200  stop
201 #endif
202 
203 #if (CALCTHK==1)
204  call calc_thk_sia_expl(time, dtime, dxi, deta, z_mar)
205 #elif (CALCTHK==2 || CALCTHK==3)
206  call calc_thk_sia_impl(time, dtime, dxi, deta, z_mar, mean_accum, ii, jj, nn)
207 #elif (CALCTHK==4)
208  call calc_thk_expl(time, dtime, dxi, deta, z_mar)
209 #elif (CALCTHK==5 || CALCTHK==6)
210  call calc_thk_impl(time, dtime, dxi, deta, z_mar, mean_accum, ii, jj, nn)
211 #else
212  stop ' >>> sico_main_loop: Parameter CALCTHK must be between 1 and 6!'
213 #endif
214 
215 #if (MARGIN==3) /* coupled SIA/SSA or SIA/SStA/SSA dynamics */
216  call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 3_i2b)
217 #elif (DYNAMICS==2) /* hybrid SIA/SStA dynamics */
218  call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 2_i2b)
219 #else /* SIA-only dynamics */
220 #if (CALCTHK==1 || CALCTHK==2 || CALCTHK==3)
221  call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 1_i2b)
222 #elif (CALCTHK==4 || CALCTHK==5 || CALCTHK==6)
223  call calc_thk_mask_update(time, dtime, dxi, deta, z_sl, z_mar, 2_i2b)
224 #endif
225 #endif
226 
227  ! ------ New values --> old values
228 
229  zs = zs_neu
230  zm = zm_neu
231  zb = zb_neu
232  zl = zl_neu
233  h_c= h_c_neu
234  h_t= h_t_neu
235 
236  !-------- Melting temperature --------
237 
238  call calc_temp_melt()
239 
240  !-------- Basal temperature --------
241 
242  call calc_temp_bas()
243 
244  !-------- Basal melting rate --------
245 
246  call calc_qbm(time, z_sl, dzeta_c, dzeta_r)
247 
248  !-------- Effective thickness of subglacial water --------
249 
250  call calc_thk_water_bas(z_sl)
251 
252  !-------- Data output --------
253 #ifndef ALLOW_OPENAD
254 #if (OUTPUT==1)
255 
256  if ( mod(itercount, iter_out) == 0 ) then
257 
258 #if (ERGDAT==0)
259  flag_3d_output = .false.
260 #elif (ERGDAT==1)
261  flag_3d_output = .true.
262 #endif
263 
264  call output1(runname, time, delta_ts, glac_index, z_sl, &
265  flag_3d_output, ndat2d, ndat3d)
266 
267  end if
268 
269 #elif (OUTPUT==2)
270 
271  do n=1, n_output
272 
273  if (itercount == iter_output(n)) then
274 
275 #if (ERGDAT==0)
276  flag_3d_output = .false.
277 #elif (ERGDAT==1)
278  flag_3d_output = .true.
279 #endif
280 
281  call output1(runname, time, delta_ts, glac_index, z_sl, &
282  flag_3d_output, ndat2d, ndat3d)
283 
284  end if
285 
286  end do
287 
288 #elif (OUTPUT==3)
289 
290  if ( mod(itercount, iter_out) == 0 ) then
291 
292  flag_3d_output = .false.
293 
294  call output1(runname, time, delta_ts, glac_index, z_sl, &
295  flag_3d_output, ndat2d, ndat3d)
296 
297  end if
298 
299  do n=1, n_output
300 
301  if (itercount == iter_output(n)) then
302 
303  flag_3d_output = .true.
304 
305  call output1(runname, time, delta_ts, glac_index, z_sl, &
306  flag_3d_output, ndat2d, ndat3d)
307 
308  end if
309 
310  end do
311 
312 #endif
313 
314  if ( mod(itercount, iter_ser) == 0 ) then
315  call output2(time, dxi, deta, delta_ts, glac_index, z_sl)
316  call output4(time, dxi, deta, delta_ts, glac_index, z_sl)
317 #if (defined(ASF) && WRITE_SER_FILE_STAKES>0)
318  call output5(time, dxi, deta, delta_ts, glac_index, z_sl)
319 #endif
320 
321  end if
322 #endif
323 
324  end subroutine sico_main_loop_iter
325 
326 end module sico_main_loop_iter_m
subroutine, public calc_thk_sia_impl(time, dtime, dxi, deta, z_mar, mean_accum, ii, jj, nn)
Over-implicit solver for the diffusive SIA ice surface equation.
Definition: calc_thk_m.F90:251
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
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
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
integer(i2b), 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...
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
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
subroutine, public output1(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in native binary or NetCDF format.
Definition: output_m.F90:59
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(i2b), 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 ...
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:173
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...
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:611
subroutine, public calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
subroutine, public output4(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data of the deep ice cores on file in ASCII format (and optionally in NetCDF f...
Definition: output_m.F90:4877
integer(i2b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
subroutine, public calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
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:35
subroutine, public calc_temp_melt()
Computation of the melting temperatures.
subroutine, public calc_thk_impl(time, dtime, dxi, deta, z_mar, mean_accum, ii, jj, nn)
Over-implicit solver for the general ice thickness equation.
Definition: calc_thk_m.F90:753
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
integer(i2b), dimension(0:jmax, 0:imax) kc_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
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.
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:54
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
integer(i2b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
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
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, ii, jj, nn, ndat2d, ndat3d, n_output, runname, itercount, iter_temp, iter_wss, iter_ser, iter_out, iter_output)
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 (...
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:482
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:57
subroutine, public output2(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data on file in ASCII format (and optionally in NetCDF format).
Definition: output_m.F90:3521
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:56
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
Ice discharge parameterization for the Greenland ice sheet.
integer(i2b), dimension(0:jmax, 0:imax) n_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
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
Writing of output data on files.
Definition: output_m.F90:36
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:53
Declarations of global variables for SICOPOLIS.