forked from altMITgcm/MITgcm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
salt_plume_volfrac.F
189 lines (171 loc) · 6.02 KB
/
salt_plume_volfrac.F
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
C $Header: /u/gcmpack/MITgcm/pkg/salt_plume/salt_plume_volfrac.F,v 1.2 2014/05/21 13:59:39 atn Exp $
C $Name: $
#include "SALT_PLUME_OPTIONS.h"
CBOP
C !ROUTINE: SALT_PLUME_VOLFRAC
C !INTERFACE:
SUBROUTINE SALT_PLUME_VOLFRAC(
I bi, bj, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE SALT_PLUME_VOLFRAC
C | o Compute saltplume penetration.
C *==========================================================*
C | Compute fraction of volume flux associated with saltplume
C | flux penetrating through the entire water columns due to
C | rejected salt during freezing.
C |
C | For example, if surface value is Saltplume0,
C | and each level gets equal fraction 1/5 down to SPDepth=5,
C | SALT_PLUME_VOLFRAC will report
C | dSPvolkLev2Above[2to1,3to2,4to3,5to4,6to5] = [4/5,3/5,2/5,1/5, 0]
C | dSPvolSurf2kLev [1to1,1to2,1to3,1to4,1to5] = [1/5,1/5,1/5,1/5,1/5]
C | sum [into5] = 1to5 + 6to5 - 5to4 = 1/5 + 0 - 1/5 = 0
C | [into4] = 1to4 + 5to4 - 4to3 = 1/5 + 1/5 - 2/5 = 0
C | [into3] = 1to3 + 4to3 - 3to2 = 1/5 + 2/5 - 3/5 = 0
C | [into2] = 1to2 + 3to2 - 2to1 = 1/5 + 3/5 - 4/5 = 0
C | [into1] = 1to1 + 2to1 - 1to[1,2,3,4,5] = 1/5 + 4/5 - 5/5 = 0
C | NOTE: volume will always be conserved.
C | =====
C | Written by : ATN (based on SALT_PLUME_FRAC)
C | Date : Apr 14, 2014
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
#include "SIZE.h"
#include "GRID.h"
#include "SALT_PLUME.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
C !INPUT/OUTPUT PARAMETERS:
C input arguments
C SPDpeth :: corresponding SaltPlumeDepth(i,j) at this grid point
C myTime :: Current time in simulation
C myIter :: Current iteration number in simulation
C myThid :: My Thread Id. number
INTEGER bi,bj
_RL myTime
INTEGER myIter
INTEGER myThid
C input/output arguments
C CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP
#ifdef ALLOW_SALT_PLUME
#ifdef SALT_PLUME_VOLUME
C !LOCAL VARIABLES:
_RL dMbdt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL temp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dplumek
INTEGER SPkBottom (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
INTEGER i,j,k,kp1,Nlev,Nrp1
INTEGER imt
parameter( imt=(sNx+2*OLx)*(sNy+2*OLy) )
C initialize at every time step
Nrp1=Nr+1
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
dSPvolSurf2kLev(i,j,k,bi,bj) = 0. _d 0
dSPvolkLev2Above(i,j,k,bi,bj) = 0. _d 0
SPplumek(i,j,k,bi,bj) = 1. _d 0
ENDDO
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
SPplumek(i,j,Nrp1,bi,bj) = 1. _d 0
SPbrineVolFlux(i,j,bi,bj) = 0. _d 0
SPkBottom(i,j) = 0
ENDDO
ENDDO
C call salt_plume_frac to fill in SPplumek and SPkBottom
C use dMbdt+temp as a temporary arrays here to save memory:
DO k = Nrp1,1,-1
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
temp(i,j)=SaltPlumeDepth(i,j,bi,bj)
dMbdt(i,j)=abs(rF(k))
ENDDO
ENDDO
CALL SALT_PLUME_FRAC(
I imt,oneRS,temp,
#ifdef SALT_PLUME_SPLIT_BASIN
I XC(1-Olx,1-Oly,bi,bj),YC(1-Olx,1-Oly,bi,bj),
#endif
U dMbdt,
I myTime, 1, myThid )
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
SPplumek(i,j,k,bi,bj)=dMbdt(i,j)
IF(SPplumek(i,j,k,bi,bj).GT. 0.9999999) THEN
SPkBottom(i,j)=k
ENDIF
ENDDO
ENDDO
ENDDO
C reinitialize dMbdt = 0
DO j=1-Oly,sNy+Oly
DO i=1-Olx,sNx+Olx
dMbdt(i,j)=0. _d 0
ENDDO
ENDDO
C Now calculating dplumek, dSPvolumeUp, dSPvolSurf2kLev
C units:
C Sbrine=dsb/dt*dt/(rhoConst*SPalpha*drF)[psu kg/m2/s*s/(kg/m3*m)]=[psu]
C SPplumek : fraction : unitless
C SaltPlumeFlux: dsb/dt [psu.kg/m^2/s = g/m^2/s]
C brine_mass_flux dMb/dt = dsb/dt / Sbrine [kg/m2/s]
C = dsb/dt / (dsb/dt*dt/(rhoConst*SPalpha*drF))
C = rhoConst*SPalpha*drF/dt [kg/m3 m/s]=[kg/m2/s]
C dVbrine/dt = dMb/dt 1/rhoConst [m/s]
C has 2 ways to define brine properties: either provide
C (A) SPalpha: vol frac or (B) SPbrineSalt: brine salinity.
C (A) SPalpha: can calc SPbrineSalt as fxn of dhice/dt,
C constrained by SPbrineSaltmax:
C SPbrineSalt=SaltPlumeFlux/rhoConst/SPalpha/drF(1)*dt
C SPbrineSalt=min(SPbrineSalt,SPbrineSaltmax)
C dMbdt = saltPlumeFlux / SPbrineSalt
C = rhoConst*SPalpha*drF(1)/dt <-- a function of SPalpha
C (B) SPbrinesalt provided
C dMbdt = saltPlumeFlux / SPbrineSalt <-- fxn of SPbrineSalt
C Assuming we go with (B) here:
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
C brine mass and volume at surface:
dMbdt(i,j)=saltPlumeFlux(i,j,bi,bj)/SPbrineSconst
SPbrineVolFlux(i,j,bi,bj)=dMbdt(i,j)*mass2rUnit
ENDDO
ENDDO
C Distributing down: this is always from level 1 to depth
DO k=Nr,1,-1
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
dplumek=SPplumek(i,j,k+1,bi,bj)-SPplumek(i,j,k,bi,bj)
dSPvolSurf2kLev(i,j,k,bi,bj)=dplumek*SPbrineVolFlux(i,j,bi,bj)
ENDDO
ENDDO
ENDDO
C Now volume up: need to scan from bottom of SPDepth
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
Nlev=SPkBottom(i,j)
IF(Nlev.GE.1 .AND. Nlev.LE.Nr) THEN
DO k=Nlev,1,-1
kp1=k+1
dSPvolkLev2Above(i,j,k,bi,bj)=dSPvolkLev2Above(i,j,kp1,bi,bj)
& - dSPvolSurf2kLev(i,j,k,bi,bj)
ENDDO
ENDIF
ENDDO
ENDDO
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(
& SPplumek,'PLUMEKB1',0,Nr,1,bi,bj,myThid )
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
#endif /* SALT_PLUME_VOLUME */
#endif /* ALLOW_SALT_PLUME */
RETURN
END