/[MITgcm]/MITgcm_contrib/ksnow/press_release/code_expt/shelfice_thermodynamics.F
ViewVC logotype

Annotation of /MITgcm_contrib/ksnow/press_release/code_expt/shelfice_thermodynamics.F

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

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