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