/[MITgcm]/MITgcm_contrib/dcarroll/highres_darwin/code/shelfice_solve4fluxes.F
ViewVC logotype

Contents of /MITgcm_contrib/dcarroll/highres_darwin/code/shelfice_solve4fluxes.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Sun Sep 22 21:23:47 2019 UTC (5 years, 10 months ago) by dcarroll
Branch: MAIN
CVS Tags: HEAD
Error occurred while calculating annotation data.
Initial check in of high resolution Darwin simulation code

1 #include "SHELFICE_OPTIONS.h"
2 #ifdef ALLOW_AUTODIFF
3 # include "AUTODIFF_OPTIONS.h"
4 #endif
5 #ifdef ALLOW_CTRL
6 # include "CTRL_OPTIONS.h"
7 #endif
8
9 CBOP
10 C !ROUTINE: SHELFICE_SOLVE4FLUXES
11 C !INTERFACE:
12 SUBROUTINE SHELFICE_SOLVE4FLUXES(
13 I tLoc, sLoc, pLoc,
14 I gammaT, gammaS,
15 I iceConductionDistance, thetaIceConduction,
16 O heatFlux, fwFlux,
17 O forcingT, forcingS,
18 I bi, bj, myTime, myIter, myThid )
19
20 C !DESCRIPTION: \bv
21 C *==========================================================*
22 C | SUBROUTINE SOLVE4FLUXES
23 C | o Calculate
24 C *==========================================================*
25 C \ev
26
27 C !USES:
28 IMPLICIT NONE
29
30 C === Global variables ===
31 #include "SIZE.h"
32 #include "EEPARAMS.h"
33 #include "PARAMS.h"
34 #include "GRID.h"
35 #include "DYNVARS.h"
36 #include "FFIELDS.h"
37 #include "SHELFICE.h"
38 #include "SHELFICE_COST.h"
39 #ifdef ALLOW_AUTODIFF
40 # include "CTRL_SIZE.h"
41 # include "ctrl.h"
42 # include "ctrl_dummy.h"
43 #endif /* ALLOW_AUTODIFF */
44 #ifdef ALLOW_AUTODIFF_TAMC
45 # ifdef SHI_ALLOW_GAMMAFRICT
46 # include "tamc.h"
47 # include "tamc_keys.h"
48 # endif /* SHI_ALLOW_GAMMAFRICT */
49 #endif /* ALLOW_AUTODIFF_TAMC */
50
51 C !INPUT PARAMETERS:
52 C tLoc ::
53 C sLoc ::
54 C pLoc ::
55 C insitutLoc ::
56 C gammaT ::
57 C gammaS ::
58 C bi,bj :: tile indices
59 C myTime :: current time in simulation
60 C myIter :: iteration number in simulation
61 C myThid :: my Thread Id number
62
63 C !OUTPUT PARAMETERS:
64 C heatFlux ::
65 C fwFlux ::
66 C forcingT ::
67 C forcingS ::
68 C----------
69
70 _RL tLoc, sLoc, pLoc, insitutLoc
71 _RL gammaT, gammaS
72 _RL iceConductionDistance, thetaIceConduction
73 _RL heatFlux, fwFlux, forcingS, forcingT
74 INTEGER i, j, bi, bj
75 _RL myTime
76 INTEGER myIter, myThid
77 character*200 msgBuf
78 CEOP
79
80 #ifndef ALLOW_OPENAD
81 _RL SW_TEMP
82 EXTERNAL SW_TEMP
83 #endif
84
85
86 C !LOCAL VARIABLES:
87 C === Local variables ===
88 _RL thetaFreeze, saltFreeze
89 _RL eps1, eps2, eps3, eps4, eps5, eps6, eps7, eps8
90 _RL aqe, bqe, cqe, discrim, recip_aqe
91 _RL a0, a1, a2, b0, c0
92 _RL w_I
93
94 C === Useful Units ===
95 C-- gammaT, m s^-1
96 C-- gammaS, m s^-1
97 C-- rUnit2mass (rhoConst), kg m^-3
98 C-- mass2rUnit (recip_rhoConst), m^3 kg^-1
99 C-- eps3, W K^-1 m^-2
100 C-- fwFlux, kg m^-2 s^-1
101 C-- heatFlux, kg m^-2 s^-1
102 C-- forcing T, K m/s
103 C-- forcing S, psu m/s
104 C-- SHELFICEkappa, m^2/s
105
106 C-- eps1, W K^-1 m^-2 : kg/m^3 * J/kg/K * m/s
107 C-- eps2, W m^-2 : kg/m^3 * J/kg * m/s
108 C-- eps3, W K^-1 m^-2 : kg m^-3 * J kg^-1 K^-1 * m^2 s^-1 * m^-1
109
110 C-- fwFlux : fresh water flux due to melting (kg m^-2 s^-1)
111
112 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
113 C linear dependence of freezing point on salinity
114 a0 = -0.0575 _d 0
115 c0 = 0.0901 _d 0
116 b0 = -7.61 _d -4
117
118 C-- convert potential T into in-situ T relative to surface
119 insitutLoc = SW_TEMP(sLoc,tLoc,pLoc, zeroRL)
120
121 C-- DEFINE SOME CONSTANTS
122 eps1 = rUnit2mass*HeatCapacity_Cp * gammaT
123 eps2 = rUnit2mass*SHELFICElatentHeat * gammaS
124 eps3 = rhoShelfIce*SHELFICEheatCapacity_Cp * SHELFICEkappa
125 & /iceConductionDistance;
126 eps4 = b0*pLoc + c0
127 eps6 = eps4 - insitutLoc
128 eps7 = eps4 - thetaIceConduction
129
130 IF ( debugLevel.GE.debLevE ) THEN
131 WRITE(msgBuf,'(A25,9E16.8)')
132 &'ZZZ7 r2mass, Cp, gmaT,SIlh,gmaS, rhoSI, SI_Cp, Kap, ICondDis ',
133 & rUnit2mass,HeatCapacity_Cp,gammaT, SHELFICElatentHeat,gammaS,
134 & rhoShelfIce, SHELFICEheatCapacity_Cp, SHELFICEkappa,
135 & iceConductionDistance
136
137 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
138 & SQUEEZE_RIGHT , 1)
139 ENDIF
140
141 C-- Default thermodynamics specify a linear T gradient
142 C-- through the ice (Holland and Jenkins, 1999, Section 2.
143 aqe = a0*(eps1 + eps3)
144 bqe = eps1*eps6 + eps3*eps7 - eps2 - SHELFICEsalinity*aqe
145 cqe = eps2 * sLoc - SHELFICEsalinity*(eps1*eps6 + eps3*eps7)
146
147 C-- Alterantively, we can have a nonlinear T gradient
148 C-- through the ice (Holland and Jenkins, 1999, Section 3.
149 C-- This demands a different set of constants
150 IF (SHELFICEadvDiffHeatFlux .EQV. .TRUE.) THEN
151
152 IF ( debugLevel.GE.debLevE ) THEN
153 WRITE(msgBuf,'(A)')
154 & '1useSHELFICEadvDiffHeatFlux '
155
156 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
157 & SQUEEZE_RIGHT , 1)
158 ENDIF
159
160 eps8 = rUnit2mass * gammaS * SHELFICEheatCapacity_Cp
161 aqe = a0 *(eps1 - eps8)
162 bqe = eps1*eps6 + sLoc*eps8*a0 - eps8*eps7 - eps2 -
163 & SHELFICEsalinity*eps1*a0
164 cqe = sLoc*(eps8*eps7 + eps2) - SHELFICEsalinity*eps1
165 ENDIF
166
167 C solve quadratic equation for salinity at shelfice-ocean interface
168 recip_aqe = 0. _d 0
169 IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe
170
171 C-- Sb = \frac{-bqe \pm \sqrt(bqe^2 - 4 aqc cqe)}{2 aqe}
172 discrim = bqe*bqe - 4. _d 0*aqe*cqe
173
174 C-- Try the negative root (- SQRT(discrim)) of the quadratic eq.
175 saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
176
177 C--- If the negative root yields a negative salinity, then use the
178 C-- positive root (+ SQRT(discrim))
179 IF ( saltFreeze .LT. 0. _d 0 ) THEN
180 saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
181 ENDIF
182
183 C-- in situ seawater freezing point using linearization
184 thetaFreeze = a0*saltFreeze + eps4
185
186 C-- Calculate the upward heat and fresh water fluxes;
187 C-- MITgcm sign conventions: downward (negative) fresh water flux
188 C-- implies melting and due to upward (positive) heat flux
189
190 C-- This formulation of fwflux, derived from the heat balance equation
191 C-- instead of the salt balance equation, can handle the case when the
192 C-- salinity of the ocean, boundary layer, and ice are identical.
193 fwFlux = 1/SHELFICElatentHeat*(
194 & eps3*(thetaFreeze - thetaIceConduction) -
195 & eps1*(insitutLoc - thetaFreeze) )
196
197 C-- If a nonlinear local ice T gradient near the ice-ocean interface
198 C-- is allowed and fwflux is positive (ice growth) then
199 C-- we must solve the quadratic equation using a different set of
200 C-- coeffs (Holland Jenkins, 1999).
201 IF ((SHELFICEadvDiffHeatFlux .EQV. .TRUE.) .AND.
202 & (fwFlux .GT. zeroRL)) THEN
203 IF ( debugLevel.GE.debLevE ) THEN
204 WRITE(msgBuf,'(A)')
205 & '2useSHELFICEadvDiffHeatFlux '
206
207 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
208 & SQUEEZE_RIGHT , 1)
209 ENDIF
210 aqe = a0 *(eps1)
211 bqe = eps1*eps6 - eps2 - SHELFICEsalinity*eps1*a0
212 cqe = sLoc*(eps2) - SHELFICEsalinity*eps1
213
214 recip_aqe = 0. _d 0
215 IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe
216
217 discrim = bqe*bqe - 4. _d 0*aqe*cqe
218 saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
219 IF ( saltFreeze .LT. 0. _d 0 ) THEN
220 saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
221 ENDIF
222
223 thetaFreeze = a0*saltFreeze + eps4
224
225 fwFlux = 1/SHELFICElatentHeat* (
226 & eps3*(thetaFreeze - thetaIceConduction) -
227 & eps1*(insitutLoc - thetaFreeze) )
228 ENDIF
229
230 C meltwater advective flux velocity at ice-ocean interface (m/s)
231 C * negative corresponds to downward flux (melting)
232 w_I = fwFlux * mass2rUnit
233
234 C-- Calculate the upward heat fluxes:
235 C-- melting requires upward (positive) heat flux from ocean to ice.
236
237 C-- The heatFlux variable corresponds with the change of energy in the
238 C-- ocean grid cell volume. In the conservative case (J2001),
239 C-- advective heat fluxes change the energy of the volume whereas in
240 C-- the non-conservative case there are no advective heat fluxes
241 C-- melting or freezing have no associated advective heat fluxes.
242 IF (SHELFICEconserve) THEN
243
244 C-- In the conservative case (J2001) there are two cases, fixed and
245 C-- non-fixed ocean volume.
246
247 IF (useRealFreshWaterFlux ) THEN
248
249 C-- If the ocean volume can change (realFWFlux=true) then advection of
250 C-- meltwater does not displace water at T=insitutLoc in the cell and the
251 C-- heat flux correpsonding to the total energy flux of the volume
252 C-- consists of only two terms: turbulent fluxes (positive out)
253 C-- and advective meltwater fluxes (negative in).
254 heatFlux = rUnit2mass*HeatCapacity_Cp * (
255 & gammaT * (insitutLoc - thetaFreeze)
256 & + w_I * thetaFreeze )
257 ELSE
258
259 C-- If the volume is fixed (realFWFlux=false) then the advection of
260 C-- meltwater does displace ambient water at T=insitutLoc in the cell.
261 C-- Displacement reduction volume energy by w_I * insitutLoc (positive)
262 heatFlux = rUnit2mass*HeatCapacity_Cp * (
263 & gammaT * (insitutLoc - thetaFreeze)
264 & + w_I * (thetaFreeze - insitutLoc) )
265 ENDIF
266
267 ELSE
268
269 C-- In the non-conservative form, only fluxes are turbulent fluxes
270 heatFlux = rUnit2mass*HeatCapacity_Cp *
271 & gammaT * (insitutLoc - thetaFreeze)
272 ENDIF
273
274 C-- Calculate the T and S tendency terms. T tendency term is
275 C-- not necessarily proportional to the heat flux term above because
276 C-- the heat flux term corresponds to total energy change in the grid
277 C-- cell and not the change of energy per unit volume.
278
279 IF (SHELFICEconserve) THEN
280 C-- In the conservative case, meltwater advection contributes (J2001)
281 C-- to T and S tendencies
282 C-- * forcing T (K m/s)
283 forcingT =
284 & (gammaT - w_I)*(thetaFreeze - insitutLoc)
285
286 C-- * forcing S (psu m/s)
287 forcingS =
288 & (gammaS - w_I)*(saltFreeze - sLoc)
289 ELSE
290 C-- Otherwise, the only fluxes out of the ocean that change T and S
291 C-- are the turbulent fluxes.
292 forcingT = gammaT * ( thetaFreeze - insitutLoc )
293 forcingS = gammaS * ( saltFreeze - sLoc )
294 ENDIF
295
296 IF ( debugLevel.GE.debLevE ) THEN
297 WRITE(msgBuf,'(A25,7E16.8)')
298 & 'ZZZ6 aqe, bqe, ceq ',
299 & aqe,bqe,cqe
300
301 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
302 & SQUEEZE_RIGHT , 1)
303
304 WRITE(msgBuf,'(A25,7E16.8)')
305 & 'ZZZ2 T,S,P,t,TFrz,SFrz,w_I ',
306 & tLoc,sLoc,pLoc, insitutLoc, thetaFreeze, saltFreeze, w_I
307
308 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
309 & SQUEEZE_RIGHT , 1)
310 ENDIF
311
312
313 RETURN
314 END
315

  ViewVC Help
Powered by ViewVC 1.1.22