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

Annotation of /MITgcm_contrib/verification_other/shelfice_remeshing/code/shelfice_thermodynamics.F

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

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