SICOPOLIS V5-dev  Revision 1368
calc_thk_water_bas_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : c a l c _ t h k _ w a t e r _ b a s _ m
4 !
5 !> @file
6 !!
7 !! Computation of the thickness of the water column under the ice base.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2018 Ralf Greve
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 !> Computation of the thickness of the water column under the ice base.
34 !<------------------------------------------------------------------------------
36 
37  use sico_types_m
39  use sico_vars_m
40 
41 #if (BASAL_HYDROLOGY==1)
42  use hydro_m
43 #endif
44 
45  implicit none
46 
47  private
48  public :: calc_thk_water_bas
49 
50 contains
51 
52 !-------------------------------------------------------------------------------
53 !> Main subroutine of calc_thk_water_bas_m:
54 !! Computation of the thickness of the water column under the ice base.
55 !<------------------------------------------------------------------------------
56  subroutine calc_thk_water_bas(z_sl)
57 
58  implicit none
59 
60  real(dp), intent(in) :: z_sl
61 
62 #ifdef ALLOW_OPENAD /* OpenAD */
63  integer(i4b) :: i, j
64 #endif /* OpenAD */
65 
66  logical, save :: firstcall = .true.
67 
68 #if (BASAL_HYDROLOGY==1)
69  real(dp), save :: rho_rho_w_ratio
70  integer , dimension(0:IMAX,0:JMAX) :: hydro_icemask
71  real(dp), dimension(0:IMAX,0:JMAX) :: hydro_topg, hydro_thk, &
72  hydro_temppabase, hydro_supply, &
73  hydro_sflux, hydro_vflux, &
74  hydro_vfluxX, hydro_vfluxY, &
75  hydro_bwat
76 
77 #ifdef ALLOW_OPENAD /* OpenAD */
78  integer(i1b), dimension(0:IMAX,0:JMAX) :: t_maske
79  real(dp) , dimension(0:IMAX,0:JMAX) :: t_H_c, t_H_t, t_Q_b_tot
80 #endif /* OpenAD */
81 
82  type(hydro_t), save :: hydro
83  !!! RG: Does this need a save attribute?
84 #endif
85 
86 !-------- Water column --------
87 
88 #ifndef ALLOW_OPENAD /* Normal */
89 
90 #if (BASAL_HYDROLOGY==1)
91 
92  if (firstcall) then
93 
94  rho_rho_w_ratio = rho/rho_w
95 
96  call hydro_init(hydro, xi, eta)
97  call hydro_gen_conf(hydro, &
98  & method='quinn', &
99  & avoid_frz=.false., &
100  & filter_len=0.0_dp, &
101  & rho_seawater=rho_sw, &
102  & rho_freshwater=rho_w, &
103  & rho_ice=rho)
104 
105  end if
106 
107  hydro_topg = transpose(zl)-z_sl
108  hydro_temppabase = transpose(temph_b)
109 
110  where (transpose(maske)==0_i1b) ! grounded ice
111  hydro_icemask = 1
112  hydro_thk = transpose(h_c+h_t)
113  hydro_supply = rho_rho_w_ratio*transpose(q_b_tot)
114  elsewhere
115  hydro_icemask = 0
116  hydro_thk = 0.0_dp
117  hydro_supply = 0.0_dp
118  end where
119 
120  call hydro_set_topg(hydro, hydro_topg)
121  call hydro_set_thk(hydro, hydro_thk)
122  call hydro_set_temppabase(hydro, hydro_temppabase)
123  call hydro_set_supply(hydro, hydro_supply)
124  call hydro_set_mask(hydro, hydro_icemask)
125 
126  call hydro_update(hydro)
127 
128  call hydro_get_sflux(hydro, hydro_sflux)
129  call hydro_get_vflux(hydro, hydro_vflux, hydro_vfluxx, hydro_vfluxy)
130  call hydro_get_bwat(hydro, hydro_bwat)
131 
132  q_w = transpose(hydro_vflux)
133  q_w_x = transpose(hydro_vfluxx)
134  q_w_y = transpose(hydro_vfluxy)
135  h_w = transpose(hydro_bwat)
136 
137 #else
138 
139  where (maske==0_i1b) ! grounded ice
140  q_w = 0.0_dp
141  q_w_x = 0.0_dp
142  q_w_y = 0.0_dp
143  h_w = 0.0_dp
144  end where
145 
146 #endif
147 
148  where (maske==2_i1b) ! ocean
149  q_w = 0.0_dp
150  q_w_x = 0.0_dp
151  q_w_y = 0.0_dp
152  h_w = z_sl-zl
153  elsewhere (maske==3_i1b) ! floating ice
154  q_w = 0.0_dp
155  q_w_x = 0.0_dp
156  q_w_y = 0.0_dp
157  h_w = zb-zl
158  elsewhere (maske==1_i1b) ! ice-free land
159  q_w = 0.0_dp
160  q_w_x = 0.0_dp
161  q_w_y = 0.0_dp
162  h_w = 0.0_dp
163  end where
164 
165  if (firstcall) firstcall = .false.
166 
167 #else /* OpenAD */
168 
169 #if (BASAL_HYDROLOGY==1)
170 
171  if (firstcall) then
172 
173  rho_rho_w_ratio = rho/rho_w
174 
175  call hydro_init(hydro, xi, eta)
176  call hydro_gen_conf(hydro, &
177  & method='quinn', &
178  & avoid_frz=.false., &
179  & filter_len=0.0_dp, &
180  & rho_seawater=rho_sw, &
181  & rho_freshwater=rho_w, &
182  & rho_ice=rho)
183 
184  end if
185 
186  ! hard-coding transpose
187  do j=0,jmax
188  do i=0,imax
189  hydro_topg(i,j) = zl(j,i) - z_sl
190  hydro_temppabase(i,j) = temph_b(j,i)
191  ! transpose these arrays for easy searching below
192  t_maske(i,j) = maske(j,i)
193  t_h_c(i,j) = h_c(j,i)
194  t_h_t(i,j) = h_t(j,i)
195  t_q_b_tot(i,j) = q_b_tot(j,i)
196  end do
197  end do
198 
199  do j=0,jmax
200  do i=0,imax
201  if (t_maske(i,j)==0_i1b) then
202  hydro_icemask(i,j) = 1
203  hydro_thk(i,j) = t_h_c(i,j) + t_h_t(i,j)
204  hydro_supply(i,j) = rho_rho_w_ratio*t_q_b_tot(i,j)
205  else
206  hydro_icemask(i,j) = 0
207  hydro_thk(i,j) = 0.0_dp
208  hydro_supply(i,j) = 0.0_dp
209  end if
210  end do
211  end do
212 
213  call hydro_set_topg(hydro, hydro_topg)
214  call hydro_set_thk(hydro, hydro_thk)
215  call hydro_set_temppabase(hydro, hydro_temppabase)
216  call hydro_set_supply(hydro, hydro_supply)
217  call hydro_set_mask(hydro, hydro_icemask)
218 
219  call hydro_update(hydro)
220 
221  call hydro_get_sflux(hydro, hydro_sflux)
222  call hydro_get_vflux(hydro, hydro_vflux, hydro_vfluxx, hydro_vfluxy)
223  call hydro_get_bwat(hydro, hydro_bwat)
224 
225  do i=0,imax
226  do j=0,jmax
227  q_w(j,i) = hydro_vflux(i,j)
228  q_w_x(j,i) = hydro_vfluxx(i,j)
229  q_w_y(j,i) = hydro_vfluxy(i,j)
230  h_w(j,i) = hydro_bwat(i,j)
231  end do
232  end do
233 
234 #else
235 
236  do i=0,imax
237  do j=0,jmax
238  if ( maske(j,i)==0_i1b ) then ! grounded ice
239  q_w(j,i) = 0.0_dp
240  q_w_x(j,i) = 0.0_dp
241  q_w_y(j,i) = 0.0_dp
242  h_w(j,i) = 0.0_dp
243  end if
244  end do
245  end do
246 
247 #endif
248 
249  do i=0,imax
250  do j=0,jmax
251 
252  if ( maske(j,i)==2_i1b ) then
253  q_w(j,i) = 0.0_dp
254  q_w_x(j,i) = 0.0_dp
255  q_w_y(j,i) = 0.0_dp
256  h_w(j,i) = z_sl-zl(j,i)
257  else if ( maske(j,i)==3_i1b ) then
258  q_w(j,i) = 0.0_dp
259  q_w_x(j,i) = 0.0_dp
260  q_w_y(j,i) = 0.0_dp
261  h_w(j,i) = zb(j,i)-zl(j,i)
262  else if ( maske(j,i)==1_i1b ) then
263  q_w(j,i) = 0.0_dp
264  q_w_x(j,i) = 0.0_dp
265  q_w_y(j,i) = 0.0_dp
266  h_w(j,i) = 0.0_dp
267  end if
268 
269  end do
270  end do
271 
272  if (firstcall) firstcall = .false.
273 
274 #endif /* Normal vs. OpenAD */
275 
276  end subroutine calc_thk_water_bas
277 
278 !-------------------------------------------------------------------------------
279 
280 end module calc_thk_water_bas_m
281 !
subroutine, public calc_thk_water_bas(z_sl)
Main subroutine of calc_thk_water_bas_m: Computation of the thickness of the water column under the i...
real(dp) rho_w
RHO_W: Density of pure water.
Computation of the thickness of the water column under the ice base.
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
Declarations of kind types for SICOPOLIS.
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...
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 ...
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_w
q_w(j,i): Scalar volume flux of the basal water
real(dp), dimension(0:jmax, 0:imax) q_w_y
q_w_y(j,i): Scalar volume flux of the basal water in y-direction
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) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
real(dp) rho_sw
RHO_SW: Density of sea water.
real(dp), dimension(0:jmax, 0:imax) q_w_x
q_w_x(j,i): Scalar volume flux of the basal water in x-direction
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...
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.