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

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

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

Revision 1.3 - (show annotations) (download)
Thu Sep 10 14:41:57 2015 UTC (9 years, 10 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +1 -1 lines
FILE REMOVED 
*** empty log message ***
1 C $Header: /u/gcmpack/MITgcm_contrib/shelfice_remeshing/code/shelfice_thermodynamics.F,v 1.2 2015/08/07 10:35:32 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
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
99 _RL aqe, bqe, cqe, discrim, recip_aqe
100 _RL drKp1, 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 C with "real fresh water flux" (affecting ETAN),
177 C there is more to modify
178 rFac = 1. _d 0
179 IF ( SHELFICEconserve .AND. useRealFreshWaterFlux ) rFac = 0. _d 0
180 C heat flux into the ice shelf, default is diffusive flux
181 C (Holland and Jenkins, 1999, eq.21)
182 dFac = 0. _d 0
183 IF ( SHELFICEadvDiffHeatFlux ) dFac = 1. _d 0
184 fwflxFac = 0. _d 0
185 C linear dependence of freezing point on salinity
186 a0 = -0.0575 _d 0
187 a1 = 0.0 _d -0
188 a2 = 0.0 _d -0
189 c0 = 0.0901 _d 0
190 b = -7.61 _d -4
191 #ifdef ALLOW_ISOMIP_TD
192 IF ( useISOMIPTD ) THEN
193 C non-linear dependence of freezing point on salinity
194 a0 = -0.0575 _d 0
195 a1 = 1.710523 _d -3
196 a2 = -2.154996 _d -4
197 b = -7.53 _d -4
198 c0 = 0. _d 0
199 ENDIF
200 convertFW2SaltLoc = convertFW2Salt
201 C hardcoding this value here is OK because it only applies to ISOMIP
202 C where this value is part of the protocol
203 IF ( convertFW2SaltLoc .EQ. -1. ) convertFW2SaltLoc = 33.4 _d 0
204 #endif /* ALLOW_ISOMIP_TD */
205
206 DO bj = myByLo(myThid), myByHi(myThid)
207 DO bi = myBxLo(myThid), myBxHi(myThid)
208 DO J = 1-OLy,sNy+OLy
209 DO I = 1-OLx,sNx+OLx
210 shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
211 shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
212 shelficeForcingT (I,J,bi,bj) = 0. _d 0
213 shelficeForcingS (I,J,bi,bj) = 0. _d 0
214 #if (defined SHI_ALLOW_GAMMAFRICT && defined ALLOW_DIAGNOSTICS)
215 uStarDiag (I,J,bi,bj) = 0. _d 0
216 #endif /* SHI_ALLOW_GAMMAFRICT and ALLOW_DIAGNOSTICS */
217 ENDDO
218 ENDDO
219 ENDDO
220 ENDDO
221 #ifdef ALLOW_SHIFWFLX_CONTROL
222 DO bj = myByLo(myThid), myByHi(myThid)
223 DO bi = myBxLo(myThid), myBxHi(myThid)
224 DO J = 1-OLy,sNy+OLy
225 DO I = 1-OLx,sNx+OLx
226 xx_shifwflx_loc(I,J,bi,bj) = 0. _d 0
227 ENDDO
228 ENDDO
229 ENDDO
230 ENDDO
231 #ifdef ALLOW_CTRL
232 if (useCTRL) CALL CTRL_GET_GEN (
233 & xx_shifwflx_file, xx_shifwflxstartdate, xx_shifwflxperiod,
234 & maskSHI, xx_shifwflx_loc, xx_shifwflx0, xx_shifwflx1,
235 & xx_shifwflx_dummy,
236 & xx_shifwflx_remo_intercept, xx_shifwflx_remo_slope,
237 & wshifwflx,
238 & myTime, myIter, myThid )
239 #endif
240 #endif /* ALLOW_SHIFWFLX_CONTROL */
241 DO bj = myByLo(myThid), myByHi(myThid)
242 DO bi = myBxLo(myThid), myBxHi(myThid)
243
244 IF ( SHELFICEBoundaryLayer ) THEN
245 C-- average over boundary layer width
246 DO J = 1, sNy+1
247 DO I = 1, sNx+1
248 u_topdr(I,J,bi,bj) = 0.0
249 v_topdr(I,J,bi,bj) = 0.0
250 ENDDO
251 ENDDO
252 ENDIF
253
254 #ifdef ALLOW_AUTODIFF_TAMC
255 # ifdef SHI_ALLOW_GAMMAFRICT
256 act1 = bi - myBxLo(myThid)
257 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
258 act2 = bj - myByLo(myThid)
259 max2 = myByHi(myThid) - myByLo(myThid) + 1
260 act3 = myThid - 1
261 max3 = nTx*nTy
262 act4 = ikey_dynamics - 1
263 ikey = (act1 + 1) + act2*max1
264 & + act3*max1*max2
265 & + act4*max1*max2*max3
266 # endif /* SHI_ALLOW_GAMMAFRICT */
267 #endif /* ALLOW_AUTODIFF_TAMC */
268 DO J = 1, sNy
269 DO I = 1, sNx
270 C-- make local copies of temperature, salinity and depth (pressure in deci-bar)
271 C-- underneath the ice
272 K = MAX(1,kTopC(I,J,bi,bj))
273 pLoc(I,J) = ABS(R_shelfIce(I,J,bi,bj))
274 c pLoc(I,J) = shelficeMass(I,J,bi,bj)*gravity*1. _d -4
275 tLoc(I,J) = theta(I,J,K,bi,bj)
276 sLoc(I,J) = MAX(salt(I,J,K,bi,bj), zeroRL)
277 c print *, 'JJHA' , hFacC(I,J,K,bi,bj)
278 IF ( .not.SHELFICEBoundaryLayer ) THEN
279 uLoc(I,J) = recip_hFacC(I,J,K,bi,bj) *
280 & ( uVel(I, J,K,bi,bj) * _hFacW(I, J,K,bi,bj)
281 & + uVel(I+1,J,K,bi,bj) * _hFacW(I+1,J,K,bi,bj) )
282 vLoc(I,J) = recip_hFacC(I,J,K,bi,bj) *
283 & ( vVel(I, J,K,bi,bj) * _hFacS(I, J,K,bi,bj)
284 & + vVel(I,J+1,K,bi,bj) * _hFacS(I,J+1,K,bi,bj) )
285 ENDIF
286 ! tmpfac =
287 ! & ( _hFacS(I,J, K,bi,bj) +
288 ! & _hFacS(I,J+1,K,bi,bj) ) / 2.0
289 ! vLoc(I,J) =
290 ! & ( vVel(I,J, K,bi,bj) * _hFacS(I,J, K,bi,bj)
291 ! & + vVel(I,J+1,K,bi,bj) * _hFacS(I,J+1,K,bi,bj) )
292 ! if (tmpfac .gt. 0.0) then
293 ! if (j.lt.100) then
294 ! print *, 1./tmpfac, recip_hFacC(I,J,K,bi,bj)
295 ! endif
296 ! vLoc(I,J) = vLoc(I,J) / tmpfac
297 ! else
298 ! vLoc(I,J) = 0.0
299 ! endif
300 ENDDO
301 ENDDO
302
303 IF ( SHELFICEBoundaryLayer ) THEN
304 DO J = 1, sNy+1
305 DO I = 1, sNx+1
306 K = ksurfW(I,J,bi,bj)
307 Kp1 = K+1
308 IF (K.lt.Nr) then
309 drKp1 = drF(K)*(1. _d 0-_hFacW(I,J,K,bi,bj))
310 drKp1 = max (drKp1, 0. _d 0)
311 recip_drLoc = 1.0 /
312 & (drF(K)*_hFacW(I,J,K,bi,bj)+drKp1)
313 u_topdr(I,J,bi,bj) =
314 & (drF(K)*_hFacW(I,J,K,bi,bj)*uVel(I,J,K,bi,bj) +
315 & drKp1*uVel(I,J,Kp1,bi,bj))
316 & * recip_drLoc
317 ELSE
318 u_topdr(I,J,bi,bj) = 0. _d 0
319 ENDIF
320
321 K = ksurfS(I,J,bi,bj)
322 Kp1 = K+1
323 IF (K.lt.Nr) then
324 drKp1 = drF(K)*(1. _d 0-_hFacS(I,J,K,bi,bj))
325 drKp1 = max (drKp1, 0. _d 0)
326 recip_drLoc = 1.0 /
327 & (drF(K)*_hFacS(I,J,K,bi,bj)+drKp1)
328 v_topdr(I,J,bi,bj) =
329 & (drF(K)*_hFacS(I,J,K,bi,bj)*vVel(I,J,K,bi,bj) +
330 & drKp1*vVel(I,J,Kp1,bi,bj))
331 & * recip_drLoc
332 ELSE
333 v_topdr(I,J,bi,bj) = 0. _d 0
334 ENDIF
335
336 ENDDO
337 ENDDO
338 ENDIF
339
340 IF ( SHELFICEBoundaryLayer ) THEN
341 C-- average over boundary layer width
342 DO J = 1, sNy
343 DO I = 1, sNx
344 K = kTopC(I,J,bi,bj)
345 IF ( K .NE. 0 .AND. K .LT. Nr ) THEN
346 Kp1 = MIN(Nr,K+1)
347 C-- overlap into lower cell
348 drKp1 = drF(K)*( 1. _d 0 - _hFacC(I,J,K,bi,bj) )
349 C-- Dans fix
350 drKp1 = MAX(drKp1, 0.)
351 C-- lower cell may not be as thick as required
352 drKp1 = MIN( drKp1, drF(Kp1) * _hFacC(I,J,Kp1,bi,bj) )
353 recip_drLoc = 1. _d 0 /
354 & ( drF(K)*_hFacC(I,J,K,bi,bj) + drKp1 )
355 tLoc(I,J) = ( tLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
356 & + theta(I,J,Kp1,bi,bj) *drKp1 )
357 & * recip_drLoc
358 sLoc(I,J) = ( sLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
359 & + MAX(salt(I,J,Kp1,bi,bj), zeroRL) * drKp1 )
360 & * recip_drLoc
361
362 ! uLoc(I,J) = ( uLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
363 ! & + drKp1 * recip_hFacC(I,J,Kp1,bi,bj) *
364 ! & ( uVel(I, J,Kp1,bi,bj) * _hFacW(I, J,Kp1,bi,bj)
365 ! & + uVel(I+1,J,Kp1,bi,bj) * _hFacW(I+1,J,Kp1,bi,bj) )
366 ! & ) * recip_drLoc
367 ! vLoc(I,J) = ( vLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
368 ! & + drKp1 * recip_hFacC(I,J,Kp1,bi,bj) *
369 ! & ( vVel(I,J, Kp1,bi,bj) * _hFacS(I,J, Kp1,bi,bj)
370 ! & + vVel(I,J+1,Kp1,bi,bj) * _hFacS(I,J+1,Kp1,bi,bj) )
371 ! & ) * recip_drLoc
372 ENDIF
373 ENDDO
374 ENDDO
375 ENDIF
376
377 IF ( SHELFICEBoundaryLayer ) THEN
378 DO J = 1, sNy
379 DO I = 1, sNx
380 uLoc(I,J) =
381 & u_topdr(I,J,bi,bj) + u_topdr(I+1,J,bi,bj)
382 vLoc(I,J) =
383 & v_topdr(I,J,bi,bj) + v_topdr(I,J+1,bi,bj)
384 ENDDO
385 ENDDO
386 ENDIF
387
388 C-- turn potential temperature into in-situ temperature relative
389 C-- to the surface
390 DO J = 1, sNy
391 DO I = 1, sNx
392 #ifndef ALLOW_OPENAD
393 tLoc(I,J) = SW_TEMP(sLoc(I,J),tLoc(I,J),pLoc(I,J),zeroRL)
394 #else
395 CALL SW_TEMP(sLoc(I,J),tLoc(I,J),pLoc(I,J),zeroRL,tLoc(I,J))
396 #endif
397 ENDDO
398 ENDDO
399
400 #ifdef SHI_ALLOW_GAMMAFRICT
401 IF ( SHELFICEuseGammaFrict ) THEN
402 DO J = 1, sNy
403 DO I = 1, sNx
404 K = kTopC(I,J,bi,bj)
405 IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN
406 ustarSq = shiCdrag * MAX( 1.D-6,
407 & 0.25 _d 0 *(uLoc(I,J)*uLoc(I,J)+vLoc(I,J)*vLoc(I,J)) )
408 ustar = SQRT(ustarSq)
409 #ifdef ALLOW_DIAGNOSTICS
410 uStarDiag(I,J,bi,bj) = ustar
411 #endif /* ALLOW_DIAGNOSTICS */
412 C instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc))
413 C etastar = 1. _d 0
414 C gammaTurbConst = 1. _d 0 / (2. _d 0 * shiZetaN*etastar)
415 C & - recip_shiKarman
416 IF ( fCori(I,J,bi,bj) .NE. 0. _d 0 ) THEN
417 gammaTurb = LOG( ustarSq * shiZetaN * etastar**2
418 & / ABS(fCori(I,J,bi,bj) * 5.0 _d 0 * shiKinVisc))
419 & * recip_shiKarman
420 & + gammaTurbConst
421 C Do we need to catch the unlikely case of very small ustar
422 C that can lead to negative gammaTurb?
423 C gammaTurb = MAX(0.D0, gammaTurb)
424 ELSE
425 gammaTurb = gammaTurbConst
426 ENDIF
427 shiTransCoeffT(i,j,bi,bj) = MAX( zeroRL,
428 & ustar/(gammaTurb + gammaTmoleT) )
429 shiTransCoeffS(i,j,bi,bj) = MAX( zeroRL,
430 & ustar/(gammaTurb + gammaTmoleS) )
431 ENDIF
432 ENDDO
433 ENDDO
434 ENDIF
435 #endif /* SHI_ALLOW_GAMMAFRICT */
436
437 #ifdef ALLOW_AUTODIFF_TAMC
438 # ifdef SHI_ALLOW_GAMMAFRICT
439 CADJ STORE shiTransCoeffS(:,:,bi,bj) = comlev1_bibj,
440 CADJ & key=ikey, byte=isbyte
441 CADJ STORE shiTransCoeffT(:,:,bi,bj) = comlev1_bibj,
442 CADJ & key=ikey, byte=isbyte
443 # endif /* SHI_ALLOW_GAMMAFRICT */
444 #endif /* ALLOW_AUTODIFF_TAMC */
445 #ifdef ALLOW_ISOMIP_TD
446 IF ( useISOMIPTD ) THEN
447 DO J = 1, sNy
448 DO I = 1, sNx
449 K = kTopC(I,J,bi,bj)
450 IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN
451 C-- Calculate freezing temperature as a function of salinity and pressure
452 thetaFreeze =
453 & sLoc(I,J) * ( a0 + a1*sqrt(sLoc(I,J)) + a2*sLoc(I,J) )
454 & + b*pLoc(I,J) + c0
455 C-- Calculate the upward heat and fresh water fluxes
456 shelfIceHeatFlux(I,J,bi,bj) = maskC(I,J,K,bi,bj)
457 & * shiTransCoeffT(i,j,bi,bj)
458 & * ( tLoc(I,J) - thetaFreeze )
459 & * HeatCapacity_Cp*rUnit2mass
460 #ifdef ALLOW_SHIFWFLX_CONTROL
461 & - xx_shifwflx_loc(I,J,bi,bj)*SHELFICElatentHeat
462 #endif /* ALLOW_SHIFWFLX_CONTROL */
463 C upward heat flux into the shelf-ice implies basal melting,
464 C thus a downward (negative upward) fresh water flux (as a mass flux),
465 C and vice versa
466 shelfIceFreshWaterFlux(I,J,bi,bj) =
467 & - shelfIceHeatFlux(I,J,bi,bj)
468 & *recip_latentHeat
469 C-- compute surface tendencies
470 shelficeForcingT(i,j,bi,bj) =
471 & - shelfIceHeatFlux(I,J,bi,bj)
472 & *recip_Cp*mass2rUnit
473 & - cFac * shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit
474 & * ( thetaFreeze - tLoc(I,J) )
475 shelficeForcingS(i,j,bi,bj) =
476 & shelfIceFreshWaterFlux(I,J,bi,bj) * mass2rUnit
477 & * ( cFac*sLoc(I,J) + (1. _d 0-cFac)*convertFW2SaltLoc )
478 C-- stress at the ice/water interface is computed in separate
479 C routines that are called from mom_fluxform/mom_vecinv
480 ELSE
481 shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
482 shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
483 shelficeForcingT (I,J,bi,bj) = 0. _d 0
484 shelficeForcingS (I,J,bi,bj) = 0. _d 0
485 ENDIF
486 ENDDO
487 ENDDO
488 ELSE
489 #else
490 IF ( .TRUE. ) THEN
491 #endif /* ALLOW_ISOMIP_TD */
492 C use BRIOS thermodynamics, following Hellmers PhD thesis:
493 C Hellmer, H., 1989, A two-dimensional model for the thermohaline
494 C circulation under an ice shelf, Reports on Polar Research, No. 60
495 C (in German).
496
497 DO J = 1, sNy
498 DO I = 1, sNx
499 K = kTopC(I,J,bi,bj)
500 IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN
501 C heat flux into the ice shelf, default is diffusive flux
502 C (Holland and Jenkins, 1999, eq.21)
503 thetaFreeze = a0*sLoc(I,J)+c0+b*pLoc(I,J)
504 fwflxFac = 0. _d 0
505 IF ( tLoc(I,J) .GT. thetaFreeze ) fwflxFac = dFac
506 C a few abbreviations
507 eps1 = rUnit2mass*HeatCapacity_Cp
508 & *shiTransCoeffT(i,j,bi,bj)
509 eps2 = rUnit2mass*SHELFICElatentHeat
510 & *shiTransCoeffS(i,j,bi,bj)
511 eps5 = rUnit2mass*HeatCapacity_Cp
512 & *shiTransCoeffS(i,j,bi,bj)
513
514 C solve quadratic equation for salinity at shelfice-ocean interface
515 C note: this part of the code is not very intuitive as it involves
516 C many arbitrary abbreviations that were introduced to derive the
517 C correct form of the quadratic equation for salinity. The abbreviations
518 C only make sense in connection with my notes on this (M.Losch)
519 C
520 C eps3a was introduced as a constant variant of eps3 to avoid AD of
521 C code of typ (pLoc-const)/pLoc
522 eps3a = rhoShelfIce*SHELFICEheatCapacity_Cp
523 & * SHELFICEkappa * ( 1. _d 0 - dFac )
524 eps3 = eps3a/pLoc(I,J)
525 eps4 = b*pLoc(I,J) + c0
526 eps6 = eps4 - tLoc(I,J)
527 eps7 = eps4 - SHELFICEthetaSurface
528 eps8 = rUnit2mass*SHELFICEheatCapacity_Cp
529 & *shiTransCoeffS(i,j,bi,bj) * fwflxFac
530 aqe = a0 *(eps1+eps3-eps8)
531 recip_aqe = 0. _d 0
532 IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe
533 c bqe = eps1*eps6 + eps3*eps7 - eps2
534 bqe = eps1*eps6
535 & + eps3a*( b
536 & + ( c0 - SHELFICEthetaSurface )/pLoc(I,J) )
537 & - eps2
538 & + eps8*( a0*sLoc(I,J) - eps7 )
539 cqe = ( eps2 + eps8*eps7 )*sLoc(I,J)
540 discrim = bqe*bqe - 4. _d 0*aqe*cqe
541 #undef ALLOW_SHELFICE_DEBUG
542 #ifdef ALLOW_SHELFICE_DEBUG
543 IF ( discrim .LT. 0. _d 0 ) THEN
544 print *, 'ml-shelfice: discrim = ', discrim,aqe,bqe,cqe
545 print *, 'ml-shelfice: pLoc = ', pLoc(I,J)
546 print *, 'ml-shelfice: tLoc = ', tLoc(I,J)
547 print *, 'ml-shelfice: sLoc = ', sLoc(I,J)
548 print *, 'ml-shelfice: tsurface= ',
549 & SHELFICEthetaSurface
550 print *, 'ml-shelfice: eps1 = ', eps1
551 print *, 'ml-shelfice: eps2 = ', eps2
552 print *, 'ml-shelfice: eps3 = ', eps3
553 print *, 'ml-shelfice: eps4 = ', eps4
554 print *, 'ml-shelfice: eps5 = ', eps5
555 print *, 'ml-shelfice: eps6 = ', eps6
556 print *, 'ml-shelfice: eps7 = ', eps7
557 print *, 'ml-shelfice: eps8 = ', eps8
558 print *, 'ml-shelfice: rU2mass = ', rUnit2mass
559 print *, 'ml-shelfice: rhoIce = ', rhoShelfIce
560 print *, 'ml-shelfice: cFac = ', cFac
561 print *, 'ml-shelfice: Cp_W = ', HeatCapacity_Cp
562 print *, 'ml-shelfice: Cp_I = ',
563 & SHELFICEHeatCapacity_Cp
564 print *, 'ml-shelfice: gammaT = ',
565 & SHELFICEheatTransCoeff
566 print *, 'ml-shelfice: gammaS = ',
567 & SHELFICEsaltTransCoeff
568 print *, 'ml-shelfice: lat.heat= ',
569 & SHELFICElatentHeat
570 STOP 'ABNORMAL END in S/R SHELFICE_THERMODYNAMICS'
571 ENDIF
572 #endif /* ALLOW_SHELFICE_DEBUG */
573 saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
574 IF ( saltFreeze .LT. 0. _d 0 )
575 & saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
576 thetaFreeze = a0*saltFreeze + eps4
577 C-- upward fresh water flux due to melting (in kg/m^2/s)
578 cph change to identical form
579 cph freshWaterFlux = rUnit2mass
580 cph & * shiTransCoeffS(i,j,bi,bj)
581 cph & * ( saltFreeze - sLoc(I,J) ) / saltFreeze
582 freshWaterFlux = rUnit2mass
583 & * shiTransCoeffS(i,j,bi,bj)
584 & * ( 1. _d 0 - sLoc(I,J) / saltFreeze )
585 #ifdef ALLOW_SHIFWFLX_CONTROL
586 & + xx_shifwflx_loc(I,J,bi,bj)
587 #endif /* ALLOW_SHIFWFLX_CONTROL */
588 C-- Calculate the upward heat and fresh water fluxes;
589 C-- MITgcm sign conventions: downward (negative) fresh water flux
590 C-- implies melting and due to upward (positive) heat flux
591 shelfIceHeatFlux(I,J,bi,bj) =
592 & ( eps3
593 & - freshWaterFlux*SHELFICEheatCapacity_Cp*fwflxFac )
594 & * ( thetaFreeze - SHELFICEthetaSurface )
595 & - cFac*freshWaterFlux*( SHELFICElatentHeat
596 & - HeatCapacity_Cp*( thetaFreeze - rFac*tLoc(I,J) ) )
597 shelfIceFreshWaterFlux(I,J,bi,bj) = freshWaterFlux
598 C-- compute surface tendencies
599 shelficeForcingT(i,j,bi,bj) =
600 & ( shiTransCoeffT(i,j,bi,bj)
601 & - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )
602 & * ( thetaFreeze - tLoc(I,J) )
603 shelficeForcingS(i,j,bi,bj) =
604 & ( shiTransCoeffS(i,j,bi,bj)
605 & - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )
606 & * ( saltFreeze - sLoc(I,J) )
607 ELSE
608 shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
609 shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
610 shelficeForcingT (I,J,bi,bj) = 0. _d 0
611 shelficeForcingS (I,J,bi,bj) = 0. _d 0
612 ENDIF
613 ENDDO
614 ENDDO
615 ENDIF
616 C endif (not) useISOMIPTD
617 ENDDO
618 ENDDO
619
620 IF (SHELFICEMassStepping) THEN
621 CALL SHELFICE_STEP_ICEMASS( myTime, myIter, myThid )
622 ENDIF
623
624 C-- Calculate new loading anomaly (in case the ice-shelf mass was updated)
625 #ifndef ALLOW_AUTODIFF
626 c IF ( SHELFICEloadAnomalyFile .EQ. ' ' ) THEN
627 DO bj = myByLo(myThid), myByHi(myThid)
628 DO bi = myBxLo(myThid), myBxHi(myThid)
629 DO j = 1-OLy, sNy+OLy
630 DO i = 1-OLx, sNx+OLx
631 shelficeLoadAnomaly(i,j,bi,bj) = gravity
632 & *( shelficeMass(i,j,bi,bj) + rhoConst*Ro_surf(i,j,bi,bj) )
633 ENDDO
634 ENDDO
635 ENDDO
636 ENDDO
637 c ENDIF
638 #endif /* ndef ALLOW_AUTODIFF */
639
640 #ifdef ALLOW_DIAGNOSTICS
641 IF ( useDiagnostics ) THEN
642 CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx',
643 & 0,1,0,1,1,myThid)
644 CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux, 'SHIhtFlx',
645 & 0,1,0,1,1,myThid)
646 C SHIForcT (Ice shelf forcing for theta [W/m2], >0 increases theta)
647 tmpFac = HeatCapacity_Cp*rUnit2mass
648 CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingT,tmpFac,1,
649 & 'SHIForcT',0,1,0,1,1,myThid)
650 C SHIForcS (Ice shelf forcing for salt [g/m2/s], >0 increases salt)
651 tmpFac = rUnit2mass
652 CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingS,tmpFac,1,
653 & 'SHIForcS',0,1,0,1,1,myThid)
654 C Transfer coefficients
655 CALL DIAGNOSTICS_FILL(shiTransCoeffT,'SHIgammT',
656 & 0,1,0,1,1,myThid)
657 CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS',
658 & 0,1,0,1,1,myThid)
659 C Friction velocity
660 #ifdef SHI_ALLOW_GAMMAFRICT
661 IF ( SHELFICEuseGammaFrict )
662 & CALL DIAGNOSTICS_FILL(uStarDiag,'SHIuStar',0,1,0,1,1,myThid)
663 #endif /* SHI_ALLOW_GAMMAFRICT */
664 ENDIF
665 #endif /* ALLOW_DIAGNOSTICS */
666
667 #endif /* ALLOW_SHELFICE */
668
669 RETURN
670
671 END
  ViewVC Help
Powered by ViewVC 1.1.22