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