1 |
C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/shelfice/shelfice_thermodynamics.F,v 1.47 2015/12/17 01:52:05 jmc Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "SHELFICE_OPTIONS.h" |
5 |
#ifdef ALLOW_AUTODIFF |
6 |
# include "AUTODIFF_OPTIONS.h" |
7 |
#endif |
8 |
#ifdef ALLOW_CTRL |
9 |
# include "CTRL_OPTIONS.h" |
10 |
#endif |
11 |
|
12 |
CBOP |
13 |
C !ROUTINE: SHELFICE_THERMODYNAMICS |
14 |
C !INTERFACE: |
15 |
SUBROUTINE SHELFICE_THERMODYNAMICS( |
16 |
I myTime, myIter, myThid ) |
17 |
C !DESCRIPTION: \bv |
18 |
C *=============================================================* |
19 |
C | S/R SHELFICE_THERMODYNAMICS |
20 |
C | o shelf-ice main routine. |
21 |
C | compute temperature and (virtual) salt flux at the |
22 |
C | shelf-ice ocean interface |
23 |
C | |
24 |
C | stresses at the ice/water interface are computed in separate |
25 |
C | routines that are called from mom_fluxform/mom_vecinv |
26 |
|
27 |
CIGF | ASSUMES |
28 |
C--- | * SHELFICEconserve = true |
29 |
C *=============================================================* |
30 |
C \ev |
31 |
|
32 |
C !USES: |
33 |
IMPLICIT NONE |
34 |
|
35 |
C === Global variables === |
36 |
#include "SIZE.h" |
37 |
#include "EEPARAMS.h" |
38 |
#include "PARAMS.h" |
39 |
#include "GRID.h" |
40 |
#include "DYNVARS.h" |
41 |
#include "FFIELDS.h" |
42 |
#include "SHELFICE.h" |
43 |
#include "SHELFICE_COST.h" |
44 |
#ifdef ALLOW_AUTODIFF |
45 |
# include "CTRL_SIZE.h" |
46 |
# include "ctrl.h" |
47 |
# include "ctrl_dummy.h" |
48 |
#endif /* ALLOW_AUTODIFF */ |
49 |
#ifdef ALLOW_AUTODIFF_TAMC |
50 |
# ifdef SHI_ALLOW_GAMMAFRICT |
51 |
# include "tamc.h" |
52 |
# include "tamc_keys.h" |
53 |
# endif /* SHI_ALLOW_GAMMAFRICT */ |
54 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
55 |
|
56 |
C !INPUT/OUTPUT PARAMETERS: |
57 |
C === Routine arguments === |
58 |
C myIter :: iteration counter for this thread |
59 |
C myTime :: time counter for this thread |
60 |
C myThid :: thread number for this instance of the routine. |
61 |
_RL myTime |
62 |
INTEGER myIter |
63 |
INTEGER myThid |
64 |
|
65 |
#ifdef ALLOW_SHELFICE |
66 |
C !LOCAL VARIABLES : |
67 |
C === Local variables === |
68 |
C I,J,K,Kp1,bi,bj :: loop counters |
69 |
C tLoc, sLoc, pLoc :: local potential temperature, salinity, pressure |
70 |
C theta/saltFreeze :: temperature and salinity of water at the |
71 |
C ice-ocean interface (at the freezing point) |
72 |
C freshWaterFlux :: local variable for fresh water melt flux due |
73 |
C to melting in kg/m^2/s |
74 |
C (negative density x melt rate) |
75 |
C iceFrontCellThickness :: the ratio of the horizontal length |
76 |
C of the ice front in each model grid cell |
77 |
C divided by the grid cell area. The "thickness" |
78 |
C of the colum perpendicular to the front |
79 |
C iceFrontWidth :: the width of the ice front. |
80 |
|
81 |
INTEGER I,J,K,Kp1 |
82 |
INTEGER bi,bj |
83 |
INTEGER CURI, CURJ, FRONT_K |
84 |
|
85 |
_RL tLoc |
86 |
_RL sLoc |
87 |
_RL pLoc |
88 |
_RL uLoc |
89 |
_RL vLoc |
90 |
_RL wLoc |
91 |
_RL speedLoc |
92 |
|
93 |
C#ifndef SHI_USTAR_WETPOINT |
94 |
C _RL uLoc(1-olx:snx+olx,1-oly:sny+oly) |
95 |
C _RL vLoc(1-olx:snx+olx,1-oly:sny+oly) |
96 |
C#endif |
97 |
C _RL velSq(1-olx:snx+olx,1-oly:sny+oly) |
98 |
|
99 |
_RL freshWaterFlux |
100 |
|
101 |
_RL ice_bottom_Z_C, seafloor_N |
102 |
_RL wet_top_Z_N, wet_bottom_Z_N |
103 |
_RL iceFrontWetContact_Z_max, iceFrontContact_Z_min |
104 |
_RL iceFrontContact_H |
105 |
_RL iceFrontVertContactFrac, iceFrontCellThickness |
106 |
_RL iceFrontWidth, iceFrontFaceArea |
107 |
_RL thermalConductionDistance, thermalConductionTemp |
108 |
_RL tmpHeatFlux, tmpFWFLX |
109 |
_RL tmpForcingT, tmpForcingS |
110 |
_RL tmpFac, icfgridareaFrac |
111 |
INTEGER SI |
112 |
|
113 |
#ifdef ALLOW_DIAGNOSTICS |
114 |
_RL uStarDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
115 |
#endif /* ALLOW_DIAGNOSTICS */ |
116 |
|
117 |
_RL epsilon_H |
118 |
|
119 |
#ifdef ALLOW_SHIFWFLX_CONTROL |
120 |
_RL xx_shifwflx_loc(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) |
121 |
#endif |
122 |
|
123 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
124 |
|
125 |
C-- minimum fraction of a cell adjacent to an ice front that must be |
126 |
C-- wet for exchange to happen |
127 |
epsilon_H = 1. _d -03 |
128 |
|
129 |
C-- hard coded for now. |
130 |
thermalConductionDistance = 100.0 _d 0 |
131 |
thermalConductionTemp = -20.0 _d 0 |
132 |
icfgridareaFrac = 1.0 _d 0 |
133 |
|
134 |
C heat flux into the ice shelf, default is diffusive flux |
135 |
C (Holland and Jenkins, 1999, eq.21) |
136 |
|
137 |
DO bj = myByLo(myThid), myByHi(myThid) |
138 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
139 |
DO J = 1-OLy,sNy+OLy |
140 |
DO I = 1-OLx,sNx+OLx |
141 |
shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0 |
142 |
shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0 |
143 |
SHIICFHeatFlux (I,J,bi,bj) = 0. _d 0 |
144 |
SHIICFFreshWaterFlux(I,J,bi,bj) = 0. _d 0 |
145 |
shelficeForcingT (I,J,bi,bj) = 0. _d 0 |
146 |
shelficeForcingS (I,J,bi,bj) = 0. _d 0 |
147 |
shelficeForcingTR (I,J,bi,bj) = 0. _d 0 |
148 |
|
149 |
DO K = 1, NR |
150 |
iceFrontHeatFlux(I,J,K,bi,bj) = 0. _d 0 |
151 |
iceFrontFreshWaterFlux(I,J,K,bi,bj) = 0. _d 0 |
152 |
iceFrontForcingT(I,J,K,bi,bj) = 0. _d 0 |
153 |
iceFrontForcingS(I,J,K,bi,bj) = 0. _d 0 |
154 |
iceFrontForcingTR(I,J,K,bi,bj) = 0. _d 0 |
155 |
ENDDO /* K */ |
156 |
|
157 |
ENDDO /* I */ |
158 |
ENDDO /* J */ |
159 |
|
160 |
C 3-D velocity-dependent turbulent transfer coefficients |
161 |
|
162 |
DO J = 1-OLy,sNy+OLy |
163 |
DO I = 1-OLx,sNx+OLx |
164 |
DO K = 1, NR |
165 |
|
166 |
C make local copies of velocity |
167 |
uLoc = 0.5*(uVel(I,J,K,bi,bj)+uVel(I+1,J,K,bi,bj)) |
168 |
vLoc = 0.5*(vVel(I,J,K,bi,bj)+vVel(I,J+1,K,bi,bj)) |
169 |
wLoc = 0.5*(wVel(I,J,K,bi,bj)+wVel(I,J,K+1,bi,bj)) |
170 |
speedLoc = SQRT(uLoc*uLoc + vLoc*vLoc + wLoc*wLoc) |
171 |
|
172 |
#ifdef ALLOW_shiTransCoeff_3d |
173 |
shiTransCoeffT(I,J,K,bi,bj) = 1.1 _d -3 * |
174 |
& abs(speedLoc) + 1. _d -4 |
175 |
|
176 |
shiTransCoeffS(I,J,K,bi,bj) = 3.1 _d -5 * |
177 |
& abs(speedLoc) + 5.05 _d -7 |
178 |
#endif |
179 |
ENDDO /* K */ |
180 |
ENDDO /* I */ |
181 |
ENDDO /* J */ |
182 |
|
183 |
C-- First ice front then ice shelf. Loop through each i,j point |
184 |
C-- process ice fronts in k, then process ice shelf. |
185 |
DO J = 1, sNy |
186 |
DO I = 1, sNx |
187 |
|
188 |
C-- The K index where the ice front ends (0 if no ice front) |
189 |
FRONT_K = K_icefront(I,J,bi,bj) |
190 |
|
191 |
C-- If there is an ice front at this (I,J) continue |
192 |
IF (FRONT_K .GT. 0) THEN |
193 |
|
194 |
C-- Loop through all depths where the ice front is fround |
195 |
DO K = 1, FRONT_K |
196 |
C-- Loop around the four laterally neighboring cells of the ice front. |
197 |
C-- If any neighboring points has wet volume in contact with the ice |
198 |
C-- front at (I,J) then calculate ice-ocean exchanges. |
199 |
C-- The four laterally neighboring point are at (CURI,CURJ) |
200 |
DO SI = 1,4 |
201 |
IF (SI .EQ. 1) THEN |
202 |
C-- Looking to right |
203 |
CURI = I+1 |
204 |
CURJ = J |
205 |
|
206 |
iceFrontWidth = dyG(I+1,J,bi,bj) |
207 |
|
208 |
ELSEIF (SI .EQ. 2) THEN |
209 |
C-- Looking to LEFT |
210 |
CURI = I-1 |
211 |
CURJ = J |
212 |
|
213 |
iceFrontWidth = dyG(I,J,bi,bj) |
214 |
ELSEIF (SI .EQ. 3) THEN |
215 |
C-- Looking to NORTH |
216 |
CURI = I |
217 |
CURJ = J+1 |
218 |
|
219 |
iceFrontWidth = dxG(I,J+1,bi,bj) |
220 |
ELSEIF (SI .EQ. 4) THEN |
221 |
C-- Looking to south |
222 |
CURI = I |
223 |
CURJ = J-1 |
224 |
|
225 |
iceFrontWidth = dxG(I,J,bi,bj) |
226 |
endif |
227 |
|
228 |
C-- cell depth describes the average distance |
229 |
C-- perpendicular to the ice front fact |
230 |
|
231 |
iceFrontCellThickness = 0. _d 0 |
232 |
IF(iceFrontWidth.NE.0. _d 0) |
233 |
& iceFrontCellThickness = RA(CURI,CURJ,bi,bj) |
234 |
& /iceFrontWidth |
235 |
iceFrontFaceArea = DRF(K)*iceFrontWidth |
236 |
|
237 |
C-- First, make sure the adjacent point has at least some water in it. |
238 |
IF (_hFacC(CURI,CURJ,K,bi,bj) .GT. zeroRL) THEN |
239 |
|
240 |
C-- we need to determine how much of the ice front is in contact with |
241 |
C-- water in the neighboring grid cell at this depth level. |
242 |
|
243 |
C-- 1. Determine the top depth with water in the current cell |
244 |
C-- 2. Determine the top depth with water in the neighbor cell |
245 |
C-- 3. Determine the depth where water gap between (1) and (2). |
246 |
C-- 4. If there is a gap then ice front is in contact with water in |
247 |
C-- the neighboring cell |
248 |
|
249 |
C-- ice_bottom_Z_C: the depth (m) of the bottom of the ice in the |
250 |
C-- current cell. Bounded between rF(K) and rF(K+1). |
251 |
C-- * If the ice extends past the bottom of the cell then |
252 |
C-- ice_bottom_Z_C = rF(K+1) |
253 |
C-- [rF(k) >= ice_bottom_Z_C >= rF(K+1)] (rF is negative) |
254 |
ice_bottom_Z_C = max(rF(K+1), |
255 |
& min(Ro_surf(I,J, bi,bj), rF(K))) |
256 |
|
257 |
C-- wet_top_Z_N: the depth (m) of the bottom of the ice in the |
258 |
C-- neighboring grid. If the neighboring cell has ice in |
259 |
C-- (in the form of a shelf or front) then wet_top_Z_N is |
260 |
C-- the depth of this neighboring ice. |
261 |
C-- |
262 |
C-- * If neighbor cell has no ice, then Ro_surf = 0 and |
263 |
C-- wet_top_Z_N = rF(K) |
264 |
C-- [rF(k) >= wet_top_Z_N >= rF(K+1)] (rF is negative) |
265 |
|
266 |
wet_top_Z_N = max(rF(K+1), |
267 |
& min(Ro_surf(CURI,CURJ, bi,bj), rF(K))) |
268 |
|
269 |
C-- wet_bottom_Z_N: the depth (m) of the bottom of the wet part of the |
270 |
C-- neighboring cell. If the seafloor reaches into |
271 |
C-- the grid cell then the bottom of the wet part of the |
272 |
C-- grid cell is at the seafloor. |
273 |
C-- |
274 |
C-- * If the seafloor is deeper than this grid cell then |
275 |
C-- wet_bottom_Z = rF(K+1) |
276 |
C-- * If the seafloor is shallower than this grid cell then |
277 |
C-- wet_bottom_Z = rF(K) |
278 |
C-- * If the seafloor reaches partly into this grid cell |
279 |
C-- then wet_bottom_Z = R_low |
280 |
|
281 |
C-- [rF(k) >= wet_bottom_Z >= rF(K+1)] (rF is negative) |
282 |
|
283 |
wet_bottom_Z_N = min(rF(K), |
284 |
& max(R_low(CURI,CURJ, bi,bj), rF(K+1))) |
285 |
|
286 |
C-- iceFrontWetContact_Z_max: The deepest point where the |
287 |
C-- the ice front at (I,J) is in contact with water |
288 |
C-- in the neighboring cell. The shallower of |
289 |
C-- wet_bottom_Z_N (seafloor depth of neighboring point) and |
290 |
C-- ice_bottom_Z_C (bottom of ice front in this center cell). |
291 |
|
292 |
C-- * wet_bottom_Z_N if the seafloor of the neighboring |
293 |
C-- cell is shallower than the ice draft at (I,J). |
294 |
C-- * ice_bottom_Z_C if the ice draft at (I,J) is shallower |
295 |
C-- than the seafloor of the neighboring cell. |
296 |
|
297 |
IF (ice_bottom_Z_C .GT. wet_bottom_Z_N) THEN |
298 |
iceFrontWetContact_Z_max = ice_bottom_Z_C |
299 |
ELSE |
300 |
iceFrontWetContact_Z_max = wet_bottom_Z_N |
301 |
ENDIF |
302 |
|
303 |
C-- The shallowest depth where the ice front at (I,J) is in contact |
304 |
C-- with water in the neighboring cell. If the neighboring cell has |
305 |
C-- no ice draft then wet_top_Z_N = rF(k), the top of the cell. |
306 |
C-- Otherwise, the shallowest depth where the ice front at (I,J) can |
307 |
C-- be in in contact with water (not ice) in (CURI, CURJ) |
308 |
C-- is wet_top_Z_N. |
309 |
|
310 |
C-- the fraction of the grid cell height that has ice draft in contact |
311 |
C-- with water in the neighboring cell. |
312 |
iceFrontVertContactFrac = |
313 |
& (wet_top_Z_N - iceFrontWetContact_Z_max)/ DRF(K) |
314 |
|
315 |
|
316 |
C-- Only proceed if iceFrontVertContactFrac is > 0, the |
317 |
C-- ice draft at (I,J) |
318 |
C-- is in contact with some water in the neighboring grid cell. |
319 |
IF (iceFrontVertContactFrac .GT. epsilon_H) THEN |
320 |
tLoc = theta(CURI,CURJ,K,bi,bj) |
321 |
sLoc = MAX(salt(CURI,CURJ,K,bi,bj), zeroRL) |
322 |
|
323 |
C-- use pressure at the halfway point between the top and bottom of |
324 |
C-- points of the ice front where the ice front is in contact with |
325 |
C-- open water. |
326 |
pLoc = 0.5 _d 0 * ABS(wet_top_Z_N + |
327 |
& iceFrontWetContact_Z_max) |
328 |
|
329 |
CALL SHELFICE_SOLVE4FLUXES( |
330 |
I tLoc, sLoc, pLoc, |
331 |
#ifndef ALLOW_shiTransCoeff_3d |
332 |
I shiTransCoeffT(CURI,CURJ,bi,bj), |
333 |
I shiTransCoeffS(CURI,CURJ,bi,bj), |
334 |
#else |
335 |
I shiTransCoeffT(CURI,CURJ,K,bi,bj), |
336 |
I shiTransCoeffS(CURI,CURJ,K,bi,bj), |
337 |
#endif |
338 |
I thermalConductionDistance, |
339 |
I thermalConductionTemp, |
340 |
O tmpHeatFlux, tmpFWFLX, |
341 |
O tmpForcingT, tmpForcingS, |
342 |
I bi, bj, myTime, myIter, myThid ) |
343 |
|
344 |
C-- fluxes and forcing must be scaled by iceFrontVertContactFract and |
345 |
C-- iceFrontContactFrac some fraction of the heigth and width of the |
346 |
C-- grid cell face may not ice in contact with water. |
347 |
|
348 |
C tmpHeatFlux and tmpFWFLX come as W/m^2 and kg/m^2/s respectively |
349 |
C-- but these rates only apply to the |
350 |
C-- fraction of the grid cell that has ice in contact with seawater. |
351 |
C-- we must scale by iceFrontVertContactFrac to get to the average |
352 |
C-- fluxes in this grid cell. |
353 |
|
354 |
C-- In units W/m^2 |
355 |
iceFrontHeatFlux(CURI,CURJ,K,bi,bj) = |
356 |
& iceFrontHeatFlux(CURI,CURJ,K,bi,bj) + |
357 |
& tmpHeatFlux*iceFrontVertContactFrac |
358 |
|
359 |
C In units of kg/m^2/s |
360 |
iceFrontFreshWaterFlux(CURI,CURJ,K,bi,bj) = |
361 |
& iceFrontFreshWaterFlux(CURI,CURJ,K,bi,bj) + |
362 |
& tmpFWFLX*iceFrontVertContactFrac |
363 |
|
364 |
iceFrontForcingTR(CURI,CURJ,K,bi,bj) = |
365 |
& iceFrontFreshWaterFlux(CURI,CURJ,K,bi,bj) * |
366 |
& mass2rUnit |
367 |
|
368 |
C ow - 06/29/2018 |
369 |
C ow - Verticallly sum up the 3D icefront heat and freshwater fluxes to |
370 |
C ow - compute the total flux for the water column. The shelfice fluxes, |
371 |
C ow - which are 2D, will be added later. NOTE that only |
372 |
C ow - ice-front melts below shelf-ice are included to be consistent |
373 |
C ow - with Rignot's data |
374 |
if(k.GE.kTopC(I,J,bi,bj))then |
375 |
if(RA(CURI,CURJ,bi,bj).NE.0. _d 0)then |
376 |
icfgridareaFrac = |
377 |
& iceFrontFaceArea/RA(CURI,CURJ,bi,bj) |
378 |
SHIICFHeatFlux(CURI,CURJ,bi,bj) = |
379 |
& SHIICFHeatFlux(CURI,CURJ,bi,bj) + |
380 |
& iceFrontHeatFlux(CURI,CURJ,K,bi,bj) |
381 |
& * icfgridareaFrac |
382 |
SHIICFFreshWaterFlux(CURI,CURJ,bi,bj) = |
383 |
& SHIICFFreshWaterFlux(CURI,CURJ,bi,bj) + |
384 |
& iceFrontFreshWaterFlux(CURI,CURJ,K,bi,bj) |
385 |
& * icfgridareaFrac |
386 |
endif |
387 |
endif |
388 |
C iceFrontForcing[T,S] X m/s but these rates only apply to the |
389 |
C-- fraction of the grid cell that has ice in contact with seawater. |
390 |
C-- we must scale by iceFrontVertContactFrac to get to the average |
391 |
C-- fluxes in this grid cell. We must also divide the by the length |
392 |
C-- of the grid cell perpendicular to the face. |
393 |
|
394 |
IF (iceFrontCellThickness .NE. 0. _d 0) THEN |
395 |
C In units of K / s |
396 |
iceFrontForcingT(CURI,CURJ,K,bi,bj) = |
397 |
& iceFrontForcingT(CURI,CURJ,K,bi,bj) + |
398 |
& tmpForcingT/iceFrontCellThickness* |
399 |
& iceFrontVertContactFrac |
400 |
|
401 |
C In units of psu /s |
402 |
iceFrontForcingS(CURI,CURJ,K,bi,bj) = |
403 |
& iceFrontForcingS(CURI,CURJ,K,bi,bj) + |
404 |
& tmpForcingS/iceFrontCellThickness* |
405 |
& iceFrontVertContactFrac |
406 |
|
407 |
ENDIF /* iceFrontCellThickness */ |
408 |
C In units of kg /s |
409 |
addMass(CURI,CURJ,K,bi,bj) = |
410 |
& addMass(CURI,CURJ,K,bi,bj) - |
411 |
& tmpFWFLX*iceFrontFaceArea* |
412 |
& iceFrontVertContactFrac |
413 |
ENDIF /* iceFrontVertContactFrac */ |
414 |
ENDIF /* hFacC(CURI,CURJ,K,bi,bj) */ |
415 |
ENDDO /* SI loop for adjacent cells */ |
416 |
ENDDO /* K LOOP */ |
417 |
ENDIF /* FRONT K */ |
418 |
|
419 |
C-- ice shelf |
420 |
K = kTopC(I,J,bi,bj) |
421 |
|
422 |
C-- If there is an ice front at this (I,J) continue |
423 |
C-- I am assuming K is only .GT. when there is at least some |
424 |
C-- nonzero wet point below the shelf in the grid cell. |
425 |
IF (K .GT. 0) THEN |
426 |
C-- Initialize these values to zero |
427 |
pLoc = 0 _d 0 |
428 |
tLoc = 0 _d 0 |
429 |
sLoc = 0 _d 0 |
430 |
|
431 |
C-- make local copies of temperature, salinity and depth |
432 |
C-- (pressure in deci-bar) underneath the ice |
433 |
C-- for the ice shelf case we use hydrostatic pressure at the ice |
434 |
C-- base of the ice shelf, top of the cavity. |
435 |
|
436 |
pLoc = ABS(R_shelfIce(I,J,bi,bj)) |
437 |
tLoc = theta(I,J,K,bi,bj) |
438 |
sLoc = MAX(salt(I,J,K,bi,bj), zeroRL) |
439 |
|
440 |
CALL SHELFICE_SOLVE4FLUXES( |
441 |
I tLoc, sLoc, pLoc, |
442 |
#ifndef ALLOW_shiTransCoeff_3d |
443 |
I shiTransCoeffT(I,J,bi,bj), |
444 |
I shiTransCoeffS(I,J,bi,bj), |
445 |
#else |
446 |
I shiTransCoeffT(I,J,K,bi,bj), |
447 |
I shiTransCoeffS(I,J,K,bi,bj), |
448 |
#endif |
449 |
I pLoc, thermalConductionTemp, |
450 |
O tmpHeatFlux, tmpFWFLX, |
451 |
O tmpForcingT, tmpForcingS, |
452 |
I bi, bj, myTime, myIter, myThid ) |
453 |
|
454 |
C In units of W/m^2 |
455 |
shelficeHeatFlux(I,J,bi,bj) = tmpHeatFlux |
456 |
C In units of kg/m^2/s |
457 |
shelfIceFreshWaterFlux(I,J,bi,bj) = tmpFWFLX |
458 |
|
459 |
shelficeForcingTR(I,J,bi,bj) = |
460 |
& shelfIceFreshWaterFlux(I,J,bi,bj) * mass2rUnit |
461 |
|
462 |
C ow - 06/29/2018 |
463 |
C ow - Now add shelfice heat and freshwater fluxes |
464 |
SHIICFHeatFlux(i,j,bi,bj) = |
465 |
& SHIICFHeatFlux(i,j,bi,bj) + |
466 |
& shelficeHeatFlux(i,j,bi,bj) |
467 |
SHIICFFreshWaterFlux(i,j,bi,bj) = |
468 |
& SHIICFFreshWaterFlux(i,j,bi,bj) + |
469 |
& shelfIceFreshWaterFlux(i,j,bi,bj) |
470 |
C In units of K/s -- division by drF required first |
471 |
shelficeForcingT(I,J,bi,bj) = tmpForcingT* |
472 |
& recip_drF(K)* _recip_hFacC(i,j,K,bi,bj) |
473 |
C In units of psu/s -- division by drF required first |
474 |
shelficeForcingS(I,J,bi,bj) = tmpForcingS* |
475 |
& recip_drF(K)* _recip_hFacC(i,j,K,bi,bj) |
476 |
C In units of kg/s -- multiplication of area required first |
477 |
addMass(I,J,K, bi,bj) = addMass(I,J,K, bi,bj) - |
478 |
& tmpFWFLX*RA(I,J,bi,bj) |
479 |
ENDIF /* SHELF K > 0 */ |
480 |
ENDDO /* i */ |
481 |
ENDDO /* j */ |
482 |
ENDDO /* bi */ |
483 |
ENDDO /* bj */ |
484 |
|
485 |
|
486 |
C-- Calculate new loading anomaly (in case the ice-shelf mass was updated) |
487 |
#ifndef ALLOW_AUTODIFF |
488 |
c IF ( SHELFICEloadAnomalyFile .EQ. ' ' ) THEN |
489 |
DO bj = myByLo(myThid), myByHi(myThid) |
490 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
491 |
DO j = 1-OLy, sNy+OLy |
492 |
DO i = 1-OLx, sNx+OLx |
493 |
shelficeLoadAnomaly(i,j,bi,bj) = gravity |
494 |
& *( shelficeMass(i,j,bi,bj) + rhoConst*Ro_surf(i,j,bi,bj) ) |
495 |
ENDDO |
496 |
ENDDO |
497 |
ENDDO |
498 |
ENDDO |
499 |
c ENDIF |
500 |
#endif /* ndef ALLOW_AUTODIFF */ |
501 |
|
502 |
#ifdef ALLOW_DIAGNOSTICS |
503 |
IF ( useDiagnostics ) THEN |
504 |
CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx', |
505 |
& 0,1,0,1,1,myThid) |
506 |
CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux, 'SHIhtFlx', |
507 |
& 0,1,0,1,1,myThid) |
508 |
|
509 |
CALL DIAGNOSTICS_FILL_RS(SHIICFFreshWaterFlux,'SHIICFfwFlx', |
510 |
& 0,1,0,1,1,myThid) |
511 |
CALL DIAGNOSTICS_FILL_RS(SHIICFHeatFlux, 'SHIICFhtFlx', |
512 |
& 0,1,0,1,1,myThid) |
513 |
|
514 |
CALL DIAGNOSTICS_FILL(iceFrontFreshWaterFlux, 'ICFfwFlx', |
515 |
& 0,Nr,0,1,1,myThid) |
516 |
CALL DIAGNOSTICS_FILL(iceFrontHeatFlux, 'ICFhtFlx', |
517 |
& 0,Nr,0,1,1,myThid) |
518 |
|
519 |
CALL DIAGNOSTICS_FILL(shelfIceForcingTR, 'SHITR ', |
520 |
& 0,Nr,0,1,1,myThid) |
521 |
|
522 |
CALL DIAGNOSTICS_FILL(iceFrontForcingTR, 'ICFTR ', |
523 |
& 0,Nr,0,1,1,myThid) |
524 |
|
525 |
C SHIForcT (Ice shelf forcing for theta [W/m2], >0 increases theta) |
526 |
tmpFac = HeatCapacity_Cp*rUnit2mass |
527 |
CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingT,tmpFac,1, |
528 |
& 'SHIForcT',0,1,0,1,1,myThid) |
529 |
C SHIForcS (Ice shelf forcing for salt [g/m2/s], >0 increases salt) |
530 |
tmpFac = rUnit2mass |
531 |
CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingS,tmpFac,1, |
532 |
& 'SHIForcS',0,1,0,1,1,myThid) |
533 |
|
534 |
C ICFForcT (Ice front forcing for theta [W/m2], >0 increases theta) |
535 |
tmpFac = HeatCapacity_Cp*rUnit2mass |
536 |
CALL DIAGNOSTICS_SCALE_FILL(iceFrontForcingT,tmpFac,1, |
537 |
& 'ICFForcT',0,Nr,0,1,1,myThid) |
538 |
C ICFForcS (Ice front forcing for salt [g/m2/s], >0 increases salt) |
539 |
tmpFac = rUnit2mass |
540 |
CALL DIAGNOSTICS_SCALE_FILL(iceFrontForcingS,tmpFac,1, |
541 |
& 'ICFForcS',0,Nr,0,1,1,myThid) |
542 |
|
543 |
C Transfer coefficients |
544 |
#ifndef ALLOW_shiTransCoeff_3d |
545 |
CALL DIAGNOSTICS_FILL(shiTransCoeffT,'SHIgammT', |
546 |
& 0,1,0,1,1,myThid) |
547 |
CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS', |
548 |
& 0,1,0,1,1,myThid) |
549 |
#else |
550 |
CALL DIAGNOSTICS_FILL(shiTransCoeffT,'SHIgammT', |
551 |
& 0,Nr,0,1,1,myThid) |
552 |
CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS', |
553 |
& 0,Nr,0,1,1,myThid) |
554 |
#endif |
555 |
C Friction velocity |
556 |
#ifdef SHI_ALLOW_GAMMAFRICT |
557 |
IF ( SHELFICEuseGammaFrict ) |
558 |
& CALL DIAGNOSTICS_FILL(uStarDiag,'SHIuStar',0,1,0,1,1,myThid) |
559 |
#endif /* SHI_ALLOW_GAMMAFRICT */ |
560 |
ENDIF |
561 |
#endif |
562 |
|
563 |
#endif /* ALLOW_SHELFICE */ |
564 |
RETURN |
565 |
END |