SICOPOLIS V5-dev  Revision 1368
openad_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : o p e n a d _ 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 
32 !-------------------------------------------------------------------------------
33 !> Module for all openad-related subroutines
34 !<------------------------------------------------------------------------------
35 
36 
37 module openad_m
38 
39  implicit none
40 
41  private
42  public :: adjoint_master
43 
44 #if (defined (ALLOW_GRDCHK))
45  public :: grdchk_main
46 #elif (defined (ALLOW_OPENAD))
47 ! public :: adjoint_init
48 ! public :: tape_mode
49 ! public :: adjoint_mode
50  private :: direct_substitution
51  private :: print_output
52  private :: print_var
53 #endif
54 
55 contains
56 
57 !-------------------------------------------------------------------------------
58 !> Adjoint master is the main tool by which sicopolis.F90 invokes the adjoint
59 !! code. Its job is to figure out what mode of the adjoint code is being invoked
60 !! and run the appropriate subroutine.
61 !<------------------------------------------------------------------------------
62  subroutine adjoint_master
63 
64  implicit none
65 
66 #ifdef ALLOW_OPENAD
67  integer :: mode
68  character(len=32) :: arg
69  logical :: ISPLAIN, ISTAPE, ISADJOINT
70 #endif
71 
72 #ifdef ALLOW_GRDCHK
73  call grdchk_main
74 #endif
75 
76 #ifdef ALLOW_OPENAD
77  call get_ad_mode(mode,arg,isplain, istape, isadjoint)
78  call direct_substitution(isplain, istape, isadjoint)
79 #endif
80 
81  end subroutine adjoint_master
82 
83 #ifdef ALLOW_GRDCHK
84 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
85 ! Subroutine : g r d c h k _ m a i n
86 ! Purpose : Gradient check top level routine
87 ! Compares the gradients calculated by the adjoint model
88 ! to the gradients calculated by finite difference
89 ! approximations
90 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
91 
92  subroutine grdchk_main
93 
94  use sico_types_m
96  use sico_vars_m
97 
98  use sico_init_m
100  use sico_end_m
101 
103  cost_final
104 
105  implicit none
106 
107  integer(i4b) :: ii((imax+1)*(jmax+1)), &
108  jj((imax+1)*(jmax+1)), &
109  nn(0:jmax,0:imax)
110  integer(i4b) :: ndat2d, ndat3d
111  integer(i4b) :: n_output
112  real(dp) :: delta_ts, glac_index
113  real(dp) :: mean_accum
114  real(dp) :: dtime, dtime_temp, dtime_wss, &
115  dtime_out, dtime_ser
116  real(dp) :: time, time_init, time_end, time_output(100)
117  real(dp) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
118  real(dp) :: z_sl, dzsl_dtau, z_mar
119  character(len=100) :: runname
120 
121  !-------- Variable declarations needed for this routine specifically
122  ! Should we be perturbing these
123  ! values by basically 5% of ori-
124  ! ginal?
125  !real(dp) :: orig_val, perturb_val = 1e-6
126  real(dp) :: orig_val, perturb_val = 5e-2
127  real(dp), dimension(3) :: fc_collected
128  real(dp), dimension(3) :: direction
129  real(dp) :: gfd0,gfd
130  integer(i4b), parameter :: points = 10
131  integer(i4b), dimension(points) :: ipoints, jpoints
132  integer(i4b) :: i, j, p, d
133  character(len=100) :: fname
134 
135  !-------- This array holds the direction of perturbations to follow:
136  direction(1) = 0
137  direction(2) = 1
138  direction(3) = -1
139 
140 #if (!defined(GRL) && !defined(ANT))
141  print *, ">>> Adjoint only available for GRL and ANT right now; kill code."
142 #endif
143 
144  do p = 1, points
145 #if (defined(GRL))
146  ipoints(p) = int(real(imax/2))
147  jpoints(p) = int(real(jmax/5)) + (p-1) * points
148 #elif (defined(ANT))
149 ! ipoints(p) = int(real(IMAX/3)) + int(real((.85-.33)*IMAX/points)) * (p - 1)
150 ! jpoints(p) = int(real(JMAX/2))
151  ipoints = (/94, 240/) ! tsurf = temp_c(:,:,KCMAX)
152  jpoints = (/187, 66/) ! tsurf = temp_c(:,:,KCMAX)
153 #endif
154  end do
155 
156  !-------- Priming output file for gradients
157  open(99, file='GradientVals_'//trim(runname)//'.dat',&
158  form="FORMATTED", status="REPLACE")
159  open(98, file='CostVals_'//trim(runname)//'.dat',&
160  form="FORMATTED", status="REPLACE")
161 
162  !-------- Loop over points
163  do p = 1, points
164  i = ipoints(p)
165  j = jpoints(p)
166 
167  !-------- Loop over perturbation direction (0, +, -)
168  do d = 1, 3 ! (4-th loop)
169 
170  !-------- Let yourself know where you are:
171  print *, ' point (p, i, j), direction (d) [ ', p , ', ', i, ', ', j, ', ', d, ' ] '
172 
173  ! A complete sico loop
174  call deldirs
176 
177  call sico_init(delta_ts, glac_index, &
178  mean_accum, &
179  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
180  time, time_init, time_end, time_output, &
181  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
182  z_sl, dzsl_dtau, z_mar, &
183  ii, jj, nn, &
184  ndat2d, ndat3d, n_output, &
185  runname)
186 
187  ! store original value that will be perturbed
188  ! and then perturb it (note: direction(1) = 0)
189  ! -- H_c
190  !orig_val = H_c(j,i)
191  !H_c(j,i) = orig_val * (1 + direction(d)*perturb_val)
192 
193  ! -- temp_c (KCMAX==SURFACE, 0==BASE)
194  !orig_val = temp_c(KCMAX,j,i)
195  !temp_c(KCMAX,j,i) = temp_c(KCMAX,j,i)*(1 +direction(d)*perturb_val)
196 
197  ! -- mean annual temp
198  !orig_val = temp_ma_present(j,i)
199  !temp_ma_present(j,i) = orig_val*(1 +direction(d)*perturb_val)
200 
201  ! -- mean annual temp
202  !orig_val = q_geo(j,i)
203  !q_geo(j,i) = orig_val*(1 +direction(d)*perturb_val)
204 
205  ! -- mean annual temp
206  orig_val = vx_c(kcmax,j,i)
207  vx_c(kcmax,j,i) = orig_val*(1 +direction(d)*perturb_val)
208 
209  ! -- mean annual temp
210  !orig_val = c_slide(j,i)
211  !c_slide(j,i) = orig_val*(1 +direction(d)*perturb_val)
212 
213  ! -- mean annual temp
214  !orig_val = c_drag(j,i)
215  !c_drag(j,i) = orig_val*(1 +direction(d)*perturb_val)
216 
217  ! -- precip_present
218  !orig_val = precip_present(j,i,1)
219  !precip_present(j,i,1) = orig_val*(1 +direction(d)*perturb_val)
220 
221  ! -- experimental field acc_fact
222  !orig_val = acc_fact(j,i)
223  !acc_fact(j,i) = orig_val*(1 +direction(d)*perturb_val)
224 
225  ! -- sanity check
226  write(6,fmt='(a,f40.20)') "orig_val = ", orig_val
227  write(6,fmt='(a,f40.20)') "pert_val = ", orig_val*(1+direction(d)*perturb_val)
228 
229  call sico_main_loop(delta_ts, glac_index, &
230  mean_accum, &
231  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
232  time, time_init, time_end, time_output, &
233  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
234  z_sl, dzsl_dtau, z_mar, &
235  ii, jj, nn, &
236  ndat2d, ndat3d, n_output, &
237  runname)
238 
239  call cost_final(runname)
240  call sico_end
241 
242  ! store cost
243  fc_collected(d) = fc
244 
245  end do ! (close perturb loop)
246 
247  ! --------- calculate simple 2-sided finite difference due to
248  ! (+) and (-) perturbation
249  gfd = (fc_collected(2) - fc_collected(3))/(2.d0*perturb_val*orig_val)
250 
251  write(6, fmt='(a,f40.20)') "Finite difference is = ", gfd
252 
253  ! --------- write these values to output file
254  write(99, fmt='(f40.20)') gfd
255  write(98, fmt='(f40.20)') fc_collected(1)
256  write(98, fmt='(f40.20)') fc_collected(2)
257  write(98, fmt='(f40.20)') fc_collected(3)
258  write(98, fmt='(a)') '----------------------------------'
259 
260  end do ! (close point loop)
261 
262  close(unit=99)
263  close(unit=98)
264 
265  end subroutine grdchk_main
266 #endif
267 
268 !!-------------------------------------------------------------------------------
269 !!> Checks to see if output dir exists. If so, deletes it.
270 !!<------------------------------------------------------------------------------
271  subroutine deldirs
272 
273  implicit none
274 
275  character(len=256) :: shell_command
276 
277  !-------- deleting directories
278 
279  shell_command = 'if [ -d'
280  shell_command = trim(shell_command)//' '//outpath
281  shell_command = trim(shell_command)//' '//'] ; then rm -rf'
282  shell_command = trim(shell_command)//' '//outpath
283  shell_command = trim(shell_command)//' '//'; fi'
284 
285  call system(trim(shell_command))
286 
287  end subroutine deldirs
288 
289 #ifdef ALLOW_OPENAD
290 !-------------------------------------------------------------------------------
291 !> This is the top-level routine for the adjoint mode sweep. It takes arguments
292 !! from the executable {-p, -t, -a} supplied at run time for the plain, tape, and
293 !! adjoint modes of execution.
294 !<------------------------------------------------------------------------------
295  subroutine direct_substitution(ISPLAIN, ISTAPE, ISADJOINT)
296 
297  use oad_rev
298  use oad_tape
299  use sico_variables_m
300 #ifdef GRL
302 #endif
305  use sico_init_m
306  use sico_end_m
307 
308  !integer :: mode
309  !character(len=32) :: arg
310  logical :: ISPLAIN, ISTAPE, ISADJOINT
311 
312  integer(i4b) :: ndat2d, ndat3d
313  integer(i4b) :: n_output
314  integer(i4b), dimension((IMAX+1)*(JMAX+1)) :: ii, jj
315  integer(i4b), dimension(0:JMAX,0:IMAX) :: nn
316  real(dp) :: delta_ts, glac_index
317  real(dp) :: mean_accum
318  real(dp) :: dtime, dtime_temp, &
319  dtime_wss, dtime_out, dtime_ser
320  real(dp) :: time, time_init, time_end
321  real(dp), dimension(100) :: time_output
322  real(dp) :: dxi, deta, dzeta_c, &
323  dzeta_t, dzeta_r
324  real(dp) :: z_sl, dzsl_dtau, z_mar
325  character(len=100) :: runname
326 
327  our_rev_mode%arg_store = .false.
328  our_rev_mode%arg_restore= .false.
329 
330  our_rev_mode%res_store = .false.
331  our_rev_mode%res_restore= .false.
332 
333  our_rev_mode%plain = .false.
334  our_rev_mode%tape = .false.
335  our_rev_mode%adjoint = .false.
336 
338 
339  call sico_init(delta_ts, glac_index, &
340  mean_accum, &
341  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
342  time, time_init, time_end, time_output, &
343  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
344  z_sl, dzsl_dtau, z_mar, &
345  ii, jj, nn, &
346  ndat2d, ndat3d, n_output, &
347  runname)
348 
349  call var_transfer&
381  rf,rf_imq,kappa_imq,c_imq,year_zero,&
391  kei_r_incr,fc,objf_test,&
392  mult_test,target_topo_tau,&
393  c_int_table,c_int_inv_table,c,n_temp_min,n_temp_max,n_enth_min,&
394  n_enth_max,l_inv,l_eto,&
395  rho_i_imq,rho_c_imq,kappa_c_imq,c_c_imq,&
398 #if (defined(INITMIP_SMB_ANOM_FILE))
399  as_anom_initmip,&
400 #endif
401 #if (defined(INITMIP_BMB_ANOM_FILE))
402  ab_anom_initmip,&
403 #endif
404 #if (ICE_STREAM==2)
405  maske_sedi,&
406 #endif
407 #if (DISC==2)
408  glann_time_min,glann_time_stp,glann_time_max,&
409  ndata_glann,dt_glann_climber,&
410 #endif
411 #ifdef GRL
412  disc_dw,n_discharge_call_dw,iter_mar_coa_dw,c_dis_0_dw,s_dis_dw,&
413  c_dis_fac_dw,t_sub_pd_dw,alpha_sub_dw,alpha_o_dw,m_h_dw,&
414  m_d_dw,r_mar_eff_dw,t_sea_freeze_dw,c_dis_dw,&
415 #elif (defined(ANT))
416  n_sector,&
417 #endif
419 #if (defined(AGE_COST))
420  age_data,age_unc, &
421 #endif
422  n_temp_min_imq,n_temp_max_imq,&
424 
425  our_rev_mode%plain=isplain
426  our_rev_mode%tape=istape
427 
428  if(istape.eqv..true.) call oad_tape_init()
429  call sicopolis_openad(delta_ts, glac_index, &
430  mean_accum, &
431  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
432  time, time_init, time_end, time_output, &
433  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
434  z_sl, dzsl_dtau, z_mar, &
435  ii, jj, nn, &
436  ndat2d, ndat3d, n_output, &
437  runname)
438 
439  if(isadjoint) then
440  our_rev_mode%tape=.false.
441  our_rev_mode%adjoint=.true.
442  call sicopolis_openad(delta_ts, glac_index, &
443  mean_accum, &
444  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
445  time, time_init, time_end, time_output, &
446  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
447  z_sl, dzsl_dtau, z_mar, &
448  ii, jj, nn, &
449  ndat2d, ndat3d, n_output, &
450  runname)
451  call print_output(runname)
452  call oad_tape_delete()
453  else
454  if(istape.eqv..true.) call oad_tape_delete()
455  end if
456 
457  call sico_end
458 
459  end subroutine direct_substitution
460 
461 !-------------------------------------------------------------------------------
462 !> This is where the final cost fc is initialized. See ctrl_m.F90 (subroutine
463 !! final_cost) to see its structure.
464 !<------------------------------------------------------------------------------
465  subroutine oad_independent_init()
466  use oad_active
467  use oad_sico_variables_m
468 
469  implicit none
470 
471  fc%d = 1.0_dp
472 
473  end subroutine oad_independent_init
474 
475 !-------------------------------------------------------------------------------
476 !> Take arguments passed to executable called driver in OAD mode.
477 !! Only options are (should be)
478 !! "-h help"
479 !! "-p PLAIN MODE"
480 !! "-t TAPE MODE"
481 !! "-a ADJOINT MODE"
482 !!
483 !<------------------------------------------------------------------------------
484  subroutine get_ad_mode(mode,arg,ISPLAIN, ISTAPE, ISADJOINT)
485 
486  implicit none
487 
488  integer :: mode
489  character(len=32) :: arg
490 
491  logical :: ISPLAIN, ISTAPE, ISADJOINT
492 
493  isplain = .false.
494  istape = .false.
495  isadjoint = .false.
496 
497  do mode=1,command_argument_count()
498  ! grabbing the input to driver
499  call get_command_argument(mode,arg)
500 
501  ! going through possible arguments to driver
502  select case(adjustl(arg))
503  case("-p")
504  isplain = .true.
505  case("-t")
506  istape = .true.
507  case("-a")
508  isadjoint = .true.
509  case("-h")
510  print *, "Command line options:"
511  print *, "-h help"
512  print *, "-p PLAIN MODE"
513  print *, "-t TAPE MODE"
514  print *, "-a ADJOINT MODE"
515  stop
516  case default
517  print *, "Unknown command line option ", arg
518  print *, "Command line options are:"
519  print *, "-h help"
520  print *, "-p PLAIN MODE"
521  print *, "-t TAPE MODE"
522  print *, "-a ADJOINT MODE"
523  stop
524  end select
525  end do
526 
527  if ( (isplain .NEQV. .true.) &
528  .AND. (istape .NEQV. .true.) &
529  .AND. (isadjoint .NEQV. .true.)) then
530  print *, " No valid option specified."
531  print *, "Command line options are:"
532  print *, "-h help"
533  print *, "-p PLAIN MODE"
534  print *, "-t TAPE MODE"
535  print *, "-a ADJOINT MODE"
536  stop
537  end if
538 
539  if(isadjoint .EQV. .true.) then
540  istape = .true.
541  end if
542 
543  if((isplain .EQV. .true.) .AND. (istape .EQV. .true.)) then
544  print *, " Cannot specify both -p along with any other option"
545  stop
546  end if
547 
548  end subroutine get_ad_mode
549 
550 !-------------------------------------------------------------------------------
551 !> This is an absolutely essential routine that takes all of the variables
552 !! foobar used in the forward (tape) sweep and assigns them to a_foobar
553 !! variables *needed* to perform the reverse (adjoint) sweep. If your
554 !! compilation FAILS with the error message complaining that something was not
555 !! ACTIVE or was ACTIVE and needed not be, the solution is to be found inside
556 !! this routine, by taking variable foobar and making it either active by
557 !! assigning:
558 !! foobar%v = a_foobar (to make foobar ACTIVE) or
559 !! foobar = a_foobar (to make foobar NOT ACTIVE). But that's not all!
560 !! If your regression script suddenly is doing quite poorly (FDs vs. ADs are not
561 !! within tolerance anymore) is it likely that some new variable as been
562 !! introduced in the trunk and needs to be *transfered* here. In other words, in
563 !! the commit history of sicopolis, someone introduced new_foobar, and
564 !! new_foobar is completely missing in the below *and* in the calling (in
565 !! subroutine direct_substitution()) of var_transfer above. This needs to be
566 !! amended before proceeding. If you feel unsure, email
567 !! liz.curry.logan@gmail.com.
568 !<------------------------------------------------------------------------------
569  subroutine var_transfer&
570  (a_maske,a_maske_old,a_maske_neu,a_n_cts,a_n_cts_neu,a_kc_cts,&
571  a_kc_cts_neu,a_flag_inner_point,&
572  a_flag_grounding_line_1,a_flag_grounding_line_2,a_flag_calving_front_1,&
573  a_flag_calving_front_2,a_flag_shelfy_stream_x,a_flag_shelfy_stream_y,&
574  a_flag_shelfy_stream,a_xi,a_eta,a_zeta_c,a_zeta_t,a_zeta_r,a_aa,a_flag_aa_nonzero,&
575  a_ea,a_eaz_c,a_eaz_c_quotient,a_lambda,a_phi,a_area,a_sq_g11_g,a_sq_g22_g,&
576  a_insq_g11_g,a_insq_g22_g,a_sq_g11_sgx,a_sq_g11_sgy,a_sq_g22_sgx,a_sq_g22_sgy,&
577  a_insq_g11_sgx,a_insq_g22_sgy,a_zs,a_zm,a_zb,a_zl,a_zl0,a_wss,a_flex_rig_lith,&
578  a_time_lag_asth,a_h_c,a_h_t,a_dzs_dxi,a_dzm_dxi,a_dzb_dxi,a_dh_c_dxi,a_dh_t_dxi,&
579  a_dzs_deta,a_dzm_deta,a_dzb_deta,a_dh_c_deta,a_dh_t_deta,a_dzs_dxi_g,&
580  a_dzm_dxi_g,a_dzb_dxi_g,a_dh_c_dxi_g,a_dh_t_dxi_g,a_dzs_deta_g,a_dzm_deta_g,&
581  a_dzb_deta_g,a_dh_c_deta_g,a_dh_t_deta_g,a_dzs_dtau,a_dzm_dtau,a_dzb_dtau,&
582  a_dzl_dtau,a_dh_c_dtau,a_dh_t_dtau,a_p_weert,a_q_weert,a_p_weert_inv,&
583  a_c_slide,a_d_help_b,a_c_drag,a_p_b_w,a_vx_b,a_vy_b,a_vx_m,a_vy_m,a_ratio_sl_x,&
584  a_ratio_sl_y,a_vx_b_g,a_vy_b_g,a_vz_b,a_vz_m,a_vx_s_g,a_vy_s_g,a_vz_s,&
585  a_flui_ave_sia,a_h_diff,a_qx,a_qy,a_q_gl_g,a_q_geo,a_temp_b,a_temph_b,a_q_bm,a_q_b_apl,&
586  a_q_tld,a_q_b_tot,a_h_w,a_accum,a_evap,a_runoff,a_runoff_apl,a_as_perp,a_temp_s,a_am_perp,&
587  a_am_perp_st,a_zs_neu,a_zm_neu,a_zb_neu,a_zl_neu,a_h_c_neu,a_h_t_neu,a_zs_ref,&
588  a_accum_present,a_precip_ma_present,a_precip_ma_lgm_anom,&
589  a_temp_ma_present,a_temp_mj_present,a_temp_ma_lgm_anom,a_temp_mj_lgm_anom,&
590  a_dist_dxdy,a_acc_fact,a_precip_present,a_precip_lgm_anom,a_gamma_precip_lgm_anom,&
591  a_temp_mm_present,a_temp_mm_lgm_anom,a_d_help_c,a_vx_c,a_vy_c,a_vz_c,a_temp_c,&
592  a_temp_c_neu,a_temp_c_m,a_age_c,a_age_c_neu,a_txz_c,a_tyz_c,a_sigma_c,a_enh_c,&
593  a_de_ssa,a_vis_int_g,a_vx_g,a_vy_g,a_d_help_t,a_vx_t,a_vy_t,a_vz_t,a_omega_t,&
594  a_omega_t_neu,a_temp_t_m,a_age_t,a_age_t_neu,a_txz_t,a_tyz_t,a_sigma_t,a_enh_t,&
595  a_temp_r,a_temp_r_neu,a_enth_c,a_enth_c_neu,a_omega_c,a_omega_c_neu,a_enth_t,&
596  a_enth_t_neu,a_dxx_c,a_dyy_c,a_dxy_c,a_dxz_c,a_dyz_c,a_de_c,a_lambda_shear_c,&
597  a_dxx_t,a_dyy_t,a_dxy_t,a_dxz_t,a_dyz_t,a_de_t,a_lambda_shear_t,a_rho,a_rho_w,&
598  a_rho_sw,a_l,a_g,a_nue,a_beta,a_delta_tm_sw,a_omega_max,a_h_r,a_rho_c_r,a_kappa,a_kappa_r,&
599  a_rho_a,a_r_t,a_r_t_imq,a_r,a_a,a_b,a_lambda0,a_phi0,a_s_stat_0,a_beta1_0,a_beta1_lt_0,&
600  a_beta1_ht_0,a_beta2_0,a_beta2_lt_0,a_beta2_ht_0,a_phi_sep_0,a_pmax_0,a_mu_0,&
601  a_rf,a_rf_imq,a_kappa_imq,a_c_imq,a_year_zero,&
602  a_forcing_flag,a_n_core,a_lambda_core,a_phi_core,a_x_core,&
603  a_y_core,a_grip_time_min,a_grip_time_stp,a_grip_time_max,a_ndata_grip,&
604  a_griptemp,a_gi_time_min,a_gi_time_stp,a_gi_time_max,a_ndata_gi,&
605  a_glacial_index,a_specmap_time_min,a_specmap_time_stp,a_specmap_time_max,&
606  a_ndata_specmap,a_specmap_zsl,a_time_target_topo_init,&
607  a_time_target_topo_final,a_maske_target,a_zs_target,a_zb_target,&
608  a_zl_target,a_h_target,a_maske_maxextent,a_mask_ablation_type,a_ncid_temp_precip,&
609  a_ndata_temp_precip,a_temp_precip_time_min,a_temp_precip_time_stp,&
610  a_temp_precip_time_max,a_temp_precip_time,a_kei,a_kei_r_max,&
611  a_kei_r_incr,a_fc,a_objf_test,&
612  a_mult_test,a_target_topo_tau,&
613  a_c_int_table,a_c_int_inv_table,a_c,a_n_temp_min,a_n_temp_max,a_n_enth_min,&
614  a_n_enth_max,a_l_inv,a_l_eto,&
615  a_rho_i_imq,a_rho_c_imq,a_kappa_c_imq,a_c_c_imq, &
616  a_as_perp_apl,a_mb_source_apl, &
617  a_calv_grounded,a_calv_grounded_apl,a_ncid_ser,a_ncid_core,a_n_data_kei,&
618 #if (defined(INITMIP_SMB_ANOM_FILE))
619  a_as_anom_initmip,&
620 #endif
621 #if (defined(INITMIP_BMB_ANOM_FILE))
622  a_ab_anom_initmip,&
623 #endif
624 #if (ICE_STREAM==2)
625  a_maske_sedi,&
626 #endif
627 #if (DISC==2)
628  a_glann_time_min,a_glann_time_stp,a_glann_time_max,&
629  a_ndata_glann,a_dt_glann_climber,&
630 #endif
631 #ifdef GRL
632  a_disc_dw,a_n_discharge_call_dw,a_iter_mar_coa_dw,a_c_dis_0_dw,a_s_dis_dw,&
633  a_c_dis_fac_dw,a_t_sub_pd_dw,a_alpha_sub_dw,a_alpha_o_dw,a_m_h_dw,&
634  a_m_d_dw,a_r_mar_eff_dw,a_t_sea_freeze_dw,a_c_dis_dw,&
635 #elif (defined(ANT))
636  a_n_sector,&
637 #endif
638  a_et,a_melt,a_melt_star,a_rainfall,a_snowfall,&
639 #if (defined(AGE_COST))
640  a_age_data,a_age_unc, &
641 #endif
642  a_n_temp_min_imq,a_n_temp_max_imq,&
643  a_smb_corr,a_smb_corr_prescribed)
644 
645  use oad_active
646  use oad_sico_variables_m
647  use oad_sico_vars_m
648  use oad_enth_temp_omega_m
649  use oad_ice_material_properties_m
650 #ifdef GRL
651  use oad_discharge_workers_m
652 #endif
653 
654  implicit none
655 
656  integer(i1b), dimension(0:JMAX,0:IMAX) :: a_maske
657  integer(i1b), dimension(0:JMAX,0:IMAX) :: a_maske_old
658  integer(i1b), dimension(0:JMAX,0:IMAX) :: a_maske_neu
659  integer(i1b), dimension(0:JMAX,0:IMAX) :: a_n_cts
660  integer(i1b), dimension(0:JMAX,0:IMAX) :: a_n_cts_neu
661  integer(i4b), dimension(0:JMAX,0:IMAX) :: a_kc_cts
662  integer(i4b), dimension(0:JMAX,0:IMAX) :: a_kc_cts_neu
663  logical, dimension(0:JMAX,0:IMAX) :: a_flag_inner_point
664  logical, dimension(0:JMAX,0:IMAX) :: a_flag_grounding_line_1
665  logical, dimension(0:JMAX,0:IMAX) :: a_flag_grounding_line_2
666  logical, dimension(0:JMAX,0:IMAX) :: a_flag_calving_front_1
667  logical, dimension(0:JMAX,0:IMAX) :: a_flag_calving_front_2
668  logical, dimension(0:JMAX,0:IMAX) :: a_flag_shelfy_stream_x
669  logical, dimension(0:JMAX,0:IMAX) :: a_flag_shelfy_stream_y
670  logical, dimension(0:JMAX,0:IMAX) :: a_flag_shelfy_stream
671  real(dp), dimension(0:IMAX) :: a_xi
672  real(dp), dimension(0:JMAX) :: a_eta
673  real(dp), dimension(0:KCMAX) :: a_zeta_c
674  real(dp), dimension(0:KTMAX) :: a_zeta_t
675  real(dp), dimension(0:KRMAX) :: a_zeta_r
676  real(dp) :: a_aa
677  logical :: a_flag_aa_nonzero
678  real(dp) :: a_ea
679  real(dp), dimension(0:KCMAX) :: a_eaz_c
680  real(dp), dimension(0:KCMAX) :: a_eaz_c_quotient
681  real(dp), dimension(0:JMAX,0:IMAX) :: a_lambda
682  real(dp), dimension(0:JMAX,0:IMAX) :: a_phi
683  real(dp), dimension(0:JMAX,0:IMAX) :: a_area
684  real(dp), dimension(0:JMAX,0:IMAX) :: a_sq_g11_g
685  real(dp), dimension(0:JMAX,0:IMAX) :: a_sq_g22_g
686  real(dp), dimension(0:JMAX,0:IMAX) :: a_insq_g11_g
687  real(dp), dimension(0:JMAX,0:IMAX) :: a_insq_g22_g
688  real(dp), dimension(0:JMAX,0:IMAX) :: a_sq_g11_sgx
689  real(dp), dimension(0:JMAX,0:IMAX) :: a_sq_g11_sgy
690  real(dp), dimension(0:JMAX,0:IMAX) :: a_sq_g22_sgx
691  real(dp), dimension(0:JMAX,0:IMAX) :: a_sq_g22_sgy
692  real(dp), dimension(0:JMAX,0:IMAX) :: a_insq_g11_sgx
693  real(dp), dimension(0:JMAX,0:IMAX) :: a_insq_g22_sgy
694  real(dp), dimension(0:JMAX,0:IMAX) :: a_zs
695  real(dp), dimension(0:JMAX,0:IMAX) :: a_zm
696  real(dp), dimension(0:JMAX,0:IMAX) :: a_zb
697  real(dp), dimension(0:JMAX,0:IMAX) :: a_zl
698  real(dp), dimension(0:JMAX,0:IMAX) :: a_zl0
699  real(dp), dimension(0:JMAX,0:IMAX) :: a_wss
700  real(dp), dimension(0:JMAX,0:IMAX) :: a_flex_rig_lith
701  real(dp), dimension(0:JMAX,0:IMAX) :: a_time_lag_asth
702  real(dp), dimension(0:JMAX,0:IMAX) :: a_H_c
703  real(dp), dimension(0:JMAX,0:IMAX) :: a_H_t
704  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzs_dxi
705  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzm_dxi
706  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzb_dxi
707  real(dp), dimension(0:JMAX,0:IMAX) :: a_dH_c_dxi
708  real(dp), dimension(0:JMAX,0:IMAX) :: a_dH_t_dxi
709  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzs_deta
710  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzm_deta
711  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzb_deta
712  real(dp), dimension(0:JMAX,0:IMAX) :: a_dH_c_deta
713  real(dp), dimension(0:JMAX,0:IMAX) :: a_dH_t_deta
714  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzs_dxi_g
715  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzm_dxi_g
716  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzb_dxi_g
717  real(dp), dimension(0:JMAX,0:IMAX) :: a_dH_c_dxi_g
718  real(dp), dimension(0:JMAX,0:IMAX) :: a_dH_t_dxi_g
719  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzs_deta_g
720  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzm_deta_g
721  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzb_deta_g
722  real(dp), dimension(0:JMAX,0:IMAX) :: a_dH_c_deta_g
723  real(dp), dimension(0:JMAX,0:IMAX) :: a_dH_t_deta_g
724  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzs_dtau
725  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzm_dtau
726  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzb_dtau
727  real(dp), dimension(0:JMAX,0:IMAX) :: a_dzl_dtau
728  real(dp), dimension(0:JMAX,0:IMAX) :: a_dH_c_dtau
729  real(dp), dimension(0:JMAX,0:IMAX) :: a_dH_t_dtau
730  integer(i4b), dimension(0:JMAX,0:IMAX) :: a_p_weert
731  integer(i4b), dimension(0:JMAX,0:IMAX) :: a_q_weert
732  real(dp), dimension(0:JMAX,0:IMAX) :: a_p_weert_inv
733  real(dp), dimension(0:JMAX,0:IMAX) :: a_c_slide
734  real(dp), dimension(0:JMAX,0:IMAX) :: a_d_help_b
735  real(dp), dimension(0:JMAX,0:IMAX) :: a_c_drag
736  real(dp), dimension(0:JMAX,0:IMAX) :: a_p_b_w
737  real(dp), dimension(0:JMAX,0:IMAX) :: a_vx_b
738  real(dp), dimension(0:JMAX,0:IMAX) :: a_vy_b
739  real(dp), dimension(0:JMAX,0:IMAX) :: a_vx_m
740  real(dp), dimension(0:JMAX,0:IMAX) :: a_vy_m
741  real(dp), dimension(0:JMAX,0:IMAX) :: a_ratio_sl_x
742  real(dp), dimension(0:JMAX,0:IMAX) :: a_ratio_sl_y
743  real(dp), dimension(0:JMAX,0:IMAX) :: a_vx_b_g
744  real(dp), dimension(0:JMAX,0:IMAX) :: a_vy_b_g
745  real(dp), dimension(0:JMAX,0:IMAX) :: a_vz_b
746  real(dp), dimension(0:JMAX,0:IMAX) :: a_vz_m
747  real(dp), dimension(0:JMAX,0:IMAX) :: a_vx_s_g
748  real(dp), dimension(0:JMAX,0:IMAX) :: a_vy_s_g
749  real(dp), dimension(0:JMAX,0:IMAX) :: a_vz_s
750  real(dp), dimension(0:JMAX,0:IMAX) :: a_flui_ave_sia
751  real(dp), dimension(0:JMAX,0:IMAX) :: a_h_diff
752  real(dp), dimension(0:JMAX,0:IMAX) :: a_qx
753  real(dp), dimension(0:JMAX,0:IMAX) :: a_qy
754  real(dp), dimension(0:JMAX,0:IMAX) :: a_q_gl_g
755  real(dp), dimension(0:JMAX,0:IMAX) :: a_q_geo
756  real(dp), dimension(0:JMAX,0:IMAX) :: a_temp_b
757  real(dp), dimension(0:JMAX,0:IMAX) :: a_temph_b
758  real(dp), dimension(0:JMAX,0:IMAX) :: a_Q_bm
759  real(dp), dimension(0:JMAX,0:IMAX) :: a_Q_b_apl
760  real(dp), dimension(0:JMAX,0:IMAX) :: a_Q_tld
761  real(dp), dimension(0:JMAX,0:IMAX) :: a_Q_b_tot
762  real(dp), dimension(0:JMAX,0:IMAX) :: a_H_w
763  real(dp), dimension(0:JMAX,0:IMAX) :: a_accum
764  real(dp), dimension(0:JMAX,0:IMAX) :: a_evap
765  real(dp), dimension(0:JMAX,0:IMAX) :: a_runoff
766  real(dp), dimension(0:JMAX,0:IMAX) :: a_runoff_apl
767  real(dp), dimension(0:JMAX,0:IMAX) :: a_as_perp
768  real(dp), dimension(0:JMAX,0:IMAX) :: a_as_perp_apl
769  real(dp), dimension(0:JMAX,0:IMAX) :: a_mb_source_apl
770  real(dp), dimension(0:JMAX,0:IMAX) :: a_calv_grounded
771  real(dp), dimension(0:JMAX,0:IMAX) :: a_calv_grounded_apl
772  real(dp), dimension(0:JMAX,0:IMAX) :: a_temp_s
773  real(dp), dimension(0:JMAX,0:IMAX) :: a_am_perp
774  real(dp), dimension(0:JMAX,0:IMAX) :: a_am_perp_st
775  real(dp), dimension(0:JMAX,0:IMAX) :: a_zs_neu
776  real(dp), dimension(0:JMAX,0:IMAX) :: a_zm_neu
777  real(dp), dimension(0:JMAX,0:IMAX) :: a_zb_neu
778  real(dp), dimension(0:JMAX,0:IMAX) :: a_zl_neu
779  real(dp), dimension(0:JMAX,0:IMAX) :: a_H_c_neu
780  real(dp), dimension(0:JMAX,0:IMAX) :: a_H_t_neu
781  real(dp), dimension(0:JMAX,0:IMAX) :: a_zs_ref
782  real(dp), dimension(0:JMAX,0:IMAX) :: a_accum_present
783  real(dp), dimension(0:JMAX,0:IMAX) :: a_precip_ma_present
784  real(dp), dimension(0:JMAX,0:IMAX) :: a_precip_ma_lgm_anom
785  real(dp), dimension(0:JMAX,0:IMAX) :: a_temp_ma_present
786  real(dp), dimension(0:JMAX,0:IMAX) :: a_temp_mj_present
787  real(dp), dimension(0:JMAX,0:IMAX) :: a_temp_ma_lgm_anom
788  real(dp), dimension(0:JMAX,0:IMAX) :: a_temp_mj_lgm_anom
789  real(dp), dimension(0:JMAX,0:IMAX) :: a_q_imp
790  real(dp), dimension(-JMAX:JMAX,-IMAX:IMAX) :: a_dist_dxdy
791  real(dp), dimension(0:JMAX,0:IMAX,12) :: a_precip_present
792  real(dp), dimension(0:JMAX,0:IMAX) :: a_acc_fact
793  real(dp), dimension(0:JMAX,0:IMAX,12) :: a_precip_lgm_anom
794  real(dp), dimension(0:JMAX,0:IMAX,12) :: a_gamma_precip_lgm_anom
795  real(dp), dimension(0:JMAX,0:IMAX,12) :: a_temp_mm_present
796  real(dp), dimension(0:JMAX,0:IMAX,12) :: a_temp_mm_lgm_anom
797  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_d_help_c
798  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_vx_c
799  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_vy_c
800  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_vz_c
801  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_temp_c
802  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_temp_c_neu
803  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_temp_c_m
804  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_age_c
805  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_age_c_neu
806  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_txz_c
807  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_tyz_c
808  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_sigma_c
809  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_enh_c
810  real(dp), dimension(0:JMAX,0:IMAX) :: a_de_ssa
811  real(dp), dimension(0:JMAX,0:IMAX) :: a_vis_int_g
812  real(dp), dimension(0:JMAX,0:IMAX) :: a_vx_g
813  real(dp), dimension(0:JMAX,0:IMAX) :: a_vy_g
814  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_d_help_t
815  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_vx_t
816  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_vy_t
817  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_vz_t
818  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_omega_t
819  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_omega_t_neu
820  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_temp_t_m
821  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_age_t
822  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_age_t_neu
823  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_txz_t
824  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_tyz_t
825  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_sigma_t
826  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_enh_t
827  real(dp), dimension(0:KRMAX,0:JMAX,0:IMAX) :: a_temp_r
828  real(dp), dimension(0:KRMAX,0:JMAX,0:IMAX) :: a_temp_r_neu
829  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_enth_c
830  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_enth_c_neu
831  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_omega_c
832  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_omega_c_neu
833  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_enth_t
834  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_enth_t_neu
835  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_dxx_c
836  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_dyy_c
837  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_dxy_c
838  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_dxz_c
839  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_dyz_c
840  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_de_c
841  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_lambda_shear_c
842  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_dxx_t
843  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_dyy_t
844  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_dxy_t
845  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_dxz_t
846  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_dyz_t
847  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_de_t
848  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: a_lambda_shear_t
849  real(dp) :: a_RHO
850  real(dp) :: a_RHO_W
851  real(dp) :: a_RHO_SW
852  real(dp) :: a_L
853  real(dp) :: a_G
854  real(dp) :: a_NUE
855  real(dp) :: a_BETA
856  real(dp) :: a_DELTA_TM_SW
857  real(dp) :: a_OMEGA_MAX
858  real(dp) :: a_H_R
859  real(dp) :: a_RHO_C_R
860  real(dp) :: a_KAPPA_R
861  real(dp) :: a_RHO_A
862  real(dp) :: a_R_T
863  real(dp) :: a_R_T_IMQ
864  real(dp) :: a_R
865  real(dp) :: a_A
866  real(dp) :: a_B
867  real(dp) :: a_LAMBDA0
868  real(dp) :: a_PHI0
869  real(dp) :: a_S_STAT_0
870  real(dp) :: a_BETA1_0
871  real(dp) :: a_BETA1_LT_0
872  real(dp) :: a_BETA1_HT_0
873  real(dp) :: a_BETA2_0
874  real(dp) :: a_BETA2_LT_0
875  real(dp) :: a_BETA2_HT_0
876  real(dp) :: a_PHI_SEP_0
877  real(dp) :: a_PMAX_0
878  real(dp) :: a_MU_0
879  real(dp), dimension(-190:10) :: a_RF
880  real(dp), dimension(-256:255) :: a_RF_IMQ
881  real(dp), dimension(-256:255) :: a_KAPPA_IMQ
882  real(dp), dimension(-256:255) :: a_C_IMQ
883  real(dp) :: a_year_zero
884  integer(i1b) :: a_forcing_flag
885  integer(i4b) :: a_n_core
886  real(dp), dimension(a_n_core):: a_lambda_core
887  real(dp), dimension(a_n_core) :: a_phi_core
888  real(dp), dimension(a_n_core) :: a_x_core
889  real(dp), dimension(a_n_core) :: a_y_core
890  integer(i4b) :: a_grip_time_min
891  integer(i4b) :: a_grip_time_stp
892  integer(i4b) :: a_grip_time_max
893  integer(i4b) :: a_ndata_grip
894  real(dp), dimension(0:a_ndata_grip) :: a_griptemp
895  integer(i4b) :: a_gi_time_min
896  integer(i4b) :: a_gi_time_stp
897  integer(i4b) :: a_gi_time_max
898  integer(i4b) :: a_ndata_gi
899  real(dp), dimension(0:a_ndata_gi) :: a_glacial_index
900  integer(i4b) :: a_specmap_time_min
901  integer(i4b) :: a_specmap_time_stp
902  integer(i4b) :: a_specmap_time_max
903  integer(i4b) :: a_ndata_specmap
904  real(dp), dimension(0:a_ndata_specmap) :: a_specmap_zsl
905  real(dp) :: a_target_topo_tau
906  real(dp) :: a_time_target_topo_init
907  real(dp) :: a_time_target_topo_final
908  integer(i1b), dimension(0:JMAX,0:IMAX) :: a_maske_target
909  real(dp), dimension(0:JMAX,0:IMAX) :: a_zs_target
910  real(dp), dimension(0:JMAX,0:IMAX) :: a_zb_target
911  real(dp), dimension(0:JMAX,0:IMAX) :: a_zl_target
912  real(dp), dimension(0:JMAX,0:IMAX) :: a_H_target
913  integer(i1b), dimension(0:JMAX,0:IMAX) :: a_maske_maxextent
914  integer(i1b), dimension(0:JMAX,0:IMAX) :: a_mask_ablation_type
915  integer(i4b) :: a_ncid_ser
916  integer(i4b) :: a_ncid_core
917  integer(i4b) :: a_ncid_temp_precip
918  integer(i4b) :: a_ndata_temp_precip
919  real(dp) :: a_temp_precip_time_min
920  real(dp) :: a_temp_precip_time_stp
921  real(dp) :: a_temp_precip_time_max
922  real(dp), dimension(0:a_ndata_temp_precip) :: a_temp_precip_time
923  real(dp), dimension(-10000:10000) :: a_kei
924  integer(i4b) :: a_n_data_kei
925  real(dp) :: a_kei_r_max
926  real(dp) :: a_kei_r_incr
927  real(dp) :: a_fc
928  real(dp) :: a_objf_test
929  real(dp) :: a_mult_test
930  real(dp), dimension(-256:255) :: a_c_int_table
931  real(dp), dimension(-524288:524287) :: a_c_int_inv_table
932  real(dp), dimension(-190:10) :: a_c
933  integer(i4b) :: a_n_temp_min
934  integer(i4b) :: a_n_temp_max
935  integer(i4b) :: a_n_enth_min
936  integer(i4b) :: a_n_enth_max
937  real(dp) :: a_L_inv
938  real(dp) :: a_L_eto
939  real(dp) :: a_RHO_I_IMQ
940  real(dp) :: a_RHO_C_IMQ
941  real(dp), dimension(-190:10) :: a_KAPPA
942  real(dp) :: a_KAPPA_C_IMQ
943  real(dp) :: a_C_C_IMQ
944 #if (defined(INITMIP_SMB_ANOM_FILE))
945  real(dp), dimension(0:JMAX,0:IMAX) :: a_as_anom_initmip
946 #endif
947 #if (defined(INITMIP_BMB_ANOM_FILE))
948  real(dp), dimension(0:JMAX,0:IMAX) :: a_ab_anom_initmip
949 #endif
950 #if (ICE_STREAM==2) /* with ice streams */
951  integer(i1b), dimension(0:JMAX,0:IMAX) :: a_maske_sedi
952 #endif
953 #if (DISC==2)
954  integer(i4b) :: a_glann_time_min
955  integer(i4b) :: a_glann_time_stp
956  integer(i4b) :: a_glann_time_max
957  integer(i4b) :: a_ndata_glann
958  real(dp), dimension(:), allocatable :: a_dT_glann_CLIMBER
959 #endif
960 #ifdef GRL
961  integer(i4b) :: a_disc_DW
962  integer(i4b) :: a_n_discharge_call_DW
963  integer(i4b) :: a_iter_mar_coa_DW
964  real(dp) :: a_c_dis_0_DW
965  real(dp) :: a_s_dis_DW
966  real(dp) :: a_c_dis_fac_DW
967  real(dp) :: a_T_sub_PD_DW
968  real(dp) :: a_alpha_sub_DW
969  real(dp) :: a_alpha_o_DW
970  real(dp) :: a_m_H_DW
971  real(dp) :: a_m_D_DW
972  real(dp) :: a_r_mar_eff_DW
973  real(dp) :: a_T_sea_freeze_DW
974  real(dp), dimension(0:JMAX,0:IMAX) :: a_c_dis_DW
975 #endif
976 #if (defined(AGE_COST))
977 #if (CALCMOD==1)
978  real(dp), dimension(0:KTMAX+KCMAX+1,0:JMAX,0:IMAX) :: a_age_data
979  real(dp), dimension(0:KTMAX+KCMAX+1,0:JMAX,0:IMAX) :: a_age_unc
980 #else
981  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_age_data
982  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: a_age_unc
983 #endif
984 #endif /* No age cost used */
985  integer(i4b) :: a_n_temp_min_IMQ
986  integer(i4b) :: a_n_temp_max_IMQ
987 #ifdef ANT
988  integer(i4b), dimension(0:JMAX,0:IMAX) :: a_n_sector
989 #endif
990  real(dp), dimension(0:JMAX,0:IMAX) :: a_smb_corr_prescribed
991  real(dp), dimension(0:JMAX,0:IMAX) :: a_smb_corr
992  real(dp), dimension(0:JMAX,0:IMAX) :: a_melt
993  real(dp), dimension(0:JMAX,0:IMAX) :: a_melt_star
994  real(dp), dimension(0:JMAX,0:IMAX) :: a_rainfall
995  real(dp), dimension(0:JMAX,0:IMAX) :: a_snowfall
996  real(dp), dimension(0:JMAX,0:IMAX) :: a_et
997 
998  a = a_a
999  aa = a_aa
1000 #if (defined(INITMIP_BMB_ANOM_FILE))
1001  ab_anom_initmip = a_ab_anom_initmip
1002 #endif
1003  accum%v = a_accum
1004  accum_present = a_accum_present
1005  acc_fact = a_acc_fact
1006  age_c%v = a_age_c
1007  age_c_neu%v = a_age_c_neu
1008 #if (defined(AGE_COST))
1009  age_data = a_age_data
1010  age_unc = a_age_unc
1011 #endif
1012 #if (defined(AGE_COST) && CALCMOD==1)
1013  age_t%v = a_age_t
1014  age_t_neu%v = a_age_t_neu
1015 #else
1016  age_t = a_age_t
1017  age_t_neu = a_age_t_neu
1018 #endif
1019 #ifdef GRL
1020  alpha_o_dw = a_alpha_o_dw !ns
1021  alpha_sub_dw = a_alpha_sub_dw !ns
1022 #endif
1023  am_perp = a_am_perp
1024  am_perp_st = a_am_perp_st
1025  area = a_area
1026  as_perp%v = a_as_perp
1027  as_perp_apl = a_as_perp_apl
1028 #if (defined(INITMIP_SMB_ANOM_FILE))
1029  as_anom_initmip = a_as_anom_initmip
1030 #endif
1031  b = a_b
1032  beta = a_beta
1033  beta1_0 = a_beta1_0
1034  beta1_ht_0 = a_beta1_ht_0
1035  beta1_lt_0 = a_beta1_lt_0
1036  beta2_0 = a_beta2_0
1037  beta2_ht_0 = a_beta2_ht_0
1038  beta2_lt_0 = a_beta2_lt_0
1039  c = a_c
1040  calv_grounded%v = a_calv_grounded
1041  calv_grounded_apl = a_calv_grounded_apl
1042  c_imq = a_c_imq !ns
1043  c_c_imq = a_c_c_imq !ns
1044  c_int_table = a_c_int_table !ns
1045  c_int_inv_table = a_c_int_inv_table !ns
1046 #ifdef GRL
1047  c_dis_dw = a_c_dis_dw !ns
1048  c_dis_0_dw = a_c_dis_0_dw !ns
1049 #endif
1050 #if (DYNAMICS==2 || MARGIN==3)
1051  c_drag%v = a_c_drag
1052 #else
1053  c_drag = a_c_drag
1054 #endif
1055 #ifdef GRL
1056  c_slide%v = a_c_slide
1057 #elif (defined(ANT))
1058  c_slide = a_c_slide
1059 #endif
1060  delta_tm_sw = a_delta_tm_sw
1061 #if (ENHMOD == 5 || DYNAMICS==2)
1062  de_c%v = a_de_c
1063 #else
1064  de_c = a_de_c
1065 #endif
1066 #if (MARGIN==3 || DYNAMICS==2)
1067  de_ssa%v = a_de_ssa
1068 #else
1069  de_ssa = a_de_ssa
1070 #endif
1071  de_t = a_de_t
1072  dh_c_deta = a_dh_c_deta
1073  dh_c_deta_g%v = a_dh_c_deta_g
1074  dh_c_dtau%v = a_dh_c_dtau
1075  dh_c_dxi = a_dh_c_dxi
1076  dh_c_dxi_g%v = a_dh_c_dxi_g
1077  dh_t_deta = a_dh_t_deta
1078  dh_t_deta_g%v = a_dh_t_deta_g
1079  dh_t_dtau%v = a_dh_t_dtau
1080  dh_t_dxi = a_dh_t_dxi
1081  dh_t_dxi_g%v = a_dh_t_dxi_g
1082 #ifdef GRL
1083  disc_dw = a_disc_dw !ns
1084 #endif
1085  dist_dxdy = a_dist_dxdy
1086 #if (CALCMOD == 5 || ENHMOD == 5 || DYNAMICS==2)
1087  dxx_c%v = a_dxx_c
1088  dxy_c%v = a_dxy_c
1089  dxz_c%v = a_dxz_c
1090  dyy_c%v = a_dyy_c
1091  dyz_c%v = a_dyz_c
1092 #else
1093  dxx_c = a_dxx_c
1094  dxy_c = a_dxy_c
1095  dxz_c = a_dxz_c
1096  dyy_c = a_dyy_c
1097  dyz_c = a_dyz_c
1098 #endif
1099  dxx_t = a_dxx_t
1100  dxy_t = a_dxy_t
1101  dxz_t = a_dxz_t
1102  dyy_t = a_dyy_t
1103  dyz_t = a_dyz_t
1104  dzb_deta = a_dzb_deta
1105  dzb_deta_g%v = a_dzb_deta_g
1106  dzb_dtau%v = a_dzb_dtau
1107  dzb_dxi = a_dzb_dxi
1108  dzb_dxi_g%v = a_dzb_dxi_g
1109 #if (REBOUND >= 1)
1110  dzl_dtau%v = a_dzl_dtau
1111 #else
1112  dzl_dtau = a_dzl_dtau
1113 #endif
1114  dzm_deta = a_dzm_deta
1115  dzm_deta_g%v = a_dzm_deta_g
1116  dzm_dtau%v = a_dzm_dtau
1117  dzm_dxi = a_dzm_dxi
1118  dzm_dxi_g%v = a_dzm_dxi_g
1119  dzs_deta%v = a_dzs_deta
1120  dzs_deta_g%v = a_dzs_deta_g
1121  dzs_dtau%v = a_dzs_dtau
1122  dzs_dxi%v = a_dzs_dxi
1123  dzs_dxi_g%v = a_dzs_dxi_g
1124  d_help_b%v = a_d_help_b
1125  d_help_c%v = a_d_help_c
1126  d_help_t%v = a_d_help_t
1127  ea = a_ea
1128  eaz_c = a_eaz_c
1129  eaz_c_quotient = a_eaz_c_quotient
1130 #if (ENHMOD == 5)
1131  enh_c%v = a_enh_c
1132  enh_t%v = a_enh_t
1133 #else
1134  enh_c = a_enh_c
1135  enh_t = a_enh_t
1136 #endif
1137 #if (CALCMOD >= 2)
1138  enth_c%v = a_enth_c
1139  enth_c_neu%v = a_enth_c_neu
1140 #else
1141  enth_c = a_enth_c
1142  enth_c_neu = a_enth_c_neu
1143 #endif
1144  enth_t = a_enth_t
1145  enth_t_neu = a_enth_t_neu
1146  et%v = a_et
1147  eta = a_eta
1148  evap = a_evap
1149  fc%v = a_fc
1150  flag_aa_nonzero = a_flag_aa_nonzero
1151  flag_calving_front_1 = a_flag_calving_front_1
1152  flag_calving_front_2 = a_flag_calving_front_2
1153  flag_grounding_line_1 = a_flag_grounding_line_1
1154  flag_grounding_line_2 = a_flag_grounding_line_2
1155  flag_inner_point = a_flag_inner_point
1156  flag_shelfy_stream = a_flag_shelfy_stream
1157  flag_shelfy_stream_x = a_flag_shelfy_stream_x
1158  flag_shelfy_stream_y = a_flag_shelfy_stream_y
1159  flex_rig_lith = a_flex_rig_lith
1160 #if (DYNAMICS==2 || MARGIN==3)
1161  flui_ave_sia%v = a_flui_ave_sia
1162 #else
1163  flui_ave_sia = a_flui_ave_sia
1164 #endif
1165  forcing_flag = a_forcing_flag
1166  g = a_g
1167  gamma_precip_lgm_anom = a_gamma_precip_lgm_anom
1168  gi_time_min = a_gi_time_min
1169  gi_time_max = a_gi_time_max
1170  gi_time_stp = a_gi_time_stp
1171 #if TSURFACE==5
1172  allocate(glacial_index(0:a_ndata_gi))
1173  glacial_index = a_glacial_index
1174 #endif
1175 #if TSURFACE==4
1176  allocate(griptemp(0:a_ndata_grip))
1177  griptemp = a_griptemp
1178 #endif
1179 #if (DISC==2)
1180  glann_time_min = a_glann_time_min
1181  glann_time_stp = a_glann_time_stp
1182  glann_time_max = a_glann_time_max
1183  ndata_glann = a_ndata_glann
1184  dt_glann_climber = a_dt_glann_climber
1185 #endif
1186  h_c%v = a_h_c
1187  h_c_neu%v = a_h_c_neu
1188  h_diff%v = a_h_diff
1189  h_r = a_h_r
1190  h_t%v = a_h_t
1191  h_target = a_h_target
1192  h_t_neu%v = a_h_t_neu
1193  h_w = a_h_w
1194  insq_g11_g = a_insq_g11_g
1195  insq_g11_sgx = a_insq_g11_sgx
1196  insq_g22_g = a_insq_g22_g
1197  insq_g22_sgy = a_insq_g22_sgy
1198 #ifdef GRL
1199  iter_mar_coa_dw = a_iter_mar_coa_dw !ns
1200 #endif
1201  kappa = a_kappa
1202  kappa_r = a_kappa_r
1203  kappa_c_imq = a_kappa_c_imq
1204  kappa_imq = a_kappa_imq !ns? called KAPPA, no _IMQ
1205  kc_cts = a_kc_cts
1206  kc_cts_neu = a_kc_cts_neu
1207  kei = a_kei
1208  kei_r_incr = a_kei_r_incr
1209  kei_r_max = a_kei_r_max
1210  l = a_l
1211  l_inv = a_l_inv !ns
1212  l_eto = a_l_eto !ns
1213  lambda = a_lambda
1214  lambda0 = a_lambda0
1215 !#if OUTSER==3
1216  allocate(lambda_core(a_n_core))
1217  lambda_core = a_lambda_core ! there but is it working?
1218 !#endif
1219 #if (ENHMOD == 5)
1220  lambda_shear_c%v = a_lambda_shear_c
1221  lambda_shear_t%v = a_lambda_shear_t
1222 #else
1223  lambda_shear_c = a_lambda_shear_c
1224  lambda_shear_t = a_lambda_shear_t
1225 #endif
1226  maske = a_maske
1227  maske_maxextent = a_maske_maxextent
1228  maske_neu = a_maske_neu
1229  maske_old = a_maske_old
1230  maske_target = a_maske_target
1231  mask_ablation_type = a_mask_ablation_type
1232 #if (ICE_STREAM==2) /* with ice streams */
1233  maske_sedi = a_maske_sedi
1234 #endif
1235  mb_source_apl = a_mb_source_apl
1236  melt%v = a_melt
1237  melt_star%v = a_melt_star
1238  mult_test = a_mult_test
1239  mu_0 = a_mu_0
1240 #ifdef GRL
1241  m_d_dw = a_m_d_dw !ns
1242  m_h_dw = a_m_h_dw !ns
1243 #endif
1244  ncid_core = a_ncid_core
1245  ncid_ser = a_ncid_ser
1246  ncid_temp_precip = a_ncid_temp_precip
1247  ndata_gi = a_ndata_gi
1248  ndata_grip = a_ndata_grip
1249  ndata_specmap = a_ndata_specmap
1250  ndata_temp_precip = a_ndata_temp_precip
1251  nue = a_nue
1252  n_core = a_n_core
1253  n_cts = a_n_cts
1254  n_cts_neu = a_n_cts_neu
1255 #ifdef GRL
1256  n_discharge_call_dw = a_n_discharge_call_dw !ns
1257 #endif
1258  n_data_kei = a_n_data_kei !ns
1259  n_enth_min = a_n_enth_min !ns
1260  n_enth_max = a_n_enth_max !ns
1261 #ifdef ANT
1262  n_sector = a_n_sector
1263 #endif
1264  n_temp_min_imq = a_n_temp_min_imq !ns
1265  n_temp_max_imq = a_n_temp_max_imq !ns
1266  objf_test%v = a_objf_test
1267 #if (CALCMOD >= 2)
1268  omega_c%v = a_omega_c
1269  omega_c_neu%v = a_omega_c_neu
1270 #elif (CALCMOD == 0 || CALCMOD == 1)
1271  omega_c = a_omega_c
1272  omega_c_neu = a_omega_c_neu
1273 #endif
1274  omega_max = a_omega_max
1275 #if (CALCMOD >= 1)
1276  omega_t%v = a_omega_t
1277  omega_t_neu%v = a_omega_t_neu
1278 #else
1279  omega_t = a_omega_t
1280  omega_t_neu = a_omega_t_neu
1281 #endif
1282  phi = a_phi
1283  phi0 = a_phi0
1284 !#if OUTSER==3
1285  allocate(phi_core(a_n_core))
1286  phi_core = a_phi_core
1287 !#endif
1288  phi_sep_0 = a_phi_sep_0
1289  pmax_0 = a_pmax_0
1290  precip_lgm_anom = a_precip_lgm_anom
1291  precip_ma_lgm_anom = a_precip_ma_lgm_anom
1292  precip_ma_present = a_precip_ma_present
1293  precip_present%v = a_precip_present
1294 #if (defined(ANT) && SLIDE_LAW==2)
1295  p_b_w%v = a_p_b_w
1296 #elif (SLIDE_LAW==1)
1297  p_b_w = a_p_b_w
1298 #endif
1299  p_weert = a_p_weert
1300  p_weert_inv = a_p_weert_inv
1301 #if (ENHMOD == 5 || DYNAMICS==2 || MARGIN==3)
1302  qx%v = a_qx
1303  qy%v = a_qy
1304 #else
1305  qx = a_qx
1306  qy = a_qy
1307 #endif
1308  q_bm%v = a_q_bm
1309  q_b_apl = a_q_b_apl
1310  q_b_tot%v = a_q_b_tot
1311  q_geo%v = a_q_geo
1312  q_gl_g = a_q_gl_g
1313 #if (CALCMOD >= 0)
1314  q_tld%v = a_q_tld
1315 #else
1316  q_tld = a_q_tld
1317 #endif
1318  q_weert = a_q_weert
1319  r = a_r
1320  rainfall%v = a_rainfall
1321 #if (DYNAMICS==2 || MARGIN==3)
1322  ratio_sl_x%v = a_ratio_sl_x
1323  ratio_sl_y%v = a_ratio_sl_y
1324 #else
1325  ratio_sl_x = a_ratio_sl_x
1326  ratio_sl_y = a_ratio_sl_y
1327 #endif
1328  rf = a_rf
1329  rf_imq = a_rf_imq !just called RF in numcore
1330  rho = a_rho
1331  rho_a = a_rho_a
1332  rho_c_imq = a_rho_c_imq
1333  rho_c_r = a_rho_c_r
1334  rho_i_imq = a_rho_i_imq !ns
1335  rho_sw = a_rho_sw
1336  rho_w = a_rho_w
1337  runoff%v = a_runoff
1338  runoff_apl = a_runoff_apl
1339 #ifdef GRL
1340  r_mar_eff_dw = a_r_mar_eff_dw !ns
1341 #endif
1342  r_t = a_r_t ! just called R_T
1343  r_t_imq = a_r_t_imq ! just called R_T
1344  sigma_c%v = a_sigma_c
1345  sigma_t%v = a_sigma_t
1346  smb_corr = a_smb_corr
1347  smb_corr_prescribed = a_smb_corr_prescribed
1348  snowfall%v = a_snowfall
1349  specmap_time_max = a_specmap_time_max
1350  specmap_time_min = a_specmap_time_min
1351  specmap_time_stp = a_specmap_time_stp
1352 #if (SEA_LEVEL==3)
1353  allocate(specmap_zsl(0:a_ndata_specmap))
1354  specmap_zsl = a_specmap_zsl !present but does header use this?
1355 #endif
1356  sq_g11_g = a_sq_g11_g
1357  sq_g11_sgx = a_sq_g11_sgx
1358  sq_g11_sgy = a_sq_g11_sgy
1359  sq_g22_g = a_sq_g22_g
1360  sq_g22_sgx = a_sq_g22_sgx
1361  sq_g22_sgy = a_sq_g22_sgy
1362 #ifdef GRL
1363  s_dis_dw = a_s_dis_dw !ns
1364  c_dis_fac_dw = a_c_dis_fac_dw !ns
1365 #endif
1366  s_stat_0 = a_s_stat_0
1367  target_topo_tau = a_target_topo_tau
1368  temph_b = a_temph_b
1369  temp_b = a_temp_b
1370  temp_c%v = a_temp_c
1371  temp_c_m%v = a_temp_c_m
1372  temp_c_neu%v = a_temp_c_neu
1373  temp_ma_lgm_anom = a_temp_ma_lgm_anom
1374  temp_ma_present%v = a_temp_ma_present
1375  temp_mj_lgm_anom = a_temp_mj_lgm_anom
1376  temp_mj_present%v = a_temp_mj_present
1377  temp_mm_lgm_anom = a_temp_mm_lgm_anom
1378  temp_mm_present = a_temp_mm_present
1379 #if ( (TSURFACE==6) && (ACCSURFACE==6) )
1380  allocate(temp_precip_time(0:a_ndata_temp_precip))
1381  temp_precip_time = a_temp_precip_time ! in there but does header file nee
1382 #endif
1383  temp_precip_time_max = a_temp_precip_time_max
1384  temp_precip_time_min = a_temp_precip_time_min
1385  temp_precip_time_stp = a_temp_precip_time_stp
1386  temp_r%v = a_temp_r
1387  temp_r_neu%v = a_temp_r_neu
1388  temp_s%v = a_temp_s
1389  temp_t_m%v = a_temp_t_m
1390  time_lag_asth = a_time_lag_asth
1391  time_target_topo_final = a_time_target_topo_final
1392  time_target_topo_init = a_time_target_topo_init
1393 #ifdef GRL
1394  t_sub_pd_dw = a_t_sub_pd_dw !ns
1395  t_sea_freeze_dw = a_t_sea_freeze_dw !ns
1396 #endif
1397  txz_c = a_txz_c
1398  txz_t = a_txz_t
1399  tyz_c = a_tyz_c
1400  tyz_t = a_tyz_t
1401 #if (DYNAMICS==2 || MARGIN==3)
1402  vis_int_g%v = a_vis_int_g
1403 #else
1404  vis_int_g = a_vis_int_g
1405 #endif
1406 #if (DYNAMICS==2 || MARGIN==3)
1407  vx_b%v = a_vx_b
1408  vx_b_g%v = a_vx_b_g
1409 #else
1410  vx_b = a_vx_b
1411  vx_b_g = a_vx_b_g
1412 #endif
1413  vx_c%v = a_vx_c
1414  vx_g = a_vx_g
1415 #if (ENHMOD == 5 || DYNAMICS==2 || MARGIN==3)
1416  vx_m%v = a_vx_m
1417  vy_m%v = a_vy_m
1418 #else
1419  vx_m = a_vx_m
1420  vy_m = a_vy_m
1421 #endif
1422 #if (DYNAMICS==2 || MARGIN==3)
1423  vx_s_g%v = a_vx_s_g
1424 #else
1425  vx_s_g = a_vx_s_g
1426 #endif
1427  vx_t%v = a_vx_t
1428 #if (DYNAMICS==2 || MARGIN==3)
1429  vy_b%v = a_vy_b
1430  vy_b_g%v = a_vy_b_g
1431 #else
1432  vy_b = a_vy_b
1433  vy_b_g = a_vy_b_g
1434 #endif
1435  vy_c%v = a_vy_c
1436  vy_g = a_vy_g
1437 #if (DYNAMICS==2 || MARGIN==3)
1438  vy_s_g%v = a_vy_s_g
1439 #else
1440  vy_s_g = a_vy_s_g
1441 #endif
1442  vy_t%v = a_vy_t
1443  vz_b%v = a_vz_b
1444  vz_c%v = a_vz_c
1445  vz_m%v = a_vz_m
1446 #if (ENHMOD == 5 || DYNAMICS==2)
1447  vz_s%v = a_vz_s
1448 #else
1449  vz_s = a_vz_s
1450 #endif
1451  vz_t%v = a_vz_t
1452 #if (REBOUND == 2)
1453  wss%v = a_wss
1454 #else
1455  wss = a_wss
1456 #endif
1457  xi = a_xi
1458 #if OUTSER==3
1459  allocate(x_core(a_n_core))
1460  allocate(y_core(a_n_core))
1461  x_core = a_x_core
1462  y_core = a_y_core
1463 #endif
1464  year_zero = a_year_zero
1465  zb%v = a_zb
1466  zb_neu%v = a_zb_neu
1467  zb_target = a_zb_target
1468  zeta_c = a_zeta_c
1469  zeta_r = a_zeta_r
1470  zeta_t = a_zeta_t
1471 #if (REBOUND >= 1)
1472  zl%v = a_zl
1473  zl_neu%v = a_zl_neu
1474 #else
1475  zl = a_zl
1476  zl_neu = a_zl_neu
1477 #endif
1478  zl0 = a_zl0
1479  zl_target = a_zl_target
1480  zm%v = a_zm
1481  zm_neu%v = a_zm_neu
1482  zs%v = a_zs
1483  zs_neu%v = a_zs_neu
1484  zs_ref = a_zs_ref
1485  zs_target = a_zs_target
1486 
1487  end subroutine var_transfer
1488 
1489 !!-------------------------------------------------------------------------------
1490 !!> This is another absolutely essential subroutine that is responsible for the
1491 !! writing of adjoint sensitivities to files for comparison against the gradient
1492 !! check as well as for other more science-targeted uses. As yet, it is rather
1493 !! ad-hoc and is structured for user-supplied controls. Every output file below is
1494 !! supplied and named directly by the user, and if for some reason the user runs
1495 !! an adjoint simulation where one of the controls is no longer ACTIVE, then
1496 !! errors may occur if syntaxes below are not proper.
1497 !! For instance, if foobar is a control that is ACTIVE:
1498 !! in simulation (1) but not (2), then the proper syntax
1499 !! for the outputting of sensitivities (ignoring that outputting sensitivities for
1500 !! non-ACTIVE variables is a moot exercise) is:
1501 !! Simulation (1)
1502 !! "print stuff as below" foobar(j,i,k)%d
1503 !! Simulation (2)
1504 !! "print stuff as below" foobar(j,i,k)
1505 !! Or, just comment out foobar from the writing to file entirely since all that
1506 !! will be printed is zeros anyway. And again, write to
1507 !! liz.curry.logan@gmail.com with error messages.
1508 !!<------------------------------------------------------------------------------
1509  subroutine print_output(runname)
1510  use oad_active
1511  use oad_sico_variables_m
1512 
1513  implicit none
1514 
1515  integer(i4b), parameter :: points = 10
1516  integer(i4b), dimension(points) :: ipoints, jpoints
1517  integer :: i, j, p
1518  character(len=100), intent(out) :: runname
1519 
1520  !-------- Open output files:
1521  ! this file is for comparison with the gradient check vals:
1522  open(unit=96,file='AD_Vals_for_grdchk.dat',status='replace')
1523  open(unit=95,file='AD_Vals_H_c.dat',status='replace')
1524  open(unit=94,file='AD_Vals_temp_c_surf.dat',status='replace')
1525  open(unit=93,file='AD_Vals_temp_c_base.dat',status='replace')
1526  open(unit=90,file='AD_Vals_q_geo.dat',status='replace')
1527  open(unit=88,file='AD_Vals_c_drag.dat',status='replace')
1528  open(unit=87,file='AD_Vals_precip_present_july.dat',status='replace')
1529  !open(unit=87,file='AD_Vals_precip_present_ma.dat',status='replace')
1530  open(unit=79,file='AD_Vals_acc_fact.dat',status='replace')
1531  open(unit=86,file='AD_Vals_vis_int_g.dat',status='replace')
1532  open(unit=85,file='AD_Vals_dzs_dxi_g.dat',status='replace')
1533  open(unit=84,file='AD_Vals_dzs_deta_g.dat',status='replace')
1534  open(unit=83,file='AD_Vals_calv_grounded.dat',status='replace')
1535  open(unit=82,file='AD_Vals_basal_stress.dat',status='replace')
1536  open(unit=81,file='AD_Vals_Q_bm.dat',status='replace')
1537  open(unit=80,file='AD_Vals_Q_tld.dat',status='replace')
1538 
1539  !-------- Outputting ad values for gradient-check-comparison:
1540  write(96, *) " (j,i) control(j,i)"
1541 
1542  !-------- For Greenland we're hand-picking the test points:
1543 #if (defined(GRL))
1544  do p = 1, points
1545  ipoints(p) = int(real(imax/2))
1546  jpoints(p) = int(real(jmax/5)) + (p-1) * points
1547  !write(96, '(f40.20)') c_drag(jpoints(p),ipoints(p))%d
1548  write(96, '(f40.20)') vx_c(kcmax,jpoints(p),ipoints(p))%d
1549  !write(96, '(f40.20)') c_slide(jpoints(p),ipoints(p))%d
1550  !write(96, '(f40.20)') q_geo(jpoints(p),ipoints(p))%d
1551  !write(96, '(f40.20)') temp_ma_present(jpoints(p),ipoints(p))%d
1552  !write(96, '(f40.20)') temp_c(KCMAX,jpoints(p),ipoints(p))%d
1553  !write(96, '(f40.20)') H_c(jpoints(p),ipoints(p))%d
1554  !write(96, '(f40.20)') precip_present(jpoints(p),ipoints(p),1)%d
1555  !write(96, '(f40.20)') acc_fact(jpoints(p),ipoints(p))%d
1556  end do
1557  !-------- For Antarctica it's a line halfway through the domain:
1558 #elif (defined(ANT))
1559  ipoints = (/94, 240/) ! tsurf = temp_c(:,:,KCMAX)
1560  jpoints = (/187, 66/) ! tsurf = temp_c(:,:,KCMAX)
1561  do p = 1, points
1562  !ipoints(p) = int(real(IMAX/3)) + int(real((.85-.33)*IMAX/points)) * (p - 1)
1563  !jpoints(p) = int(real(JMAX/2))
1564  write(96, '(f40.20)') h_c(jpoints(p),ipoints(p))%d
1565  end do
1566 #else
1567  print *, "Adjoint mode only for Antarctica and Greenland right now."
1568 #endif
1569 
1570  !-------- Outputting ALL sensitivities to file (entire 2D fields):
1571  do j=0,jmax
1572  do i=0,imax
1573 
1574  if ( isnan(h_c(j,i)%d) ) then
1575  write(unit=95,fmt='(t1,en18.6e4)',advance='no') -66.6
1576  else
1577  write(unit=95,fmt='(t1,en18.6e4)',advance='no') h_c(j,i)%d
1578  end if
1579 
1580  if ( isnan(temp_c(kcmax,j,i)%d) ) then
1581  write(unit=94,fmt='(t1,en18.6e4)',advance='no') -66.6
1582  else
1583  write(unit=94,fmt='(t1,en18.6e4)',advance='no') temp_c(kcmax,j,i)%d
1584  end if
1585 
1586  if ( isnan(temp_c(0,j,i)%d) ) then
1587  write(unit=93,fmt='(t1,en18.6e4)',advance='no') -66.6
1588  else
1589  write(unit=93,fmt='(t1,en18.6e4)',advance='no') temp_c(0,j,i)%d
1590  end if
1591 
1592  if ( isnan(q_geo(j,i)%d) ) then
1593  write(unit=90,fmt='(t1,en18.6e4)',advance='no') -66.6
1594  else
1595  write(unit=90,fmt='(t1,en18.6e4)',advance='no') q_geo(j,i)%d
1596  end if
1597 
1598 #if (defined(ANT))
1599  if ( isnan(c_drag(j,i)%d) ) then
1600  write(unit=88,fmt='(t1,en18.6e4)',advance='no') -66.6
1601  else
1602  write(unit=88,fmt='(t1,en18.6e4)',advance='no') c_drag(j,i)%d
1603  end if
1604 #endif
1605 
1606  if ( isnan(precip_present(j,i,7)%d) ) then
1607  write(unit=87,fmt='(t1,en18.6e4)',advance='no') -66.6
1608  else
1609  write(unit=87,fmt='(t1,en18.6e4)',advance='no') precip_present(j,i,7)%d
1610  end if
1611 
1612  if ( isnan(temp_ma_present(j,i)%d) ) then
1613  write(unit=93,fmt='(t1,en18.6e4)',advance='no') -66.6
1614  else
1615  write(unit=93,fmt='(t1,en18.6e4)',advance='no') temp_ma_present(j,i)%d
1616  end if
1617 
1618  ! if ( isnan(acc_fact(j,i)%d) ) then
1619  ! write(unit=79,fmt='(t1,en18.6e4)',advance='no') -66.6
1620  ! else
1621  ! write(unit=79,fmt='(t1,en18.6e4)',advance='no') acc_fact(j,i)%d
1622  ! end if
1623 
1624 
1625 #if (defined(ANT))
1626  if ( isnan(vis_int_g(j,i)%d) ) then
1627  write(unit=86,fmt='(t1,en18.6e4)',advance='no') -66.6
1628  else
1629  write(unit=86,fmt='(t1,en18.6e4)',advance='no') vis_int_g(j,i)%d
1630  end if
1631 #endif
1632 
1633  if ( isnan(dzs_dxi_g(j,i)%d) ) then
1634  write(unit=85,fmt='(t1,en18.6e4)',advance='no') -66.6
1635  else
1636  write(unit=85,fmt='(t1,en18.6e4)',advance='no') dzs_dxi_g(j,i)%d
1637  end if
1638 
1639  if ( isnan(dzs_deta_g(j,i)%d) ) then
1640  write(unit=84,fmt='(t1,en18.6e4)',advance='no') -66.6
1641  else
1642  write(unit=84,fmt='(t1,en18.6e4)',advance='no') dzs_deta_g(j,i)%d
1643  end if
1644 
1645 #ifdef GRL
1646  if ( isnan(calv_grounded(j,i)%d) ) then
1647  write(unit=83,fmt='(t1,en18.6e4)',advance='no') -66.6
1648  else
1649  write(unit=83,fmt='(t1,en18.6e4)',advance='no') calv_grounded(j,i)%d
1650  end if
1651 #endif
1652 
1653  if ( isnan(sigma_c(kcmax,j,i)%d) ) then
1654  write(unit=82,fmt='(t1,en18.6e4)',advance='no') -66
1655  else
1656  write(unit=82,fmt='(t1,en18.6e4)',advance='no') sigma_c(kcmax,j,i)%d
1657  end if
1658 
1659  if ( isnan(q_bm(j,i)%d) ) then
1660  write(unit=81,fmt='(t1,en18.6e4)',advance='no') -66
1661  else
1662  write(unit=81,fmt='(t1,en18.6e4)',advance='no') q_bm(j,i)%d
1663  end if
1664 
1665  if ( isnan(q_tld(j,i)%d) ) then
1666  write(unit=80,fmt='(t1,en18.6e4)',advance='no') -66
1667  else
1668  write(unit=80,fmt='(t1,en18.6e4)',advance='no') q_tld(j,i)%d
1669  end if
1670 
1671  end do
1672  write(unit=95,fmt='(t1)')
1673  write(unit=94,fmt='(t1)')
1674  write(unit=93,fmt='(t1)')
1675  write(unit=90,fmt='(t1)')
1676  write(unit=88,fmt='(t1)')
1677  write(unit=87,fmt='(t1)')
1678  write(unit=86,fmt='(t1)')
1679  write(unit=85,fmt='(t1)')
1680  write(unit=84,fmt='(t1)')
1681  write(unit=83,fmt='(t1)')
1682  write(unit=82,fmt='(t1)')
1683  write(unit=81,fmt='(t1)')
1684  write(unit=80,fmt='(t1)')
1685  end do
1686 
1687  close(unit=96)
1688  close(unit=95)
1689  close(unit=94)
1690  close(unit=93)
1691  close(unit=90)
1692  close(unit=88)
1693  close(unit=87)
1694  close(unit=86)
1695  close(unit=85)
1696  close(unit=84)
1697  close(unit=83)
1698  close(unit=82)
1699  close(unit=81)
1700  close(unit=80)
1701 
1702  end subroutine print_output
1703 
1704  subroutine print_var(tmp)
1705  use oad_active
1706 
1707  implicit none
1708 
1709  real*8 :: tmp
1710 
1711  write(6,fmt='(f40.20)') tmp
1712 
1713  end subroutine print_var
1714 
1715 #endif
1716 
1717 end module openad_m
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vx_c
vx_c(kc,j,i): Velocity in x-direction in the upper (kc) ice domain (at (i+1/2,j,kc)) ...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp) kappa_r
KAPPA_R: Heat conductivity of the lithosphere.
real(dp), dimension(0:jmax, 0:imax) c_slide
c_slide(j,i): Basal sliding coefficient
real(dp) delta_tm_sw
DELTA_TM_SW: Melting point depression of sea water due to its average salinity.
real(dp), dimension(0:jmax, 0:imax) calv_grounded_apl
calv_grounded_apl(j,i): Applied calving rate of grounded ice
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
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) de_c
de_c(kc,j,i): Full effective strain rate in the upper (kc) ice domain
real(dp) phi_sep_0
PHI_SEP_0: Separation latitude for the computation of the degree-day factors beta1 and beta2: Equator...
real(dp), dimension(:), allocatable phi_core
phi_core(n): Geographical latitude of the prescribed borehole positions
logical, dimension(0:jmax, 0:imax) flag_inner_point
flag_inner_point(j,i): Inner-point flag. .true.: inner point, .false.: margin point ...
integer(i4b) grip_time_stp
grip_time_stp: Time step of the data values for the surface temperature anomaly
integer(i4b) ncid_ser
ncid_ser: ID of the NetCDF time-series output file
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) lambda_shear_t
lambda_shear_t(kt,j,i): Shear fraction in the lower (kt) ice domain
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) d_help_c
d_help_c(kc,j,i): Auxiliary quantity for the computation of vx, vy und zs
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) dxy_t
dxy_t(kt,j,i): Strain rate dxy in the lower (kt) ice domain
real(dp), dimension(0:jmax, 0:imax) runoff_apl
runoff_apl(j,i): Applied runoff rate at the ice surface
integer(i4b) ndata_specmap
ndata_specmap: Number of data values for the sea level
real(dp), dimension(0:jmax, 0:imax) insq_g22_sgy
insq_g22_sgy(j,i): Inverse square root of g22, at (i,j+1/2)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) de_t
de_t(kt,j,i): Full effective strain rate in the lower (kt) ice domain
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) dyy_t
dyy_t(kt,j,i): Strain rate dyy in the lower (kt) ice domain
real(dp), dimension(0:jmax, 0:imax) vz_m
vz_m(j,i): Velocity in z-direction at the position z=zm (interface between the upper (kc) and the low...
real(dp) beta2_ht_0
BETA2_HT_0: Degree-day factor for ice at high summer temperatures.
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
real(dp), dimension(0:jmax, 0:imax, 12) temp_mm_present
temp_mm_present(j,i,n): Present-day mean monthly surface temperature
real(dp), dimension(0:jmax, 0:imax) dh_c_deta
dH_c_deta(j,i): Derivative of H_c by eta (at i,j+1/2)
Definition: ctrl_m.F90:9
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
real(dp) beta1_ht_0
BETA1_HT_0: Degree-day factor for snow at high summer temperatures.
real(dp), dimension(0:jmax, 0:imax) vx_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
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:kcmax, 0:jmax, 0:imax) enh_c
enh_c(kc,j,i): Flow enhancement factor in the upper (kc) ice domain
real(dp) rho_w
RHO_W: Density of pure water.
real(dp), dimension(0:jmax, 0:imax) vx_m
vx_m(j,i): Mean (depth-averaged) velocity in x-direction, at (i+1/2,j)
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
real(dp), dimension(0:jmax, 0:imax) mb_source_apl
mb_source_apl(j,i): Applied mass balance source (SMB, BMB, calving)
real(dp) h_r
H_R: Thickness of the modelled lithosphere layer.
real(dp) a
A: Semi-major axis of the planet.
real(dp), dimension(0:jmax, 0:imax) vz_b
vz_b(j,i): Velocity in z-direction at the ice base, at (i,j)
integer(i4b) specmap_time_max
specmap_time_max: Maximum time of the data values for the sea level
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp) temp_precip_time_stp
temp_precip_time_stp: Time step of the surface-temperature and precipitation data ...
real(dp) mu_0
MU_0: Firn-warming correction.
real(dp), dimension(0:jmax, 0:imax) vy_g
vy_g(j,i): Velocity in y-direction of the SSA, at (i,j)
real(dp), dimension(0:jmax, 0:imax) vx_s_g
vx_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
real(dp) temp_precip_time_max
temp_precip_time_max: Maximum time of the surface-temperature and precipitation data ...
real(dp) omega_max
OMEGA_MAX: Threshold value for the water content of temperate ice.
real(dp), dimension(0:jmax, 0:imax) calv_grounded
calv_grounded(j,i): Calving rate of grounded ice
real(dp), dimension(0:jmax, 0:imax) rainfall
rainfall(j,i): Rainfall rate at the ice surface
logical, dimension(0:jmax, 0:imax) flag_grounding_line_1
flag_grounding_line_1(j,i): Grounding line flag. .true.: grounding line point (grounded ice point wit...
real(dp), dimension(0:jmax, 0:imax) insq_g11_sgx
insq_g11_sgx(j,i): Inverse square root of g11, at (i+1/2,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t
enth_t(kt,j,i): Enthalpy in the lower (kt) ice domain
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
real(dp), dimension(0:jmax, 0:imax) vz_s
vz_s(j,i): Velocity in z-direction at the ice surface, at (i,j)
real(dp), dimension(0:jmax, 0:imax) temp_mj_present
temp_mj_present(j,i): Present-day mean summer (northern hemisphere: July, southern hemisphere: Januar...
real(dp), dimension(0:jmax, 0:imax) insq_g22_g
insq_g22_g(j,i): Inverse square root of g22 on grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) lambda
lambda(j,i): Geographic longitude of grid point (i,j)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) dyz_c
dyz_c(kc,j,i): Strain rate dyz in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
real(dp), dimension(0:jmax, 0:imax) p_b_w
p_b_w(j,i): Basal water pressure
Main loop of SICOPOLIS.
integer(i4b) specmap_time_min
specmap_time_min: Minimum time of the data values for the sea level
real(dp), dimension(0:jmax, 0:imax) sq_g22_sgx
sq_g22_sgx(j,i): Square root of g22, at (i+1/2,j)
integer(i1b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
real(dp), dimension(0:jmax, 0:imax) dh_c_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) vy_b
vy_b(j,i): Velocity in y-direction at the ice base, at (i,j+1/2)
real(dp), dimension(0:jmax, 0:imax) dh_c_dxi
dH_c_dxi(j,i): Derivative of H_c by xi (at i+1/2,j)
integer(i1b), dimension(0:jmax, 0:imax) maske_neu
maske_neu(j,i): New value of maske computed during an integration step
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dzs_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(:), allocatable glacial_index
glacial_index(n): Data values for the glacial index
real(dp) s_stat_0
S_STAT_0: Standard deviation of the air termperature for the degree-day model.
real(dp), dimension(0:jmax, 0:imax) zl_target
zl_target(j,i): Target topography (lithosphere surface)
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream_x
flag_shelfy_stream_x(j,i): Shelfy stream flag in x-direction, at (i+1/2,j). .true.: shelfy stream point .false.: otherwise
real(dp) beta2_lt_0
BETA2_LT_0: Degree-day factor for ice at low summer temperatures.
real(dp), dimension(0:jmax, 0:imax) dzm_deta
dzm_deta(j,i): Derivative of zm by eta (at i,j+1/2)
real(dp), dimension(0:jmax, 0:imax) runoff
runoff(j,i): Runoff rate at the ice surface
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) temp_t_m
temp_t_m(kt,j,i): Melting temperature in the lower (kt) ice domain
real(dp), dimension(0:jmax, 0:imax) de_ssa
de_ssa(j,i): Effective strain rate of the SSA, at (i,j)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) txz_c
txz_c(kc,j,i): Shear stress txz in the upper (kc) ice domain (at (i+1/2,j,kc))
real(dp), dimension(0:jmax, 0:imax) dzs_dxi
dzs_dxi(j,i): Derivative of zs by xi (at i+1/2,j)
real(dp), dimension(0:jmax, 0:imax) melt_star
melt_star(j,i): Superimposed ice formation rate at the ice surface
real(dp), dimension(0:jmax, 0:imax) h_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) zs_ref
zs_ref(j,i): Reference elevation for precip_present, temp_ma_present and temp_mj_present ...
real(dp), dimension(0:jmax, 0:imax) sq_g22_sgy
sq_g22_sgy(j,i): Square root of g22, at (i,j+1/2)
real(dp), dimension(0:jmax, 0:imax) snowfall
snowfall(j,i): Snowfall rate at the ice surface
real(dp) target_topo_tau
target_topo_tau: Relaxation time for target-topography adjustment
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.
real(dp), dimension(0:jmax, 0:imax) zs_target
zs_target(j,i): Target topography (ice surface)
real(dp) b
B: Semi-minor axis of the planet.
real(dp), dimension(0:jmax, 0:imax) temp_s
temp_s(j,i): Ice surface temperature
real(dp), dimension(0:jmax, 0:imax) melt
melt(j,i): Melting rate at the ice surface
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vz_t
vz_t(kt,j,i): Velocity in z-direction in the lower (kt) ice domain (at (i,j,kt+1/2)) ...
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
real(dp) g
G: Acceleration due to gravity.
real(dp), dimension(0:jmax, 0:imax) dzb_dxi
dzb_dxi(j,i): Derivative of zb by xi (at i+1/2,j)
real(dp) r
R: Radius of the planet.
integer(i1b), dimension(0:jmax, 0:imax) maske_target
maske_target(j,i): Target topography (ice-land-ocean mask)
real(dp), dimension(0:jmax, 0:imax) vy_m
vy_m(j,i): Mean (depth-averaged) velocity in y-direction, at (i,j+1/2)
subroutine sico_end()
Main routine of sico_end_m: Ending of SICOPOLIS.
Definition: sico_end_m.F90:51
integer(i4b) grip_time_max
grip_time_max: Maximum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
real(dp), dimension(0:jmax, 0:imax) zl_neu
(.)_neu: New value of quantity (.) computed during an integration step
integer(i1b), dimension(0:jmax, 0:imax) mask_ablation_type
mask_ablation_type(j,i): Mask indicating ablation type. 2: visible (ocean, for later developments)...
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
real(dp), dimension(0:jmax, 0:imax) vis_int_g
vis_int_g(j,i): Depth-integrated viscosity of the SSA, at (i,j)
Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:35
real(dp), dimension(0:jmax, 0:imax) dh_t_deta
dH_t_deta(j,i): Derivative of H_t by eta (at i,j+1/2)
integer(i4b) ndata_gi
ndata_gi: Number of data values for the glacial index
real(dp), dimension(0:jmax, 0:imax) ratio_sl_y
ratio_sl_y(j,i): Ratio of basal to surface velocity (slip ratio) in y-direction, at (i...
real(dp), dimension(0:jmax, 0:imax) temp_b
temp_b(j,i): Basal temperature
integer(i4b) ndata_grip
ndata_grip: Number of data values for the surface temperature anomaly
integer(i4b), dimension(0:jmax, 0:imax) q_weert
q_weert(j,i): Weertman exponent for the basal pressure
real(dp) beta1_0
BETA1_0: Degree-day factor for snow.
real(dp), dimension(0:jmax, 0:imax) d_help_b
d_help_b(j,i): Auxiliary quantity for the computation of vx_b and vy_b
real(dp), dimension(0:jmax, 0:imax) et
ET(j,i): Temperature excess at the ice surface (positive degree days divided by time) ...
real(dp) ea
ea: Abbreviation for exp(aa)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) lambda_shear_c
lambda_shear_c(kc,j,i): Shear fraction in the upper (kc) ice domain
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
real(dp), dimension(0:jmax, 0:imax) h_diff
h_diff(j,i): Diffusivity of the SIA evolution equation of the ice surface
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) dyz_t
dyz_t(kt,j,i): Strain rate dyz in the lower (kt) ice domain
Module for all openad-related subroutines.
Definition: openad_m.F90:37
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) d_help_t
d_help_t(kt,j,i): Auxiliary quantity for the computation of vx, vy und zs
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
real(dp) rho_c_r
RHO_C_R: Density times specific heat of the lithosphere.
real(dp), dimension(0:jmax, 0:imax) wss
wss(j,i): Isostatic steady-state displacement of the lithosphere
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) dxz_t
dxz_t(kt,j,i): Strain rate dxz in the lower (kt) ice domain
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
real(dp) pmax_0
PMAX_0: Saturation factor for the formation of superimposed ice.
real(dp), dimension(0:jmax, 0:imax) evap
evap(j,i): Evaporation rate at the ice surface
real(dp), dimension(0:jmax, 0:imax) flui_ave_sia
flui_ave_sia(j,i): Depth-averaged fluidity of the SIA
real(dp), dimension(0:jmax, 0:imax) accum
accum(j,i): Accumulation rate at the ice surface (includes liquid precipitation = rainfall!) ...
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:ktmax, 0:jmax, 0:imax) txz_t
txz_t(kt,j,i): Shear stress txz in the lower (kt) ice domain (at (i+1/2,j,kt))
real(dp), dimension(0:jmax, 0:imax) dh_c_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enh_t
enh_t(kt,j,i): Flow enhancement factor in the lower (kt) ice domain
real(dp), dimension(0:jmax, 0:imax) precip_ma_present
precip_ma_present(j,i): Present-day mean annual precipitation rate at the ice surface ...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
integer(i4b) gi_time_min
gi_time_min: Minimum time of the data values for the glacial index
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...
real(dp), dimension(0:jmax, 0:imax) zb_target
zb_target(j,i): Target topography (ice base)
real(dp) temp_precip_time_min
temp_precip_time_min: Minimum time of the surface-temperature and precipitation data ...
real(dp), dimension(0:jmax, 0:imax) smb_corr
smb_corr(j,i): Diagnosed SMB correction
real(dp) kei_r_max
kei_r_max: Maximum value of the argument r of the tabulated kei function
real(dp), dimension(0:jmax, 0:imax) q_b_apl
Q_b_apl(j,i): Applied basal melting rate.
real(dp) year_zero
year_zero: SICOPOLIS year zero in astronomical year numbering [ = signed year CE (AD) ] ...
real(dp), dimension(:), allocatable specmap_zsl
specmap_zsl(n): Data values for the sea level
real(dp), dimension(0:jmax, 0:imax) vy_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) dxx_c
dxx_c(kc,j,i): Strain rate dxx in the upper (kc) ice domain
integer(i4b) ncid_core
ncid_core: ID of the NetCDF time-series output file for the deep ice cores
subroutine oad_independent_init()
subroutine, public cost_independent_init()
The independent variable is perturbed by the xxVar at all places during the adjoint mode...
Definition: ctrl_m.F90:102
real(dp) beta2_0
BETA2_0: Degree-day factor for ice.
integer(i1b), dimension(0:jmax, 0:imax) n_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
integer(i4b) gi_time_stp
gi_time_stp: Time step of the data values for the glacial index
integer(i4b) ndata_temp_precip
ndata_temp_precip: Number of surface-temperature and precipitation data
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:krmax) zeta_r
zeta_r(kr): Sigma coordinate zeta_r of grid point kr
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(:), allocatable y_core
y_core(n): Coordinate eta (= y) of the prescribed borehole positions
real(dp), dimension(:), allocatable temp_precip_time
temp_precip_time(n): Times of the surface-temperature and precipitation data
real(dp), dimension(0:jmax, 0:imax) sq_g11_sgy
sq_g11_sgy(j,i): Square root of g11, at (i,j+1/2)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_m
temp_c_m(kc,j,i): Melting temperature in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) qx
qx(j,i): Volume flux in x-direction (depth-integrated vx, at (i+1/2,j))
real(dp), dimension(0:jmax, 0:imax) accum_present
accum_present(j,i): Present-day accumulation rate at the ice surface (for EISMINT, ISMIP HEINO and the north and south polar caps of Mars)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) dxx_t
dxx_t(kt,j,i): Strain rate dxx in the lower (kt) ice domain
real(dp), dimension(0:jmax, 0:imax) dzb_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
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) temp_ma_present
temp_ma_present(j,i): Present-day mean annual surface temperature
subroutine, public cost_dependent_init()
The dependent variable for the cost routine. This must be a scalar variable, although option for a su...
Definition: ctrl_m.F90:173
real(dp) nue
NUE: Water diffusivity in ice.
real(dp), dimension(0:jmax, 0:imax) precip_ma_lgm_anom
precip_ma_lgm_anom(j,i): LGM anomaly (ratio LGM/present) of the mean annual precipitation rate at the...
real(dp), dimension(0:jmax, 0:imax) zs_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) zl0
zl0(j,i): zl for isostatically relaxed ice-free conditions
real(dp), dimension(0:jmax, 0:imax) p_weert_inv
p_weert_inv(j,i): Inverse of p_weert
integer(i4b) gi_time_max
gi_time_max: Maximum time of the data values for the glacial index
real(dp) beta1_lt_0
BETA1_LT_0: Degree-day factor for snow at low summer temperatures.
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream
flag_shelfy_stream(j,i): Shelfy stream flag on the main grid. .true.: grounded ice, and at least one neighbour on the staggered grid is a shelfy stream point .false.: otherwise
real(dp), dimension(0:jmax, 0:imax) dzs_deta
dzs_deta(j,i): Derivative of zs by eta (at i,j+1/2)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) tyz_c
tyz_c(kc,j,i): Shear stress tyz in the upper (kc) ice domain (at (i,j+1/2,kc))
real(dp), dimension(0:jmax, 0:imax) q_gl_g
q_gl_g(j,i): Volume flux across the grounding line, at (i,j)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vz_c
vz_c(kc,j,i): Velocity in z-direction in the upper (kc) ice domain (at (i,j,kc+1/2)) ...
real(dp) phi0
PHI0: Standard parallel of the stereographic projection.
real(dp), dimension(0:jmax, 0:imax) vy_s_g
vy_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
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) sq_g11_sgx
sq_g11_sgx(j,i): Square root of g11, at (i+1/2,j)
real(dp), dimension(0:jmax, 0:imax, 12) gamma_precip_lgm_anom
gamma_precip_lgm_anom(j,i,n): negative natural logarithm of precip_lgm_anom(j,i,n) ...
integer(i1b), dimension(0:jmax, 0:imax) maske_maxextent
maske_maxextent(j,i): Maximum ice extent mask. 0: not allowed to glaciate, 1: allowed to glaciate...
real(dp) time_target_topo_final
time_target_topo_final: Final time for target-topography adjustment
real(dp), dimension(0:jmax, 0:imax) h_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
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 (...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) dyy_c
dyy_c(kc,j,i): Strain rate dyy in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) temp_mj_lgm_anom
temp_mj_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean summer (northern hemisphere...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) dxz_c
dxz_c(kc,j,i): Strain rate dxz in the upper (kc) 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 ...
real(dp), dimension(0:jmax, 0:imax) q_b_tot
Q_b_tot(j,i): Sum of Q_bm and Q_tld.
real(dp), dimension(0:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) dxy_c
dxy_c(kc,j,i): Strain rate dxy in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) vx_g
vx_g(j,i): Velocity in x-direction of the SSA, at (i,j)
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
real(dp), dimension(-10000:10000) kei
kei(n): Tabulated values of the kei function (Kelvin function of zero order)
real(dp), dimension(0:jmax, 0:imax) dzb_deta
dzb_deta(j,i): Derivative of zb by eta (at i,j+1/2)
real(dp), dimension(0:jmax, 0:imax) dh_t_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp) time_target_topo_init
time_target_topo_init: Initial time for target-topography adjustment
real(dp), dimension(0:jmax, 0:imax) dh_t_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dzm_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
integer(i4b) specmap_time_stp
specmap_time_stp: Time step of the data values for the sea level
real(dp), dimension(0:jmax, 0:imax) dzm_dxi
dzm_dxi(j,i): Derivative of zm by xi (at i+1/2,j)
real(dp), dimension(0:jmax, 0:imax) temp_ma_lgm_anom
temp_ma_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean annual surface temperature ...
integer(i1b) forcing_flag
forcing_flag: Flag for the forcing type. 1: forcing by a spatially constant surface temperature anoma...
real(dp), dimension(0:jmax, 0:imax) temph_b
temph_b(j,i): Basal temperature relative to the pressure melting point
real(dp), dimension(0:jmax, 0:imax) dzs_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vy_t
vy_t(kt,j,i): Velocity in y-direction in the lower (kt) ice domain (at (i,j+1/2,kt)) ...
real(dp), dimension(0:jmax, 0:imax) am_perp
am_perp(j,i): Ice volume flux across the z=zm interface
real(dp) beta
BETA: Clausius-Clapeyron gradient of ice.
real(dp), dimension(-jmax:jmax,-imax:imax) dist_dxdy
dist_dxdy(jr,ir): Distance between grid points with delta_i=ir, delta_j=jr
subroutine, public ctrl_init()
Initialiation of control variable. Recognized by OpenAD with the prefix xxVar, where Var is the varia...
Definition: ctrl_m.F90:33
real(dp), dimension(0:jmax, 0:imax) dzm_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
integer(i4b) ncid_temp_precip
ncid_temp_precip: ID of the NetCDF file containing the surface-temperature and precipitation data as ...
real(dp), dimension(0:jmax, 0:imax) sq_g11_g
sq_g11_g(j,i): Square root of the coefficient g11 of the metric tensor on grid point (i...
real(dp), dimension(0:jmax, 0:imax) as_perp_apl
as_perp_apl(j,i): Applied accumulation-ablation function (SMB)
logical, dimension(0:jmax, 0:imax) flag_calving_front_1
flag_calving_front_1(j,i): Calving front flag. .true.: calving front point (floating ice point with a...
real(dp), dimension(0:jmax, 0:imax) ratio_sl_x
ratio_sl_x(j,i): Ratio of basal to surface velocity (slip ratio) in x-direction, at (i+1/2...
integer(i4b) grip_time_min
grip_time_min: Minimum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) time_lag_asth
time_lag_asth(j,i): Time lag of the relaxing asthenosphere
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp), dimension(:), allocatable x_core
x_core(n): Coordinate xi (= x) of the prescribed borehole positions
Ending of SICOPOLIS.
Definition: sico_end_m.F90:35
real(dp), dimension(0:jmax, 0:imax) insq_g11_g
insq_g11_g(j,i): Inverse square root of g11 on grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) c_drag
c_drag(j,i): Auxiliary quantity for the computation of the basal drag
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vy_c
vy_c(kc,j,i): Velocity in y-direction in the upper (kc) ice domain (at (i,j+1/2,kc)) ...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
real(dp) kei_r_incr
kei_r_incr: Increment of the argument r of the tabulated kei function
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
real(dp) rho_sw
RHO_SW: Density of sea water.
logical, dimension(0:jmax, 0:imax) flag_calving_front_2
flag_calving_front_2(j,i): Calving front flag. .true.: calving front point (ocean point with at least...
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream_y
flag_shelfy_stream_y(j,i): Shelfy stream flag in y-direction, at (i,j+1/2). .true.: shelfy stream point .false.: otherwise
real(dp), dimension(0:jmax, 0:imax) zb_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) zm_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) dh_t_dxi
dH_t_dxi(j,i): Derivative of H_t by xi (at i+1/2,j)
integer(i4b) n_data_kei
n_data_kei: Number of tabulated values of the kei function
Ice discharge parameterization for the Greenland ice sheet.
real(dp) rho
RHO: Density of ice.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) sigma_c
sigma_c(kc,j,i): Effective stress in the upper (kc) ice domain
subroutine, public cost_final(runname)
This is the final cost calculation. The cost function structure is defined here. Currently is a "obse...
Definition: ctrl_m.F90:196
real(dp), dimension(0:jmax, 0:imax) vx_b
vx_b(j,i): Velocity in x-direction at the ice base, at (i+1/2,j)
real(dp), dimension(0:jmax, 0:imax) smb_corr_prescribed
smb_corr_prescribed(j,i): Prescribed SMB correction
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) sigma_t
sigma_t(kt,j,i): Effective stress in the lower (kt) ice domain
real(dp) rho_a
RHO_A: Density of the asthenosphere.
integer(i4b), dimension(0:jmax, 0:imax) p_weert
p_weert(j,i): Weertman exponent for the basal shear stress
real(dp), dimension(0:jmax, 0:imax) h_target
H_target(j,i): Target topography (ice thickness)
real(dp), dimension(:), allocatable lambda_core
lambda_core(n): Geographical longitude of the prescribed borehole positions
integer(i4b) n_core
n_core: Number of positions to be considered in the time-series file for deep boreholes ...
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
logical, dimension(0:jmax, 0:imax) flag_grounding_line_2
flag_grounding_line_2(j,i): Grounding line flag. .true.: grounding line point (floating ice point wit...
real(dp), dimension(0:jmax, 0:imax, 12) precip_present
precip_present(j,i,n): Present-day mean monthly precipitation rate at the ice surface ...
real(dp), dimension(0:jmax, 0:imax) am_perp_st
am_perp_st(j,i): Steady-state part of am_perp (without contribution of dzm_dtau)
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
real(dp), dimension(0:kcmax) zeta_c
zeta_c(kc): Sigma coordinate zeta_c of grid point kc
real(dp) lambda0
LAMBDA0: Reference longitude (central meridian) of the stereographic projection.
real(dp), dimension(0:jmax, 0:imax, 12) precip_lgm_anom
precip_lgm_anom(j,i,n): LGM anomaly (ratio LGM/present) of the mean monthly precipitation rate at the...
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...
real(dp), dimension(0:jmax, 0:imax, 12) temp_mm_lgm_anom
temp_mm_lgm_anom(j,i,n): LGM anomaly (difference LGM - present) of the mean monthly surface temperatu...
real(dp), dimension(0:jmax, 0:imax) h_w
H_w(j,i): Thickness of the water column under the ice base.
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) dzs_dtau
dzs_dtau(j,i): Derivative of zs by tau (time)
real(dp), dimension(0:jmax, 0:imax) flex_rig_lith
flex_rig_lith(j,i): Flexural rigidity of the lithosphere
real(dp), dimension(:), allocatable griptemp
griptemp(n): Data values for the surface temperature anomaly
real(dp), dimension(0:jmax, 0:imax) dzl_dtau
dzl_dtau(j,i): Derivative of zl by tau (time)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) tyz_t
tyz_t(kt,j,i): Shear stress tyz in the lower (kt) ice domain (at (i,j+1/2,kt))
real(dp), dimension(0:jmax, 0:imax) sq_g22_g
sq_g22_g(j,i): Square root of the coefficient g22 of the metric tensor on grid point (i...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vx_t
vx_t(kt,j,i): Velocity in x-direction in the lower (kt) ice domain (at (i+1/2,j,kt)) ...
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
real(dp), dimension(0:jmax, 0:imax) dzb_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) qy
qy(j,i): Volume flux in y-direction (depth-integrated vy, at (i,j+1/2))