/[MITgcm]/MITgcm_contrib/shelfice_remeshing/AUTO/code/shelfice_thermodynamics.F
ViewVC logotype

Contents of /MITgcm_contrib/shelfice_remeshing/AUTO/code/shelfice_thermodynamics.F

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

Revision 1.3 - (show annotations) (download)
Thu Oct 22 12:26:45 2015 UTC (10 years ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +5 -1 lines 
*** empty log message ***
1 C $Header: /u/gcmpack/MITgcm_contrib/shelfice_remeshing/AUTO/code/shelfice_thermodynamics.F,v 1.2 2015/10/06 10:38:53 dgoldberg 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 C *=============================================================*
27 C \ev
28
29 C !USES:
30 IMPLICIT NONE
31
32 C === Global variables ===
33 #include "SIZE.h"
34 #include "EEPARAMS.h"
35 #include "PARAMS.h"
36 #include "GRID.h"
37 #include "DYNVARS.h"
38 #include "FFIELDS.h"
39 #include "SHELFICE.h"
40 #include "SHELFICE_COST.h"
41 #ifdef ALLOW_AUTODIFF
42 # include "CTRL_SIZE.h"
43 # include "ctrl.h"
44 # include "ctrl_dummy.h"
45 #endif /* ALLOW_AUTODIFF */
46 #ifdef ALLOW_AUTODIFF_TAMC
47 # ifdef SHI_ALLOW_GAMMAFRICT
48 # include "tamc.h"
49 # include "tamc_keys.h"
50 # endif /* SHI_ALLOW_GAMMAFRICT */
51 #endif /* ALLOW_AUTODIFF_TAMC */
52
53 C !INPUT/OUTPUT PARAMETERS:
54 C === Routine arguments ===
55 C myIter :: iteration counter for this thread
56 C myTime :: time counter for this thread
57 C myThid :: thread number for this instance of the routine.
58 _RL myTime
59 INTEGER myIter
60 INTEGER myThid
61
62 #ifdef ALLOW_SHELFICE
63 C !LOCAL VARIABLES :
64 C === Local variables ===
65 C I,J,K,Kp1,bi,bj :: loop counters
66 C tLoc, sLoc, pLoc :: local in-situ temperature, salinity, pressure
67 C theta/saltFreeze :: temperature and salinity of water at the
68 C ice-ocean interface (at the freezing point)
69 C freshWaterFlux :: local variable for fresh water melt flux due
70 C to melting in kg/m^2/s
71 C (negative density x melt rate)
72 C convertFW2SaltLoc:: local copy of convertFW2Salt
73 C cFac :: 1 for conservative form, 0, otherwise
74 C rFac :: realFreshWaterFlux factor
75 C dFac :: 0 for diffusive heat flux (Holland and Jenkins, 1999,
76 C eq21)
77 C 1 for advective and diffusive heat flux (eq22, 26, 31)
78 C fwflxFac :: only effective for dFac=1, 1 if we expect a melting
79 C fresh water flux, 0 otherwise
80 C auxiliary variables and abbreviations:
81 C a0, a1, a2, b, c0
82 C eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8
83 C aqe, bqe, cqe, discrim, recip_aqe
84 C drKp1, recip_drLoc
85 INTEGER I,J,K,Kp1,kp2
86 INTEGER bi,bj
87 _RL tLoc(1:sNx,1:sNy)
88 _RL sLoc(1:sNx,1:sNy)
89 _RL pLoc(1:sNx,1:sNy)
90 _RL uLoc(1:sNx,1:sNy)
91 _RL vLoc(1:sNx,1:sNy)
92 _RL u_topdr(1:sNx+1,1:sNy+1,nSx,nSy)
93 _RL v_topdr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
94 _RL thetaFreeze, saltFreeze, recip_Cp
95 _RL freshWaterFlux, convertFW2SaltLoc
96 _RL a0, a1, a2, b, c0
97 _RL eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8
98 _RL cFac, rFac, dFac, fwflxFac, realfwFac
99 _RL aqe, bqe, cqe, discrim, recip_aqe
100 _RL drKp1, drKp2, recip_drLoc
101 _RL recip_latentHeat
102 _RL tmpFac
103
104 #ifdef SHI_ALLOW_GAMMAFRICT
105 _RL shiPr, shiSc, shiLo, recip_shiKarman, shiTwoThirds
106 _RL gammaTmoleT, gammaTmoleS, gammaTurb, gammaTurbConst
107 _RL ustar, ustarSq, etastar
108 PARAMETER ( shiTwoThirds = 0.66666666666666666666666666667D0 )
109 #ifdef ALLOW_DIAGNOSTICS
110 _RL uStarDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
111 #endif /* ALLOW_DIAGNOSTICS */
112 #endif
113
114 #ifndef ALLOW_OPENAD
115 _RL SW_TEMP
116 EXTERNAL SW_TEMP
117 #endif
118
119 #ifdef ALLOW_SHIFWFLX_CONTROL
120 _RL xx_shifwflx_loc(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
121 #endif
122 CEOP
123 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
124
125 #ifdef SHI_ALLOW_GAMMAFRICT
126 #ifdef ALLOW_AUTODIFF
127 C re-initialize here again, curtesy to TAF
128 DO bj = myByLo(myThid), myByHi(myThid)
129 DO bi = myBxLo(myThid), myBxHi(myThid)
130 DO J = 1-OLy,sNy+OLy
131 DO I = 1-OLx,sNx+OLx
132 shiTransCoeffT(i,j,bi,bj) = SHELFICEheatTransCoeff
133 shiTransCoeffS(i,j,bi,bj) = SHELFICEsaltTransCoeff
134 ENDDO
135 ENDDO
136 ENDDO
137 ENDDO
138 #endif /* ALLOW_AUTODIFF */
139 IF ( SHELFICEuseGammaFrict ) THEN
140 C Implement friction velocity-dependent transfer coefficient
141 C of Holland and Jenkins, JPO, 1999
142 recip_shiKarman= 1. _d 0 / 0.4 _d 0
143 shiLo = 0. _d 0
144 shiPr = shiPrandtl**shiTwoThirds
145 shiSc = shiSchmidt**shiTwoThirds
146 cph shiPr = (viscArNr(1)/diffKrNrT(1))**shiTwoThirds
147 cph shiSc = (viscArNr(1)/diffKrNrS(1))**shiTwoThirds
148 gammaTmoleT = 12.5 _d 0 * shiPr - 6. _d 0
149 gammaTmoleS = 12.5 _d 0 * shiSc - 6. _d 0
150 C instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc))
151 etastar = 1. _d 0
152 gammaTurbConst = 1. _d 0 / (2. _d 0 * shiZetaN*etastar)
153 & - recip_shiKarman
154 #ifdef ALLOW_AUTODIFF
155 DO bj = myByLo(myThid), myByHi(myThid)
156 DO bi = myBxLo(myThid), myBxHi(myThid)
157 DO J = 1-OLy,sNy+OLy
158 DO I = 1-OLx,sNx+OLx
159 shiTransCoeffT(i,j,bi,bj) = 0. _d 0
160 shiTransCoeffS(i,j,bi,bj) = 0. _d 0
161 ENDDO
162 ENDDO
163 ENDDO
164 ENDDO
165 #endif /* ALLOW_AUTODIFF */
166 ENDIF
167 #endif /* SHI_ALLOW_GAMMAFRICT */
168
169 recip_latentHeat = 0. _d 0
170 IF ( SHELFICElatentHeat .NE. 0. _d 0 )
171 & recip_latentHeat = 1. _d 0/SHELFICElatentHeat
172 C are we doing the conservative form of Jenkins et al. (2001)?
173 recip_Cp = 1. _d 0 / HeatCapacity_Cp
174 cFac = 0. _d 0
175 IF ( SHELFICEconserve ) cFac = 1. _d 0
176
177 realFWfac = 0. _d 0
178 IF ( SHELFICErealFWflux ) realFWfac = 1. _d 0
179 C with "real fresh water flux" (affecting ETAN),
180 C there is more to modify
181 rFac = 1. _d 0
182 IF ( SHELFICEconserve .AND. useRealFreshWaterFlux ) rFac = 0. _d 0
183 C heat flux into the ice shelf, default is diffusive flux
184 C (Holland and Jenkins, 1999, eq.21)
185 dFac = 0. _d 0
186 IF ( SHELFICEadvDiffHeatFlux ) dFac = 1. _d 0
187 fwflxFac = 0. _d 0
188 C linear dependence of freezing point on salinity
189 a0 = -0.0575 _d 0
190 a1 = 0.0 _d -0
191 a2 = 0.0 _d -0
192 c0 = 0.0901 _d 0
193 b = -7.61 _d -4
194 #ifdef ALLOW_ISOMIP_TD
195 IF ( useISOMIPTD ) THEN
196 C non-linear dependence of freezing point on salinity
197 a0 = -0.0575 _d 0
198 a1 = 1.710523 _d -3
199 a2 = -2.154996 _d -4
200 b = -7.53 _d -4
201 c0 = 0. _d 0
202 ENDIF
203 convertFW2SaltLoc = convertFW2Salt
204 C hardcoding this value here is OK because it only applies to ISOMIP
205 C where this value is part of the protocol
206 IF ( convertFW2SaltLoc .EQ. -1. ) convertFW2SaltLoc = 33.4 _d 0
207 #endif /* ALLOW_ISOMIP_TD */
208
209 DO bj = myByLo(myThid), myByHi(myThid)
210 DO bi = myBxLo(myThid), myBxHi(myThid)
211 DO J = 1-OLy,sNy+OLy
212 DO I = 1-OLx,sNx+OLx
213 shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
214 shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
215 shelficeForcingT (I,J,bi,bj) = 0. _d 0
216 shelficeForcingS (I,J,bi,bj) = 0. _d 0
217 #if (defined SHI_ALLOW_GAMMAFRICT && defined ALLOW_DIAGNOSTICS)
218 uStarDiag (I,J,bi,bj) = 0. _d 0
219 #endif /* SHI_ALLOW_GAMMAFRICT and ALLOW_DIAGNOSTICS */
220 ENDDO
221 ENDDO
222 ENDDO
223 ENDDO
224 #ifdef ALLOW_SHIFWFLX_CONTROL
225 DO bj = myByLo(myThid), myByHi(myThid)
226 DO bi = myBxLo(myThid), myBxHi(myThid)
227 DO J = 1-OLy,sNy+OLy
228 DO I = 1-OLx,sNx+OLx
229 xx_shifwflx_loc(I,J,bi,bj) = 0. _d 0
230 ENDDO
231 ENDDO
232 ENDDO
233 ENDDO
234 #ifdef ALLOW_CTRL
235 if (useCTRL) CALL CTRL_GET_GEN (
236 & xx_shifwflx_file, xx_shifwflxstartdate, xx_shifwflxperiod,
237 & maskSHI, xx_shifwflx_loc, xx_shifwflx0, xx_shifwflx1,
238 & xx_shifwflx_dummy,
239 & xx_shifwflx_remo_intercept, xx_shifwflx_remo_slope,
240 & wshifwflx,
241 & myTime, myIter, myThid )
242 #endif
243 #endif /* ALLOW_SHIFWFLX_CONTROL */
244 DO bj = myByLo(myThid), myByHi(myThid)
245 DO bi = myBxLo(myThid), myBxHi(myThid)
246
247 IF ( SHELFICEBoundaryLayer ) THEN
248 C-- average over boundary layer width
249 DO J = 1, sNy+1
250 DO I = 1, sNx+1
251 u_topdr(I,J,bi,bj) = 0.0
252 v_topdr(I,J,bi,bj) = 0.0
253 ENDDO
254 ENDDO
255 ENDIF
256
257 #ifdef ALLOW_AUTODIFF_TAMC
258 # ifdef SHI_ALLOW_GAMMAFRICT
259 act1 = bi - myBxLo(myThid)
260 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
261 act2 = bj - myByLo(myThid)
262 max2 = myByHi(myThid) - myByLo(myThid) + 1
263 act3 = myThid - 1
264 max3 = nTx*nTy
265 act4 = ikey_dynamics - 1
266 ikey = (act1 + 1) + act2*max1
267 & + act3*max1*max2
268 & + act4*max1*max2*max3
269 # endif /* SHI_ALLOW_GAMMAFRICT */
270 #endif /* ALLOW_AUTODIFF_TAMC */
271 DO J = 1, sNy
272 DO I = 1, sNx
273 C-- make local copies of temperature, salinity and depth (pressure in deci-bar)
274 C-- underneath the ice
275 K = MAX(1,kTopC(I,J,bi,bj))
276 pLoc(I,J) = ABS(R_shelfIce(I,J,bi,bj))
277 c pLoc(I,J) = shelficeMass(I,J,bi,bj)*gravity*1. _d -4
278 tLoc(I,J) = theta(I,J,K,bi,bj)
279 sLoc(I,J) = MAX(salt(I,J,K,bi,bj), zeroRL)
280 IF ( .not.SHELFICEBoundaryLayer ) THEN
281 uLoc(I,J) = recip_hFacC(I,J,K,bi,bj) *
282 & ( uVel(I, J,K,bi,bj) * _hFacW(I, J,K,bi,bj)
283 & + uVel(I+1,J,K,bi,bj) * _hFacW(I+1,J,K,bi,bj) )
284 vLoc(I,J) = recip_hFacC(I,J,K,bi,bj) *
285 & ( vVel(I, J,K,bi,bj) * _hFacS(I, J,K,bi,bj)
286 & + vVel(I,J+1,K,bi,bj) * _hFacS(I,J+1,K,bi,bj) )
287 ENDIF
288 ENDDO
289 ENDDO
290
291 ! IF ( SHELFICEBoundaryLayer ) THEN
292 ! DO J = 1, sNy+1
293 ! DO I = 1, sNx+1
294 !
295 ! K = ksurfW(I,J,bi,bj)
296 ! Kp1 = K+1
297 ! Kp2 = K+2
298 !
299 ! IF (ShelficeThickBoundaryLayer .and.
300 ! & (K.ne.0.and.K.LT.Nr-1)) THEN
301 !
302 ! drKp1 = drF(K)*( 1.5 - _hFacW(I,J,K,bi,bj) )
303 ! drKp2 = drKp1 - drF(kp1)*_hFacW(I,J,kp1,bi,bj)
304 ! drKp2 = MAX( drKp2, 0. _d 0)
305 ! drKp2 = MIN( drKp2,
306 ! & drF(kp2)*_hFacW(I,J,kp2,bi,bj))
307 ! drKp1 = drKp1 - drKp2
308 ! drKp1 = MAX( drKp1, 0. _d 0)
309 ! recip_drLoc = 1. _d 0 /
310 ! & (drF(K)*_hFacW(I,J,K,bi,bj)+drKp1+drKp2)
311 ! u_topdr(I,J,bi,bj) =
312 ! & (drF(K)*_hFacW(I,J,K,bi,bj)*uVel(I,J,K,bi,bj) +
313 ! & drKp1*uVel(I,J,Kp1,bi,bj)) * recip_drLoc
314 ! u_topdr(I,J,bi,bj) = u_topdr(I,J,bi,bj) +
315 ! & drKp2 * uVel(I,J,Kp2,bi,bj) * recip_drLoc
316 !
317 ! ELSEIF ( (K .NE. 0 .AND. K.EQ.Nr-1) .OR.
318 ! & (.not.SHELFICEthickboundarylayer.AND.
319 ! & (K .NE. 0 .AND. K .LT. Nr) ) ) THEN
320 !
321 ! drKp1 = drF(K)*(1. _d 0-_hFacW(I,J,K,bi,bj))
322 ! drKp1 = max (drKp1, 0. _d 0)
323 ! recip_drLoc = 1.0 /
324 ! & (drF(K)*_hFacW(I,J,K,bi,bj)+drKp1)
325 ! u_topdr(I,J,bi,bj) =
326 ! & (drF(K)*_hFacW(I,J,K,bi,bj)*uVel(I,J,K,bi,bj) +
327 ! & drKp1*uVel(I,J,Kp1,bi,bj))
328 ! & * recip_drLoc
329 !
330 ! ELSE
331 !
332 ! u_topdr(I,J,bi,bj) = 0. _d 0
333 !
334 ! ENDIF
335 !
336 ! K = ksurfS(I,J,bi,bj)
337 ! Kp1 = K+1
338 ! Kp2 = K+2
339 !
340 ! IF (ShelficeThickBoundaryLayer .and.
341 ! & (K.ne.0.and.K.LT.Nr-1)) THEN
342 !
343 ! drKp1 = drF(K)*( 1.5 - _hFacS(I,J,K,bi,bj) )
344 ! drKp2 = drKp1 - drF(kp1)*_hFacS(I,J,kp1,bi,bj)
345 ! drKp2 = MAX( drKp2, 0. _d 0)
346 ! drKp2 = MIN( drKp2,
347 ! & drF(kp2)*_hFacS(I,J,kp2,bi,bj))
348 ! drKp1 = drKp1 - drKp2
349 ! drKp1 = MAX( drKp1, 0. _d 0)
350 ! recip_drLoc = 1. _d 0 /
351 ! & (drF(K)*_hFacS(I,J,K,bi,bj)+drKp1+drKp2)
352 ! v_topdr(I,J,bi,bj) =
353 ! & (drF(K)*_hFacS(I,J,K,bi,bj)*vVel(I,J,K,bi,bj) +
354 ! & drKp1*vVel(I,J,Kp1,bi,bj)) * recip_drLoc
355 ! v_topdr(I,J,bi,bj) = v_topdr(I,J,bi,bj) +
356 ! & drKp2 * vVel(I,J,Kp2,bi,bj) * recip_drLoc
357 !
358 ! ELSEIF ( (K .NE. 0 .AND. K.EQ.Nr-1) .OR.
359 ! & ((.NOT.SHELFICEthickboundarylayer).AND.
360 ! & (K .NE. 0 .AND. K .LT. Nr) ) ) THEN
361 !
362 ! drKp1 = drF(K)*(1. _d 0-_hFacS(I,J,K,bi,bj))
363 ! drKp1 = max (drKp1, 0. _d 0)
364 ! recip_drLoc = 1.0 /
365 ! & (drF(K)*_hFacS(I,J,K,bi,bj)+drKp1)
366 ! v_topdr(I,J,bi,bj) =
367 ! & (drF(K)*_hFacS(I,J,K,bi,bj)*vVel(I,J,K,bi,bj) +
368 ! & drKp1*vVel(I,J,Kp1,bi,bj))
369 ! & * recip_drLoc
370 !
371 ! ELSE
372 !
373 ! v_topdr(I,J,bi,bj) = 0. _d 0
374 !
375 ! ENDIF
376 !
377 ! ENDDO
378 ! ENDDO
379 ! ENDIF
380
381 IF ( SHELFICEBoundaryLayer ) THEN
382 DO J = 1, sNy+1
383 DO I = 1, sNx+1
384 K = ksurfW(I,J,bi,bj)
385 Kp1 = K+1
386 IF (K.lt.Nr) then
387 drKp1 = drF(K)*(1. _d 0-_hFacW(I,J,K,bi,bj))
388 drKp1 = max (drKp1, 0. _d 0)
389 recip_drLoc = 1.0 /
390 & (drF(K)*_hFacW(I,J,K,bi,bj)+drKp1)
391 u_topdr(I,J,bi,bj) =
392 & (drF(K)*_hFacW(I,J,K,bi,bj)*uVel(I,J,K,bi,bj) +
393 & drKp1*uVel(I,J,Kp1,bi,bj))
394 & * recip_drLoc
395 ELSE
396 u_topdr(I,J,bi,bj) = 0. _d 0
397 ENDIF
398
399 K = ksurfS(I,J,bi,bj)
400 Kp1 = K+1
401 IF (K.lt.Nr) then
402 drKp1 = drF(K)*(1. _d 0-_hFacS(I,J,K,bi,bj))
403 drKp1 = max (drKp1, 0. _d 0)
404 recip_drLoc = 1.0 /
405 & (drF(K)*_hFacS(I,J,K,bi,bj)+drKp1)
406 v_topdr(I,J,bi,bj) =
407 & (drF(K)*_hFacS(I,J,K,bi,bj)*vVel(I,J,K,bi,bj) +
408 & drKp1*vVel(I,J,Kp1,bi,bj))
409 & * recip_drLoc
410 ELSE
411 v_topdr(I,J,bi,bj) = 0. _d 0
412 ENDIF
413
414 ENDDO
415 ENDDO
416 ENDIF
417
418 IF ( SHELFICEBoundaryLayer ) THEN
419 C-- average over boundary layer width
420 DO J = 1, sNy
421 DO I = 1, sNx
422 K = kTopC(I,J,bi,bj)
423 IF ( K .NE. 0 .AND. K .LT. Nr ) THEN
424 Kp1 = MIN(Nr,K+1)
425 C-- overlap into lower cell
426 drKp1 = drF(K)*( 1. _d 0 - _hFacC(I,J,K,bi,bj) )
427 C-- Dans fix
428 drKp1 = MAX(drKp1, 0.)
429 C-- lower cell may not be as thick as required
430 drKp1 = MIN( drKp1, drF(Kp1) * _hFacC(I,J,Kp1,bi,bj) )
431 recip_drLoc = 1. _d 0 /
432 & ( drF(K)*_hFacC(I,J,K,bi,bj) + drKp1 )
433 tLoc(I,J) = ( tLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
434 & + theta(I,J,Kp1,bi,bj) *drKp1 )
435 & * recip_drLoc
436 sLoc(I,J) = ( sLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
437 & + MAX(salt(I,J,Kp1,bi,bj), zeroRL) * drKp1 )
438 & * recip_drLoc
439
440 ! uLoc(I,J) = ( uLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
441 ! & + drKp1 * recip_hFacC(I,J,Kp1,bi,bj) *
442 ! & ( uVel(I, J,Kp1,bi,bj) * _hFacW(I, J,Kp1,bi,bj)
443 ! & + uVel(I+1,J,Kp1,bi,bj) * _hFacW(I+1,J,Kp1,bi,bj) )
444 ! & ) * recip_drLoc
445 ! vLoc(I,J) = ( vLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
446 ! & + drKp1 * recip_hFacC(I,J,Kp1,bi,bj) *
447 ! & ( vVel(I,J, Kp1,bi,bj) * _hFacS(I,J, Kp1,bi,bj)
448 ! & + vVel(I,J+1,Kp1,bi,bj) * _hFacS(I,J+1,Kp1,bi,bj) )
449 ! & ) * recip_drLoc
450 ENDIF
451 ENDDO
452 ENDDO
453 ENDIF
454
455
456 IF ( SHELFICEBoundaryLayer ) THEN
457 DO J = 1, sNy
458 DO I = 1, sNx
459 uLoc(I,J) =
460 & u_topdr(I,J,bi,bj) + u_topdr(I+1,J,bi,bj)
461 vLoc(I,J) =
462 & v_topdr(I,J,bi,bj) + v_topdr(I,J+1,bi,bj)
463 ENDDO
464 ENDDO
465 ENDIF
466
467 C-- turn potential temperature into in-situ temperature relative
468 C-- to the surface
469 DO J = 1, sNy
470 DO I = 1, sNx
471 #ifndef ALLOW_OPENAD
472 tLoc(I,J) = SW_TEMP(sLoc(I,J),tLoc(I,J),pLoc(I,J),zeroRL)
473 #else
474 CALL SW_TEMP(sLoc(I,J),tLoc(I,J),pLoc(I,J),zeroRL,tLoc(I,J))
475 #endif
476 ENDDO
477 ENDDO
478
479 #ifdef SHI_ALLOW_GAMMAFRICT
480 IF ( SHELFICEuseGammaFrict ) THEN
481 DO J = 1, sNy
482 DO I = 1, sNx
483 K = kTopC(I,J,bi,bj)
484 IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN
485 ustarSq = shiCdrag * MAX( 1.D-6,
486 & 0.25 _d 0 *(uLoc(I,J)*uLoc(I,J)+vLoc(I,J)*vLoc(I,J)) )
487 ustar = SQRT(ustarSq)
488 #ifdef ALLOW_DIAGNOSTICS
489 uStarDiag(I,J,bi,bj) = ustar
490 #endif /* ALLOW_DIAGNOSTICS */
491 C instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc))
492 C etastar = 1. _d 0
493 C gammaTurbConst = 1. _d 0 / (2. _d 0 * shiZetaN*etastar)
494 C & - recip_shiKarman
495 IF ( fCori(I,J,bi,bj) .NE. 0. _d 0 ) THEN
496 gammaTurb = LOG( ustarSq * shiZetaN * etastar**2
497 & / ABS(fCori(I,J,bi,bj) * 5.0 _d 0 * shiKinVisc))
498 & * recip_shiKarman
499 & + gammaTurbConst
500 C Do we need to catch the unlikely case of very small ustar
501 C that can lead to negative gammaTurb?
502 C gammaTurb = MAX(0.D0, gammaTurb)
503 ELSE
504 gammaTurb = gammaTurbConst
505 ENDIF
506 shiTransCoeffT(i,j,bi,bj) = MAX( zeroRL,
507 & ustar/(gammaTurb + gammaTmoleT) )
508 shiTransCoeffS(i,j,bi,bj) = MAX( zeroRL,
509 & ustar/(gammaTurb + gammaTmoleS) )
510 ENDIF
511 ENDDO
512 ENDDO
513 ENDIF
514 #endif /* SHI_ALLOW_GAMMAFRICT */
515
516 #ifdef ALLOW_AUTODIFF_TAMC
517 # ifdef SHI_ALLOW_GAMMAFRICT
518 CADJ STORE shiTransCoeffS(:,:,bi,bj) = comlev1_bibj,
519 CADJ & key=ikey, byte=isbyte
520 CADJ STORE shiTransCoeffT(:,:,bi,bj) = comlev1_bibj,
521 CADJ & key=ikey, byte=isbyte
522 # endif /* SHI_ALLOW_GAMMAFRICT */
523 #endif /* ALLOW_AUTODIFF_TAMC */
524 #ifdef ALLOW_ISOMIP_TD
525 IF ( useISOMIPTD ) THEN
526 DO J = 1, sNy
527 DO I = 1, sNx
528 K = kTopC(I,J,bi,bj)
529 IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN
530 C-- Calculate freezing temperature as a function of salinity and pressure
531 thetaFreeze =
532 & sLoc(I,J) * ( a0 + a1*sqrt(sLoc(I,J)) + a2*sLoc(I,J) )
533 & + b*pLoc(I,J) + c0
534 C-- Calculate the upward heat and fresh water fluxes
535 shelfIceHeatFlux(I,J,bi,bj) = maskC(I,J,K,bi,bj)
536 & * shiTransCoeffT(i,j,bi,bj)
537 & * ( tLoc(I,J) - thetaFreeze )
538 & * HeatCapacity_Cp*rUnit2mass
539 #ifdef ALLOW_SHIFWFLX_CONTROL
540 & - xx_shifwflx_loc(I,J,bi,bj)*SHELFICElatentHeat
541 #endif /* ALLOW_SHIFWFLX_CONTROL */
542 C upward heat flux into the shelf-ice implies basal melting,
543 C thus a downward (negative upward) fresh water flux (as a mass flux),
544 C and vice versa
545 shelfIceFreshWaterFlux(I,J,bi,bj) =
546 & - shelfIceHeatFlux(I,J,bi,bj)
547 & *recip_latentHeat
548 C-- compute surface tendencies
549 shelficeForcingT(i,j,bi,bj) =
550 & - shelfIceHeatFlux(I,J,bi,bj)
551 & *recip_Cp*mass2rUnit
552 & - cFac * shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit
553 & * ( thetaFreeze - tLoc(I,J) )
554 shelficeForcingS(i,j,bi,bj) =
555 & shelfIceFreshWaterFlux(I,J,bi,bj) * mass2rUnit
556 & * ( cFac*sLoc(I,J) + (1. _d 0-cFac)*convertFW2SaltLoc )
557 C-- stress at the ice/water interface is computed in separate
558 C routines that are called from mom_fluxform/mom_vecinv
559 ELSE
560 shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
561 shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
562 shelficeForcingT (I,J,bi,bj) = 0. _d 0
563 shelficeForcingS (I,J,bi,bj) = 0. _d 0
564 ENDIF
565 ENDDO
566 ENDDO
567 ELSE
568 #else
569 IF ( .TRUE. ) THEN
570 #endif /* ALLOW_ISOMIP_TD */
571 C use BRIOS thermodynamics, following Hellmers PhD thesis:
572 C Hellmer, H., 1989, A two-dimensional model for the thermohaline
573 C circulation under an ice shelf, Reports on Polar Research, No. 60
574 C (in German).
575
576 DO J = 1, sNy
577 DO I = 1, sNx
578 K = kTopC(I,J,bi,bj)
579 IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN
580 C heat flux into the ice shelf, default is diffusive flux
581 C (Holland and Jenkins, 1999, eq.21)
582 thetaFreeze = a0*sLoc(I,J)+c0+b*pLoc(I,J)
583 fwflxFac = 0. _d 0
584 IF ( tLoc(I,J) .GT. thetaFreeze ) fwflxFac = dFac
585 C a few abbreviations
586 eps1 = rUnit2mass*HeatCapacity_Cp
587 & *shiTransCoeffT(i,j,bi,bj)
588 eps2 = rUnit2mass*SHELFICElatentHeat
589 & *shiTransCoeffS(i,j,bi,bj)
590 eps5 = rUnit2mass*HeatCapacity_Cp
591 & *shiTransCoeffS(i,j,bi,bj)
592
593 C solve quadratic equation for salinity at shelfice-ocean interface
594 C note: this part of the code is not very intuitive as it involves
595 C many arbitrary abbreviations that were introduced to derive the
596 C correct form of the quadratic equation for salinity. The abbreviations
597 C only make sense in connection with my notes on this (M.Losch)
598 C
599 C eps3a was introduced as a constant variant of eps3 to avoid AD of
600 C code of typ (pLoc-const)/pLoc
601 eps3a = rhoShelfIce*SHELFICEheatCapacity_Cp
602 & * SHELFICEkappa * ( 1. _d 0 - dFac )
603 eps3 = eps3a/pLoc(I,J)
604 eps4 = b*pLoc(I,J) + c0
605 eps6 = eps4 - tLoc(I,J)
606 eps7 = eps4 - SHELFICEthetaSurface
607 eps8 = rUnit2mass*SHELFICEheatCapacity_Cp
608 & *shiTransCoeffS(i,j,bi,bj) * fwflxFac
609 aqe = a0 *(eps1+eps3-eps8)
610 recip_aqe = 0. _d 0
611 IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe
612 c bqe = eps1*eps6 + eps3*eps7 - eps2
613 bqe = eps1*eps6
614 & + eps3a*( b
615 & + ( c0 - SHELFICEthetaSurface )/pLoc(I,J) )
616 & - eps2
617 & + eps8*( a0*sLoc(I,J) - eps7 )
618 cqe = ( eps2 + eps8*eps7 )*sLoc(I,J)
619 discrim = bqe*bqe - 4. _d 0*aqe*cqe
620 #undef ALLOW_SHELFICE_DEBUG
621 #ifdef ALLOW_SHELFICE_DEBUG
622 IF ( discrim .LT. 0. _d 0 ) THEN
623 print *, 'ml-shelfice: discrim = ', discrim,aqe,bqe,cqe
624 print *, 'ml-shelfice: pLoc = ', pLoc(I,J)
625 print *, 'ml-shelfice: tLoc = ', tLoc(I,J)
626 print *, 'ml-shelfice: sLoc = ', sLoc(I,J)
627 print *, 'ml-shelfice: tsurface= ',
628 & SHELFICEthetaSurface
629 print *, 'ml-shelfice: eps1 = ', eps1
630 print *, 'ml-shelfice: eps2 = ', eps2
631 print *, 'ml-shelfice: eps3 = ', eps3
632 print *, 'ml-shelfice: eps4 = ', eps4
633 print *, 'ml-shelfice: eps5 = ', eps5
634 print *, 'ml-shelfice: eps6 = ', eps6
635 print *, 'ml-shelfice: eps7 = ', eps7
636 print *, 'ml-shelfice: eps8 = ', eps8
637 print *, 'ml-shelfice: rU2mass = ', rUnit2mass
638 print *, 'ml-shelfice: rhoIce = ', rhoShelfIce
639 print *, 'ml-shelfice: cFac = ', cFac
640 print *, 'ml-shelfice: Cp_W = ', HeatCapacity_Cp
641 print *, 'ml-shelfice: Cp_I = ',
642 & SHELFICEHeatCapacity_Cp
643 print *, 'ml-shelfice: gammaT = ',
644 & SHELFICEheatTransCoeff
645 print *, 'ml-shelfice: gammaS = ',
646 & SHELFICEsaltTransCoeff
647 print *, 'ml-shelfice: lat.heat= ',
648 & SHELFICElatentHeat
649 STOP 'ABNORMAL END in S/R SHELFICE_THERMODYNAMICS'
650 ENDIF
651 #endif /* ALLOW_SHELFICE_DEBUG */
652 saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
653 IF ( saltFreeze .LT. 0. _d 0 )
654 & saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
655 thetaFreeze = a0*saltFreeze + eps4
656 C-- upward fresh water flux due to melting (in kg/m^2/s)
657 cph change to identical form
658 cph freshWaterFlux = rUnit2mass
659 cph & * shiTransCoeffS(i,j,bi,bj)
660 cph & * ( saltFreeze - sLoc(I,J) ) / saltFreeze
661 freshWaterFlux = rUnit2mass
662 & * shiTransCoeffS(i,j,bi,bj)
663 & * ( 1. _d 0 - sLoc(I,J) / saltFreeze )
664 #ifdef ALLOW_SHIFWFLX_CONTROL
665 & + xx_shifwflx_loc(I,J,bi,bj)
666 #endif /* ALLOW_SHIFWFLX_CONTROL */
667 C-- Calculate the upward heat and fresh water fluxes;
668 C-- MITgcm sign conventions: downward (negative) fresh water flux
669 C-- implies melting and due to upward (positive) heat flux
670 shelfIceHeatFlux(I,J,bi,bj) =
671 & ( eps3
672 & - freshWaterFlux*SHELFICEheatCapacity_Cp*fwflxFac )
673 & * ( thetaFreeze - SHELFICEthetaSurface )
674 & - cFac*freshWaterFlux*( SHELFICElatentHeat
675 & - HeatCapacity_Cp*( thetaFreeze - rFac*tLoc(I,J) ) )
676 shelfIceFreshWaterFlux(I,J,bi,bj) = freshWaterFlux
677 C-- compute surface tendencies
678 shelficeForcingT(i,j,bi,bj) =
679 & ( shiTransCoeffT(i,j,bi,bj)
680 & - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )
681 & * ( thetaFreeze - tLoc(I,J) )
682 & - realFWfac*shelfIceFreshWaterFlux(I,J,bi,bj)*
683 & mass2rUnit*
684 & ( tLoc(I,J) - theta(I,J,K,bi,bj) )
685 shelficeForcingS(i,j,bi,bj) =
686 & ( shiTransCoeffS(i,j,bi,bj)
687 & - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )
688 & * ( saltFreeze - sLoc(I,J) )
689 & - realFWfac*shelfIceFreshWaterFlux(I,J,bi,bj)*
690 & mass2rUnit*
691 & ( sLoc(I,J) - salt(I,J,K,bi,bj) )
692 ELSE
693 shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
694 shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
695 shelficeForcingT (I,J,bi,bj) = 0. _d 0
696 shelficeForcingS (I,J,bi,bj) = 0. _d 0
697 ENDIF
698 ENDDO
699 ENDDO
700 ENDIF
701 C endif (not) useISOMIPTD
702 ENDDO
703 ENDDO
704
705 IF (SHELFICEMassStepping) THEN
706 ! CALL SHELFICE_STEP_ICEMASS( myTime, myIter, myThid )
707 ENDIF
708
709 C-- Calculate new loading anomaly (in case the ice-shelf mass was updated)
710 #ifndef ALLOW_AUTODIFF
711 c IF ( SHELFICEloadAnomalyFile .EQ. ' ' ) THEN
712 DO bj = myByLo(myThid), myByHi(myThid)
713 DO bi = myBxLo(myThid), myBxHi(myThid)
714 DO j = 1-OLy, sNy+OLy
715 DO i = 1-OLx, sNx+OLx
716 shelficeLoadAnomaly(i,j,bi,bj) = gravity
717 & *( shelficeMass(i,j,bi,bj) + rhoConst*Ro_surf(i,j,bi,bj) )
718 ENDDO
719 ENDDO
720 ENDDO
721 ENDDO
722 c ENDIF
723 #endif /* ndef ALLOW_AUTODIFF */
724
725 #ifdef ALLOW_DIAGNOSTICS
726 IF ( useDiagnostics ) THEN
727 CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx',
728 & 0,1,0,1,1,myThid)
729 CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux, 'SHIhtFlx',
730 & 0,1,0,1,1,myThid)
731 C SHIForcT (Ice shelf forcing for theta [W/m2], >0 increases theta)
732 tmpFac = HeatCapacity_Cp*rUnit2mass
733 CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingT,tmpFac,1,
734 & 'SHIForcT',0,1,0,1,1,myThid)
735 C SHIForcS (Ice shelf forcing for salt [g/m2/s], >0 increases salt)
736 tmpFac = rUnit2mass
737 CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingS,tmpFac,1,
738 & 'SHIForcS',0,1,0,1,1,myThid)
739 C Transfer coefficients
740 CALL DIAGNOSTICS_FILL(shiTransCoeffT,'SHIgammT',
741 & 0,1,0,1,1,myThid)
742 CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS',
743 & 0,1,0,1,1,myThid)
744 C Friction velocity
745 #ifdef SHI_ALLOW_GAMMAFRICT
746 IF ( SHELFICEuseGammaFrict )
747 & CALL DIAGNOSTICS_FILL(uStarDiag,'SHIuStar',0,1,0,1,1,myThid)
748 #endif /* SHI_ALLOW_GAMMAFRICT */
749 ENDIF
750 CALL DIAGNOSTICS_FILL(R_shelfice,'SHI_Rshelfice',
751 & 0,1,0,1,1,myThid)
752
753
754 #endif /* ALLOW_DIAGNOSTICS */
755
756 #endif /* ALLOW_SHELFICE */
757 RETURN
758 END
  ViewVC Help
Powered by ViewVC 1.1.22