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