SICOPOLIS V5-dev  Revision 1173
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-2017 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 
49  use compare_float_m
50 
51  implicit none
52 
53  integer(i4b), private :: disc
54  integer(i4b), private :: n_discharge_call, iter_mar_coa
55  real(dp), private :: c_dis_0, s_dis, c_dis_fac
56  real(dp), private :: t_sub_pd, alpha_sub, alpha_o, m_h, m_d, r_mar_eff
57  real(dp), private :: t_sea_freeze
58  real(dp), public :: dt_glann, dt_sub
59 
60  integer(i2b), dimension(0:JMAX,0:IMAX), public :: mask_mar
61  real(dp), dimension(0:JMAX,0:IMAX), private :: c_dis
62  real(dp), dimension(0:JMAX,0:IMAX), public :: cst_dist, cos_grad_tc
63  real(dp), dimension(0:JMAX,0:IMAX), public :: dis_perp
64 
65  private
67 
68 contains
69 
70 !-------------------------------------------------------------------------------
71 !> Ice discharge parameters (Greenland).
72 !! [Assign ice discharge parameters.]
73 !<------------------------------------------------------------------------------
74  subroutine disc_param(dtime)
75 
76  ! Author: Reinhard Calov
77  ! Institution: Potsdam Institute for Climate Impact Research
78  ! Date: 10.6.16
79 
80  ! Purpose: Calculate discharge parameters.
81 
82  implicit none
83 
84  real(dp), intent(in) :: dtime
85 
86  real(dp) :: dtime_mar_coa
87 
88  disc = disc
89 
90  c_dis_0 = c_dis_0
91  c_dis_fac = c_dis_fac
92 
93  m_h = m_h
94  m_d = m_d
95 
96  r_mar_eff = r_mar_eff *1000.0_dp ! km -> m
97 
98 #if (defined(S_DIS))
99  s_dis = s_dis
100 #else
101  s_dis = 1.0_dp ! default value
102 #endif
103 
104 #if (defined(ALPHA_SUB))
105  alpha_sub = alpha_sub
106 #else
107  alpha_sub = 0.5_dp ! default value
108 #endif
109 
110 #if (defined(ALPHA_O))
111  alpha_o = alpha_o
112 #else
113  alpha_o = 1.0_dp ! default value
114 #endif
115 
116  t_sea_freeze = 0.0_dp - delta_tm_sw ! freezing temperature of sea water
117 
118  n_discharge_call = -1
119 
120 #if (defined(DTIME_MAR_COA0))
121  dtime_mar_coa = dtime_mar_coa0*year_sec ! a -> s
122 #else
123  stop ' >>> disc_param: DTIME_MAR_COA0 not defined in header file!'
124 #endif
125 
126  if (.not.approx_integer_multiple(dtime_mar_coa, dtime, eps_sp_dp)) &
127  stop ' >>> disc_param: dtime_mar_coa must be a multiple of dtime!'
128 
129  iter_mar_coa = nint(dtime_mar_coa/dtime)
130 
131 #if (GRID > 1)
132  stop ' >>> disc_param: GRID==2 not allowed for this application!'
133 #endif
134 
135 #if (ANF_DAT != 1 && defined(EXEC_MAKE_C_DIS_0))
136  stop ' >>> disc_param: EXEC_MAKE_C_DIS_0 defined requires ANF_DAT==1!'
137 #endif
138 
139  end subroutine disc_param
140 
141 !-------------------------------------------------------------------------------
142 !> Constant in ice discharge parameterization (Greenland).
143 !! [Determine (amount of magnitude of) constant in ice discharge
144 !! parameterization.]
145 !<------------------------------------------------------------------------------
146  subroutine calc_c_dis_0(dxi, deta)
148  ! Author: Reinhard Calov
149  ! Institution: Potsdam Institute for Climate Impact Research
150  ! Date: 8.6.16
151  ! Purpose: Compute c_dis_0 indirectly from present-day topography
152  ! from the condition that the parameterized total discharge
153  ! equals the discharge give by dis_target. Note: This c_dis_0 could
154  ! differ from the optimal one yielded from the siumualted ice thickness.
155  !
156  ! This is very useful, because c_dis_0 can vary a lot for different powers
157  ! m_D and m_H. This way we easier yield the approximately right value
158  ! for c_dis_0 for given powers m_D and m_H.
159 
160  ! This routine is not to be used regularly, and it is only executed if the
161  ! parameter EXEC_MAKE_C_DIS_0 is defined in the header file.
162 
163  implicit none
164 
165  real(dp), intent(in) :: dxi, deta
166 
167  integer(i4b) :: i, j
168  real(dp) :: disc_target, disc_tot, dT_sub
169 
170  ! targeted total ice discharge
171  ! 350 Gt/yr Calov et al. (2015) !!! 400 Gt/yr van den Broeke et al (2016)
172  disc_target = 350.0_dp ! in Gt/yr
173 
174  h_c=zs-zb; h_t=0.0_dp
175 
176  c_dis_0 = 1.0_dp
177  c_dis_fac = 1.0_dp
178  s_dis = 1.0_dp
179  call disc_fields()
180 
181  ! ensure that we get present-day discharge here (disc=1).
182  dt_glann = 0.0_dp
183  n_discharge_call = -1
184  call discharge(dxi, deta)
185 
186  ! disc_tot=0.0_dp
187  ! do i=0, IMAX
188  ! do j=0, JMAX
189  ! if(maske(j,i).eq.0.or.maske(j,i).eq.3) then
190  ! disc_tot=disc_tot+dis_perp(j,i)*area(j,i)
191  ! end if
192  ! end do
193  ! end do
194  !
195  ! disc_tot = disc_tot*YEAR_SEC ! m^3/s -> m^3/a
196 
197  disc_tot = 0.0_dp
198 
199  do i=1, imax-1
200  do j=1, jmax-1
201 
202  if (mask_mar(j,i) == 1) then
203  disc_tot = disc_tot + dis_perp(j,i)*area(j,i)
204  end if
205 
206  end do
207  end do
208 
209  disc_tot = disc_tot *year_sec *(rho/rho_w)
210  ! m^3/s ice equiv. -> m^3/a water equiv.
211 
212  disc_target = disc_target*1.0e+12_dp/rho_w ! Gt/a -> m^3/a water equiv.
213 
214  write(6, fmt='(a)') ' '
215  write(6, fmt='(a)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
216 
217  write(6, fmt='(a)') ' '
218  write(6, fmt='(3x,a,es12.4)') 'c_dis_0_init = ', c_dis_0
219 
220  c_dis_0 = c_dis_0 * disc_target/disc_tot
221 
222  write(6, fmt='(a)') ' '
223  write(6, fmt='(3x,a,es12.4)') 'disc_target = ', disc_target
224  write(6, fmt='(3x,a,es12.4)') 'disc_tot = ', disc_tot
225 
226  write(6, fmt='(a)') ' '
227  write(6, fmt='(3x,a,es12.4)') '--> c_dis_0 = ', c_dis_0
228 
229  end subroutine calc_c_dis_0
230 
231 !-------------------------------------------------------------------------------
232 !> Dependence of ice discharge coefficient on latitude (Greenland).
233 !! [Determine dependence of ice discharge coefficient on latitude.
234 !! This can be improved. For now I recommend s_dis=1.]
235 !<------------------------------------------------------------------------------
236  subroutine disc_fields()
238  ! Author: Reinhard Calov
239  ! Institution: Potsdam Institute for Climate Impact Research
240  ! Date: 8.6.16
241 
242  ! Purpose: Assign fields relevant for ice discharge parameterization
243  ! 1. Discharge coefficient field using linear dependence on
244  ! latitude.
245  ! 2. Subsurface temperature dependence on longitude and latitude.
246 
247  implicit none
248 
249  integer(i4b) :: i, j
250  real(dp) :: delta_phi=25.0_dp ! approximately the phi
251  ! spanning the middle of domain
252 
253  c_dis = c_dis_0*c_dis_fac &
254  *(1.0_dp-(1.0_dp-s_dis)*(phi*pi_180_inv-60.0_dp)/delta_phi)
255 
256  c_dis = max(c_dis, 0.0_dp)
257 
258  t_sub_pd = 3.0_dp ! Can depend on lambda, phi later on
259 
260  end subroutine disc_fields
261 
262 !-------------------------------------------------------------------------------
263 !> Ice discharge parameterization main formula, controler (general).
264 !! [Compute ice discharge via a parameterization using distance of ice
265 !! margin to coast and ice thickness as parameters.]
266 !<------------------------------------------------------------------------------
267  subroutine discharge(dxi, deta)
269  ! Authors: Reinhard Calov, Andrey Ganopolski
270  ! Institution: Potsdam Institute for Climate Impact Research
271  ! Date: 11.6.16
272 
273  implicit none
274 
275  real(dp), intent(in) :: dxi, deta
276 
277  real(dp), parameter :: alpha_tc=60.0_dp ! maximal allowed angle
278  real(dp) :: cos_tc
279 
280  n_discharge_call = n_discharge_call + 1
281 
282  cos_tc=dcos(alpha_tc*pi_180)
283 
284  if ((mod(n_discharge_call, iter_mar_coa)==0).or.(n_discharge_call==1)) then
285  write(6, fmt='(10x,a)') 'Computation of mask_mar, cst_dist, cos_grad_tc'
286  call marginal_ring(dxi, deta)
287  call coastal_distance(dxi, deta)
288  end if
289 
290  if(disc.ge.1) then !------------------ disc >= 1
291  where(mask_mar.eq.1.and.cos_grad_tc.ge.cos_tc)
292  dis_perp=c_dis*(h_c+h_t)**m_h/cst_dist**m_d
293  elsewhere
294  dis_perp=0.0_dp
295  end where
296  if(disc.eq.2) then !----------------- disc = 2
297  call calc_sub_oc_dt(t_sub_pd, dt_glann, dt_sub)
298  dis_perp=dis_perp*(1.0_dp+alpha_sub*dt_sub) ! actual discharge present-day one corrected
299  ! via sub-ocean temperature anomaly
300  end if
301  dis_perp=max(0.0_dp, dis_perp) ! ensure positive values
302  else if(disc.eq.0) then !-------------- disc = 0
303  dis_perp=0.0_dp
304  end if
305 
306  end subroutine discharge
307 
308 !-------------------------------------------------------------------------------
309 !> Distance to the coast (general).
310 !! [Compute distance to the coast for every land point.]
311 !<------------------------------------------------------------------------------
312  subroutine coastal_distance(dxi, deta)
313 
314  ! Author: Reinhard Calov
315  ! Institution: Potsdam Institute for Climate Impact Research
316  ! Date: 2.11.11
317  ! Methode: Options are:
318  ! i_search=1: Brute force.
319  ! i_search=2: Defining a square with two grid distances side length around the
320  ! actual grid point and enlanging the square successively by
321  ! one grid. Because this should have been a circle, the
322  ! square is enlarged again if the coast line was found; this time
323  ! enlargement is by a rectangles a the respective four side of
324  ! the square. Smaller side lenght equals the rectangle root of 2 of the
325  ! larger side.
326 
327  implicit none
328 
329  real(dp), intent(in) :: dxi, deta
330 
331  ! --------------------------------------------------
332  real(dp), parameter :: grid_dist=5000.0_dp
333  real(dp), parameter :: c_smooth=0.2_dp
334  integer(i4b), parameter :: i_search=2
335  ! --------------------------------------------------
336  integer(i4b) :: i, j, i_pos, j_pos, l, l_e, d_l
337  real(dp) :: dxi_inv, deta_inv
338  real(dp) :: cst_dist_tmp, qrt_1, qrt_2
339  logical :: leave_loop
340 
341  real(dp), allocatable, dimension(:,:) :: zs_tmp, grad_zs_x, grad_zs_y, &
342  grad_dist_x, grad_dist_y
343 
344  allocate(zs_tmp(0:jmax,0:imax), grad_zs_x(0:jmax,0:imax), grad_zs_y(0:jmax,0:imax), &
345  grad_dist_x(0:jmax,0:imax), grad_dist_y(0:jmax,0:imax))
346 
347  select case (i_search)
348 
349  case(1)
350  ! brute force computation of minimal distance to coast for every land point
351 
352  do i_pos=0, imax
353  do j_pos=0, jmax
354  if(maske(j_pos,i_pos).ne.2) then
355  cst_dist(j_pos,i_pos)=1.d+20
356  do i=0, imax
357  do j=0, jmax
358  if(maske(j,i).eq.2) then
359  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
360  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
361  cst_dist(j_pos,i_pos)=cst_dist_tmp
362  end if
363  end if
364  end do
365  end do
366  else
367  cst_dist(j_pos,i_pos)=grid_dist
368  end if
369  end do
370  end do
371 
372  case(2)
373 
374  do i_pos=0, imax
375  do j_pos=0, jmax
376  if(maske(j_pos,i_pos).ne.2) then
377  leave_loop=.false.
378  cst_dist(j_pos,i_pos)=1.d+20
379  do l=1, max(imax,jmax)
380  do i=max(i_pos-l,0), min(i_pos+l,imax)
381  j=max(j_pos-l,0); j=min(j,jmax)
382  if(maske(j,i).eq.2) then
383  leave_loop=.true.
384  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
385  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
386  cst_dist(j_pos,i_pos)=cst_dist_tmp
387  end if
388  end if
389  j=min(j_pos+l,jmax)
390  if(maske(j,i).eq.2) then
391  leave_loop=.true.
392  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
393  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
394  cst_dist(j_pos,i_pos)=cst_dist_tmp
395  end if
396  end if
397  end do
398  do j=max(j_pos-l+1,0), min(j_pos+l-1,jmax)
399  i=max(i_pos-l,0); i=min(i,imax)
400  if(maske(j,i).eq.2) then
401  leave_loop=.true.
402  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
403  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
404  cst_dist(j_pos,i_pos)=cst_dist_tmp
405  end if
406  end if
407  i=min(i_pos+l,imax)
408  if(maske(j,i).eq.2) then
409  leave_loop=.true.
410  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
411  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
412  cst_dist(j_pos,i_pos)=cst_dist_tmp
413  end if
414  end if
415  end do
416  if(leave_loop) then
417  l_e=l
418  d_l=ceiling(l_e*sqrt(2.0_dp)+0.001_dp)
419  exit ! leave loop in l
420  end if
421  end do
422  ! left
423  do i=max(i_pos-(l_e+d_l),0),min(i_pos-(l_e-1),imax)
424  do j=max(j_pos-l_e,0),min(j_pos+l_e,jmax)
425  if(maske(j,i).eq.2) then
426  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
427  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
428  cst_dist(j_pos,i_pos)=cst_dist_tmp
429  end if
430  end if
431  end do
432  end do
433  ! right
434  do i=max(i_pos+l_e+1,0),min(i_pos+l_e+d_l,imax)
435  do j=max(j_pos-l_e,0),min(j_pos+l_e,jmax)
436  if(maske(j,i).eq.2) then
437  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
438  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
439  cst_dist(j_pos,i_pos)=cst_dist_tmp
440  end if
441  end if
442  end do
443  end do
444  ! lower
445  do i=max(i_pos-l_e,0),min(i_pos+l_e,imax)
446  do j=max(j_pos-(l_e+d_l),0),min(j_pos-(l_e-1),jmax)
447  if(maske(j,i).eq.2) then
448  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
449  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
450  cst_dist(j_pos,i_pos)=cst_dist_tmp
451  end if
452  end if
453  end do
454  end do
455  ! upper
456  do i=max(i_pos-l_e,0),min(i_pos+l_e,imax)
457  do j=max(j_pos+l_e+1,0),min(j_pos+l_e+d_l,jmax)
458  if(maske(j,i).eq.2) then
459  cst_dist_tmp=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
460  if(cst_dist_tmp.le.cst_dist(j_pos,i_pos)) then
461  cst_dist(j_pos,i_pos)=cst_dist_tmp
462  end if
463  end if
464  end do
465  end do
466  else
467  cst_dist(j_pos,i_pos)=grid_dist
468  end if
469  end do
470  end do
471 
472  end select
473 
474  ! test GIS only
475 
476  if(.false.) then
477  do i=0, imax
478  do j=0, jmax
479  if(xi(i).le.-350000.0_dp &
480  .and.-2800000.0_dp.le.eta(j) &
481  .and.eta(j).le.-2250000.0_dp &
482  .or. &
483  xi(i).ge.100000.0_dp &
484  .and.-2280000.0_dp.le.eta(j) &
485  .and.eta(j).le.-1880000.0_dp &
486  .or. &
487  eta(j).ge.-1250000.0_dp &
488  .and.0.0_dp.le.xi(i) &
489  .and.xi(i).le.230000.0_dp) then
490  cst_dist(j,i)=0.3_dp*cst_dist(j,i)
491  end if
492  end do
493  end do
494  end if
495 
496  ! forbid zero, something like grid distance is reasonable here.
497  ! Dependents of coordinate system, for spherical coordinate system:
498  ! take length on Earth's surface. For now, simply just 5 km are taken.
499 
500  cst_dist=max(grid_dist, cst_dist)
501 
502  ! forbit to large values
503  cst_dist=min(1.0e10_dp, cst_dist)
504 
505  zs_tmp=zs
506 
507  ! smoothing the surface topography
508  do i=1, imax-1
509  do j=1, jmax-1
510  zs_tmp(j,i) = c_smooth*(zs(j,i-1)+zs(j,i+1)+zs(j-1,i)+zs(j+1,i)) &
511  +(1.0_dp-4.0_dp*c_smooth)*zs(j,i)
512  end do
513  end do
514 
515  ! gradients (2. order for now ...)
516 
517  dxi_inv = 1.0_dp/dxi
518  deta_inv = 1.0_dp/deta
519 
520 ! ------ x-derivatives
521 
522  do i=1, imax-1
523  do j=0, jmax
524  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)
525  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)
526  end do
527  end do
528 
529  do j=0, jmax
530  grad_zs_x(j,0) = (zs_tmp(j,1)-zs_tmp(j,0))*dxi_inv*insq_g11_g(j,0)
531  grad_zs_x(j,imax) = (zs_tmp(j,imax)-zs_tmp(j,imax-1))*dxi_inv*insq_g11_g(j,imax)
532  grad_dist_x(j,0) = (cst_dist(j,1)-cst_dist(j,0))*dxi_inv*insq_g11_g(j,0)
533  grad_dist_x(j,imax) = (cst_dist(j,imax)-cst_dist(j,imax-1))*dxi_inv*insq_g11_g(j,imax)
534  end do
535 
536 ! ------ y-derivatives
537 
538  do i=0, imax
539  do j=1, jmax-1
540  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)
541  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)
542  end do
543  end do
544 
545  do i=0, imax
546  grad_zs_y(0,i) = (zs_tmp(1,i)-zs_tmp(0,i))*deta_inv*insq_g22_g(0,i)
547  grad_zs_y(jmax,i) = (zs_tmp(jmax,i)-zs_tmp(jmax-1,i))*deta_inv*insq_g22_g(jmax,i)
548  grad_dist_y(0,i) = (cst_dist(1,i)-cst_dist(0,i))*deta_inv*insq_g22_g(0,i)
549  grad_dist_y(jmax,i) = (cst_dist(jmax,i)-cst_dist(jmax-1,i))*deta_inv*insq_g22_g(jmax,i)
550  end do
551 
552  ! angle between topographical gradient and coastal distance gradient
553 
554  do i=0, imax
555  do j=0, jmax
556  qrt_1=grad_zs_x(j,i)*grad_zs_x(j,i)+grad_zs_y(j,i)*grad_zs_y(j,i)
557  qrt_2=grad_dist_x(j,i)*grad_dist_x(j,i)+grad_dist_y(j,i)*grad_dist_y(j,i)
558  if(qrt_1.ne.0.0_dp.and.qrt_2.ne.0.0_dp) then
559  !!! 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))/ &
560  !!! (sqrt(qrt_1)*sqrt(qrt_2)) )
561 
562  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))/ &
563  (sqrt(qrt_1)*sqrt(qrt_2))
564  else
565  cos_grad_tc(j,i)=-1.0_dp
566  end if
567  end do
568  end do
569 
570  deallocate(zs_tmp, grad_zs_x, grad_zs_y, grad_dist_x, grad_dist_y)
571 
572  end subroutine coastal_distance
573 
574 !-------------------------------------------------------------------------------
575 !> Ring along an ice sheet margin (general).
576 !! [Compute marginal ring.]
577 !<------------------------------------------------------------------------------
578  subroutine marginal_ring(dxi, deta)
579 
580  ! Author: Reinhard Calov
581  ! Institution: Potsdam Institute for Climate Impact Research
582  ! Date: 28.10.11
583 
584  ! Purpose: Determine an r_max_eff wide ring arong ice margins
585  ! towards the interior of the ice sheet. Small stripes of land are
586  ! not considered as land. For the logics De Morgan is used often.
587 
588  ! Methode: Two staggered loops in i,j. The inner i,j loop acts inside
589  ! a rectangle defined by r_mar_eff. r_mar_eff should be small compared
590  ! to the domain size.
591 
592  implicit none
593 
594  real(dp), intent(in) :: dxi, deta
595 
596  integer(i4b) :: i, j, i_pos, j_pos, di_eff, dj_eff
597  integer(i4b) :: count_tmp
598  real(dp) :: r_p
599 
600  if(r_mar_eff.le.1.0e6_dp) then
601  mask_mar=0
602  do i_pos=1, imax-1
603  do j_pos=1, jmax-1
604  if((maske(j_pos,i_pos).eq.1.or.maske(j_pos,i_pos).eq.2).and. &
605  .not.(maske(j_pos,i_pos).eq.1.and.maske(j_pos,i_pos-1).eq.0 &
606  .and.maske(j_pos,i_pos+1).eq.0.or. &
607  maske(j_pos,i_pos).eq.1.and.maske(j_pos-1,i_pos).eq.0 &
608  .and.maske(j_pos+1,i_pos).eq.0.or. &
609  maske(j_pos,i_pos).eq.1.and.maske(j_pos,i_pos-1).eq.3 &
610  .and.maske(j_pos,i_pos+1).eq.3.or. &
611  maske(j_pos,i_pos).eq.1.and.maske(j_pos-1,i_pos).eq.3 &
612  .and.maske(j_pos+1,i_pos).eq.3.or. &
613  maske(j_pos,i_pos).eq.1.and.maske(j_pos,i_pos-1).eq.3 &
614  .and.maske(j_pos,i_pos+1).eq.0.or. &
615  maske(j_pos,i_pos).eq.1.and.maske(j_pos-1,i_pos).eq.0 &
616  .and.maske(j_pos+1,i_pos).eq.3)) then ! outside ice sheet, exclude isolated land stripes
617  di_eff=int(r_mar_eff/dxi)+1; dj_eff=int(r_mar_eff/deta)+1 ! only for grid=0, 1 yet!
618  do i=max(i_pos-di_eff,0), min(i_pos+di_eff,imax)
619  do j=max(j_pos-dj_eff,0), min(j_pos+dj_eff,jmax)
620  r_p=sqrt((xi(i_pos)-xi(i))**2+(eta(j_pos)-eta(j))**2)
621  if(r_p.le.r_mar_eff.and.(maske(j,i).eq.0.or.maske(j,i).eq.3)) then
622  mask_mar(j,i)=1
623  end if
624  end do
625  end do
626  end if
627  end do
628  end do
629  where(maske.ge.1)
630  mask_mar=1
631  end where
632  else ! the ring encompassed entire Greenland
633  mask_mar=1
634  end if
635 
636  end subroutine marginal_ring
637 
638 !-------------------------------------------------------------------------------
639 !> Anomaly of subsurface temperature (general).
640 !! [Compute anomaly of subsurface temperature with a simple parameterization.]
641 !<------------------------------------------------------------------------------
642  subroutine calc_sub_oc_dt(T_sub_PD, dT_glann, dT_sub)
643 
644  ! Author: Reinhard Calov
645  ! Institution: Potsdam Institute for Climate Impact Research
646  ! Purpose: Compute anomaly of subsurface temperature
647  ! with a simple parameterization
648  ! using global annual temperature
649  ! anomaly (e.g. from CLIMBER)
650 
651  implicit none
652 
653  real(dp), intent(in) :: T_sub_PD, dT_glann
654  real(dp), intent(out) :: dT_sub
655 
656  real(dp) :: T_sub
657 
658  dt_sub=alpha_o*dt_glann ! Sub-ocean temperature anomaly
659  t_sub=max(t_sea_freeze, t_sub_pd)+dt_sub ! Actual sub-ocean temperature
660  t_sub=max(t_sea_freeze, t_sub) ! Ocean temperature should stay
661  ! above freezing point
662  dt_sub=t_sub-t_sub_pd
663 
664  end subroutine calc_sub_oc_dt
665 
666 !-------------------------------------------------------------------------------
667 
668 end module discharge_workers_m
669 !
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
real(dp) rho_w
RHO_W: Density of pure water.
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...
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
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.