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