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

Annotation of /MITgcm_contrib/dcarroll/highres_darwin/code/shelfice_thermodynamics.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 C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/shelfice/shelfice_thermodynamics.F,v 1.47 2015/12/17 01:52:05 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    
27     CIGF | ASSUMES
28     C--- | * SHELFICEconserve = true
29     C *=============================================================*
30     C \ev
31    
32     C !USES:
33     IMPLICIT NONE
34    
35     C === Global variables ===
36     #include "SIZE.h"
37     #include "EEPARAMS.h"
38     #include "PARAMS.h"
39     #include "GRID.h"
40     #include "DYNVARS.h"
41     #include "FFIELDS.h"
42     #include "SHELFICE.h"
43     #include "SHELFICE_COST.h"
44     #ifdef ALLOW_AUTODIFF
45     # include "CTRL_SIZE.h"
46     # include "ctrl.h"
47     # include "ctrl_dummy.h"
48     #endif /* ALLOW_AUTODIFF */
49     #ifdef ALLOW_AUTODIFF_TAMC
50     # ifdef SHI_ALLOW_GAMMAFRICT
51     # include "tamc.h"
52     # include "tamc_keys.h"
53     # endif /* SHI_ALLOW_GAMMAFRICT */
54     #endif /* ALLOW_AUTODIFF_TAMC */
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 potential 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 iceFrontCellThickness :: the ratio of the horizontal length
76     C of the ice front in each model grid cell
77     C divided by the grid cell area. The "thickness"
78     C of the colum perpendicular to the front
79     C iceFrontWidth :: the width of the ice front.
80    
81     INTEGER I,J,K,Kp1
82     INTEGER bi,bj
83     INTEGER CURI, CURJ, FRONT_K
84    
85     _RL tLoc
86     _RL sLoc
87     _RL pLoc
88     _RL uLoc
89     _RL vLoc
90     _RL wLoc
91     _RL speedLoc
92    
93     C#ifndef SHI_USTAR_WETPOINT
94     C _RL uLoc(1-olx:snx+olx,1-oly:sny+oly)
95     C _RL vLoc(1-olx:snx+olx,1-oly:sny+oly)
96     C#endif
97     C _RL velSq(1-olx:snx+olx,1-oly:sny+oly)
98    
99     _RL freshWaterFlux
100    
101     _RL ice_bottom_Z_C, seafloor_N
102     _RL wet_top_Z_N, wet_bottom_Z_N
103     _RL iceFrontWetContact_Z_max, iceFrontContact_Z_min
104     _RL iceFrontContact_H
105     _RL iceFrontVertContactFrac, iceFrontCellThickness
106     _RL iceFrontWidth, iceFrontFaceArea
107     _RL thermalConductionDistance, thermalConductionTemp
108     _RL tmpHeatFlux, tmpFWFLX
109     _RL tmpForcingT, tmpForcingS
110     _RL tmpFac, icfgridareaFrac
111     INTEGER SI
112    
113     #ifdef ALLOW_DIAGNOSTICS
114     _RL uStarDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
115     #endif /* ALLOW_DIAGNOSTICS */
116    
117     _RL epsilon_H
118    
119     #ifdef ALLOW_SHIFWFLX_CONTROL
120     _RL xx_shifwflx_loc(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
121     #endif
122    
123     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
124    
125     C-- minimum fraction of a cell adjacent to an ice front that must be
126     C-- wet for exchange to happen
127     epsilon_H = 1. _d -03
128    
129     C-- hard coded for now.
130     thermalConductionDistance = 100.0 _d 0
131     thermalConductionTemp = -20.0 _d 0
132     icfgridareaFrac = 1.0 _d 0
133    
134     C heat flux into the ice shelf, default is diffusive flux
135     C (Holland and Jenkins, 1999, eq.21)
136    
137     DO bj = myByLo(myThid), myByHi(myThid)
138     DO bi = myBxLo(myThid), myBxHi(myThid)
139     DO J = 1-OLy,sNy+OLy
140     DO I = 1-OLx,sNx+OLx
141     shelfIceHeatFlux (I,J,bi,bj) = 0. _d 0
142     shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
143     SHIICFHeatFlux (I,J,bi,bj) = 0. _d 0
144     SHIICFFreshWaterFlux(I,J,bi,bj) = 0. _d 0
145     shelficeForcingT (I,J,bi,bj) = 0. _d 0
146     shelficeForcingS (I,J,bi,bj) = 0. _d 0
147     shelficeForcingTR (I,J,bi,bj) = 0. _d 0
148    
149     DO K = 1, NR
150     iceFrontHeatFlux(I,J,K,bi,bj) = 0. _d 0
151     iceFrontFreshWaterFlux(I,J,K,bi,bj) = 0. _d 0
152     iceFrontForcingT(I,J,K,bi,bj) = 0. _d 0
153     iceFrontForcingS(I,J,K,bi,bj) = 0. _d 0
154     iceFrontForcingTR(I,J,K,bi,bj) = 0. _d 0
155     ENDDO /* K */
156    
157     ENDDO /* I */
158     ENDDO /* J */
159    
160     C 3-D velocity-dependent turbulent transfer coefficients
161    
162     DO J = 1-OLy,sNy+OLy
163     DO I = 1-OLx,sNx+OLx
164     DO K = 1, NR
165    
166     C make local copies of velocity
167     uLoc = 0.5*(uVel(I,J,K,bi,bj)+uVel(I+1,J,K,bi,bj))
168     vLoc = 0.5*(vVel(I,J,K,bi,bj)+vVel(I,J+1,K,bi,bj))
169     wLoc = 0.5*(wVel(I,J,K,bi,bj)+wVel(I,J,K+1,bi,bj))
170     speedLoc = SQRT(uLoc*uLoc + vLoc*vLoc + wLoc*wLoc)
171    
172     #ifdef ALLOW_shiTransCoeff_3d
173     shiTransCoeffT(I,J,K,bi,bj) = 1.1 _d -3 *
174     & abs(speedLoc) + 1. _d -4
175    
176     shiTransCoeffS(I,J,K,bi,bj) = 3.1 _d -5 *
177     & abs(speedLoc) + 5.05 _d -7
178     #endif
179     ENDDO /* K */
180     ENDDO /* I */
181     ENDDO /* J */
182    
183     C-- First ice front then ice shelf. Loop through each i,j point
184     C-- process ice fronts in k, then process ice shelf.
185     DO J = 1, sNy
186     DO I = 1, sNx
187    
188     C-- The K index where the ice front ends (0 if no ice front)
189     FRONT_K = K_icefront(I,J,bi,bj)
190    
191     C-- If there is an ice front at this (I,J) continue
192     IF (FRONT_K .GT. 0) THEN
193    
194     C-- Loop through all depths where the ice front is fround
195     DO K = 1, FRONT_K
196     C-- Loop around the four laterally neighboring cells of the ice front.
197     C-- If any neighboring points has wet volume in contact with the ice
198     C-- front at (I,J) then calculate ice-ocean exchanges.
199     C-- The four laterally neighboring point are at (CURI,CURJ)
200     DO SI = 1,4
201     IF (SI .EQ. 1) THEN
202     C-- Looking to right
203     CURI = I+1
204     CURJ = J
205    
206     iceFrontWidth = dyG(I+1,J,bi,bj)
207    
208     ELSEIF (SI .EQ. 2) THEN
209     C-- Looking to LEFT
210     CURI = I-1
211     CURJ = J
212    
213     iceFrontWidth = dyG(I,J,bi,bj)
214     ELSEIF (SI .EQ. 3) THEN
215     C-- Looking to NORTH
216     CURI = I
217     CURJ = J+1
218    
219     iceFrontWidth = dxG(I,J+1,bi,bj)
220     ELSEIF (SI .EQ. 4) THEN
221     C-- Looking to south
222     CURI = I
223     CURJ = J-1
224    
225     iceFrontWidth = dxG(I,J,bi,bj)
226     endif
227    
228     C-- cell depth describes the average distance
229     C-- perpendicular to the ice front fact
230    
231     iceFrontCellThickness = 0. _d 0
232     IF(iceFrontWidth.NE.0. _d 0)
233     & iceFrontCellThickness = RA(CURI,CURJ,bi,bj)
234     & /iceFrontWidth
235     iceFrontFaceArea = DRF(K)*iceFrontWidth
236    
237     C-- First, make sure the adjacent point has at least some water in it.
238     IF (_hFacC(CURI,CURJ,K,bi,bj) .GT. zeroRL) THEN
239    
240     C-- we need to determine how much of the ice front is in contact with
241     C-- water in the neighboring grid cell at this depth level.
242    
243     C-- 1. Determine the top depth with water in the current cell
244     C-- 2. Determine the top depth with water in the neighbor cell
245     C-- 3. Determine the depth where water gap between (1) and (2).
246     C-- 4. If there is a gap then ice front is in contact with water in
247     C-- the neighboring cell
248    
249     C-- ice_bottom_Z_C: the depth (m) of the bottom of the ice in the
250     C-- current cell. Bounded between rF(K) and rF(K+1).
251     C-- * If the ice extends past the bottom of the cell then
252     C-- ice_bottom_Z_C = rF(K+1)
253     C-- [rF(k) >= ice_bottom_Z_C >= rF(K+1)] (rF is negative)
254     ice_bottom_Z_C = max(rF(K+1),
255     & min(Ro_surf(I,J, bi,bj), rF(K)))
256    
257     C-- wet_top_Z_N: the depth (m) of the bottom of the ice in the
258     C-- neighboring grid. If the neighboring cell has ice in
259     C-- (in the form of a shelf or front) then wet_top_Z_N is
260     C-- the depth of this neighboring ice.
261     C--
262     C-- * If neighbor cell has no ice, then Ro_surf = 0 and
263     C-- wet_top_Z_N = rF(K)
264     C-- [rF(k) >= wet_top_Z_N >= rF(K+1)] (rF is negative)
265    
266     wet_top_Z_N = max(rF(K+1),
267     & min(Ro_surf(CURI,CURJ, bi,bj), rF(K)))
268    
269     C-- wet_bottom_Z_N: the depth (m) of the bottom of the wet part of the
270     C-- neighboring cell. If the seafloor reaches into
271     C-- the grid cell then the bottom of the wet part of the
272     C-- grid cell is at the seafloor.
273     C--
274     C-- * If the seafloor is deeper than this grid cell then
275     C-- wet_bottom_Z = rF(K+1)
276     C-- * If the seafloor is shallower than this grid cell then
277     C-- wet_bottom_Z = rF(K)
278     C-- * If the seafloor reaches partly into this grid cell
279     C-- then wet_bottom_Z = R_low
280    
281     C-- [rF(k) >= wet_bottom_Z >= rF(K+1)] (rF is negative)
282    
283     wet_bottom_Z_N = min(rF(K),
284     & max(R_low(CURI,CURJ, bi,bj), rF(K+1)))
285    
286     C-- iceFrontWetContact_Z_max: The deepest point where the
287     C-- the ice front at (I,J) is in contact with water
288     C-- in the neighboring cell. The shallower of
289     C-- wet_bottom_Z_N (seafloor depth of neighboring point) and
290     C-- ice_bottom_Z_C (bottom of ice front in this center cell).
291    
292     C-- * wet_bottom_Z_N if the seafloor of the neighboring
293     C-- cell is shallower than the ice draft at (I,J).
294     C-- * ice_bottom_Z_C if the ice draft at (I,J) is shallower
295     C-- than the seafloor of the neighboring cell.
296    
297     IF (ice_bottom_Z_C .GT. wet_bottom_Z_N) THEN
298     iceFrontWetContact_Z_max = ice_bottom_Z_C
299     ELSE
300     iceFrontWetContact_Z_max = wet_bottom_Z_N
301     ENDIF
302    
303     C-- The shallowest depth where the ice front at (I,J) is in contact
304     C-- with water in the neighboring cell. If the neighboring cell has
305     C-- no ice draft then wet_top_Z_N = rF(k), the top of the cell.
306     C-- Otherwise, the shallowest depth where the ice front at (I,J) can
307     C-- be in in contact with water (not ice) in (CURI, CURJ)
308     C-- is wet_top_Z_N.
309    
310     C-- the fraction of the grid cell height that has ice draft in contact
311     C-- with water in the neighboring cell.
312     iceFrontVertContactFrac =
313     & (wet_top_Z_N - iceFrontWetContact_Z_max)/ DRF(K)
314    
315    
316     C-- Only proceed if iceFrontVertContactFrac is > 0, the
317     C-- ice draft at (I,J)
318     C-- is in contact with some water in the neighboring grid cell.
319     IF (iceFrontVertContactFrac .GT. epsilon_H) THEN
320     tLoc = theta(CURI,CURJ,K,bi,bj)
321     sLoc = MAX(salt(CURI,CURJ,K,bi,bj), zeroRL)
322    
323     C-- use pressure at the halfway point between the top and bottom of
324     C-- points of the ice front where the ice front is in contact with
325     C-- open water.
326     pLoc = 0.5 _d 0 * ABS(wet_top_Z_N +
327     & iceFrontWetContact_Z_max)
328    
329     CALL SHELFICE_SOLVE4FLUXES(
330     I tLoc, sLoc, pLoc,
331     #ifndef ALLOW_shiTransCoeff_3d
332     I shiTransCoeffT(CURI,CURJ,bi,bj),
333     I shiTransCoeffS(CURI,CURJ,bi,bj),
334     #else
335     I shiTransCoeffT(CURI,CURJ,K,bi,bj),
336     I shiTransCoeffS(CURI,CURJ,K,bi,bj),
337     #endif
338     I thermalConductionDistance,
339     I thermalConductionTemp,
340     O tmpHeatFlux, tmpFWFLX,
341     O tmpForcingT, tmpForcingS,
342     I bi, bj, myTime, myIter, myThid )
343    
344     C-- fluxes and forcing must be scaled by iceFrontVertContactFract and
345     C-- iceFrontContactFrac some fraction of the heigth and width of the
346     C-- grid cell face may not ice in contact with water.
347    
348     C tmpHeatFlux and tmpFWFLX come as W/m^2 and kg/m^2/s respectively
349     C-- but these rates only apply to the
350     C-- fraction of the grid cell that has ice in contact with seawater.
351     C-- we must scale by iceFrontVertContactFrac to get to the average
352     C-- fluxes in this grid cell.
353    
354     C-- In units W/m^2
355     iceFrontHeatFlux(CURI,CURJ,K,bi,bj) =
356     & iceFrontHeatFlux(CURI,CURJ,K,bi,bj) +
357     & tmpHeatFlux*iceFrontVertContactFrac
358    
359     C In units of kg/m^2/s
360     iceFrontFreshWaterFlux(CURI,CURJ,K,bi,bj) =
361     & iceFrontFreshWaterFlux(CURI,CURJ,K,bi,bj) +
362     & tmpFWFLX*iceFrontVertContactFrac
363    
364     iceFrontForcingTR(CURI,CURJ,K,bi,bj) =
365     & iceFrontFreshWaterFlux(CURI,CURJ,K,bi,bj) *
366     & mass2rUnit
367    
368     C ow - 06/29/2018
369     C ow - Verticallly sum up the 3D icefront heat and freshwater fluxes to
370     C ow - compute the total flux for the water column. The shelfice fluxes,
371     C ow - which are 2D, will be added later. NOTE that only
372     C ow - ice-front melts below shelf-ice are included to be consistent
373     C ow - with Rignot's data
374     if(k.GE.kTopC(I,J,bi,bj))then
375     if(RA(CURI,CURJ,bi,bj).NE.0. _d 0)then
376     icfgridareaFrac =
377     & iceFrontFaceArea/RA(CURI,CURJ,bi,bj)
378     SHIICFHeatFlux(CURI,CURJ,bi,bj) =
379     & SHIICFHeatFlux(CURI,CURJ,bi,bj) +
380     & iceFrontHeatFlux(CURI,CURJ,K,bi,bj)
381     & * icfgridareaFrac
382     SHIICFFreshWaterFlux(CURI,CURJ,bi,bj) =
383     & SHIICFFreshWaterFlux(CURI,CURJ,bi,bj) +
384     & iceFrontFreshWaterFlux(CURI,CURJ,K,bi,bj)
385     & * icfgridareaFrac
386     endif
387     endif
388     C iceFrontForcing[T,S] X m/s but these rates only apply to the
389     C-- fraction of the grid cell that has ice in contact with seawater.
390     C-- we must scale by iceFrontVertContactFrac to get to the average
391     C-- fluxes in this grid cell. We must also divide the by the length
392     C-- of the grid cell perpendicular to the face.
393    
394     IF (iceFrontCellThickness .NE. 0. _d 0) THEN
395     C In units of K / s
396     iceFrontForcingT(CURI,CURJ,K,bi,bj) =
397     & iceFrontForcingT(CURI,CURJ,K,bi,bj) +
398     & tmpForcingT/iceFrontCellThickness*
399     & iceFrontVertContactFrac
400    
401     C In units of psu /s
402     iceFrontForcingS(CURI,CURJ,K,bi,bj) =
403     & iceFrontForcingS(CURI,CURJ,K,bi,bj) +
404     & tmpForcingS/iceFrontCellThickness*
405     & iceFrontVertContactFrac
406    
407     ENDIF /* iceFrontCellThickness */
408     C In units of kg /s
409     addMass(CURI,CURJ,K,bi,bj) =
410     & addMass(CURI,CURJ,K,bi,bj) -
411     & tmpFWFLX*iceFrontFaceArea*
412     & iceFrontVertContactFrac
413     ENDIF /* iceFrontVertContactFrac */
414     ENDIF /* hFacC(CURI,CURJ,K,bi,bj) */
415     ENDDO /* SI loop for adjacent cells */
416     ENDDO /* K LOOP */
417     ENDIF /* FRONT K */
418    
419     C-- ice shelf
420     K = kTopC(I,J,bi,bj)
421    
422     C-- If there is an ice front at this (I,J) continue
423     C-- I am assuming K is only .GT. when there is at least some
424     C-- nonzero wet point below the shelf in the grid cell.
425     IF (K .GT. 0) THEN
426     C-- Initialize these values to zero
427     pLoc = 0 _d 0
428     tLoc = 0 _d 0
429     sLoc = 0 _d 0
430    
431     C-- make local copies of temperature, salinity and depth
432     C-- (pressure in deci-bar) underneath the ice
433     C-- for the ice shelf case we use hydrostatic pressure at the ice
434     C-- base of the ice shelf, top of the cavity.
435    
436     pLoc = ABS(R_shelfIce(I,J,bi,bj))
437     tLoc = theta(I,J,K,bi,bj)
438     sLoc = MAX(salt(I,J,K,bi,bj), zeroRL)
439    
440     CALL SHELFICE_SOLVE4FLUXES(
441     I tLoc, sLoc, pLoc,
442     #ifndef ALLOW_shiTransCoeff_3d
443     I shiTransCoeffT(I,J,bi,bj),
444     I shiTransCoeffS(I,J,bi,bj),
445     #else
446     I shiTransCoeffT(I,J,K,bi,bj),
447     I shiTransCoeffS(I,J,K,bi,bj),
448     #endif
449     I pLoc, thermalConductionTemp,
450     O tmpHeatFlux, tmpFWFLX,
451     O tmpForcingT, tmpForcingS,
452     I bi, bj, myTime, myIter, myThid )
453    
454     C In units of W/m^2
455     shelficeHeatFlux(I,J,bi,bj) = tmpHeatFlux
456     C In units of kg/m^2/s
457     shelfIceFreshWaterFlux(I,J,bi,bj) = tmpFWFLX
458    
459     shelficeForcingTR(I,J,bi,bj) =
460     & shelfIceFreshWaterFlux(I,J,bi,bj) * mass2rUnit
461    
462     C ow - 06/29/2018
463     C ow - Now add shelfice heat and freshwater fluxes
464     SHIICFHeatFlux(i,j,bi,bj) =
465     & SHIICFHeatFlux(i,j,bi,bj) +
466     & shelficeHeatFlux(i,j,bi,bj)
467     SHIICFFreshWaterFlux(i,j,bi,bj) =
468     & SHIICFFreshWaterFlux(i,j,bi,bj) +
469     & shelfIceFreshWaterFlux(i,j,bi,bj)
470     C In units of K/s -- division by drF required first
471     shelficeForcingT(I,J,bi,bj) = tmpForcingT*
472     & recip_drF(K)* _recip_hFacC(i,j,K,bi,bj)
473     C In units of psu/s -- division by drF required first
474     shelficeForcingS(I,J,bi,bj) = tmpForcingS*
475     & recip_drF(K)* _recip_hFacC(i,j,K,bi,bj)
476     C In units of kg/s -- multiplication of area required first
477     addMass(I,J,K, bi,bj) = addMass(I,J,K, bi,bj) -
478     & tmpFWFLX*RA(I,J,bi,bj)
479     ENDIF /* SHELF K > 0 */
480     ENDDO /* i */
481     ENDDO /* j */
482     ENDDO /* bi */
483     ENDDO /* bj */
484    
485    
486     C-- Calculate new loading anomaly (in case the ice-shelf mass was updated)
487     #ifndef ALLOW_AUTODIFF
488     c IF ( SHELFICEloadAnomalyFile .EQ. ' ' ) THEN
489     DO bj = myByLo(myThid), myByHi(myThid)
490     DO bi = myBxLo(myThid), myBxHi(myThid)
491     DO j = 1-OLy, sNy+OLy
492     DO i = 1-OLx, sNx+OLx
493     shelficeLoadAnomaly(i,j,bi,bj) = gravity
494     & *( shelficeMass(i,j,bi,bj) + rhoConst*Ro_surf(i,j,bi,bj) )
495     ENDDO
496     ENDDO
497     ENDDO
498     ENDDO
499     c ENDIF
500     #endif /* ndef ALLOW_AUTODIFF */
501    
502     #ifdef ALLOW_DIAGNOSTICS
503     IF ( useDiagnostics ) THEN
504     CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx',
505     & 0,1,0,1,1,myThid)
506     CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux, 'SHIhtFlx',
507     & 0,1,0,1,1,myThid)
508    
509     CALL DIAGNOSTICS_FILL_RS(SHIICFFreshWaterFlux,'SHIICFfwFlx',
510     & 0,1,0,1,1,myThid)
511     CALL DIAGNOSTICS_FILL_RS(SHIICFHeatFlux, 'SHIICFhtFlx',
512     & 0,1,0,1,1,myThid)
513    
514     CALL DIAGNOSTICS_FILL(iceFrontFreshWaterFlux, 'ICFfwFlx',
515     & 0,Nr,0,1,1,myThid)
516     CALL DIAGNOSTICS_FILL(iceFrontHeatFlux, 'ICFhtFlx',
517     & 0,Nr,0,1,1,myThid)
518    
519     CALL DIAGNOSTICS_FILL(shelfIceForcingTR, 'SHITR ',
520     & 0,Nr,0,1,1,myThid)
521    
522     CALL DIAGNOSTICS_FILL(iceFrontForcingTR, 'ICFTR ',
523     & 0,Nr,0,1,1,myThid)
524    
525     C SHIForcT (Ice shelf forcing for theta [W/m2], >0 increases theta)
526     tmpFac = HeatCapacity_Cp*rUnit2mass
527     CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingT,tmpFac,1,
528     & 'SHIForcT',0,1,0,1,1,myThid)
529     C SHIForcS (Ice shelf forcing for salt [g/m2/s], >0 increases salt)
530     tmpFac = rUnit2mass
531     CALL DIAGNOSTICS_SCALE_FILL(shelficeForcingS,tmpFac,1,
532     & 'SHIForcS',0,1,0,1,1,myThid)
533    
534     C ICFForcT (Ice front forcing for theta [W/m2], >0 increases theta)
535     tmpFac = HeatCapacity_Cp*rUnit2mass
536     CALL DIAGNOSTICS_SCALE_FILL(iceFrontForcingT,tmpFac,1,
537     & 'ICFForcT',0,Nr,0,1,1,myThid)
538     C ICFForcS (Ice front forcing for salt [g/m2/s], >0 increases salt)
539     tmpFac = rUnit2mass
540     CALL DIAGNOSTICS_SCALE_FILL(iceFrontForcingS,tmpFac,1,
541     & 'ICFForcS',0,Nr,0,1,1,myThid)
542    
543     C Transfer coefficients
544     #ifndef ALLOW_shiTransCoeff_3d
545     CALL DIAGNOSTICS_FILL(shiTransCoeffT,'SHIgammT',
546     & 0,1,0,1,1,myThid)
547     CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS',
548     & 0,1,0,1,1,myThid)
549     #else
550     CALL DIAGNOSTICS_FILL(shiTransCoeffT,'SHIgammT',
551     & 0,Nr,0,1,1,myThid)
552     CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS',
553     & 0,Nr,0,1,1,myThid)
554     #endif
555     C Friction velocity
556     #ifdef SHI_ALLOW_GAMMAFRICT
557     IF ( SHELFICEuseGammaFrict )
558     & CALL DIAGNOSTICS_FILL(uStarDiag,'SHIuStar',0,1,0,1,1,myThid)
559     #endif /* SHI_ALLOW_GAMMAFRICT */
560     ENDIF
561     #endif
562    
563     #endif /* ALLOW_SHELFICE */
564     RETURN
565     END

  ViewVC Help
Powered by ViewVC 1.1.22