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