/[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.15 - (hide annotations) (download)
Mon Apr 11 16:02:53 2016 UTC (9 years, 3 months ago) by dgoldberg
Branch: MAIN
Changes since 1.14: +105 -55 lines
turn off all fluxes under grounded ice

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

  ViewVC Help
Powered by ViewVC 1.1.22