/[MITgcm]/MITgcm_contrib/dcarroll/highres_darwin/code/shelfice_solve4fluxes.F
ViewVC logotype

Annotation of /MITgcm_contrib/dcarroll/highres_darwin/code/shelfice_solve4fluxes.F

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


Revision 1.1 - (hide annotations) (download)
Sun Sep 22 21:23:47 2019 UTC (5 years, 10 months ago) by dcarroll
Branch: MAIN
CVS Tags: HEAD
Initial check in of high resolution Darwin simulation code

1 dcarroll 1.1 #include "SHELFICE_OPTIONS.h"
2     #ifdef ALLOW_AUTODIFF
3     # include "AUTODIFF_OPTIONS.h"
4     #endif
5     #ifdef ALLOW_CTRL
6     # include "CTRL_OPTIONS.h"
7     #endif
8    
9     CBOP
10     C !ROUTINE: SHELFICE_SOLVE4FLUXES
11     C !INTERFACE:
12     SUBROUTINE SHELFICE_SOLVE4FLUXES(
13     I tLoc, sLoc, pLoc,
14     I gammaT, gammaS,
15     I iceConductionDistance, thetaIceConduction,
16     O heatFlux, fwFlux,
17     O forcingT, forcingS,
18     I bi, bj, myTime, myIter, myThid )
19    
20     C !DESCRIPTION: \bv
21     C *==========================================================*
22     C | SUBROUTINE SOLVE4FLUXES
23     C | o Calculate
24     C *==========================================================*
25     C \ev
26    
27     C !USES:
28     IMPLICIT NONE
29    
30     C === Global variables ===
31     #include "SIZE.h"
32     #include "EEPARAMS.h"
33     #include "PARAMS.h"
34     #include "GRID.h"
35     #include "DYNVARS.h"
36     #include "FFIELDS.h"
37     #include "SHELFICE.h"
38     #include "SHELFICE_COST.h"
39     #ifdef ALLOW_AUTODIFF
40     # include "CTRL_SIZE.h"
41     # include "ctrl.h"
42     # include "ctrl_dummy.h"
43     #endif /* ALLOW_AUTODIFF */
44     #ifdef ALLOW_AUTODIFF_TAMC
45     # ifdef SHI_ALLOW_GAMMAFRICT
46     # include "tamc.h"
47     # include "tamc_keys.h"
48     # endif /* SHI_ALLOW_GAMMAFRICT */
49     #endif /* ALLOW_AUTODIFF_TAMC */
50    
51     C !INPUT PARAMETERS:
52     C tLoc ::
53     C sLoc ::
54     C pLoc ::
55     C insitutLoc ::
56     C gammaT ::
57     C gammaS ::
58     C bi,bj :: tile indices
59     C myTime :: current time in simulation
60     C myIter :: iteration number in simulation
61     C myThid :: my Thread Id number
62    
63     C !OUTPUT PARAMETERS:
64     C heatFlux ::
65     C fwFlux ::
66     C forcingT ::
67     C forcingS ::
68     C----------
69    
70     _RL tLoc, sLoc, pLoc, insitutLoc
71     _RL gammaT, gammaS
72     _RL iceConductionDistance, thetaIceConduction
73     _RL heatFlux, fwFlux, forcingS, forcingT
74     INTEGER i, j, bi, bj
75     _RL myTime
76     INTEGER myIter, myThid
77     character*200 msgBuf
78     CEOP
79    
80     #ifndef ALLOW_OPENAD
81     _RL SW_TEMP
82     EXTERNAL SW_TEMP
83     #endif
84    
85    
86     C !LOCAL VARIABLES:
87     C === Local variables ===
88     _RL thetaFreeze, saltFreeze
89     _RL eps1, eps2, eps3, eps4, eps5, eps6, eps7, eps8
90     _RL aqe, bqe, cqe, discrim, recip_aqe
91     _RL a0, a1, a2, b0, c0
92     _RL w_I
93    
94     C === Useful Units ===
95     C-- gammaT, m s^-1
96     C-- gammaS, m s^-1
97     C-- rUnit2mass (rhoConst), kg m^-3
98     C-- mass2rUnit (recip_rhoConst), m^3 kg^-1
99     C-- eps3, W K^-1 m^-2
100     C-- fwFlux, kg m^-2 s^-1
101     C-- heatFlux, kg m^-2 s^-1
102     C-- forcing T, K m/s
103     C-- forcing S, psu m/s
104     C-- SHELFICEkappa, m^2/s
105    
106     C-- eps1, W K^-1 m^-2 : kg/m^3 * J/kg/K * m/s
107     C-- eps2, W m^-2 : kg/m^3 * J/kg * m/s
108     C-- eps3, W K^-1 m^-2 : kg m^-3 * J kg^-1 K^-1 * m^2 s^-1 * m^-1
109    
110     C-- fwFlux : fresh water flux due to melting (kg m^-2 s^-1)
111    
112     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
113     C linear dependence of freezing point on salinity
114     a0 = -0.0575 _d 0
115     c0 = 0.0901 _d 0
116     b0 = -7.61 _d -4
117    
118     C-- convert potential T into in-situ T relative to surface
119     insitutLoc = SW_TEMP(sLoc,tLoc,pLoc, zeroRL)
120    
121     C-- DEFINE SOME CONSTANTS
122     eps1 = rUnit2mass*HeatCapacity_Cp * gammaT
123     eps2 = rUnit2mass*SHELFICElatentHeat * gammaS
124     eps3 = rhoShelfIce*SHELFICEheatCapacity_Cp * SHELFICEkappa
125     & /iceConductionDistance;
126     eps4 = b0*pLoc + c0
127     eps6 = eps4 - insitutLoc
128     eps7 = eps4 - thetaIceConduction
129    
130     IF ( debugLevel.GE.debLevE ) THEN
131     WRITE(msgBuf,'(A25,9E16.8)')
132     &'ZZZ7 r2mass, Cp, gmaT,SIlh,gmaS, rhoSI, SI_Cp, Kap, ICondDis ',
133     & rUnit2mass,HeatCapacity_Cp,gammaT, SHELFICElatentHeat,gammaS,
134     & rhoShelfIce, SHELFICEheatCapacity_Cp, SHELFICEkappa,
135     & iceConductionDistance
136    
137     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
138     & SQUEEZE_RIGHT , 1)
139     ENDIF
140    
141     C-- Default thermodynamics specify a linear T gradient
142     C-- through the ice (Holland and Jenkins, 1999, Section 2.
143     aqe = a0*(eps1 + eps3)
144     bqe = eps1*eps6 + eps3*eps7 - eps2 - SHELFICEsalinity*aqe
145     cqe = eps2 * sLoc - SHELFICEsalinity*(eps1*eps6 + eps3*eps7)
146    
147     C-- Alterantively, we can have a nonlinear T gradient
148     C-- through the ice (Holland and Jenkins, 1999, Section 3.
149     C-- This demands a different set of constants
150     IF (SHELFICEadvDiffHeatFlux .EQV. .TRUE.) THEN
151    
152     IF ( debugLevel.GE.debLevE ) THEN
153     WRITE(msgBuf,'(A)')
154     & '1useSHELFICEadvDiffHeatFlux '
155    
156     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
157     & SQUEEZE_RIGHT , 1)
158     ENDIF
159    
160     eps8 = rUnit2mass * gammaS * SHELFICEheatCapacity_Cp
161     aqe = a0 *(eps1 - eps8)
162     bqe = eps1*eps6 + sLoc*eps8*a0 - eps8*eps7 - eps2 -
163     & SHELFICEsalinity*eps1*a0
164     cqe = sLoc*(eps8*eps7 + eps2) - SHELFICEsalinity*eps1
165     ENDIF
166    
167     C solve quadratic equation for salinity at shelfice-ocean interface
168     recip_aqe = 0. _d 0
169     IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe
170    
171     C-- Sb = \frac{-bqe \pm \sqrt(bqe^2 - 4 aqc cqe)}{2 aqe}
172     discrim = bqe*bqe - 4. _d 0*aqe*cqe
173    
174     C-- Try the negative root (- SQRT(discrim)) of the quadratic eq.
175     saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
176    
177     C--- If the negative root yields a negative salinity, then use the
178     C-- positive root (+ SQRT(discrim))
179     IF ( saltFreeze .LT. 0. _d 0 ) THEN
180     saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
181     ENDIF
182    
183     C-- in situ seawater freezing point using linearization
184     thetaFreeze = a0*saltFreeze + eps4
185    
186     C-- Calculate the upward heat and fresh water fluxes;
187     C-- MITgcm sign conventions: downward (negative) fresh water flux
188     C-- implies melting and due to upward (positive) heat flux
189    
190     C-- This formulation of fwflux, derived from the heat balance equation
191     C-- instead of the salt balance equation, can handle the case when the
192     C-- salinity of the ocean, boundary layer, and ice are identical.
193     fwFlux = 1/SHELFICElatentHeat*(
194     & eps3*(thetaFreeze - thetaIceConduction) -
195     & eps1*(insitutLoc - thetaFreeze) )
196    
197     C-- If a nonlinear local ice T gradient near the ice-ocean interface
198     C-- is allowed and fwflux is positive (ice growth) then
199     C-- we must solve the quadratic equation using a different set of
200     C-- coeffs (Holland Jenkins, 1999).
201     IF ((SHELFICEadvDiffHeatFlux .EQV. .TRUE.) .AND.
202     & (fwFlux .GT. zeroRL)) THEN
203     IF ( debugLevel.GE.debLevE ) THEN
204     WRITE(msgBuf,'(A)')
205     & '2useSHELFICEadvDiffHeatFlux '
206    
207     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
208     & SQUEEZE_RIGHT , 1)
209     ENDIF
210     aqe = a0 *(eps1)
211     bqe = eps1*eps6 - eps2 - SHELFICEsalinity*eps1*a0
212     cqe = sLoc*(eps2) - SHELFICEsalinity*eps1
213    
214     recip_aqe = 0. _d 0
215     IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe
216    
217     discrim = bqe*bqe - 4. _d 0*aqe*cqe
218     saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
219     IF ( saltFreeze .LT. 0. _d 0 ) THEN
220     saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
221     ENDIF
222    
223     thetaFreeze = a0*saltFreeze + eps4
224    
225     fwFlux = 1/SHELFICElatentHeat* (
226     & eps3*(thetaFreeze - thetaIceConduction) -
227     & eps1*(insitutLoc - thetaFreeze) )
228     ENDIF
229    
230     C meltwater advective flux velocity at ice-ocean interface (m/s)
231     C * negative corresponds to downward flux (melting)
232     w_I = fwFlux * mass2rUnit
233    
234     C-- Calculate the upward heat fluxes:
235     C-- melting requires upward (positive) heat flux from ocean to ice.
236    
237     C-- The heatFlux variable corresponds with the change of energy in the
238     C-- ocean grid cell volume. In the conservative case (J2001),
239     C-- advective heat fluxes change the energy of the volume whereas in
240     C-- the non-conservative case there are no advective heat fluxes
241     C-- melting or freezing have no associated advective heat fluxes.
242     IF (SHELFICEconserve) THEN
243    
244     C-- In the conservative case (J2001) there are two cases, fixed and
245     C-- non-fixed ocean volume.
246    
247     IF (useRealFreshWaterFlux ) THEN
248    
249     C-- If the ocean volume can change (realFWFlux=true) then advection of
250     C-- meltwater does not displace water at T=insitutLoc in the cell and the
251     C-- heat flux correpsonding to the total energy flux of the volume
252     C-- consists of only two terms: turbulent fluxes (positive out)
253     C-- and advective meltwater fluxes (negative in).
254     heatFlux = rUnit2mass*HeatCapacity_Cp * (
255     & gammaT * (insitutLoc - thetaFreeze)
256     & + w_I * thetaFreeze )
257     ELSE
258    
259     C-- If the volume is fixed (realFWFlux=false) then the advection of
260     C-- meltwater does displace ambient water at T=insitutLoc in the cell.
261     C-- Displacement reduction volume energy by w_I * insitutLoc (positive)
262     heatFlux = rUnit2mass*HeatCapacity_Cp * (
263     & gammaT * (insitutLoc - thetaFreeze)
264     & + w_I * (thetaFreeze - insitutLoc) )
265     ENDIF
266    
267     ELSE
268    
269     C-- In the non-conservative form, only fluxes are turbulent fluxes
270     heatFlux = rUnit2mass*HeatCapacity_Cp *
271     & gammaT * (insitutLoc - thetaFreeze)
272     ENDIF
273    
274     C-- Calculate the T and S tendency terms. T tendency term is
275     C-- not necessarily proportional to the heat flux term above because
276     C-- the heat flux term corresponds to total energy change in the grid
277     C-- cell and not the change of energy per unit volume.
278    
279     IF (SHELFICEconserve) THEN
280     C-- In the conservative case, meltwater advection contributes (J2001)
281     C-- to T and S tendencies
282     C-- * forcing T (K m/s)
283     forcingT =
284     & (gammaT - w_I)*(thetaFreeze - insitutLoc)
285    
286     C-- * forcing S (psu m/s)
287     forcingS =
288     & (gammaS - w_I)*(saltFreeze - sLoc)
289     ELSE
290     C-- Otherwise, the only fluxes out of the ocean that change T and S
291     C-- are the turbulent fluxes.
292     forcingT = gammaT * ( thetaFreeze - insitutLoc )
293     forcingS = gammaS * ( saltFreeze - sLoc )
294     ENDIF
295    
296     IF ( debugLevel.GE.debLevE ) THEN
297     WRITE(msgBuf,'(A25,7E16.8)')
298     & 'ZZZ6 aqe, bqe, ceq ',
299     & aqe,bqe,cqe
300    
301     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
302     & SQUEEZE_RIGHT , 1)
303    
304     WRITE(msgBuf,'(A25,7E16.8)')
305     & 'ZZZ2 T,S,P,t,TFrz,SFrz,w_I ',
306     & tLoc,sLoc,pLoc, insitutLoc, thetaFreeze, saltFreeze, w_I
307    
308     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
309     & SQUEEZE_RIGHT , 1)
310     ENDIF
311    
312    
313     RETURN
314     END
315    

  ViewVC Help
Powered by ViewVC 1.1.22