SICOPOLIS V5-dev  Revision 1420
discharge_workers_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : d i s c h a r g e _ w o r k e r s _ m
4 !
5 !> @file
6 !!
7 !! Ice discharge parameterization for the Greenland ice sheet.
8 !!
9 !! References:
10 !! @li Calov, R., A. Robinson, M. Perrette and A. Ganopolski. 2015.\n
11 !! Simulating the Greenland ice sheet under present-day and
12 !! palaeo constraints including a new discharge parameterization.\n
13 !! Cryosphere 9, 179-196.
14 !!
15 !! @section Copyright
16 !!
17 !! Copyright 2011-2019 Reinhard Calov, Andrey Ganopolski
18 !! (Potsdam Institute for Climate Impact Research)
19 !!
20 !! @section License
21 !!
22 !! This file is part of SICOPOLIS.
23 !!
24 !! SICOPOLIS is free software: you can redistribute it and/or modify
25 !! it under the terms of the GNU General Public License as published by
26 !! the Free Software Foundation, either version 3 of the License, or
27 !! (at your option) any later version.
28 !!
29 !! SICOPOLIS is distributed in the hope that it will be useful,
30 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
31 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 !! GNU General Public License for more details.
33 !!
34 !! You should have received a copy of the GNU General Public License
35 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
36 !<
37 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
38 
39 !-------------------------------------------------------------------------------
40 !> Ice discharge parameterization for the Greenland ice sheet.
41 !<------------------------------------------------------------------------------
43 
44  use sico_types_m
46  use sico_vars_m
47  use error_m
48 
49 #if !defined(ALLOW_OPENAD) /* Normal */
50  use compare_float_m
51 #endif /* Normal */
52 
53  implicit none
54 
55 #if !defined(ALLOW_OPENAD) /* Normal */
56  integer(i4b), private :: disc
57  integer(i4b), private :: n_discharge_call, iter_mar_coa
58  real(dp), private :: c_dis_0, s_dis, c_dis_fac
59  real(dp), private :: t_sub_pd, alpha_sub, alpha_o, m_h, m_d, r_mar_eff
60  real(dp), private :: t_sea_freeze
61 #else /* OpenAD */
62  integer(i4b), public :: disc_dw
63  integer(i4b), public :: n_discharge_call_dw, iter_mar_coa_dw
64  real(dp), public :: c_dis_0_dw, s_dis_dw, c_dis_fac_dw
65  real(dp), public :: t_sub_pd_dw, alpha_sub_dw, alpha_o_dw, m_h_dw, m_d_dw, r_mar_eff_dw
66  real(dp), public :: t_sea_freeze_dw
67 #endif /* Normal vs. OpenAD */
68 
69  real(dp), public :: dt_glann, dt_sub
70 
71  integer(i1b), dimension(0:JMAX,0:IMAX), public :: mask_mar
72 
73 #if !defined(ALLOW_OPENAD) /* Normal */
74  real(dp), dimension(0:JMAX,0:IMAX), private :: c_dis
75 #else /* OpenAD */
76  real(dp), dimension(0:JMAX,0:IMAX), public :: c_dis_dw
77 #endif /* Normal vs. OpenAD */
78 
79  real(dp), dimension(0:JMAX,0:IMAX), public :: cst_dist, cos_grad_tc
80  real(dp), dimension(0:JMAX,0:IMAX), public :: dis_perp
81 
82 #if !defined(ALLOW_OPENAD) /* Normal */
83  private
84 #endif /* Normal */
85 
86 #if !defined(ALLOW_OPENAD_DIFFERENTIATE) /* ALLOW_OPENAD_DIFFERENTIATE unset */
88 #else /* ALLOW_OPENAD_DIFFERENTIATE set */
90 #endif /* Unset vs. set ALLOW_OPENAD_DIFFERENTIATE */
91 
92 contains
93 
94 #if !defined(ALLOW_OPENAD_DIFFERENTIATE) /* ALLOW_OPENAD_DIFFERENTIATE unset */
95 
96 !-------------------------------------------------------------------------------
97 !> Ice discharge parameters (Greenland).
98 !! [Assign ice discharge parameters.]
99 !<------------------------------------------------------------------------------
100  subroutine disc_param(dtime)
102  ! Author: Reinhard Calov
103  ! Institution: Potsdam Institute for Climate Impact Research
104  ! Date: 10.6.16
105 
106  ! Purpose: Calculate discharge parameters.
107 
108  implicit none
109 
110  real(dp), intent(in) :: dtime
111 
112  real(dp) :: dtime_mar_coa
113 
114 #if !defined(ALLOW_OPENAD) /* Normal */
115 
116  disc = disc
117  c_dis_0 = c_dis_0
118  c_dis_fac = c_dis_fac
119 
120  m_h = m_h
121  m_d = m_d
122 
123  r_mar_eff = r_mar_eff *1000.0_dp ! km -> m
124 
125 #else /* OpenAD */
126 
127  disc_dw = disc
128  c_dis_0_dw = c_dis_0
129  c_dis_fac_dw = c_dis_fac
130 
131  m_h_dw = m_h
132  m_d_dw = m_d
133 
134  r_mar_eff_dw = r_mar_eff *1000.0_dp ! km -> m
135 
136 #endif /* Normal vs. OpenAD */
137 
138 #if !defined(ALLOW_OPENAD) /* Normal */
139 
140 #if (defined(S_DIS))
141  s_dis = s_dis
142 #else
143  s_dis = 1.0_dp ! default value
144 #endif
145 
146 #else /* OpenAD */
147 
148 #if (defined(S_DIS))
149  s_dis_dw = s_dis
150 #else
151  s_dis_dw = 1.0_dp ! default value
152 #endif
153 
154 #endif /* Normal vs. OpenAD */
155 
156 #if !defined(ALLOW_OPENAD) /* Normal */
157 
158 #if (defined(ALPHA_SUB))
159  alpha_sub = alpha_sub
160 #else
161  alpha_sub = 0.5_dp ! default value
162 #endif
163 
164 #else /* OpenAD */
165 
166 #if (defined(ALPHA_SUB))
167  alpha_sub_dw = alpha_sub
168 #else
169  alpha_sub_dw = 0.5_dp ! default value
170 #endif
171 
172 #endif /* Normal vs. OpenAD */
173 
174 #if !defined(ALLOW_OPENAD) /* Normal */
175 
176 #if (defined(ALPHA_O))
177  alpha_o = alpha_o
178 #else
179  alpha_o = 1.0_dp ! default value
180 #endif
181 
182  t_sea_freeze = 0.0_dp - delta_tm_sw ! freezing temperature of sea water
183 
184 #else /* OpenAD */
185 
186 #if (defined(ALPHA_O))
187  alpha_o_dw = alpha_o
188 #else
189  alpha_o_dw = 1.0_dp ! default value
190 #endif
191 
192  t_sea_freeze_dw = 0.0_dp - delta_tm_sw ! freezing temperature of sea water
193 
194 #endif /* Normal vs. OpenAD */
195 
196 #if !defined(ALLOW_OPENAD) /* Normal */
197  n_discharge_call = -1
198 #else /* OpenAD */
199  n_discharge_call_dw = -1
200 #endif /* Normal vs. OpenAD */
201 
202 #if (defined(DTIME_MAR_COA0))
203  dtime_mar_coa = dtime_mar_coa0*year_sec ! a -> s
204 #else
205  errormsg = ' >>> disc_param: DTIME_MAR_COA0 not defined in header file!'
206  call error(errormsg)
207 #endif
208 
209 #if !defined(ALLOW_OPENAD) /* Normal */
210  if (.not.approx_integer_multiple(dtime_mar_coa, dtime, eps_sp_dp)) then
211  errormsg = ' >>> disc_param: dtime_mar_coa must be a multiple of dtime!'
212  call error(errormsg)
213  end if
214 #else /* OpenAD */
215  print *, ' >>> disc_param: compare_float_m not used in adjoint'
216  print *, ' applications! double check '
217  print *, ' dtime_mar_coa and dtime are multiples'
218 #endif /* Normal vs. OpenAD */
219 
220 #if !defined(ALLOW_OPENAD) /* Normal */
221  iter_mar_coa = nint(dtime_mar_coa/dtime)
222 #else /* OpenAD */
223  iter_mar_coa_dw = nint(dtime_mar_coa/dtime)
224 #endif /* Normal vs. OpenAD */
225 
226 #if (GRID > 1)
227  errormsg = ' >>> disc_param: GRID==2 not allowed for this application!'
228  call error(errormsg)
229 #endif
230 
231 #if (ANF_DAT != 1 && defined(EXEC_MAKE_C_DIS_0))
232  errormsg = ' >>> disc_param: EXEC_MAKE_C_DIS_0 defined requires ANF_DAT==1!'
233  call error(errormsg)
234 #endif
235 
236  end subroutine disc_param
237 
238 #endif /* ALLOW_OPENAD_DIFFERENTIATE unset */
239 
240 !-------------------------------------------------------------------------------
241 !> Constant in ice discharge parameterization (Greenland).
242 !! [Determine (amount of magnitude of) constant in ice discharge
243 !! parameterization.]
244 !<------------------------------------------------------------------------------
245  subroutine calc_c_dis_0(dxi, deta)
247  ! Author: Reinhard Calov
248  ! Institution: Potsdam Institute for Climate Impact Research
249  ! Date: 8.6.16
250  ! Purpose: Compute c_dis_0 indirectly from present-day topography
251  ! from the condition that the parameterized total discharge
252  ! equals the discharge give by dis_target. Note: This c_dis_0 could
253  ! differ from the optimal one yielded from the siumualted ice thickness.
254  !
255  ! This is very useful, because c_dis_0 can vary a lot for different powers
256  ! m_D and m_H. This way we easier yield the approximately right value
257  ! for c_dis_0 for given powers m_D and m_H.
258 
259  ! This routine is not to be used regularly, and it is only executed if the
260  ! parameter EXEC_MAKE_C_DIS_0 is defined in the header file.
261 
262  implicit none
263 
264  real(dp), intent(in) :: dxi, deta
265 
266  integer(i4b) :: i, j
267  real(dp) :: disc_target, disc_tot, dT_sub
268 
269  ! targeted total ice discharge
270  ! 350 Gt/yr Calov et al. (2015) !!! 400 Gt/yr van den Broeke et al (2016)
271  disc_target = 350.0_dp ! in Gt/yr
272 
273  h_c=zs-zb; h_t=0.0_dp
274 
275 #if !defined(ALLOW_OPENAD) /* Normal */
276  c_dis_0 = 1.0_dp
277  c_dis_fac = 1.0_dp
278  s_dis = 1.0_dp
279 #else /* OpenAD */
280  c_dis_0_dw = 1.0_dp
281  c_dis_fac_dw = 1.0_dp
282  s_dis_dw = 1.0_dp
283 #endif /* Normal vs. OpenAD */
284 
285  call disc_fields()
286 
287  ! ensure that we get present-day discharge here (disc=1).
288  dt_glann = 0.0_dp
289 
290 #if !defined(ALLOW_OPENAD) /* Normal */
291  n_discharge_call = -1
292 #else /* OpenAD */
293  n_discharge_call_dw = -1
294 #endif /* Normal vs. OpenAD */
295 
296  call discharge(dxi, deta)
297 
298  ! disc_tot=0.0_dp
299  ! do i=0, IMAX
300  ! do j=0, JMAX
301  ! if(maske(j,i).eq.0_i1b.or.maske(j,i).eq.3_i1b) then
302  ! disc_tot=disc_tot+dis_perp(j,i)*area(j,i)
303  ! end if
304  ! end do
305  ! end do
306  !
307  ! disc_tot = disc_tot*YEAR_SEC ! m^3/s -> m^3/a
308 
309  disc_tot = 0.0_dp
310 
311  do i=1, imax-1
312  do j=1, jmax-1
313 
314  if (mask_mar(j,i) == 1_i1b) then
315  disc_tot = disc_tot + dis_perp(j,i)*area(j,i)
316  end if
317 
318  end do
319  end do
320 
321  disc_tot = disc_tot *year_sec *(rho/rho_w)
322  ! m^3/s ice equiv. -> m^3/a water equiv.
323 
324  disc_target = disc_target*1.0e+12_dp/rho_w ! Gt/a -> m^3/a water equiv.
325 
326  write(6, fmt='(a)') ' '
327  write(6, fmt='(a)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
328 
329  write(6, fmt='(a)') ' '
330 
331 #if !defined(ALLOW_OPENAD) /* Normal */
332  write(6, fmt='(3x,a,es12.4)') 'c_dis_0_init = ', c_dis_0
333 #else /* OpenAD */
334  write(6, fmt='(3x,a,es12.4)') 'c_dis_0_init = ', c_dis_0_dw
335 #endif /* Normal vs. OpenAD */
336 
337 #if !defined(ALLOW_OPENAD) /* Normal */
338  c_dis_0 = c_dis_0 * disc_target/disc_tot
339 #else /* OpenAD */
340  c_dis_0_dw = c_dis_0_dw * disc_target/disc_tot
341 #endif /* Normal vs. OpenAD */
342 
343  write(6, fmt='(a)') ' '
344  write(6, fmt='(3x,a,es12.4)') 'disc_target = ', disc_target
345  write(6, fmt='(3x,a,es12.4)') 'disc_tot = ', disc_tot
346 
347  write(6, fmt='(a)') ' '
348 
349 #if !defined(ALLOW_OPENAD) /* Normal */
350  write(6, fmt='(3x,a,es12.4)') '--> c_dis_0 = ', c_dis_0
351 #else /* OpenAD */
352  write(6, fmt='(3x,a,es12.4)') '--> c_dis_0 = ', c_dis_0_dw
353 #endif /* Normal vs. OpenAD */
354 
355  end subroutine calc_c_dis_0
356 
357 !-------------------------------------------------------------------------------
358 !> Dependence of ice discharge coefficient on latitude (Greenland).
359 !! [Determine dependence of ice discharge coefficient on latitude.
360 !! This can be improved. For now I recommend s_dis=1.]
361 !<------------------------------------------------------------------------------
362  subroutine disc_fields()
364  ! Author: Reinhard Calov
365  ! Institution: Potsdam Institute for Climate Impact Research
366  ! Date: 8.6.16
367 
368  ! Purpose: Assign fields relevant for ice discharge parameterization
369  ! 1. Discharge coefficient field using linear dependence on
370  ! latitude.
371  ! 2. Subsurface temperature dependence on longitude and latitude.
372 
373  implicit none
374 
375  integer(i4b) :: i, j
376  real(dp) :: delta_phi=25.0_dp ! approximately the phi
377  ! spanning the middle of domain
378 
379 #if !defined(ALLOW_OPENAD) /* Normal */
380 
381  c_dis = c_dis_0*c_dis_fac &
382  *(1.0_dp-(1.0_dp-s_dis)*(phi*pi_180_inv-60.0_dp)/delta_phi)
383 
384  c_dis = max(c_dis, 0.0_dp)
385 
386 #else /* OpenAD */
387 
388  c_dis_dw = c_dis_0_dw*c_dis_fac_dw &
389  *(1.0_dp-(1.0_dp-s_dis_dw)*(phi*pi_180_inv-60.0_dp)/delta_phi)
390 
391  c_dis_dw = max(c_dis_dw, 0.0_dp)
392 
393 #endif /* Normal vs. OpenAD */
394 
395 #if !defined(ALLOW_OPENAD) /* Normal */
396  t_sub_pd = 3.0_dp ! Can depend on lambda, phi later on
397 #else /* OpenAD */
398  t_sub_pd_dw = 3.0_dp ! Can depend on lambda, phi later on
399 #endif /* Normal vs. OpenAD */
400 
401  end subroutine disc_fields
402 
403 !-------------------------------------------------------------------------------
404 !> Ice discharge parameterization main formula, controler (general).
405 !! [Compute ice discharge via a parameterization using distance of ice
406 !! margin to coast and ice thickness as parameters.]
407 !<------------------------------------------------------------------------------
408  subroutine discharge(dxi, deta)
410 #ifdef ALLOW_OPENAD /* OpenAD */
411  use ctrl_m, only: myfloor
412 #endif /* OpenAD */
413 
414  ! Authors: Reinhard Calov, Andrey Ganopolski
415  ! Institution: Potsdam Institute for Climate Impact Research
416  ! Date: 11.6.16
417 
418  implicit none
419 
420  real(dp), intent(in) :: dxi, deta
421 
422  real(dp), parameter :: alpha_tc=60.0_dp ! maximal allowed angle
423  real(dp) :: cos_tc
424 
425 #ifdef ALLOW_OPENAD /* OpenAD */
426  integer(i4b) :: i, j, valmodint
427  real(dp) :: tmp, valmod
428 #endif /* OpenAD */
429 
430 #if !defined(ALLOW_OPENAD) /* Normal */
431  n_discharge_call = n_discharge_call + 1
432 #else /* OpenAD */
433  n_discharge_call_dw = n_discharge_call_dw + 1
434 #endif /* Normal vs. OpenAD */
435 
436  cos_tc=dcos(alpha_tc*pi_180)
437 
438 #if !defined(ALLOW_OPENAD) /* Normal */
439 
440  if ((mod(n_discharge_call, iter_mar_coa)==0).or.(n_discharge_call==1)) then
441  write(6, fmt='(10x,a)') 'Computation of mask_mar, cst_dist, cos_grad_tc'
442  call marginal_ring(dxi, deta)
443  call coastal_distance(dxi, deta)
444  end if
445 
446 #else /* OpenAD */
447 
448  ! this is mod broken down algorithmically:
449  tmp = real(n_discharge_call_dw)/real(iter_mar_coa_dw)
450  call myfloor(tmp, valmodint)
451  valmod = real(n_discharge_call_DW) - real(iter_mar_coa_DW) * real(valmodint)
452  if ((valmod==0) .or. (n_discharge_call_dw==1)) then
453  write(6, fmt='(10x,a)') 'Computation of mask_mar, cst_dist, cos_grad_tc'
454  call marginal_ring(dxi, deta)
455  call coastal_distance(dxi, deta)
456  end if
457 
458 #endif /* Normal vs. OpenAD */
459 
460 #if !defined(ALLOW_OPENAD) /* Normal */
461 
462  if(disc.ge.1) then !------------------ disc >= 1
463  where(mask_mar.eq.1_i1b.and.cos_grad_tc.ge.cos_tc)
464  dis_perp=c_dis*(h_c+h_t)**m_h/cst_dist**m_d
465  elsewhere
466  dis_perp=0.0_dp
467  end where
468 
469  if(disc.eq.2) then !----------------- disc = 2
470  call calc_sub_oc_dt(t_sub_pd, dt_glann, dt_sub)
471  dis_perp=dis_perp*(1.0_dp+alpha_sub*dt_sub) ! actual discharge present-day one corrected
472  ! via sub-ocean temperature anomaly
473  end if
474  dis_perp=max(0.0_dp, dis_perp) ! ensure positive values
475  else if(disc.eq.0) then !-------------- disc = 0
476  dis_perp=0.0_dp
477  end if
478 
479 #else /* OpenAD */
480 
481  if(disc_dw.ge.1) then !------------------ disc >= 1
482  do i=0, imax
483  do j=0, jmax
484  if (mask_mar(j,i).eq.1_i1b .and. cos_grad_tc(j,i).ge.cos_tc) then
485  dis_perp(j,i) = c_dis_dw(j,i)*( h_c(j,i) + h_t(j,i) )**m_h_dw/cst_dist(j,i)**m_d_dw
486  else
487  dis_perp(j,i) = 0.0_dp
488  end if
489  end do
490  end do
491  if(disc_dw.eq.2) then !----------------- disc = 2
492  call calc_sub_oc_dt(t_sub_pd_dw, dt_glann, dt_sub)
493  dis_perp=dis_perp*(1.0_dp+alpha_sub_dw*dt_sub) ! actual discharge present-day one corrected
494  ! via sub-ocean temperature anomaly
495  end if
496  dis_perp=max(0.0_dp, dis_perp) ! ensure positive values
497  else if(disc_dw.eq.0) then !-------------- disc = 0
498  dis_perp=0.0_dp
499  end if
500 
501 #endif /* Normal vs. OpenAD */
502 
503  end subroutine discharge
504 
505 !-------------------------------------------------------------------------------
506 !> Distance to the coast (general).
507 !! [Compute distance to the coast for every land point.]
508 !<------------------------------------------------------------------------------
509  subroutine coastal_distance(dxi, deta)
510 
511  ! Author: Reinhard Calov
512  ! Institution: Potsdam Institute for Climate Impact Research
513  ! Date: 2.11.11
514  ! Methode: Options are:
515  ! i_search=1: Brute force.
516  ! i_search=2: Defining a square with two grid distances side length around the
517  ! actual grid point and enlanging the square successively by
518  ! one grid. Because this should have been a circle, the
519  ! square is enlarged again if the coast line was found; this time
520  ! enlargement is by a rectangles a the respective four side of
521  ! the square. Smaller side lenght equals the rectangle root of 2 of the
522  ! larger side.
523 
524 #ifdef ALLOW_OPENAD /* OpenAD */
525  use ctrl_m, only: myceiling
526 #endif /* OpenAD */
527 
528  implicit none
529 
530  real(dp), intent(in) :: dxi, deta
531 
532  ! --------------------------------------------------
533  real(dp), parameter :: grid_dist=5000.0_dp
534  real(dp), parameter :: c_smooth=0.2_dp
535  integer(i4b), parameter :: i_search=2
536  ! --------------------------------------------------
537  integer(i4b) :: i, j, i_pos, j_pos, l, l_e, d_l
538  real(dp) :: dxi_inv, deta_inv
539  real(dp) :: cst_dist_tmp, qrt_1, qrt_2
540  logical :: leave_loop
541 
542  real(dp), allocatable, dimension(:,:) :: zs_tmp, grad_zs_x, grad_zs_y, &
543  grad_dist_x, grad_dist_y
544 
545  allocate(zs_tmp(0:jmax,0:imax), grad_zs_x(0:jmax,0:imax), grad_zs_y(0:jmax,0:imax), &
546  grad_dist_x(0:jmax,0:imax), grad_dist_y(0:jmax,0:imax))
547 
548  select case (i_search)
549 
550  case(1)
551  ! brute force computation of minimal distance to coast for every land point
552 
553  do i_pos=0, imax
554  do j_pos=0, jmax
555  if(maske(j_pos,i_pos).ne.2_i1b) then
556  cst_dist(j_pos,i_pos)=1.d+20
557  do i=0, imax
558  do j=0, jmax
559  if(maske(j,i).eq.2_i1b) then
560  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
561  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
562  cst_dist(j_pos,i_pos)=cst_dist_tmp
563  end if
564  end if
565  end do
566  end do
567  else
568  cst_dist(j_pos,i_pos)=grid_dist
569  end if
570  end do
571  end do
572 
573  case(2)
574 
575  do i_pos=0, imax
576  do j_pos=0, jmax
577  if(maske(j_pos,i_pos).ne.2_i1b) then
578  leave_loop=.false.
579  cst_dist(j_pos,i_pos)=1.d+20
580 
581 #if !defined(ALLOW_OPENAD) /* Normal */
582  do l=1, max(imax,jmax)
583 #else /* OpenAD */
584  l = 1
585  do while (l<= max(imax,jmax) .and. leave_loop.eqv..false.)
586 #endif /* Normal vs. OpenAD */
587 
588  do i=max(i_pos-l,0), min(i_pos+l,imax)
589  j=max(j_pos-l,0); j=min(j,jmax)
590  if(maske(j,i).eq.2_i1b) then
591  leave_loop=.true.
592  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
593  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
594  cst_dist(j_pos,i_pos)=cst_dist_tmp
595  end if
596  end if
597  j=min(j_pos+l,jmax)
598  if(maske(j,i).eq.2_i1b) then
599  leave_loop=.true.
600  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
601  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
602  cst_dist(j_pos,i_pos)=cst_dist_tmp
603  end if
604  end if
605  end do
606  do j=max(j_pos-l+1,0), min(j_pos+l-1,jmax)
607  i=max(i_pos-l,0); i=min(i,imax)
608  if(maske(j,i).eq.2_i1b) then
609  leave_loop=.true.
610  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
611  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
612  cst_dist(j_pos,i_pos)=cst_dist_tmp
613  end if
614  end if
615  i=min(i_pos+l,imax)
616  if(maske(j,i).eq.2_i1b) then
617  leave_loop=.true.
618  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
619  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
620  cst_dist(j_pos,i_pos)=cst_dist_tmp
621  end if
622  end if
623  end do
624 
625 #if !defined(ALLOW_OPENAD) /* Normal */
626  if(leave_loop) then
627  l_e=l
628  d_l=ceiling(l_e*sqrt(2.0_dp)+0.001_dp)
629  exit ! leave loop in l
630  end if
631 #else /* OpenAD */
632  if(leave_loop) then
633  l_e=l
634  call myceiling((l_e*sqrt(2.0_dp)+0.001_dp), d_l)
635  end if
636  l = l+1
637 #endif /* Normal vs. OpenAD */
638 
639  end do
640  ! left
641  do i=max(i_pos-(l_e+d_l),0),min(i_pos-(l_e-1),imax)
642  do j=max(j_pos-l_e,0),min(j_pos+l_e,jmax)
643  if(maske(j,i).eq.2_i1b) then
644  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
645  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
646  cst_dist(j_pos,i_pos)=cst_dist_tmp
647  end if
648  end if
649  end do
650  end do
651  ! right
652  do i=max(i_pos+l_e+1,0),min(i_pos+l_e+d_l,imax)
653  do j=max(j_pos-l_e,0),min(j_pos+l_e,jmax)
654  if(maske(j,i).eq.2_i1b) then
655  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
656  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
657  cst_dist(j_pos,i_pos)=cst_dist_tmp
658  end if
659  end if
660  end do
661  end do
662  ! lower
663  do i=max(i_pos-l_e,0),min(i_pos+l_e,imax)
664  do j=max(j_pos-(l_e+d_l),0),min(j_pos-(l_e-1),jmax)
665  if(maske(j,i).eq.2_i1b) then
666  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
667  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
668  cst_dist(j_pos,i_pos)=cst_dist_tmp
669  end if
670  end if
671  end do
672  end do
673  ! upper
674  do i=max(i_pos-l_e,0),min(i_pos+l_e,imax)
675  do j=max(j_pos+l_e+1,0),min(j_pos+l_e+d_l,jmax)
676  if(maske(j,i).eq.2_i1b) then
677  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
678  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
679  cst_dist(j_pos,i_pos)=cst_dist_tmp
680  end if
681  end if
682  end do
683  end do
684  else
685  cst_dist(j_pos,i_pos)=grid_dist
686  end if
687  end do
688  end do
689 
690  end select
691 
692  ! test GIS only
693 
694  if(.false.) then
695  do i=0, imax
696  do j=0, jmax
697  if(xi(i).le.-350000.0_dp &
698  .and.-2800000.0_dp.le.eta(j) &
699  .and.eta(j).le.-2250000.0_dp &
700  .or. &
701  xi(i).ge.100000.0_dp &
702  .and.-2280000.0_dp.le.eta(j) &
703  .and.eta(j).le.-1880000.0_dp &
704  .or. &
705  eta(j).ge.-1250000.0_dp &
706  .and.0.0_dp.le.xi(i) &
707  .and.xi(i).le.230000.0_dp) then
708  cst_dist(j,i)=0.3_dp*cst_dist(j,i)
709  end if
710  end do
711  end do
712  end if
713 
714  ! forbid zero, something like grid distance is reasonable here.
715  ! Dependents of coordinate system, for spherical coordinate system:
716  ! take length on Earth's surface. For now, simply just 5 km are taken.
717 
718  cst_dist=max(grid_dist, cst_dist)
719 
720  ! forbit to large values
721  cst_dist=min(1.0e10_dp, cst_dist)
722 
723  zs_tmp=zs
724 
725  ! smoothing the surface topography
726  do i=1, imax-1
727  do j=1, jmax-1
728  zs_tmp(j,i) = c_smooth*((zs(j,i-1)+zs(j,i+1))+(zs(j-1,i)+zs(j+1,i))) &
729  +(1.0_dp-4.0_dp*c_smooth)*zs(j,i)
730  end do
731  end do
732 
733  ! gradients (2. order for now ...)
734 
735  dxi_inv = 1.0_dp/dxi
736  deta_inv = 1.0_dp/deta
737 
738 ! ------ x-derivatives
739 
740  do i=1, imax-1
741  do j=0, jmax
742  grad_zs_x(j,i) = (zs_tmp(j,i+1)-zs_tmp(j,i-1))*0.5_dp*dxi_inv*insq_g11_g(j,i)
743  grad_dist_x(j,i) = (cst_dist(j,i+1)-cst_dist(j,i-1))*0.5_dp*dxi_inv*insq_g11_g(j,i)
744  end do
745  end do
746 
747  do j=0, jmax
748  grad_zs_x(j,0) = (zs_tmp(j,1)-zs_tmp(j,0))*dxi_inv*insq_g11_g(j,0)
749  grad_zs_x(j,imax) = (zs_tmp(j,imax)-zs_tmp(j,imax-1))*dxi_inv*insq_g11_g(j,imax)
750  grad_dist_x(j,0) = (cst_dist(j,1)-cst_dist(j,0))*dxi_inv*insq_g11_g(j,0)
751  grad_dist_x(j,imax) = (cst_dist(j,imax)-cst_dist(j,imax-1))*dxi_inv*insq_g11_g(j,imax)
752  end do
753 
754 ! ------ y-derivatives
755 
756  do i=0, imax
757  do j=1, jmax-1
758  grad_zs_y(j,i) = (zs_tmp(j+1,i)-zs_tmp(j-1,i))*0.5_dp*deta_inv*insq_g22_g(j,i)
759  grad_dist_y(j,i) = (cst_dist(j+1,i)-cst_dist(j-1,i))*0.5_dp*deta_inv*insq_g22_g(j,i)
760  end do
761  end do
762 
763  do i=0, imax
764  grad_zs_y(0,i) = (zs_tmp(1,i)-zs_tmp(0,i))*deta_inv*insq_g22_g(0,i)
765  grad_zs_y(jmax,i) = (zs_tmp(jmax,i)-zs_tmp(jmax-1,i))*deta_inv*insq_g22_g(jmax,i)
766  grad_dist_y(0,i) = (cst_dist(1,i)-cst_dist(0,i))*deta_inv*insq_g22_g(0,i)
767  grad_dist_y(jmax,i) = (cst_dist(jmax,i)-cst_dist(jmax-1,i))*deta_inv*insq_g22_g(jmax,i)
768  end do
769 
770  ! angle between topographical gradient and coastal distance gradient
771 
772  do i=0, imax
773  do j=0, jmax
774  qrt_1=grad_zs_x(j,i)*grad_zs_x(j,i)+grad_zs_y(j,i)*grad_zs_y(j,i)
775  qrt_2=grad_dist_x(j,i)*grad_dist_x(j,i)+grad_dist_y(j,i)*grad_dist_y(j,i)
776  if(qrt_1.ne.0.0_dp.and.qrt_2.ne.0.0_dp) then
777  !!! top_cst_alpha(j,i)=dacos( (grad_zs_x(j,i)*grad_dist_x(j,i)+grad_zs_y(j,i)*grad_dist_y(j,i))/ &
778  !!! (sqrt(qrt_1)*sqrt(qrt_2)) )
779 
780  cos_grad_tc(j,i)= (grad_zs_x(j,i)*grad_dist_x(j,i)+grad_zs_y(j,i)*grad_dist_y(j,i))/ &
781  (sqrt(qrt_1)*sqrt(qrt_2))
782  else
783  cos_grad_tc(j,i)=-1.0_dp
784  end if
785  end do
786  end do
787 
788  deallocate(zs_tmp, grad_zs_x, grad_zs_y, grad_dist_x, grad_dist_y)
789 
790  end subroutine coastal_distance
791 
792 !-------------------------------------------------------------------------------
793 !> Ring along an ice sheet margin (general).
794 !! [Compute marginal ring.]
795 !<------------------------------------------------------------------------------
796  subroutine marginal_ring(dxi, deta)
797 
798  ! Author: Reinhard Calov
799  ! Institution: Potsdam Institute for Climate Impact Research
800  ! Date: 28.10.11
801 
802  ! Purpose: Determine an r_max_eff wide ring arong ice margins
803  ! towards the interior of the ice sheet. Small stripes of land are
804  ! not considered as land. For the logics De Morgan is used often.
805 
806  ! Methode: Two staggered loops in i,j. The inner i,j loop acts inside
807  ! a rectangle defined by r_mar_eff. r_mar_eff should be small compared
808  ! to the domain size.
809 
810  implicit none
811 
812  real(dp), intent(in) :: dxi, deta
813 
814  integer(i4b) :: i, j, i_pos, j_pos, di_eff, dj_eff
815  integer(i4b) :: count_tmp
816  real(dp) :: r_p
817 
818 #if !defined(ALLOW_OPENAD) /* Normal */
819  if(r_mar_eff.le.1.0e6_dp) then
820 #else /* OpenAD */
821  if(r_mar_eff_dw.le.1.0e6_dp) then
822 #endif /* Normal vs. OpenAD */
823 
824  mask_mar=0_i1b
825  do i_pos=1, imax-1
826  do j_pos=1, jmax-1
827  if((maske(j_pos,i_pos).eq.1_i1b.or.maske(j_pos,i_pos).eq.2_i1b).and. &
828  .not.(maske(j_pos,i_pos).eq.1_i1b.and.maske(j_pos,i_pos-1).eq.0_i1b &
829  .and.maske(j_pos,i_pos+1).eq.0_i1b.or. &
830  maske(j_pos,i_pos).eq.1_i1b.and.maske(j_pos-1,i_pos).eq.0_i1b &
831  .and.maske(j_pos+1,i_pos).eq.0_i1b.or. &
832  maske(j_pos,i_pos).eq.1_i1b.and.maske(j_pos,i_pos-1).eq.3_i1b &
833  .and.maske(j_pos,i_pos+1).eq.3_i1b.or. &
834  maske(j_pos,i_pos).eq.1_i1b.and.maske(j_pos-1,i_pos).eq.3_i1b &
835  .and.maske(j_pos+1,i_pos).eq.3_i1b.or. &
836  maske(j_pos,i_pos).eq.1_i1b.and.maske(j_pos,i_pos-1).eq.3_i1b &
837  .and.maske(j_pos,i_pos+1).eq.0_i1b.or. &
838  maske(j_pos,i_pos).eq.1_i1b.and.maske(j_pos-1,i_pos).eq.0_i1b &
839  .and.maske(j_pos+1,i_pos).eq.3_i1b)) then ! outside ice sheet, exclude isolated land stripes
840 
841 #if !defined(ALLOW_OPENAD) /* Normal */
842  di_eff=int(r_mar_eff/dxi)+1; dj_eff=int(r_mar_eff/deta)+1 ! only for grid=0, 1 yet!
843 #else /* OpenAD */
844  di_eff=int(r_mar_eff_dw/dxi)+1; dj_eff=int(r_mar_eff_dw/deta)+1 ! only for grid=0, 1 yet!
845 #endif /* Normal vs. OpenAD */
846 
847  do i=max(i_pos-di_eff,0), min(i_pos+di_eff,imax)
848  do j=max(j_pos-dj_eff,0), min(j_pos+dj_eff,jmax)
849  r_p=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
850 
851 #if !defined(ALLOW_OPENAD) /* Normal */
852  if(r_p.le.r_mar_eff.and.(maske(j,i).eq.0_i1b.or.maske(j,i).eq.3_i1b)) then
853 #else /* OpenAD */
854  if(r_p.le.r_mar_eff_dw.and.(maske(j,i).eq.0_i1b.or.maske(j,i).eq.3_i1b)) then
855 #endif /* Normal vs. OpenAD */
856 
857  mask_mar(j,i)=1_i1b
858  end if
859  end do
860  end do
861  end if
862  end do
863  end do
864 
865 #if !defined(ALLOW_OPENAD) /* Normal */
866  where(maske.ge.1_i1b)
867  mask_mar=1_i1b
868  end where
869 #else /* OpenAD */
870  do i=0,imax
871  do j=0,jmax
872  if (maske(j,i).ge.1_i1b) then
873  mask_mar(j,i)=1_i1b
874  end if
875  end do
876  end do
877 #endif /* Normal vs. OpenAD */
878 
879  else ! the ring encompassed entire Greenland
880  mask_mar=1_i1b
881  end if
882 
883  end subroutine marginal_ring
884 
885 !-------------------------------------------------------------------------------
886 !> Anomaly of subsurface temperature (general).
887 !! [Compute anomaly of subsurface temperature with a simple parameterization.]
888 !<------------------------------------------------------------------------------
889  subroutine calc_sub_oc_dt(T_sub_PD, dT_glann, dT_sub)
890 
891  ! Author: Reinhard Calov
892  ! Institution: Potsdam Institute for Climate Impact Research
893  ! Purpose: Compute anomaly of subsurface temperature
894  ! with a simple parameterization
895  ! using global annual temperature
896  ! anomaly (e.g. from CLIMBER)
897 
898  implicit none
899 
900  real(dp), intent(in) :: T_sub_PD, dT_glann
901  real(dp), intent(out) :: dT_sub
902 
903  real(dp) :: T_sub
904 
905 #if !defined(ALLOW_OPENAD) /* Normal */
906  dt_sub=alpha_o*dt_glann ! Sub-ocean temperature anomaly
907  t_sub=max(t_sea_freeze, t_sub_pd)+dt_sub ! Actual sub-ocean temperature
908  t_sub=max(t_sea_freeze, t_sub) ! Ocean temperature should stay
909  ! above freezing point
910 #else /* OpenAD */
911  dt_sub=alpha_o_dw*dt_glann ! Sub-ocean temperature anomaly
912  t_sub=max(t_sea_freeze_dw, t_sub_pd)+dt_sub ! Actual sub-ocean temperature
913  t_sub=max(t_sea_freeze_dw, t_sub) ! Ocean temperature should stay
914  ! above freezing point
915 #endif /* Normal vs. OpenAD */
916 
917  dt_sub=t_sub-t_sub_pd
918 
919  end subroutine calc_sub_oc_dt
920 
921 !-------------------------------------------------------------------------------
922 
923 end module discharge_workers_m
924 !
real(dp) delta_tm_sw
DELTA_TM_SW: Melting point depression of sea water due to its average salinity.
integer(i1b), dimension(0:jmax, 0:imax), public mask_mar
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
real(dp), dimension(0:jmax, 0:imax), public cos_grad_tc
Definition: ctrl_m.F90:9
real(dp) rho_w
RHO_W: Density of pure water.
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
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) area
area(j,i): Area of grid cell associated with grid point (i,j)
subroutine, public calc_c_dis_0(dxi, deta)
Constant in ice discharge parameterization (Greenland). [Determine (amount of magnitude of) constant ...
real(dp), dimension(0:jmax, 0:imax), public cst_dist
real(dp), parameter pi_180_inv
pi_180_inv: 180 divided by pi (-> rad to deg)
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
subroutine, public disc_param(dtime)
Ice discharge parameters (Greenland). [Assign ice discharge parameters.].
real(dp), parameter pi_180
pi_180: pi divided by 180 (-> deg to rad)
real(dp), dimension(0:jmax, 0:imax), public dis_perp
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), 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...
subroutine, public disc_fields()
Dependence of ice discharge coefficient on latitude (Greenland). [Determine dependence of ice dischar...
Writing of error messages and stopping execution.
Definition: error_m.F90:35
Comparison of single- or double-precision floating-point numbers.
real(dp), public dt_glann
integer(i1b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
subroutine, public discharge(dxi, deta)
Ice discharge parameterization main formula, controler (general). [Compute ice discharge via a parame...
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
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:imax) xi
xi(i): Coordinate xi (= x) of grid point i
real(dp), parameter eps_sp_dp
eps_sp_dp: Small number to single-precision accuracy in double precision
character(len=256) errormsg
errormsg: Error-message string
Ice discharge parameterization for the Greenland ice sheet.
real(dp) rho
RHO: Density of ice.
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
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...
Declarations of global variables for SICOPOLIS.