/[MITgcm]/MITgcm/pkg/shelfice/shelfice_thermodynamics.F
ViewVC logotype

Contents of /MITgcm/pkg/shelfice/shelfice_thermodynamics.F

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


Revision 1.23 - (show annotations) (download)
Wed Jun 29 16:24:10 2011 UTC (12 years, 10 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint63h, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint63
Changes since 1.22: +63 -17 lines
Implement friction velocity-dependent transfer coefficients following
Holland and Jenkins, JPO, 1999
Original code by M. Losch with small modifs.

1 C $Header: /u/gcmpack/MITgcm/pkg/shelfice/shelfice_thermodynamics.F,v 1.22 2011/05/30 20:21:51 jmc Exp $
2 C $Name: $
3
4 #include "SHELFICE_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: SHELFICE_THERMODYNAMICS
8 C !INTERFACE:
9 SUBROUTINE SHELFICE_THERMODYNAMICS(
10 I myTime, myIter, myThid )
11 C !DESCRIPTION: \bv
12 C *=============================================================*
13 C | S/R SHELFICE_THERMODYNAMICS
14 C | o shelf-ice main routine.
15 C | compute temperature and (virtual) salt flux at the
16 C | shelf-ice ocean interface
17 C |
18 C | stresses at the ice/water interface are computed in separate
19 C | routines that are called from mom_fluxform/mom_vecinv
20 C *=============================================================*
21
22 C !USES:
23 IMPLICIT NONE
24
25 C === Global variables ===
26 #include "SIZE.h"
27 #include "EEPARAMS.h"
28 #include "PARAMS.h"
29 #include "GRID.h"
30 #include "DYNVARS.h"
31 #include "FFIELDS.h"
32 #include "SHELFICE.h"
33 #ifdef ALLOW_AUTODIFF
34 # include "ctrl.h"
35 # include "ctrl_dummy.h"
36 #endif /* ALLOW_AUTODIFF */
37
38 C !INPUT/OUTPUT PARAMETERS:
39 C === Routine arguments ===
40 C myIter :: iteration counter for this thread
41 C myTime :: time counter for this thread
42 C myThid :: thread number for this instance of the routine.
43 _RL myTime
44 INTEGER myIter
45 INTEGER myThid
46 CEOP
47
48 #ifdef ALLOW_SHELFICE
49 C !LOCAL VARIABLES :
50 C === Local variables ===
51 C I,J,K,Kp1,bi,bj :: loop counters
52 C tLoc, sLoc, pLoc :: local in-situ temperature, salinity, pressure
53 C theta/saltFreeze :: temperature and salinity of water at the
54 C ice-ocean interface (at the freezing point)
55 C freshWaterFlux :: local variable for fresh water melt flux due
56 C to melting in kg/m^2/s
57 C (negative density x melt rate)
58 C convertFW2SaltLoc:: local copy of convertFW2Salt
59 C cFac :: 1 for conservative form, 0, otherwise
60 C auxiliary variables and abbreviations:
61 C a0, a1, a2, b, c0
62 C eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8
63 C aqe, bqe, cqe, discrim, recip_aqe
64 C drKp1, recip_drLoc
65 INTEGER I,J,K,Kp1
66 INTEGER bi,bj
67 _RL tLoc(1:sNx,1:sNy)
68 _RL sLoc(1:sNx,1:sNy)
69 _RL pLoc(1:sNx,1:sNy)
70 _RL thetaFreeze, saltFreeze
71 _RL freshWaterFlux, convertFW2SaltLoc
72 _RL a0, a1, a2, b, c0
73 _RL eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7
74 _RL cFac, rFac
75 _RL aqe, bqe, cqe, discrim, recip_aqe
76 _RL drKp1, recip_drLoc
77 _RL tmpFac
78
79 #ifdef SHI_ALLOW_GAMMAFRICT
80 _RL shiPr, shiSc, shiLo, shiKarman
81 _RL gammaTmoleT, gammaTmoleS, gammaTurb
82 _RL ustar, ustarSq, etastar
83 #endif
84
85 _RL SW_TEMP
86 EXTERNAL SW_TEMP
87
88 #ifdef ALLOW_SHIFWFLX_CONTROL
89 _RL xx_shifwflx_loc(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
90 #endif
91 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
92
93 #ifdef SHI_ALLOW_GAMMAFRICT
94 shiKarman= 0.4 _d 0
95 shiLo = 0. _d 0
96 shiPr = (viscArNr(1)/diffKrNrT(1))**(2/3)
97 shiSc = (viscArNr(1)/diffKrNrS(1))**(2/3)
98 gammaTmoleT = 12.5 _d 0 * shiPr - 6. _d 0
99 gammaTmoleS = 12.5 _d 0 * shiSc - 6. _d 0
100 #endif
101
102 C are we doing the conservative form of Jenkins et al. (2001)?
103 cFac = 0. _d 0
104 IF ( SHELFICEconserve ) cFac = 1. _d 0
105 C with "real fresh water flux" (affecting ETAN),
106 C there is more to modify
107 rFac = 1. _d 0
108 IF ( SHELFICEconserve .AND. useRealFreshWaterFlux ) rFac = 0. _d 0
109 C linear dependence of freezing point on salinity
110 a0 = -0.0575 _d 0
111 a1 = 0.0 _d -0
112 a2 = 0.0 _d -0
113 c0 = 0.0901 _d 0
114 b = -7.61 _d -4
115 #ifdef ALLOW_ISOMIP_TD
116 IF ( useISOMIPTD ) THEN
117 C non-linear dependence of freezing point on salinity
118 a0 = -0.0575 _d 0
119 a1 = 1.710523 _d -3
120 a2 = -2.154996 _d -4
121 b = -7.53 _d -4
122 c0 = 0. _d 0
123 ENDIF
124 convertFW2SaltLoc = convertFW2Salt
125 C hardcoding this value here is OK because it only applies to ISOMIP
126 C where this value is part of the protocol
127 IF ( convertFW2SaltLoc .EQ. -1. ) convertFW2SaltLoc = 33.4 _d 0
128 #endif /* ALLOW_ISOMIP_TD */
129
130 DO bj = myByLo(myThid), myByHi(myThid)
131 DO bi = myBxLo(myThid), myBxHi(myThid)
132 DO J = 1-Oly,sNy+Oly
133 DO I = 1-Olx,sNx+Olx
134 shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
135 shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
136 shelficeForcingT (I,J,bi,bj) = 0. _d 0
137 shelficeForcingS (I,J,bi,bj) = 0. _d 0
138 ENDDO
139 ENDDO
140 ENDDO
141 ENDDO
142 #ifdef ALLOW_SHIFWFLX_CONTROL
143 DO bj = myByLo(myThid), myByHi(myThid)
144 DO bi = myBxLo(myThid), myBxHi(myThid)
145 DO J = 1-Oly,sNy+Oly
146 DO I = 1-Olx,sNx+Olx
147 xx_shifwflx_loc(I,J,bi,bj) = 0. _d 0
148 ENDDO
149 ENDDO
150 ENDDO
151 ENDDO
152 CALL CTRL_GET_GEN (
153 & xx_shifwflx_file, xx_shifwflxstartdate, xx_shifwflxperiod,
154 & maskSHI, xx_shifwflx_loc, xx_shifwflx0, xx_shifwflx1,
155 & xx_shifwflx_dummy,
156 & xx_shifwflx_remo_intercept, xx_shifwflx_remo_slope,
157 & mytime, myiter, mythid )
158 #endif /* ALLOW_SHIFWFLX_CONTROL */
159 DO bj = myByLo(myThid), myByHi(myThid)
160 DO bi = myBxLo(myThid), myBxHi(myThid)
161 DO J = 1, sNy
162 DO I = 1, sNx
163 C-- make local copies of temperature, salinity and depth (pressure)
164 C-- underneath the ice
165 K = MAX(1,kTopC(I,J,bi,bj))
166 pLoc(I,J) = ABS(R_shelfIce(I,J,bi,bj))
167 tLoc(I,J) = theta(I,J,K,bi,bj)
168 sLoc(I,J) = MAX(salt(I,J,K,bi,bj), 0. _d 0)
169 ENDDO
170 ENDDO
171 IF ( SHELFICEBoundaryLayer ) THEN
172 C-- average over boundary layer width
173 DO J = 1, sNy
174 DO I = 1, sNx
175 K = kTopC(I,J,bi,bj)
176 IF ( K .NE. 0 .AND. K .LT. Nr ) THEN
177 Kp1 = MIN(Nr,K+1)
178 C-- overlap into lower cell
179 drKp1 = drF(K)*( 1. _d 0 - _hFacC(I,J,K,bi,bj) )
180 C-- lower cell may not be as thick as required
181 drKp1 = MIN( drKp1, drF(Kp1) * _hFacC(I,J,Kp1,bi,bj) )
182 recip_drLoc = 1. _d 0 /
183 & ( drF(K)*_hFacC(I,J,K,bi,bj) + drKp1 )
184 tLoc(I,J) = ( tLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
185 & + theta(I,J,Kp1,bi,bj) *drKp1 )
186 & * recip_drLoc
187 sLoc(I,J) = ( sLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
188 & + MAX(salt(I,J,Kp1,bi,bj), 0. _d 0) * drKp1 )
189 & * recip_drLoc
190 ENDIF
191 ENDDO
192 ENDDO
193 ENDIF
194
195 C-- turn potential temperature into in-situ temperature relative
196 C-- to the surface
197 DO J = 1, sNy
198 DO I = 1, sNx
199 tLoc(I,J) = SW_TEMP(sLoc(I,J),tLoc(I,J),pLoc(I,J),0.D0)
200 ENDDO
201 ENDDO
202
203 #ifdef ALLOW_ISOMIP_TD
204 IF ( useISOMIPTD ) THEN
205 DO J = 1, sNy
206 DO I = 1, sNx
207 K = kTopC(I,J,bi,bj)
208 IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN
209 C-- Calculate freezing temperature as a function of salinity and pressure
210 thetaFreeze =
211 & sLoc(I,J) * ( a0 + a1*sqrt(sLoc(I,J)) + a2*sLoc(I,J) )
212 & + b*pLoc(I,J) + c0
213 C-- Calculate the upward heat and fresh water fluxes
214 shelfIceHeatFlux(I,J,bi,bj) = maskC(I,J,K,bi,bj)
215 & * shiTransCoeffT(i,j,bi,bj)
216 & * ( tLoc(I,J) - thetaFreeze )
217 & * HeatCapacity_Cp*rUnit2mass
218 #ifdef ALLOW_SHIFWFLX_CONTROL
219 & - xx_shifwflx_loc(I,J,bi,bj)*SHELFICElatentHeat
220 #endif /* ALLOW_SHIFWFLX_CONTROL */
221 C upward heat flux into the shelf-ice implies basal melting,
222 C thus a downward (negative upward) fresh water flux (as a mass flux),
223 C and vice versa
224 shelfIceFreshWaterFlux(I,J,bi,bj) =
225 & - shelfIceHeatFlux(I,J,bi,bj)
226 & *recip_SHELFICElatentHeat
227 C-- compute surface tendencies
228 shelficeForcingT(i,j,bi,bj) =
229 & - shelfIceHeatFlux(I,J,bi,bj)
230 & *recip_Cp*mass2rUnit
231 & - cFac * shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit
232 & * ( thetaFreeze - tLoc(I,J) )
233 shelficeForcingS(i,j,bi,bj) =
234 & shelfIceFreshWaterFlux(I,J,bi,bj) * mass2rUnit
235 & * ( cFac*sLoc(I,J) + (1. _d 0-cFac)*convertFW2SaltLoc )
236 C-- stress at the ice/water interface is computed in separate
237 C routines that are called from mom_fluxform/mom_vecinv
238 ELSE
239 shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
240 shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
241 shelficeForcingT (I,J,bi,bj) = 0. _d 0
242 shelficeForcingS (I,J,bi,bj) = 0. _d 0
243 ENDIF
244 ENDDO
245 ENDDO
246 ELSE
247 #else
248 IF ( .TRUE. ) THEN
249 #endif /* ALLOW_ISOMIP_TD */
250 C use BRIOS thermodynamics, following Hellmers PhD thesis:
251 C Hellmer, H., 1989, A two-dimensional model for the thermohaline
252 C circulation under an ice shelf, Reports on Polar Research, No. 60
253 C (in German).
254
255 DO J = 1, sNy
256 DO I = 1, sNx
257 K = kTopC(I,J,bi,bj)
258 IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN
259 #ifdef SHI_ALLOW_GAMMAFRICT
260 C Implement friction velocity-dependent transfer coefficient
261 C of Holland and Jenkins, JPO, 1999
262 IF ( SHELFICEuseGammaFrict ) THEN
263 ustarSq = shiCdrag * MAX( 1D-06,
264 & (uVel(I,J,K,bi,bj)+uVel(I+1,J,K,bi,bj))
265 & *(uVel(I,J,K,bi,bj)+uVel(I+1,J,K,bi,bj))
266 & +(vVel(I,J,K,bi,bj)+vVel(I,J+1,K,bi,bj))
267 & *(vVel(I,J,K,bi,bj)+vVel(I,J+1,K,bi,bj)) )
268 ustar = SQRT(ustarSq)
269 C instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc))
270 etastar = 1. _d 0
271 gammaTurb = LOG( ustarSq * shiZetaN * etastar**2
272 & / ABS( fCori(I,J,bi,bj) * 0.5 _d 0 * viscArNr(k) ) )
273 & / shiKarman
274 & + 1. _d 0 / (2. _d 0 * shiZetaN*etastar)
275 & - 1. _d 0 / shiKarman
276 shiTransCoeffT(i,j,bi,bj) = MAX( 0.D-10,
277 & ustar/(gammaTurb + gammaTmoleT) )
278 shiTransCoeffS(i,j,bi,bj) = MAX( 0.D-10,
279 & ustar/(gammaTurb + gammaTmoleS) )
280 ENDIF
281 #endif
282 C a few abbreviations
283 eps1 = rUnit2mass*HeatCapacity_Cp
284 & *shiTransCoeffT(i,j,bi,bj)
285 eps2 = rUnit2mass*SHELFICElatentHeat
286 & *shiTransCoeffS(i,j,bi,bj)
287 eps5 = rUnit2mass*HeatCapacity_Cp
288 & *shiTransCoeffS(i,j,bi,bj)
289
290 C solve quadratic equation for salinity at shelfice-ocean interface
291 C note: this part of the code is not very intuitive as it involves
292 C many arbitrary abbreviations that were introduced to derive the
293 C correct form of the quadratic equation for salinity. The abbreviations
294 C only make sense in connection with my notes on this (M.Losch)
295 eps3 = rhoShelfIce*SHELFICEheatCapacity_Cp
296 & * SHELFICEkappa/pLoc(I,J)
297 cph introduce a constant variant of eps3 to avoid AD of
298 cph code of typ (pLoc-const)/pLoc
299 eps3a = rhoShelfIce*SHELFICEheatCapacity_Cp
300 & * SHELFICEkappa
301 eps4 = b*pLoc(I,J) + c0
302 eps6 = eps4 - tLoc(I,J)
303 eps7 = eps4 - SHELFICEthetaSurface
304 aqe = a0 *(eps1+eps3)
305 recip_aqe = 0. _d 0
306 IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe
307 c bqe = eps1*eps6 + eps3*eps7 - eps2
308 bqe = eps1*eps6
309 & + eps3a*( b
310 & + ( c0 - SHELFICEthetaSurface )/pLoc(I,J) )
311 & - eps2
312 cqe = eps2*sLoc(I,J)
313 discrim = bqe*bqe - 4. _d 0*aqe*cqe
314 #undef ALLOW_SHELFICE_DEBUG
315 #ifdef ALLOW_SHELFICE_DEBUG
316 IF ( discrim .LT. 0. _d 0 ) THEN
317 print *, 'ml-shelfice: discrim = ', discrim,aqe,bqe,cqe
318 print *, 'ml-shelfice: pLoc = ', pLoc(I,J)
319 print *, 'ml-shelfice: tLoc = ', tLoc(I,J)
320 print *, 'ml-shelfice: sLoc = ', sLoc(I,J)
321 print *, 'ml-shelfice: tsurface= ',
322 & SHELFICEthetaSurface
323 print *, 'ml-shelfice: eps1 = ', eps1
324 print *, 'ml-shelfice: eps2 = ', eps2
325 print *, 'ml-shelfice: eps3 = ', eps3
326 print *, 'ml-shelfice: eps4 = ', eps4
327 print *, 'ml-shelfice: eps5 = ', eps5
328 print *, 'ml-shelfice: eps6 = ', eps6
329 print *, 'ml-shelfice: eps7 = ', eps7
330 print *, 'ml-shelfice: rU2mass = ', rUnit2mass
331 print *, 'ml-shelfice: rhoIce = ', rhoShelfIce
332 print *, 'ml-shelfice: cFac = ', cFac
333 print *, 'ml-shelfice: Cp_W = ', HeatCapacity_Cp
334 print *, 'ml-shelfice: Cp_I = ',
335 & SHELFICEHeatCapacity_Cp
336 print *, 'ml-shelfice: gammaT = ',
337 & SHELFICEheatTransCoeff
338 print *, 'ml-shelfice: gammaS = ',
339 & SHELFICEsaltTransCoeff
340 print *, 'ml-shelfice: lat.heat= ',
341 & SHELFICElatentHeat
342 STOP 'ABNORMAL END in S/R SHELFICE_THERMODYNAMICS'
343 ENDIF
344 #endif /* ALLOW_SHELFICE_DEBUG */
345 saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
346 IF ( saltFreeze .LT. 0. _d 0 )
347 & saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
348 thetaFreeze = a0*saltFreeze + eps4
349 C-- upward fresh water flux due to melting (in kg/m^2/s)
350 cph change to identical form
351 cph freshWaterFlux = rUnit2mass
352 cph & * shiTransCoeffS(i,j,bi,bj)
353 cph & * ( saltFreeze - sLoc(I,J) ) / saltFreeze
354 freshWaterFlux = rUnit2mass
355 & * shiTransCoeffS(i,j,bi,bj)
356 & * ( 1. _d 0 - sLoc(I,J) / saltFreeze )
357 #ifdef ALLOW_SHIFWFLX_CONTROL
358 & + xx_shifwflx_loc(I,J,bi,bj)
359 #endif /* ALLOW_SHIFWFLX_CONTROL */
360 C-- Calculate the upward heat and fresh water fluxes;
361 C-- MITgcm sign conventions: downward (negative) fresh water flux
362 C-- implies melting and due to upward (positive) heat flux
363 shelfIceHeatFlux(I,J,bi,bj) =
364 & ( eps3*( thetaFreeze - SHELFICEthetaSurface )
365 & - cFac*freshWaterFlux*( SHELFICElatentHeat
366 & - HeatCapacity_Cp*( thetaFreeze - rFac*tLoc(I,J) ) )
367 & )
368 shelfIceFreshWaterFlux(I,J,bi,bj) = freshWaterFlux
369 C-- compute surface tendencies
370 shelficeForcingT(i,j,bi,bj) =
371 & ( shiTransCoeffT(i,j,bi,bj)
372 & - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )
373 & * ( thetaFreeze - tLoc(I,J) )
374 shelficeForcingS(i,j,bi,bj) =
375 & ( shiTransCoeffS(i,j,bi,bj)
376 & - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )
377 & * ( saltFreeze - sLoc(I,J) )
378 ELSE
379 shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
380 shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
381 shelficeForcingT (I,J,bi,bj) = 0. _d 0
382 shelficeForcingS (I,J,bi,bj) = 0. _d 0
383 ENDIF
384 ENDDO
385 ENDDO
386 ENDIF
387 C endif (not) useISOMIPTD
388 ENDDO
389 ENDDO
390
391 #ifdef ALLOW_DIAGNOSTICS
392 IF ( useDiagnostics ) THEN
393 CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx',
394 & 0,1,0,1,1,myThid)
395 CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux, 'SHIhtFlx',
396 & 0,1,0,1,1,myThid)
397 C SHIForcT (Ice shelf forcing for theta [W/m2], >0 increases theta)
398 tmpFac = HeatCapacity_Cp*rUnit2mass
399 CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingT,tmpFac,1,
400 & 'SHIForcT',0, 1,0,1,1,myThid)
401 C SHIForcS (Ice shelf forcing for salt [g/m2/s], >0 increases salt)
402 tmpFac = rUnit2mass
403 CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingS,tmpFac,1,
404 & 'SHIForcS',0, 1,0,1,1,myThid)
405 ENDIF
406 #endif /* ALLOW_DIAGNOSTICS */
407
408 #endif /* ALLOW_SHELFICE */
409 RETURN
410 END

  ViewVC Help
Powered by ViewVC 1.1.22