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

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

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


Revision 1.2 - (hide annotations) (download)
Fri Aug 7 10:35:32 2015 UTC (9 years, 11 months ago) by dgoldberg
Branch: MAIN
Changes since 1.1: +0 -0 lines
*** empty log message ***

1 dgoldberg 1.1 C $Header: /u/gcmpack/MITgcm/pkg/shelfice/shelfice_thermodynamics.F,v 1.44 2015/02/15 15:46:24 jmc Exp $
2     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    
53     C !INPUT/OUTPUT PARAMETERS:
54     C === Routine arguments ===
55     C myIter :: iteration counter for this thread
56     C myTime :: time counter for this thread
57     C myThid :: thread number for this instance of the routine.
58     _RL myTime
59     INTEGER myIter
60     INTEGER myThid
61    
62     #ifdef ALLOW_SHELFICE
63     C !LOCAL VARIABLES :
64     C === Local variables ===
65     C I,J,K,Kp1,bi,bj :: loop counters
66     C tLoc, sLoc, pLoc :: local in-situ temperature, salinity, pressure
67     C theta/saltFreeze :: temperature and salinity of water at the
68     C ice-ocean interface (at the freezing point)
69     C freshWaterFlux :: local variable for fresh water melt flux due
70     C to melting in kg/m^2/s
71     C (negative density x melt rate)
72     C convertFW2SaltLoc:: local copy of convertFW2Salt
73     C cFac :: 1 for conservative form, 0, otherwise
74     C rFac :: realFreshWaterFlux factor
75     C dFac :: 0 for diffusive heat flux (Holland and Jenkins, 1999,
76     C eq21)
77     C 1 for advective and diffusive heat flux (eq22, 26, 31)
78     C fwflxFac :: only effective for dFac=1, 1 if we expect a melting
79     C fresh water flux, 0 otherwise
80     C auxiliary variables and abbreviations:
81     C a0, a1, a2, b, c0
82     C eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8
83     C aqe, bqe, cqe, discrim, recip_aqe
84     C drKp1, recip_drLoc
85     INTEGER I,J,K,Kp1
86     INTEGER bi,bj
87     _RL tLoc(1:sNx,1:sNy)
88     _RL sLoc(1:sNx,1:sNy)
89     _RL pLoc(1:sNx,1:sNy)
90     _RL uLoc(1:sNx,1:sNy)
91     _RL vLoc(1:sNx,1:sNy)
92     _RL u_topdr(1:sNx+1,1:sNy+1,nSx,nSy)
93     _RL v_topdr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
94     _RL thetaFreeze, saltFreeze, recip_Cp
95     _RL freshWaterFlux, convertFW2SaltLoc
96     _RL a0, a1, a2, b, c0
97     _RL eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8
98     _RL cFac, rFac, dFac, fwflxFac
99     _RL aqe, bqe, cqe, discrim, recip_aqe
100     _RL drKp1, recip_drLoc
101     _RL recip_latentHeat
102     _RL tmpFac
103    
104     #ifdef SHI_ALLOW_GAMMAFRICT
105     _RL shiPr, shiSc, shiLo, recip_shiKarman, shiTwoThirds
106     _RL gammaTmoleT, gammaTmoleS, gammaTurb, gammaTurbConst
107     _RL ustar, ustarSq, etastar
108     PARAMETER ( shiTwoThirds = 0.66666666666666666666666666667D0 )
109     #ifdef ALLOW_DIAGNOSTICS
110     _RL uStarDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
111     #endif /* ALLOW_DIAGNOSTICS */
112     #endif
113    
114     #ifndef ALLOW_OPENAD
115     _RL SW_TEMP
116     EXTERNAL SW_TEMP
117     #endif
118    
119     #ifdef ALLOW_SHIFWFLX_CONTROL
120     _RL xx_shifwflx_loc(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
121     #endif
122     CEOP
123     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
124    
125     #ifdef SHI_ALLOW_GAMMAFRICT
126     #ifdef ALLOW_AUTODIFF
127     C re-initialize here again, curtesy to TAF
128     DO bj = myByLo(myThid), myByHi(myThid)
129     DO bi = myBxLo(myThid), myBxHi(myThid)
130     DO J = 1-OLy,sNy+OLy
131     DO I = 1-OLx,sNx+OLx
132     shiTransCoeffT(i,j,bi,bj) = SHELFICEheatTransCoeff
133     shiTransCoeffS(i,j,bi,bj) = SHELFICEsaltTransCoeff
134     ENDDO
135     ENDDO
136     ENDDO
137     ENDDO
138     #endif /* ALLOW_AUTODIFF */
139     IF ( SHELFICEuseGammaFrict ) THEN
140     C Implement friction velocity-dependent transfer coefficient
141     C of Holland and Jenkins, JPO, 1999
142     recip_shiKarman= 1. _d 0 / 0.4 _d 0
143     shiLo = 0. _d 0
144     shiPr = shiPrandtl**shiTwoThirds
145     shiSc = shiSchmidt**shiTwoThirds
146     cph shiPr = (viscArNr(1)/diffKrNrT(1))**shiTwoThirds
147     cph shiSc = (viscArNr(1)/diffKrNrS(1))**shiTwoThirds
148     gammaTmoleT = 12.5 _d 0 * shiPr - 6. _d 0
149     gammaTmoleS = 12.5 _d 0 * shiSc - 6. _d 0
150     C instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc))
151     etastar = 1. _d 0
152     gammaTurbConst = 1. _d 0 / (2. _d 0 * shiZetaN*etastar)
153     & - recip_shiKarman
154     #ifdef ALLOW_AUTODIFF
155     DO bj = myByLo(myThid), myByHi(myThid)
156     DO bi = myBxLo(myThid), myBxHi(myThid)
157     DO J = 1-OLy,sNy+OLy
158     DO I = 1-OLx,sNx+OLx
159     shiTransCoeffT(i,j,bi,bj) = 0. _d 0
160     shiTransCoeffS(i,j,bi,bj) = 0. _d 0
161     ENDDO
162     ENDDO
163     ENDDO
164     ENDDO
165     #endif /* ALLOW_AUTODIFF */
166     ENDIF
167     #endif /* SHI_ALLOW_GAMMAFRICT */
168    
169     recip_latentHeat = 0. _d 0
170     IF ( SHELFICElatentHeat .NE. 0. _d 0 )
171     & recip_latentHeat = 1. _d 0/SHELFICElatentHeat
172     C are we doing the conservative form of Jenkins et al. (2001)?
173     recip_Cp = 1. _d 0 / HeatCapacity_Cp
174     cFac = 0. _d 0
175     IF ( SHELFICEconserve ) cFac = 1. _d 0
176     C with "real fresh water flux" (affecting ETAN),
177     C there is more to modify
178     rFac = 1. _d 0
179     IF ( SHELFICEconserve .AND. useRealFreshWaterFlux ) rFac = 0. _d 0
180     C heat flux into the ice shelf, default is diffusive flux
181     C (Holland and Jenkins, 1999, eq.21)
182     dFac = 0. _d 0
183     IF ( SHELFICEadvDiffHeatFlux ) dFac = 1. _d 0
184     fwflxFac = 0. _d 0
185     C linear dependence of freezing point on salinity
186     a0 = -0.0575 _d 0
187     a1 = 0.0 _d -0
188     a2 = 0.0 _d -0
189     c0 = 0.0901 _d 0
190     b = -7.61 _d -4
191     #ifdef ALLOW_ISOMIP_TD
192     IF ( useISOMIPTD ) THEN
193     C non-linear dependence of freezing point on salinity
194     a0 = -0.0575 _d 0
195     a1 = 1.710523 _d -3
196     a2 = -2.154996 _d -4
197     b = -7.53 _d -4
198     c0 = 0. _d 0
199     ENDIF
200     convertFW2SaltLoc = convertFW2Salt
201     C hardcoding this value here is OK because it only applies to ISOMIP
202     C where this value is part of the protocol
203     IF ( convertFW2SaltLoc .EQ. -1. ) convertFW2SaltLoc = 33.4 _d 0
204     #endif /* ALLOW_ISOMIP_TD */
205    
206     DO bj = myByLo(myThid), myByHi(myThid)
207     DO bi = myBxLo(myThid), myBxHi(myThid)
208     DO J = 1-OLy,sNy+OLy
209     DO I = 1-OLx,sNx+OLx
210     shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
211     shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
212     shelficeForcingT (I,J,bi,bj) = 0. _d 0
213     shelficeForcingS (I,J,bi,bj) = 0. _d 0
214     #if (defined SHI_ALLOW_GAMMAFRICT && defined ALLOW_DIAGNOSTICS)
215     uStarDiag (I,J,bi,bj) = 0. _d 0
216     #endif /* SHI_ALLOW_GAMMAFRICT and ALLOW_DIAGNOSTICS */
217     ENDDO
218     ENDDO
219     ENDDO
220     ENDDO
221     #ifdef ALLOW_SHIFWFLX_CONTROL
222     DO bj = myByLo(myThid), myByHi(myThid)
223     DO bi = myBxLo(myThid), myBxHi(myThid)
224     DO J = 1-OLy,sNy+OLy
225     DO I = 1-OLx,sNx+OLx
226     xx_shifwflx_loc(I,J,bi,bj) = 0. _d 0
227     ENDDO
228     ENDDO
229     ENDDO
230     ENDDO
231     #ifdef ALLOW_CTRL
232     if (useCTRL) CALL CTRL_GET_GEN (
233     & xx_shifwflx_file, xx_shifwflxstartdate, xx_shifwflxperiod,
234     & maskSHI, xx_shifwflx_loc, xx_shifwflx0, xx_shifwflx1,
235     & xx_shifwflx_dummy,
236     & xx_shifwflx_remo_intercept, xx_shifwflx_remo_slope,
237     & wshifwflx,
238     & myTime, myIter, myThid )
239     #endif
240     #endif /* ALLOW_SHIFWFLX_CONTROL */
241     DO bj = myByLo(myThid), myByHi(myThid)
242     DO bi = myBxLo(myThid), myBxHi(myThid)
243    
244     IF ( SHELFICEBoundaryLayer ) THEN
245     C-- average over boundary layer width
246     DO J = 1, sNy+1
247     DO I = 1, sNx+1
248     u_topdr(I,J,bi,bj) = 0.0
249     v_topdr(I,J,bi,bj) = 0.0
250     ENDDO
251     ENDDO
252     ENDIF
253    
254     #ifdef ALLOW_AUTODIFF_TAMC
255     # ifdef SHI_ALLOW_GAMMAFRICT
256     act1 = bi - myBxLo(myThid)
257     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
258     act2 = bj - myByLo(myThid)
259     max2 = myByHi(myThid) - myByLo(myThid) + 1
260     act3 = myThid - 1
261     max3 = nTx*nTy
262     act4 = ikey_dynamics - 1
263     ikey = (act1 + 1) + act2*max1
264     & + act3*max1*max2
265     & + act4*max1*max2*max3
266     # endif /* SHI_ALLOW_GAMMAFRICT */
267     #endif /* ALLOW_AUTODIFF_TAMC */
268     DO J = 1, sNy
269     DO I = 1, sNx
270     C-- make local copies of temperature, salinity and depth (pressure in deci-bar)
271     C-- underneath the ice
272     K = MAX(1,kTopC(I,J,bi,bj))
273     pLoc(I,J) = ABS(R_shelfIce(I,J,bi,bj))
274     c pLoc(I,J) = shelficeMass(I,J,bi,bj)*gravity*1. _d -4
275     tLoc(I,J) = theta(I,J,K,bi,bj)
276     sLoc(I,J) = MAX(salt(I,J,K,bi,bj), zeroRL)
277     c print *, 'JJHA' , hFacC(I,J,K,bi,bj)
278     IF ( .not.SHELFICEBoundaryLayer ) THEN
279     uLoc(I,J) = recip_hFacC(I,J,K,bi,bj) *
280     & ( uVel(I, J,K,bi,bj) * _hFacW(I, J,K,bi,bj)
281     & + uVel(I+1,J,K,bi,bj) * _hFacW(I+1,J,K,bi,bj) )
282     vLoc(I,J) = recip_hFacC(I,J,K,bi,bj) *
283     & ( vVel(I, J,K,bi,bj) * _hFacS(I, J,K,bi,bj)
284     & + vVel(I,J+1,K,bi,bj) * _hFacS(I,J+1,K,bi,bj) )
285     ENDIF
286     ! tmpfac =
287     ! & ( _hFacS(I,J, K,bi,bj) +
288     ! & _hFacS(I,J+1,K,bi,bj) ) / 2.0
289     ! vLoc(I,J) =
290     ! & ( vVel(I,J, K,bi,bj) * _hFacS(I,J, K,bi,bj)
291     ! & + vVel(I,J+1,K,bi,bj) * _hFacS(I,J+1,K,bi,bj) )
292     ! if (tmpfac .gt. 0.0) then
293     ! if (j.lt.100) then
294     ! print *, 1./tmpfac, recip_hFacC(I,J,K,bi,bj)
295     ! endif
296     ! vLoc(I,J) = vLoc(I,J) / tmpfac
297     ! else
298     ! vLoc(I,J) = 0.0
299     ! endif
300     ENDDO
301     ENDDO
302    
303     IF ( SHELFICEBoundaryLayer ) THEN
304     DO J = 1, sNy+1
305     DO I = 1, sNx+1
306     K = ksurfW(I,J,bi,bj)
307     Kp1 = K+1
308     IF (K.lt.Nr) then
309     drKp1 = drF(K)*(1. _d 0-_hFacW(I,J,K,bi,bj))
310     drKp1 = max (drKp1, 0. _d 0)
311     recip_drLoc = 1.0 /
312     & (drF(K)*_hFacW(I,J,K,bi,bj)+drKp1)
313     u_topdr(I,J,bi,bj) =
314     & (drF(K)*_hFacW(I,J,K,bi,bj)*uVel(I,J,K,bi,bj) +
315     & drKp1*uVel(I,J,Kp1,bi,bj))
316     & * recip_drLoc
317     ELSE
318     u_topdr(I,J,bi,bj) = 0. _d 0
319     ENDIF
320    
321     K = ksurfS(I,J,bi,bj)
322     Kp1 = K+1
323     IF (K.lt.Nr) then
324     drKp1 = drF(K)*(1. _d 0-_hFacS(I,J,K,bi,bj))
325     drKp1 = max (drKp1, 0. _d 0)
326     recip_drLoc = 1.0 /
327     & (drF(K)*_hFacS(I,J,K,bi,bj)+drKp1)
328     v_topdr(I,J,bi,bj) =
329     & (drF(K)*_hFacS(I,J,K,bi,bj)*vVel(I,J,K,bi,bj) +
330     & drKp1*vVel(I,J,Kp1,bi,bj))
331     & * recip_drLoc
332     ELSE
333     v_topdr(I,J,bi,bj) = 0. _d 0
334     ENDIF
335    
336     ENDDO
337     ENDDO
338     ENDIF
339    
340     IF ( SHELFICEBoundaryLayer ) THEN
341     C-- average over boundary layer width
342     DO J = 1, sNy
343     DO I = 1, sNx
344     K = kTopC(I,J,bi,bj)
345     IF ( K .NE. 0 .AND. K .LT. Nr ) THEN
346     Kp1 = MIN(Nr,K+1)
347     C-- overlap into lower cell
348     drKp1 = drF(K)*( 1. _d 0 - _hFacC(I,J,K,bi,bj) )
349     C-- Dans fix
350     drKp1 = MAX(drKp1, 0.)
351     C-- lower cell may not be as thick as required
352     drKp1 = MIN( drKp1, drF(Kp1) * _hFacC(I,J,Kp1,bi,bj) )
353     recip_drLoc = 1. _d 0 /
354     & ( drF(K)*_hFacC(I,J,K,bi,bj) + drKp1 )
355     tLoc(I,J) = ( tLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
356     & + theta(I,J,Kp1,bi,bj) *drKp1 )
357     & * recip_drLoc
358     sLoc(I,J) = ( sLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
359     & + MAX(salt(I,J,Kp1,bi,bj), zeroRL) * drKp1 )
360     & * recip_drLoc
361    
362     ! uLoc(I,J) = ( uLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
363     ! & + drKp1 * recip_hFacC(I,J,Kp1,bi,bj) *
364     ! & ( uVel(I, J,Kp1,bi,bj) * _hFacW(I, J,Kp1,bi,bj)
365     ! & + uVel(I+1,J,Kp1,bi,bj) * _hFacW(I+1,J,Kp1,bi,bj) )
366     ! & ) * recip_drLoc
367     ! vLoc(I,J) = ( vLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
368     ! & + drKp1 * recip_hFacC(I,J,Kp1,bi,bj) *
369     ! & ( vVel(I,J, Kp1,bi,bj) * _hFacS(I,J, Kp1,bi,bj)
370     ! & + vVel(I,J+1,Kp1,bi,bj) * _hFacS(I,J+1,Kp1,bi,bj) )
371     ! & ) * recip_drLoc
372     ENDIF
373     ENDDO
374     ENDDO
375     ENDIF
376    
377     IF ( SHELFICEBoundaryLayer ) THEN
378     DO J = 1, sNy
379     DO I = 1, sNx
380     uLoc(I,J) =
381     & u_topdr(I,J,bi,bj) + u_topdr(I+1,J,bi,bj)
382     vLoc(I,J) =
383     & v_topdr(I,J,bi,bj) + v_topdr(I,J+1,bi,bj)
384     ENDDO
385     ENDDO
386     ENDIF
387    
388     C-- turn potential temperature into in-situ temperature relative
389     C-- to the surface
390     DO J = 1, sNy
391     DO I = 1, sNx
392     #ifndef ALLOW_OPENAD
393     tLoc(I,J) = SW_TEMP(sLoc(I,J),tLoc(I,J),pLoc(I,J),zeroRL)
394     #else
395     CALL SW_TEMP(sLoc(I,J),tLoc(I,J),pLoc(I,J),zeroRL,tLoc(I,J))
396     #endif
397     ENDDO
398     ENDDO
399    
400     #ifdef SHI_ALLOW_GAMMAFRICT
401     IF ( SHELFICEuseGammaFrict ) THEN
402     DO J = 1, sNy
403     DO I = 1, sNx
404     K = kTopC(I,J,bi,bj)
405     IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN
406     ustarSq = shiCdrag * MAX( 1.D-6,
407     & 0.25 _d 0 *(uLoc(I,J)*uLoc(I,J)+vLoc(I,J)*vLoc(I,J)) )
408     ustar = SQRT(ustarSq)
409     #ifdef ALLOW_DIAGNOSTICS
410     uStarDiag(I,J,bi,bj) = ustar
411     #endif /* ALLOW_DIAGNOSTICS */
412     C instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc))
413     C etastar = 1. _d 0
414     C gammaTurbConst = 1. _d 0 / (2. _d 0 * shiZetaN*etastar)
415     C & - recip_shiKarman
416     IF ( fCori(I,J,bi,bj) .NE. 0. _d 0 ) THEN
417     gammaTurb = LOG( ustarSq * shiZetaN * etastar**2
418     & / ABS(fCori(I,J,bi,bj) * 5.0 _d 0 * shiKinVisc))
419     & * recip_shiKarman
420     & + gammaTurbConst
421     C Do we need to catch the unlikely case of very small ustar
422     C that can lead to negative gammaTurb?
423     C gammaTurb = MAX(0.D0, gammaTurb)
424     ELSE
425     gammaTurb = gammaTurbConst
426     ENDIF
427     shiTransCoeffT(i,j,bi,bj) = MAX( zeroRL,
428     & ustar/(gammaTurb + gammaTmoleT) )
429     shiTransCoeffS(i,j,bi,bj) = MAX( zeroRL,
430     & ustar/(gammaTurb + gammaTmoleS) )
431     ENDIF
432     ENDDO
433     ENDDO
434     ENDIF
435     #endif /* SHI_ALLOW_GAMMAFRICT */
436    
437     #ifdef ALLOW_AUTODIFF_TAMC
438     # ifdef SHI_ALLOW_GAMMAFRICT
439     CADJ STORE shiTransCoeffS(:,:,bi,bj) = comlev1_bibj,
440     CADJ & key=ikey, byte=isbyte
441     CADJ STORE shiTransCoeffT(:,:,bi,bj) = comlev1_bibj,
442     CADJ & key=ikey, byte=isbyte
443     # endif /* SHI_ALLOW_GAMMAFRICT */
444     #endif /* ALLOW_AUTODIFF_TAMC */
445     #ifdef ALLOW_ISOMIP_TD
446     IF ( useISOMIPTD ) THEN
447     DO J = 1, sNy
448     DO I = 1, sNx
449     K = kTopC(I,J,bi,bj)
450     IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN
451     C-- Calculate freezing temperature as a function of salinity and pressure
452     thetaFreeze =
453     & sLoc(I,J) * ( a0 + a1*sqrt(sLoc(I,J)) + a2*sLoc(I,J) )
454     & + b*pLoc(I,J) + c0
455     C-- Calculate the upward heat and fresh water fluxes
456     shelfIceHeatFlux(I,J,bi,bj) = maskC(I,J,K,bi,bj)
457     & * shiTransCoeffT(i,j,bi,bj)
458     & * ( tLoc(I,J) - thetaFreeze )
459     & * HeatCapacity_Cp*rUnit2mass
460     #ifdef ALLOW_SHIFWFLX_CONTROL
461     & - xx_shifwflx_loc(I,J,bi,bj)*SHELFICElatentHeat
462     #endif /* ALLOW_SHIFWFLX_CONTROL */
463     C upward heat flux into the shelf-ice implies basal melting,
464     C thus a downward (negative upward) fresh water flux (as a mass flux),
465     C and vice versa
466     shelfIceFreshWaterFlux(I,J,bi,bj) =
467     & - shelfIceHeatFlux(I,J,bi,bj)
468     & *recip_latentHeat
469     C-- compute surface tendencies
470     shelficeForcingT(i,j,bi,bj) =
471     & - shelfIceHeatFlux(I,J,bi,bj)
472     & *recip_Cp*mass2rUnit
473     & - cFac * shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit
474     & * ( thetaFreeze - tLoc(I,J) )
475     shelficeForcingS(i,j,bi,bj) =
476     & shelfIceFreshWaterFlux(I,J,bi,bj) * mass2rUnit
477     & * ( cFac*sLoc(I,J) + (1. _d 0-cFac)*convertFW2SaltLoc )
478     C-- stress at the ice/water interface is computed in separate
479     C routines that are called from mom_fluxform/mom_vecinv
480     ELSE
481     shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
482     shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
483     shelficeForcingT (I,J,bi,bj) = 0. _d 0
484     shelficeForcingS (I,J,bi,bj) = 0. _d 0
485     ENDIF
486     ENDDO
487     ENDDO
488     ELSE
489     #else
490     IF ( .TRUE. ) THEN
491     #endif /* ALLOW_ISOMIP_TD */
492     C use BRIOS thermodynamics, following Hellmers PhD thesis:
493     C Hellmer, H., 1989, A two-dimensional model for the thermohaline
494     C circulation under an ice shelf, Reports on Polar Research, No. 60
495     C (in German).
496    
497     DO J = 1, sNy
498     DO I = 1, sNx
499     K = kTopC(I,J,bi,bj)
500     IF ( K .NE. 0 .AND. pLoc(I,J) .GT. 0. _d 0 ) THEN
501     C heat flux into the ice shelf, default is diffusive flux
502     C (Holland and Jenkins, 1999, eq.21)
503     thetaFreeze = a0*sLoc(I,J)+c0+b*pLoc(I,J)
504     fwflxFac = 0. _d 0
505     IF ( tLoc(I,J) .GT. thetaFreeze ) fwflxFac = dFac
506     C a few abbreviations
507     eps1 = rUnit2mass*HeatCapacity_Cp
508     & *shiTransCoeffT(i,j,bi,bj)
509     eps2 = rUnit2mass*SHELFICElatentHeat
510     & *shiTransCoeffS(i,j,bi,bj)
511     eps5 = rUnit2mass*HeatCapacity_Cp
512     & *shiTransCoeffS(i,j,bi,bj)
513    
514     C solve quadratic equation for salinity at shelfice-ocean interface
515     C note: this part of the code is not very intuitive as it involves
516     C many arbitrary abbreviations that were introduced to derive the
517     C correct form of the quadratic equation for salinity. The abbreviations
518     C only make sense in connection with my notes on this (M.Losch)
519     C
520     C eps3a was introduced as a constant variant of eps3 to avoid AD of
521     C code of typ (pLoc-const)/pLoc
522     eps3a = rhoShelfIce*SHELFICEheatCapacity_Cp
523     & * SHELFICEkappa * ( 1. _d 0 - dFac )
524     eps3 = eps3a/pLoc(I,J)
525     eps4 = b*pLoc(I,J) + c0
526     eps6 = eps4 - tLoc(I,J)
527     eps7 = eps4 - SHELFICEthetaSurface
528     eps8 = rUnit2mass*SHELFICEheatCapacity_Cp
529     & *shiTransCoeffS(i,j,bi,bj) * fwflxFac
530     aqe = a0 *(eps1+eps3-eps8)
531     recip_aqe = 0. _d 0
532     IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe
533     c bqe = eps1*eps6 + eps3*eps7 - eps2
534     bqe = eps1*eps6
535     & + eps3a*( b
536     & + ( c0 - SHELFICEthetaSurface )/pLoc(I,J) )
537     & - eps2
538     & + eps8*( a0*sLoc(I,J) - eps7 )
539     cqe = ( eps2 + eps8*eps7 )*sLoc(I,J)
540     discrim = bqe*bqe - 4. _d 0*aqe*cqe
541     #undef ALLOW_SHELFICE_DEBUG
542     #ifdef ALLOW_SHELFICE_DEBUG
543     IF ( discrim .LT. 0. _d 0 ) THEN
544     print *, 'ml-shelfice: discrim = ', discrim,aqe,bqe,cqe
545     print *, 'ml-shelfice: pLoc = ', pLoc(I,J)
546     print *, 'ml-shelfice: tLoc = ', tLoc(I,J)
547     print *, 'ml-shelfice: sLoc = ', sLoc(I,J)
548     print *, 'ml-shelfice: tsurface= ',
549     & SHELFICEthetaSurface
550     print *, 'ml-shelfice: eps1 = ', eps1
551     print *, 'ml-shelfice: eps2 = ', eps2
552     print *, 'ml-shelfice: eps3 = ', eps3
553     print *, 'ml-shelfice: eps4 = ', eps4
554     print *, 'ml-shelfice: eps5 = ', eps5
555     print *, 'ml-shelfice: eps6 = ', eps6
556     print *, 'ml-shelfice: eps7 = ', eps7
557     print *, 'ml-shelfice: eps8 = ', eps8
558     print *, 'ml-shelfice: rU2mass = ', rUnit2mass
559     print *, 'ml-shelfice: rhoIce = ', rhoShelfIce
560     print *, 'ml-shelfice: cFac = ', cFac
561     print *, 'ml-shelfice: Cp_W = ', HeatCapacity_Cp
562     print *, 'ml-shelfice: Cp_I = ',
563     & SHELFICEHeatCapacity_Cp
564     print *, 'ml-shelfice: gammaT = ',
565     & SHELFICEheatTransCoeff
566     print *, 'ml-shelfice: gammaS = ',
567     & SHELFICEsaltTransCoeff
568     print *, 'ml-shelfice: lat.heat= ',
569     & SHELFICElatentHeat
570     STOP 'ABNORMAL END in S/R SHELFICE_THERMODYNAMICS'
571     ENDIF
572     #endif /* ALLOW_SHELFICE_DEBUG */
573     saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
574     IF ( saltFreeze .LT. 0. _d 0 )
575     & saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
576     thetaFreeze = a0*saltFreeze + eps4
577     C-- upward fresh water flux due to melting (in kg/m^2/s)
578     cph change to identical form
579     cph freshWaterFlux = rUnit2mass
580     cph & * shiTransCoeffS(i,j,bi,bj)
581     cph & * ( saltFreeze - sLoc(I,J) ) / saltFreeze
582     freshWaterFlux = rUnit2mass
583     & * shiTransCoeffS(i,j,bi,bj)
584     & * ( 1. _d 0 - sLoc(I,J) / saltFreeze )
585     #ifdef ALLOW_SHIFWFLX_CONTROL
586     & + xx_shifwflx_loc(I,J,bi,bj)
587     #endif /* ALLOW_SHIFWFLX_CONTROL */
588     C-- Calculate the upward heat and fresh water fluxes;
589     C-- MITgcm sign conventions: downward (negative) fresh water flux
590     C-- implies melting and due to upward (positive) heat flux
591     shelfIceHeatFlux(I,J,bi,bj) =
592     & ( eps3
593     & - freshWaterFlux*SHELFICEheatCapacity_Cp*fwflxFac )
594     & * ( thetaFreeze - SHELFICEthetaSurface )
595     & - cFac*freshWaterFlux*( SHELFICElatentHeat
596     & - HeatCapacity_Cp*( thetaFreeze - rFac*tLoc(I,J) ) )
597     shelfIceFreshWaterFlux(I,J,bi,bj) = freshWaterFlux
598     C-- compute surface tendencies
599     shelficeForcingT(i,j,bi,bj) =
600     & ( shiTransCoeffT(i,j,bi,bj)
601     & - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )
602     & * ( thetaFreeze - tLoc(I,J) )
603     shelficeForcingS(i,j,bi,bj) =
604     & ( shiTransCoeffS(i,j,bi,bj)
605     & - cFac*shelfIceFreshWaterFlux(I,J,bi,bj)*mass2rUnit )
606     & * ( saltFreeze - sLoc(I,J) )
607     ELSE
608     shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
609     shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
610     shelficeForcingT (I,J,bi,bj) = 0. _d 0
611     shelficeForcingS (I,J,bi,bj) = 0. _d 0
612     ENDIF
613     ENDDO
614     ENDDO
615     ENDIF
616     C endif (not) useISOMIPTD
617     ENDDO
618     ENDDO
619    
620     IF (SHELFICEMassStepping) THEN
621     CALL SHELFICE_STEP_ICEMASS( myTime, myIter, myThid )
622     ENDIF
623    
624     C-- Calculate new loading anomaly (in case the ice-shelf mass was updated)
625     #ifndef ALLOW_AUTODIFF
626     c IF ( SHELFICEloadAnomalyFile .EQ. ' ' ) THEN
627     DO bj = myByLo(myThid), myByHi(myThid)
628     DO bi = myBxLo(myThid), myBxHi(myThid)
629     DO j = 1-OLy, sNy+OLy
630     DO i = 1-OLx, sNx+OLx
631     shelficeLoadAnomaly(i,j,bi,bj) = gravity
632     & *( shelficeMass(i,j,bi,bj) + rhoConst*Ro_surf(i,j,bi,bj) )
633     ENDDO
634     ENDDO
635     ENDDO
636     ENDDO
637     c ENDIF
638     #endif /* ndef ALLOW_AUTODIFF */
639    
640     #ifdef ALLOW_DIAGNOSTICS
641     IF ( useDiagnostics ) THEN
642     CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx',
643     & 0,1,0,1,1,myThid)
644     CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux, 'SHIhtFlx',
645     & 0,1,0,1,1,myThid)
646     C SHIForcT (Ice shelf forcing for theta [W/m2], >0 increases theta)
647     tmpFac = HeatCapacity_Cp*rUnit2mass
648     CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingT,tmpFac,1,
649     & 'SHIForcT',0,1,0,1,1,myThid)
650     C SHIForcS (Ice shelf forcing for salt [g/m2/s], >0 increases salt)
651     tmpFac = rUnit2mass
652     CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingS,tmpFac,1,
653     & 'SHIForcS',0,1,0,1,1,myThid)
654     C Transfer coefficients
655     CALL DIAGNOSTICS_FILL(shiTransCoeffT,'SHIgammT',
656     & 0,1,0,1,1,myThid)
657     CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS',
658     & 0,1,0,1,1,myThid)
659     C Friction velocity
660     #ifdef SHI_ALLOW_GAMMAFRICT
661     IF ( SHELFICEuseGammaFrict )
662     & CALL DIAGNOSTICS_FILL(uStarDiag,'SHIuStar',0,1,0,1,1,myThid)
663     #endif /* SHI_ALLOW_GAMMAFRICT */
664     ENDIF
665     #endif /* ALLOW_DIAGNOSTICS */
666    
667     #endif /* ALLOW_SHELFICE */
668    
669     RETURN
670    
671     END

  ViewVC Help
Powered by ViewVC 1.1.22