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