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