-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmo_grav_drain.f90
executable file
·281 lines (234 loc) · 11.6 KB
/
mo_grav_drain.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
!>
!! Computes the Salt fluxes caused by gravity drainage.
!!
!!
!!
!!
!! @author Philipp Griewank
!!
!!
!! COPYRIGHT
!!
!! This file is part of SAMSIM.
!!
!! SAMSIM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
!!
!! SAMSIM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License along with SAMSIM. If not, see <http://www.gnu.org/licenses/>.
!!
!
!!
!!
!! @par Revision History
!! Injected with life by Philipp Griewank, IMPRS (<2010-08-27>)
!!
!!
MODULE mo_grav_drain
USE mo_parameters
USE mo_thermo_functions, ONLY:func_S_br
USE mo_mass
IMPLICIT NONE
PUBLIC :: fl_grav_drain,fl_grav_drain_simple
CONTAINS
!>
!! Calculates fluxes caused by gravity drainage.
!!
!! If the Rayleigh number of a layer is higher then the critical value, brine leaves the layer by a theoretical brine channel.
!! The discharged brine flows downward through all layers directly into the underlying ocean.
!! To preserve mass the same amount of water flows upwards through all lower layers.
!! In contrast to the downward flux the upward flux is assumed to be in thermal equilibrium thus moving salt and heat to each layer.
!! The upward flux is a standard upwind advection.
!! The downward flux of a layer over the timestep is = x*(Ray-Ray_crit)*dt*thick.
!! _____________________________________
!! | ->__|__<-
!! |______________| v |__________________
!! | ^ ->| |<- ^
!! |______________| |__________________
!! | ^ ^ ->|||||<- ^ ^
!! |______________|vvv|__________________
!! ^ ^ ^ ^ ^ ^
!! This superb ascii art is supposed to show how the assumed fluxes flow downward through a brine channel and upwards through the layer.
!! x and r are passed on to enable easy optimization.
!! The effect of the upward moving brine is calculated in mass_transfer.
!!
!! IMPORTANT: The height assumptions are special. The bottom of the ice edge is assumed to be at psi_s(N_active)/psi_s_min *thick_0
!!
!! The first approach assumed that brine drainage occurred between two layers but performed poorly.
!!
!! If grav_heat_flag is set to 2 the amount of heat transported out of the ice will be compensated in the lowest layer
!!
!! @par Revision History
!! created by Philipp Griewank, IMPRS (2010-08-27) \n
!! Completely revised to assume brine channels by Philipp Griewank , IMPRS (2010-11-05) \n
!! Mass_transfer is used to advect H and S by Philipp Griewank, IMPRS (2010-11-05) \n
!! Added condition S_br(k)>S_br(k+1) by Philipp Griewank. IMPRS (2011-04-29) \n
!! Added harmonic mean for permeability by Philipp Griewank (2014-01-05)
SUBROUTINE fl_grav_drain (S_br,S_bu,psi_l,psi_s,psi_g,thick,S_abs,H_abs,T,m,dt,Nlayer,N_active,ray, &
T_bottom,S_bu_bottom,grav_drain,grav_temp,grav_salt,grav_heat_flag,harmonic_flag, fl_brine_bgc)
INTEGER, INTENT(in) :: Nlayer, N_active,grav_heat_flag,harmonic_flag
REAL(wp), DIMENSION(Nlayer), INTENT(in) :: S_br,S_bu,thick,T,psi_l,psi_s,psi_g
REAL(wp), INTENT(in) :: dt,S_bu_bottom,T_bottom
REAL(wp), INTENT(inout) :: grav_drain,grav_temp,grav_salt
REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: S_abs,H_abs,m
REAL(wp), DIMENSION(Nlayer-1), INTENT(out) :: ray !< Rayleigh number
REAL(wp), DIMENSION(N_active) :: fl_up !< Upward brine flux [kg] 1 is the flux from 2 to 1, N_active flux from ocean to N_active
REAL(wp), DIMENSION(N_active) :: fl_down !< Downward brine flux [kg] 1 is the flux from 1 to N_active, N_active is from N_active to the ocean
REAL(wp), DIMENSION(Nlayer) :: perm !< Permeability
REAL(wp), DIMENSION(Nlayer) :: harmonic_perm !< Harmonic mean permeability
REAL(wp), DIMENSION(Nlayer+1) :: fl_m
REAL(wp) :: flux !< Downward brine flux [kg]
REAL(wp) :: test1
REAL(wp) :: d_S_br,height,ray_mini
REAL(wp) :: heat_loss !< Amount of heat transported from the ice [J]
INTEGER :: k,kk
REAL(wp), DIMENSION(Nlayer+1,Nlayer+1), INTENT(inout),OPTIONAL :: fl_brine_bgc
ray_mini = ray_crit!Arises from the optimization of the r
heat_loss = 0._wp
perm = 0.0_wp
perm(N_active) = 9999999._wp !lowest layer is considered fully liquid
ray = 0._wp
fl_up = 0._wp
fl_down = 0._wp
harmonic_perm = 0.0_wp
!Permeability is determined
DO k=1,N_active
perm(k)=10.0_wp**(-17.0_wp) * (1000.0_wp*ABS(psi_l(k)))**3.10_wp
END DO
!Harmonic mean permeability is determined
IF (harmonic_flag == 2) THEN
DO k=1,N_active-1
test1 = minval(perm(k:N_active-1))
IF (test1 .lt. 10._wp**(-14._wp)) THEN
harmonic_perm(k)=0.0_wp
ELSE
DO kk=k,N_active-1
harmonic_perm(k)=harmonic_perm(k)+thick(kk)/perm(kk)
END DO
!bottom layer is included in a linear fashion
harmonic_perm(k)=harmonic_perm(k)+(thick(N_active)*psi_s(N_active)/psi_s_min)/perm(N_active)
harmonic_perm(k)=(SUM(thick(k:N_active-1))+thick(N_active)*psi_s(N_active)/psi_s_min)/harmonic_perm(k)
END IF
END DO
END IF
!___________Raleigh number is calculated (see Notz 2005)__________________________
DO k=1,N_active-1
d_S_br = S_br(k)-S_br(N_active)
height = SUM(thick((k+1):N_active-1))+thick(N_active)*psi_s(N_active)/psi_s_min !Height adjustment
IF (harmonic_flag == 1) THEN
ray(k) = grav*rho_l*bbeta*d_S_br*height* MINVAL(perm(k:N_active))
ELSE IF (harmonic_flag == 2) THEN
ray(k) = grav*rho_l*bbeta*d_S_br*height* harmonic_perm(k)
END IF
ray(k) = ray(k)/(kappa_l*mu)
ray(k) = MAX(ray(k),0._wp)
END DO
grav_salt = grav_salt + SUM(S_abs(:))
DO k = 1,N_active-1
IF (ray(k)>ray_mini .AND. psi_s(k)>0.001_wp .AND. S_abs(k)/m(k)>0.1_wp .AND. S_br(k)>S_br(k+1)) THEN !if psi_s equal zero then flushing is assumed
flux = x_grav*(ray(k)-ray_mini) *dt*thick(k)
flux = MIN(flux,psi_l(k)*rho_l*thick(k)) !should limit flux to brine volume of layer
S_abs(k) = S_abs(k) -flux*S_br(k)
IF (S_abs(k)<0.0_wp) THEN
PRINT*,'Negative Salinity due to gravity drainige overdrive in layer',k
PRINT*,S_abs(k),S_abs(k)+flux*S_br(k),flux,psi_l(k)*rho_l*thick(k)
STOP 21234
END IF
grav_temp = grav_temp + flux*T(k)
H_abs(k) = H_abs(k) -flux*c_l*T(k)
heat_loss = heat_loss +flux*c_l*T(k)
fl_down(k) = flux
FORALL(kk=k:N_active)
fl_up(kk) = fl_up(kk)+flux
END FORALL
fl_up(k) = MIN(fl_up(k),psi_l(k)*rho_l*thick(k)) ! should limit flux to brine volume of layer
END IF
END DO
grav_salt = grav_salt - SUM(S_abs(:))
!Influence of summed up upward transport resulting from the downward drainage
fl_m(1)=0.0_wp
fl_m(2:N_active+1)=fl_up(1:N_active)
IF( PRESENT( fl_brine_bgc ) ) THEN
fl_brine_bgc(1:(N_active-1),N_active+1) = fl_brine_bgc(1:(N_active-1),N_active) + fl_down(1:(N_active-1))
!fl_brine_bgc(N_active,N_active+1) = fl_brine_bgc(N_active,N_active+1) + sum(fl_down(:))
DO k = 1,N_active
fl_brine_bgc(k+1,k) = fl_brine_bgc(k+1,k) + fl_up(k)
END DO
END IF
CALL mass_transfer (Nlayer,N_active,T,H_abs,S_abs,S_bu,T_bottom,S_bu_bottom,fl_m)
grav_drain = grav_drain+fl_m(N_active+1)
!Heatflux correction
IF (grav_heat_flag==2) THEN
H_abs(N_active) = H_abs(N_active) +heat_loss-fl_up(N_active)*c_l*T_bottom
END IF
IF(MINVAL(S_abs)<0.0_wp) THEN
PRINT*,'negative salinity after gravity drainige, aborted',S_abs
STOP 1337
END IF
END SUBROUTINE fl_grav_drain
!>
!! Calculates salinity to imitate the effects gravity drainage.
!!
!! Based on the assumption that super critical Rayleigh numbers are quickly reduced below the critical Rayleigh number.
!! Proposed as a very simplified parametrisation of gravity drainage.
!! Includes no fluxes of any kind, instead bulk salinity is simply reduced when ever the Rayleigh number is above the critical values.
!! The parametrization begins from the bottom layers and moves upward.
!!
!! @par Revision History
!! created by Philipp Griewank, IMPRS (2012-01-01)
SUBROUTINE fl_grav_drain_simple (psi_s,psi_l,thick,S_abs,S_br,Nlayer,N_active,ray, &
grav_drain,harmonic_flag)
INTEGER, INTENT(in) :: Nlayer, N_active, harmonic_flag
REAL(wp), DIMENSION(Nlayer), INTENT(in) :: thick,psi_l,psi_s,S_br
REAL(wp), INTENT(inout) :: grav_drain
REAL(wp), DIMENSION(Nlayer), INTENT(inout) :: S_abs
REAL(wp), DIMENSION(Nlayer-1), INTENT(out) :: ray !< Rayleigh number
REAL(wp), DIMENSION(Nlayer) :: perm !< Permeability
REAL(wp), DIMENSION(Nlayer) :: harmonic_perm !< Harmonic permeability
REAL(wp) :: d_S_br,height,ray_mini,temp
INTEGER :: k,kk
ray_mini = ray_crit
perm = 0.0_wp
perm(N_active) = 9999999._wp
ray = 0._wp
!Permeability is determined
DO k = 1,N_active
perm(k)=10.0_wp**(-17.0_wp) * (1000.0_wp*ABS(psi_l(k)))**3.10_wp
END DO
!Harmonic mean permeability is determined
IF (harmonic_flag == 2) THEN
DO k=1,N_active-1
temp = minval(perm(k:N_active-1))
IF (temp .lt. 10._wp**(-14._wp)) THEN
harmonic_perm(k)=0.0_wp
ELSE
DO kk=k,N_active-1
harmonic_perm(k)=harmonic_perm(k)+thick(kk)/perm(kk)
END DO
!bottom layer is included in a linear fashion
harmonic_perm(k)=harmonic_perm(k)+(thick(N_active)*psi_s(N_active)/psi_s_min)/perm(N_active)
harmonic_perm(k)=(SUM(thick(k:N_active-1))+thick(N_active)*psi_s(N_active)/psi_s_min)/harmonic_perm(k)
END IF
END DO
END IF
!___________Raleigh number is calculated (see Notz 2005)__________________________
DO k=1,N_active-1
d_S_br = S_br(k)-S_br(N_active)
height = SUM(thick((k+1):N_active-1))+thick(N_active)*psi_s(N_active)/psi_s_min !Height adjustment
IF (harmonic_flag == 1) THEN
ray(k) = grav*rho_l*bbeta*d_S_br*height* MINVAL(perm(k:N_active))
ELSE IF (harmonic_flag == 2) THEN
ray(k) = grav*rho_l*bbeta*d_S_br*height* harmonic_perm(k)
END IF
ray(k) = ray(k)/(kappa_l*mu)
ray(k) = MAX(ray(k),0._wp)
END DO
!##########Where the magic happens#################################
DO k = N_active-1,1,-1
IF (ray(k) > ray_mini ) THEN
S_abs(k) = S_abs(k)*0.99
END IF
END DO
grav_drain = 0._wp
END SUBROUTINE fl_grav_drain_simple
END MODULE mo_grav_drain