/[MITgcm]/MITgcm_contrib/torge/itd/code/seaice_growth.F
ViewVC logotype

Annotation of /MITgcm_contrib/torge/itd/code/seaice_growth.F

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

Revision 1.1 - (hide annotations) (download)
Fri Apr 27 22:22:17 2012 UTC (13 years, 3 months ago) by dimitri
Branch: MAIN 
check-in original code, before itd modifications
seaice_advdiff.F,v 1.60 2012/02/16 01:22:02
seaice_check_pickup.F,v 1.7 2012/03/05 15:21:44
seaice_diagnostics_init.F,v 1.33 2012/02/16 01:22:02
seaice_growth.F,v 1.162 2012/03/15 03:07:31
seaice_init_fixed.F,v 1.19 2012/03/11 13:41:38
seaice_init_varia.F,v 1.72 2012/03/14 22:55:53
seaice_readparms.F,v 1.120 2012/03/14 22:55:53
seaice_write_pickup.F,v 1.14 2012/03/05 15:21:45
seaice_read_pickup.F,v 1.16 2012/03/05 15:21:44
seaice_model.F,v 1.100 2012/03/02 18:56:06
SEAICE.h,v 1.62 2012/03/06 16:51:21
SEAICE_OPTIONS.h,v 1.63 2012/03/08 01:15:02
SEAICE_PARAMS.h,v 1.91 2012/03/11 13:41:38
SEAICE_SIZE.h,v 1.5 2012/03/06 16:51:21
SIZE.h,v 1.28 2009/05/17 21:15:07
1 dimitri 1.1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_growth.F,v 1.162 2012/03/15 03:07:31 jmc Exp $
2     C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5     #ifdef ALLOW_EXF
6     # include "EXF_OPTIONS.h"
7     #endif
8    
9     CBOP
10     C !ROUTINE: SEAICE_GROWTH
11     C !INTERFACE:
12     SUBROUTINE SEAICE_GROWTH( myTime, myIter, myThid )
13     C !DESCRIPTION: \bv
14     C *==========================================================*
15     C | SUBROUTINE seaice_growth
16     C | o Updata ice thickness and snow depth
17     C *==========================================================*
18     C \ev
19    
20     C !USES:
21     IMPLICIT NONE
22     C === Global variables ===
23     #include "SIZE.h"
24     #include "EEPARAMS.h"
25     #include "PARAMS.h"
26     #include "DYNVARS.h"
27     #include "GRID.h"
28     #include "FFIELDS.h"
29     #include "SEAICE_SIZE.h"
30     #include "SEAICE_PARAMS.h"
31     #include "SEAICE.h"
32     #include "SEAICE_TRACER.h"
33     #ifdef ALLOW_EXF
34     # include "EXF_PARAM.h"
35     # include "EXF_FIELDS.h"
36     #endif
37     #ifdef ALLOW_SALT_PLUME
38     # include "SALT_PLUME.h"
39     #endif
40     #ifdef ALLOW_AUTODIFF_TAMC
41     # include "tamc.h"
42     #endif
43    
44     C !INPUT/OUTPUT PARAMETERS:
45     C === Routine arguments ===
46     C myTime :: Simulation time
47     C myIter :: Simulation timestep number
48     C myThid :: Thread no. that called this routine.
49     _RL myTime
50     INTEGER myIter, myThid
51     CEOP
52    
53     C !FUNCTIONS:
54     #ifdef ALLOW_DIAGNOSTICS
55     LOGICAL DIAGNOSTICS_IS_ON
56     EXTERNAL DIAGNOSTICS_IS_ON
57     #endif
58    
59     C !LOCAL VARIABLES:
60     C === Local variables ===
61     C
62     C unit/sign convention:
63     C Within the thermodynamic computation all stocks, except HSNOW,
64     C are in 'effective ice meters' units, and >0 implies more ice.
65     C This holds for stocks due to ocean and atmosphere heat,
66     C at the outset of 'PART 2: determine heat fluxes/stocks'
67     C and until 'PART 7: determine ocean model forcing'
68     C This strategy minimizes the need for multiplications/divisions
69     C by ice fraction, heat capacity, etc. The only conversions that
70     C occurs are for the HSNOW (in effective snow meters) and
71     C PRECIP (fresh water m/s).
72     C
73     C HEFF is effective Hice thickness (m3/m2)
74     C HSNOW is Heffective snow thickness (m3/m2)
75     C HSALT is Heffective salt content (g/m2)
76     C AREA is the seaice cover fraction (0<=AREA<=1)
77     C Q denotes heat stocks -- converted to ice stocks (m3/m2) early on
78     C
79     C For all other stocks/increments, such as d_HEFFbyATMonOCN
80     C or a_QbyATM_cover, the naming convention is as follows:
81     C The prefix 'a_' means available, the prefix 'd_' means delta
82     C (i.e. increment), and the prefix 'r_' means residual.
83     C The suffix '_cover' denotes a value for the ice covered fraction
84     C of the grid cell, whereas '_open' is for the open water fraction.
85     C The main part of the name states what ice/snow stock is concerned
86     C (e.g. QbyATM or HEFF), and how it is affected (e.g. d_HEFFbyATMonOCN
87     C is the increment of HEFF due to the ATMosphere extracting heat from the
88     C OCeaN surface, or providing heat to the OCeaN surface).
89    
90     C i,j,bi,bj :: Loop counters
91     INTEGER i, j, bi, bj
92     C number of surface interface layer
93     INTEGER kSurface
94     C constants
95     _RL tempFrz, ICE2SNOW, SNOW2ICE
96     _RL QI, QS, recip_QI
97    
98     C-- TmixLoc :: ocean surface/mixed-layer temperature (in K)
99     _RL TmixLoc (1:sNx,1:sNy)
100    
101     C a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
102     C the atmosphere and the ocean surface - for ice covered water
103     C a_QbyATM_open :: same but for open water
104     C r_QbyATM_cover :: residual of a_QbyATM_cover after freezing/melting processes
105     C r_QbyATM_open :: same but for open water
106     _RL a_QbyATM_cover (1:sNx,1:sNy)
107     _RL a_QbyATM_open (1:sNx,1:sNy)
108     _RL r_QbyATM_cover (1:sNx,1:sNy)
109     _RL r_QbyATM_open (1:sNx,1:sNy)
110     C a_QSWbyATM_open - short wave heat flux over ocean in W/m^2
111     C a_QSWbyATM_cover - short wave heat flux under ice in W/m^2
112     _RL a_QSWbyATM_open (1:sNx,1:sNy)
113     _RL a_QSWbyATM_cover (1:sNx,1:sNy)
114     C a_QbyOCN :: available heat (in in W/m^2) due to the
115     C interaction of the ice pack and the ocean surface
116     C r_QbyOCN :: residual of a_QbyOCN after freezing/melting
117     C processes have been accounted for
118     _RL a_QbyOCN (1:sNx,1:sNy)
119     _RL r_QbyOCN (1:sNx,1:sNy)
120    
121     C conversion factors to go from Q (W/m2) to HEFF (ice meters)
122     _RL convertQ2HI, convertHI2Q
123     C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
124     _RL convertPRECIP2HI, convertHI2PRECIP
125    
126     #ifdef ALLOW_DIAGNOSTICS
127     C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
128     _RL d_AREAbyATM (1:sNx,1:sNy)
129     _RL d_AREAbyOCN (1:sNx,1:sNy)
130     _RL d_AREAbyICE (1:sNx,1:sNy)
131     #endif
132    
133     #ifdef SEAICE_ALLOW_AREA_RELAXATION
134     C ICE/SNOW stocks tendency associated with relaxation towards observation
135     _RL d_AREAbyRLX (1:sNx,1:sNy)
136     c The change of mean ice thickness due to relaxation
137     _RL d_HEFFbyRLX (1:sNx,1:sNy)
138     #endif
139    
140     c The change of mean ice thickness due to out-of-bounds values following
141     c sea ice dyhnamics
142     _RL d_HEFFbyNEG (1:sNx,1:sNy)
143    
144     c The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
145     _RL d_HEFFbyOCNonICE (1:sNx,1:sNy)
146    
147     c The sum of mean ice thickness increments due to atmospheric fluxes over the open water
148     c fraction and ice-covered fractions of the grid cell
149     _RL d_HEFFbyATMonOCN (1:sNx,1:sNy)
150     c The change of mean ice thickness due to flooding by snow
151     _RL d_HEFFbyFLOODING (1:sNx,1:sNy)
152    
153     c The mean ice thickness increments due to atmospheric fluxes over the open water
154     c fraction and ice-covered fractions of the grid cell, respectively
155     _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)
156     _RL d_HEFFbyATMonOCN_cover(1:sNx,1:sNy)
157    
158     _RL d_HSNWbyNEG (1:sNx,1:sNy)
159     _RL d_HSNWbyATMonSNW (1:sNx,1:sNy)
160     _RL d_HSNWbyOCNonSNW (1:sNx,1:sNy)
161     _RL d_HSNWbyRAIN (1:sNx,1:sNy)
162    
163     _RL d_HFRWbyRAIN (1:sNx,1:sNy)
164     C
165     C a_FWbySublim :: fresh water flux implied by latent heat of
166     C sublimation to atmosphere, same sign convention
167     C as EVAP (positive upward)
168     _RL a_FWbySublim (1:sNx,1:sNy)
169     _RL r_FWbySublim (1:sNx,1:sNy)
170     _RL d_HEFFbySublim (1:sNx,1:sNy)
171     _RL d_HSNWbySublim (1:sNx,1:sNy)
172    
173     #ifdef SEAICE_CAP_SUBLIM
174     C The latent heat flux which will sublimate all snow and ice
175     C over one time step
176     _RL latentHeatFluxMax (1:sNx,1:sNy)
177     _RL latentHeatFluxMaxMult (1:sNx,1:sNy,MULTDIM)
178     #endif
179    
180     C actual ice thickness (with upper and lower limit)
181     _RL heffActual (1:sNx,1:sNy)
182     C actual snow thickness
183     _RL hsnowActual (1:sNx,1:sNy)
184     C actual ice thickness (with lower limit only) Reciprocal
185     _RL recip_heffActual (1:sNx,1:sNy)
186     C local value (=1/HO or 1/HO_south)
187     _RL recip_HO
188     C local value (=1/ice thickness)
189     _RL recip_HH
190    
191     C AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
192     _RL AREApreTH (1:sNx,1:sNy)
193     _RL HEFFpreTH (1:sNx,1:sNy)
194     _RL HSNWpreTH (1:sNx,1:sNy)
195    
196     C wind speed
197     _RL UG (1:sNx,1:sNy)
198     #ifdef ALLOW_ATM_WIND
199     _RL SPEED_SQ
200     #endif
201    
202     C Regularization values squared
203     _RL area_reg_sq, hice_reg_sq
204    
205     C pathological cases thresholds
206     _RL heffTooHeavy
207    
208     _RL lhSublim
209    
210     C temporary variables available for the various computations
211     _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
212     _RL tmparr1 (1:sNx,1:sNy)
213    
214     #ifdef SEAICE_VARIABLE_SALINITY
215     _RL saltFluxAdjust (1:sNx,1:sNy)
216     #endif
217    
218     INTEGER ilockey
219     INTEGER it
220     _RL pFac
221     _RL ticeInMult (1:sNx,1:sNy,MULTDIM)
222     _RL ticeOutMult (1:sNx,1:sNy,MULTDIM)
223     _RL heffActualMult (1:sNx,1:sNy,MULTDIM)
224     _RL a_QbyATMmult_cover (1:sNx,1:sNy,MULTDIM)
225     _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
226     _RL a_FWbySublimMult (1:sNx,1:sNy,MULTDIM)
227     C Helper variables: reciprocal of some constants
228     _RL recip_multDim
229     _RL recip_deltaTtherm
230     _RL recip_rhoIce
231    
232     C Factor by which we increase the upper ocean friction velocity (u*) when
233     C ice is absent in a grid cell (dimensionless)
234     _RL MixedLayerTurbulenceFactor
235    
236     #ifdef ALLOW_SITRACER
237     INTEGER iTr
238     CHARACTER*8 diagName
239     #endif
240     #ifdef ALLOW_DIAGNOSTICS
241     c Helper variables for diagnostics
242     _RL DIAGarrayA (1:sNx,1:sNy)
243     _RL DIAGarrayB (1:sNx,1:sNy)
244     _RL DIAGarrayC (1:sNx,1:sNy)
245     _RL DIAGarrayD (1:sNx,1:sNy)
246     #endif
247    
248    
249     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
250    
251     C ===================================================================
252     C =================PART 0: constants and initializations=============
253     C ===================================================================
254    
255     IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
256     kSurface = Nr
257     ELSE
258     kSurface = 1
259     ENDIF
260    
261     C avoid unnecessary divisions in loops
262     recip_multDim = SEAICE_multDim
263     recip_multDim = ONE / recip_multDim
264     C above/below: double/single precision calculation of recip_multDim
265     c recip_multDim = 1./float(SEAICE_multDim)
266     recip_deltaTtherm = ONE / SEAICE_deltaTtherm
267     recip_rhoIce = ONE / SEAICE_rhoIce
268    
269     C Cutoff for iceload
270     heffTooHeavy=drF(kSurface) / 5. _d 0
271    
272     C RATIO OF SEA ICE DENSITY to SNOW DENSITY
273     ICE2SNOW = SEAICE_rhoIce/SEAICE_rhoSnow
274     SNOW2ICE = ONE / ICE2SNOW
275    
276     C HEAT OF FUSION OF ICE (J/m^3)
277     QI = SEAICE_rhoIce*SEAICE_lhFusion
278     recip_QI = ONE / QI
279     C HEAT OF FUSION OF SNOW (J/m^3)
280     QS = SEAICE_rhoSnow*SEAICE_lhFusion
281    
282     C ICE LATENT HEAT CONSTANT
283     lhSublim = SEAICE_lhEvap + SEAICE_lhFusion
284    
285     C regularization constants
286     area_reg_sq = SEAICE_area_reg * SEAICE_area_reg
287     hice_reg_sq = SEAICE_hice_reg * SEAICE_hice_reg
288    
289     C conversion factors to go from Q (W/m2) to HEFF (ice meters)
290     convertQ2HI=SEAICE_deltaTtherm/QI
291     convertHI2Q = ONE/convertQ2HI
292     C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
293     convertPRECIP2HI=SEAICE_deltaTtherm*rhoConstFresh/SEAICE_rhoIce
294     convertHI2PRECIP = ONE/convertPRECIP2HI
295    
296     DO bj=myByLo(myThid),myByHi(myThid)
297     DO bi=myBxLo(myThid),myBxHi(myThid)
298    
299     #ifdef ALLOW_AUTODIFF_TAMC
300     act1 = bi - myBxLo(myThid)
301     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
302     act2 = bj - myByLo(myThid)
303     max2 = myByHi(myThid) - myByLo(myThid) + 1
304     act3 = myThid - 1
305     max3 = nTx*nTy
306     act4 = ikey_dynamics - 1
307     iicekey = (act1 + 1) + act2*max1
308     & + act3*max1*max2
309     & + act4*max1*max2*max3
310     #endif /* ALLOW_AUTODIFF_TAMC */
311    
312    
313     C array initializations
314     C =====================
315    
316     DO J=1,sNy
317     DO I=1,sNx
318     a_QbyATM_cover (I,J) = 0.0 _d 0
319     a_QbyATM_open(I,J) = 0.0 _d 0
320     r_QbyATM_cover (I,J) = 0.0 _d 0
321     r_QbyATM_open (I,J) = 0.0 _d 0
322    
323     a_QSWbyATM_open (I,J) = 0.0 _d 0
324     a_QSWbyATM_cover (I,J) = 0.0 _d 0
325    
326     a_QbyOCN (I,J) = 0.0 _d 0
327     r_QbyOCN (I,J) = 0.0 _d 0
328    
329     #ifdef ALLOW_DIAGNOSTICS
330     d_AREAbyATM(I,J) = 0.0 _d 0
331     d_AREAbyICE(I,J) = 0.0 _d 0
332     d_AREAbyOCN(I,J) = 0.0 _d 0
333     #endif
334    
335     #ifdef SEAICE_ALLOW_AREA_RELAXATION
336     d_AREAbyRLX(I,J) = 0.0 _d 0
337     d_HEFFbyRLX(I,J) = 0.0 _d 0
338     #endif
339    
340     d_HEFFbyNEG(I,J) = 0.0 _d 0
341     d_HEFFbyOCNonICE(I,J) = 0.0 _d 0
342     d_HEFFbyATMonOCN(I,J) = 0.0 _d 0
343     d_HEFFbyFLOODING(I,J) = 0.0 _d 0
344    
345     d_HEFFbyATMonOCN_open(I,J) = 0.0 _d 0
346     d_HEFFbyATMonOCN_cover(I,J) = 0.0 _d 0
347    
348     d_HSNWbyNEG(I,J) = 0.0 _d 0
349     d_HSNWbyATMonSNW(I,J) = 0.0 _d 0
350     d_HSNWbyOCNonSNW(I,J) = 0.0 _d 0
351     d_HSNWbyRAIN(I,J) = 0.0 _d 0
352     a_FWbySublim(I,J) = 0.0 _d 0
353     r_FWbySublim(I,J) = 0.0 _d 0
354     d_HEFFbySublim(I,J) = 0.0 _d 0
355     d_HSNWbySublim(I,J) = 0.0 _d 0
356     #ifdef SEAICE_CAP_SUBLIM
357     latentHeatFluxMax(I,J) = 0.0 _d 0
358     #endif
359     c
360     d_HFRWbyRAIN(I,J) = 0.0 _d 0
361    
362     tmparr1(I,J) = 0.0 _d 0
363    
364     #ifdef SEAICE_VARIABLE_SALINITY
365     saltFluxAdjust(I,J) = 0.0 _d 0
366     #endif
367     DO IT=1,SEAICE_multDim
368     ticeInMult(I,J,IT) = 0.0 _d 0
369     ticeOutMult(I,J,IT) = 0.0 _d 0
370     a_QbyATMmult_cover(I,J,IT) = 0.0 _d 0
371     a_QSWbyATMmult_cover(I,J,IT) = 0.0 _d 0
372     a_FWbySublimMult(I,J,IT) = 0.0 _d 0
373     #ifdef SEAICE_CAP_SUBLIM
374     latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
375     #endif
376     ENDDO
377     ENDDO
378     ENDDO
379     #if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
380     DO J=1-oLy,sNy+oLy
381     DO I=1-oLx,sNx+oLx
382     frWtrAtm(I,J,bi,bj) = 0.0 _d 0
383     ENDDO
384     ENDDO
385     #endif
386    
387    
388     C =====================================================================
389     C ===========PART 1: treat pathological cases (post advdiff)===========
390     C =====================================================================
391    
392     #ifdef SEAICE_GROWTH_LEGACY
393    
394     DO J=1,sNy
395     DO I=1,sNx
396     HEFFpreTH(I,J)=HEFFNM1(I,J,bi,bj)
397     HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
398     AREApreTH(I,J)=AREANM1(I,J,bi,bj)
399     d_HEFFbyNEG(I,J)=0. _d 0
400     d_HSNWbyNEG(I,J)=0. _d 0
401     #ifdef ALLOW_DIAGNOSTICS
402     DIAGarrayA(I,J) = AREANM1(I,J,bi,bj)
403     DIAGarrayB(I,J) = AREANM1(I,J,bi,bj)
404     DIAGarrayC(I,J) = HEFFNM1(I,J,bi,bj)
405     DIAGarrayD(I,J) = HSNOW(I,J,bi,bj)
406     #endif
407     ENDDO
408     ENDDO
409    
410     #else /* SEAICE_GROWTH_LEGACY */
411    
412     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
413     Cgf no dependency through pathological cases treatment
414     IF ( SEAICEadjMODE.EQ.0 ) THEN
415     #endif
416    
417     #ifdef SEAICE_ALLOW_AREA_RELAXATION
418     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
419     CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
420     C 0) relax sea ice concentration towards observation
421     IF (SEAICE_tauAreaObsRelax .GT. 0.) THEN
422     DO J=1,sNy
423     DO I=1,sNx
424     C d_AREAbyRLX(i,j) = 0. _d 0
425     C d_HEFFbyRLX(i,j) = 0. _d 0
426     IF ( obsSIce(I,J,bi,bj).GT.AREA(I,J,bi,bj)) THEN
427     d_AREAbyRLX(i,j) =
428     & SEAICE_deltaTtherm/SEAICE_tauAreaObsRelax
429     & * (obsSIce(I,J,bi,bj) - AREA(I,J,bi,bj))
430     ENDIF
431     IF ( obsSIce(I,J,bi,bj).GT.0. _d 0 .AND.
432     & AREA(I,J,bi,bj).EQ.0. _d 0) THEN
433     C d_HEFFbyRLX(i,j) = 1. _d 1 * siEps * d_AREAbyRLX(i,j)
434     d_HEFFbyRLX(i,j) = 1. _d 1 * siEps
435     ENDIF
436     AREA(I,J,bi,bj) = AREA(I,J,bi,bj) + d_AREAbyRLX(i,j)
437     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + d_HEFFbyRLX(i,j)
438     ENDDO
439     ENDDO
440     ENDIF
441     #endif /* SEAICE_ALLOW_AREA_RELAXATION */
442    
443     C 1) treat the case of negative values:
444    
445     #ifdef ALLOW_AUTODIFF_TAMC
446     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
447     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
448     CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
449     #endif /* ALLOW_AUTODIFF_TAMC */
450     DO J=1,sNy
451     DO I=1,sNx
452     d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
453     HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
454     d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
455     HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
456     AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
457     ENDDO
458     ENDDO
459    
460     C 1.25) treat the case of very thin ice:
461    
462     #ifdef ALLOW_AUTODIFF_TAMC
463     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
464     #endif /* ALLOW_AUTODIFF_TAMC */
465     DO J=1,sNy
466     DO I=1,sNx
467     tmpscal2=0. _d 0
468     tmpscal3=0. _d 0
469     IF (HEFF(I,J,bi,bj).LE.siEps) THEN
470     tmpscal2=-HEFF(I,J,bi,bj)
471     tmpscal3=-HSNOW(I,J,bi,bj)
472     TICE(I,J,bi,bj)=celsius2K
473     DO IT=1,SEAICE_multDim
474     TICES(I,J,IT,bi,bj)=celsius2K
475     ENDDO
476     ENDIF
477     HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal2
478     d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
479     HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
480     d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
481     ENDDO
482     ENDDO
483    
484     C 1.5) treat the case of area but no ice/snow:
485    
486     #ifdef ALLOW_AUTODIFF_TAMC
487     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
488     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
489     #endif /* ALLOW_AUTODIFF_TAMC */
490     DO J=1,sNy
491     DO I=1,sNx
492     IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
493     & (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
494     ENDDO
495     ENDDO
496    
497     C 2) treat the case of very small area:
498    
499     #ifndef DISABLE_AREA_FLOOR
500     #ifdef ALLOW_AUTODIFF_TAMC
501     CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
502     #endif /* ALLOW_AUTODIFF_TAMC */
503     DO J=1,sNy
504     DO I=1,sNx
505     IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN
506     AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor)
507     ENDIF
508     ENDDO
509     ENDDO
510     #endif /* DISABLE_AREA_FLOOR */
511    
512     C 2.5) treat case of excessive ice cover, e.g., due to ridging:
513    
514     #ifdef ALLOW_AUTODIFF_TAMC
515     CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
516     #endif /* ALLOW_AUTODIFF_TAMC */
517     DO J=1,sNy
518     DO I=1,sNx
519     #ifdef ALLOW_DIAGNOSTICS
520     DIAGarrayA(I,J) = AREA(I,J,bi,bj)
521     #endif
522     #ifdef ALLOW_SITRACER
523     SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
524     #endif
525     AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)
526     ENDDO
527     ENDDO
528    
529     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
530     ENDIF
531     #endif
532    
533     C 3) store regularized values of heff, hsnow, area at the onset of thermo.
534     DO J=1,sNy
535     DO I=1,sNx
536     HEFFpreTH(I,J)=HEFF(I,J,bi,bj)
537     HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
538     AREApreTH(I,J)=AREA(I,J,bi,bj)
539     #ifdef ALLOW_DIAGNOSTICS
540     DIAGarrayB(I,J) = AREA(I,J,bi,bj)
541     DIAGarrayC(I,J) = HEFF(I,J,bi,bj)
542     DIAGarrayD(I,J) = HSNOW(I,J,bi,bj)
543     #endif
544     #ifdef ALLOW_SITRACER
545     SItrHEFF(I,J,bi,bj,1)=HEFF(I,J,bi,bj)
546     SItrAREA(I,J,bi,bj,2)=AREA(I,J,bi,bj)
547     #endif
548     ENDDO
549     ENDDO
550    
551     C 4) treat sea ice salinity pathological cases
552     #ifdef SEAICE_VARIABLE_SALINITY
553     #ifdef ALLOW_AUTODIFF_TAMC
554     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
555     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
556     #endif /* ALLOW_AUTODIFF_TAMC */
557     DO J=1,sNy
558     DO I=1,sNx
559     IF ( (HSALT(I,J,bi,bj) .LT. 0.0).OR.
560     & (HEFF(I,J,bi,bj) .EQ. 0.0) ) THEN
561     saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
562     & HSALT(I,J,bi,bj) * recip_deltaTtherm
563     HSALT(I,J,bi,bj) = 0.0 _d 0
564     ENDIF
565     ENDDO
566     ENDDO
567     #endif /* SEAICE_VARIABLE_SALINITY */
568    
569     #endif /* SEAICE_GROWTH_LEGACY */
570    
571     #ifdef ALLOW_DIAGNOSTICS
572     IF ( useDiagnostics ) THEN
573     CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIareaPR',0,1,3,bi,bj,myThid)
574     CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIareaPT',0,1,3,bi,bj,myThid)
575     CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIheffPT',0,1,3,bi,bj,myThid)
576     CALL DIAGNOSTICS_FILL(DIAGarrayD,'SIhsnoPT',0,1,3,bi,bj,myThid)
577     #ifdef ALLOW_SITRACER
578     DO iTr = 1, SItrNumInUse
579     WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'PT'
580     IF (SItrMate(iTr).EQ.'HEFF') THEN
581     CALL DIAGNOSTICS_FRACT_FILL(
582     I SItracer(1-OLx,1-OLy,bi,bj,iTr),HEFF(1-OLx,1-OLy,bi,bj),
583     I ONE, 1, diagName,0,1,2,bi,bj,myThid )
584     ELSE
585     CALL DIAGNOSTICS_FRACT_FILL(
586     I SItracer(1-OLx,1-OLy,bi,bj,iTr),AREA(1-OLx,1-OLy,bi,bj),
587     I ONE, 1, diagName,0,1,2,bi,bj,myThid )
588     ENDIF
589     ENDDO
590     #endif /* ALLOW_SITRACER */
591     ENDIF
592     #endif /* ALLOW_DIAGNOSTICS */
593    
594     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
595     Cgf no additional dependency of air-sea fluxes to ice
596     IF ( SEAICEadjMODE.GE.1 ) THEN
597     DO J=1,sNy
598     DO I=1,sNx
599     HEFFpreTH(I,J) = 0. _d 0
600     HSNWpreTH(I,J) = 0. _d 0
601     AREApreTH(I,J) = 0. _d 0
602     ENDDO
603     ENDDO
604     ENDIF
605     #endif
606    
607     #if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
608     DO J=1,sNy
609     DO I=1,sNx
610     AREAforAtmFW(I,J,bi,bj) = AREApreTH(I,J)
611     ENDDO
612     ENDDO
613     #endif
614    
615     C 4) COMPUTE ACTUAL ICE/SNOW THICKNESS; USE MIN/MAX VALUES
616     C TO REGULARIZE SEAICE_SOLVE4TEMP/d_AREA COMPUTATIONS
617    
618     #ifdef ALLOW_AUTODIFF_TAMC
619     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
620     CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
621     CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
622     #endif /* ALLOW_AUTODIFF_TAMC */
623     DO J=1,sNy
624     DO I=1,sNx
625    
626     IF (HEFFpreTH(I,J) .GT. ZERO) THEN
627     #ifdef SEAICE_GROWTH_LEGACY
628     tmpscal1 = MAX(SEAICE_area_reg,AREApreTH(I,J))
629     hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1
630     tmpscal2 = HEFFpreTH(I,J)/tmpscal1
631     heffActual(I,J) = MAX(tmpscal2,SEAICE_hice_reg)
632     #else /* SEAICE_GROWTH_LEGACY */
633     cif regularize AREA with SEAICE_area_reg
634     tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)
635     cif heffActual calculated with the regularized AREA
636     tmpscal2 = HEFFpreTH(I,J) / tmpscal1
637     cif regularize heffActual with SEAICE_hice_reg (add lower bound)
638     heffActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
639     cif hsnowActual calculated with the regularized AREA
640     hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1
641     #endif /* SEAICE_GROWTH_LEGACY */
642     cif regularize the inverse of heffActual by hice_reg
643     recip_heffActual(I,J) = AREApreTH(I,J) /
644     & sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)
645     cif Do not regularize when HEFFpreTH = 0
646     ELSE
647     heffActual(I,J) = ZERO
648     hsnowActual(I,J) = ZERO
649     recip_heffActual(I,J) = ZERO
650     ENDIF
651    
652     ENDDO
653     ENDDO
654    
655     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
656     CALL ZERO_ADJ_1D( sNx*sNy, heffActual, myThid)
657     CALL ZERO_ADJ_1D( sNx*sNy, hsnowActual, myThid)
658     CALL ZERO_ADJ_1D( sNx*sNy, recip_heffActual, myThid)
659     #endif
660    
661     #ifdef SEAICE_CAP_SUBLIM
662     C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
663     C AND SNOW THICKNESS
664     DO J=1,sNy
665     DO I=1,sNx
666     c The latent heat flux over the sea ice which
667     c will sublimate all of the snow and ice over one time
668     c step (W/m^2)
669     IF (HEFFpreTH(I,J) .GT. ZERO) THEN
670     latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *
671     & (HEFFpreTH(I,J) * SEAICE_rhoIce +
672     & HSNWpreTH(I,J) * SEAICE_rhoSnow)/AREApreTH(I,J)
673     ELSE
674     latentHeatFluxMax(I,J) = ZERO
675     ENDIF
676     ENDDO
677     ENDDO
678     #endif /* SEAICE_CAP_SUBLIM */
679    
680     C ===================================================================
681     C ================PART 2: determine heat fluxes/stocks===============
682     C ===================================================================
683    
684     C determine available heat due to the atmosphere -- for open water
685     C ================================================================
686    
687     DO j=1,sNy
688     DO i=1,sNx
689     C ocean surface/mixed layer temperature
690     TmixLoc(i,j) = theta(i,j,kSurface,bi,bj)+celsius2K
691     C wind speed from exf
692     UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))
693     ENDDO
694     ENDDO
695    
696     #ifdef ALLOW_AUTODIFF_TAMC
697     CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
698     CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
699     cCADJ STORE UG = comlev1_bibj, key = iicekey,byte=isbyte
700     cCADJ STORE TmixLoc = comlev1_bibj, key = iicekey,byte=isbyte
701     #endif /* ALLOW_AUTODIFF_TAMC */
702    
703     CALL SEAICE_BUDGET_OCEAN(
704     I UG,
705     I TmixLoc,
706     O a_QbyATM_open, a_QSWbyATM_open,
707     I bi, bj, myTime, myIter, myThid )
708    
709     C determine available heat due to the atmosphere -- for ice covered water
710     C =======================================================================
711    
712     #ifdef ALLOW_ATM_WIND
713     IF (useRelativeWind) THEN
714     C Compute relative wind speed over sea ice.
715     DO J=1,sNy
716     DO I=1,sNx
717     SPEED_SQ =
718     & (uWind(I,J,bi,bj)
719     & +0.5 _d 0*(uVel(i,j,kSurface,bi,bj)
720     & +uVel(i+1,j,kSurface,bi,bj))
721     & -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj)))**2
722     & +(vWind(I,J,bi,bj)
723     & +0.5 _d 0*(vVel(i,j,kSurface,bi,bj)
724     & +vVel(i,j+1,kSurface,bi,bj))
725     & -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj)))**2
726     IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
727     UG(I,J)=SEAICE_EPS
728     ELSE
729     UG(I,J)=SQRT(SPEED_SQ)
730     ENDIF
731     ENDDO
732     ENDDO
733     ENDIF
734     #endif /* ALLOW_ATM_WIND */
735    
736     #ifdef ALLOW_AUTODIFF_TAMC
737     CADJ STORE tice(:,:,bi,bj)
738     CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
739     CADJ STORE hsnowActual = comlev1_bibj, key = iicekey, byte = isbyte
740     CADJ STORE heffActual = comlev1_bibj, key = iicekey, byte = isbyte
741     CADJ STORE UG = comlev1_bibj, key = iicekey, byte = isbyte
742     CADJ STORE tices(:,:,:,bi,bj)
743     CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
744     CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
745     CADJ & key = iicekey, byte = isbyte
746     #endif /* ALLOW_AUTODIFF_TAMC */
747    
748     C-- Start loop over multi-categories
749     DO IT=1,SEAICE_multDim
750     c homogeneous distribution between 0 and 2 x heffActual
751     pFac = (2.0 _d 0*real(IT)-1.0 _d 0)*recip_multDim
752     DO J=1,sNy
753     DO I=1,sNx
754     heffActualMult(I,J,IT)= heffActual(I,J)*pFac
755     #ifdef SEAICE_CAP_SUBLIM
756     latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac
757     #endif
758     ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
759     ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
760     TICE(I,J,bi,bj) = ZERO
761     TICES(I,J,IT,bi,bj) = ZERO
762     ENDDO
763     ENDDO
764     ENDDO
765    
766     #ifdef ALLOW_AUTODIFF_TAMC
767     CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
768     CADJ STORE ticeInMult = comlev1_bibj, key = iicekey, byte = isbyte
769     # ifdef SEAICE_CAP_SUBLIM
770     CADJ STORE latentHeatFluxMaxMult
771     CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
772     # endif
773     CADJ STORE a_QbyATMmult_cover =
774     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
775     CADJ STORE a_QSWbyATMmult_cover =
776     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
777     CADJ STORE a_FWbySublimMult =
778     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
779     #endif /* ALLOW_AUTODIFF_TAMC */
780    
781     DO IT=1,SEAICE_multDim
782     CALL SEAICE_SOLVE4TEMP(
783     I UG, heffActualMult(1,1,IT), hsnowActual,
784     #ifdef SEAICE_CAP_SUBLIM
785     I latentHeatFluxMaxMult(1,1,IT),
786     #endif
787     U ticeInMult(1,1,IT), ticeOutMult(1,1,IT),
788     O a_QbyATMmult_cover(1,1,IT), a_QSWbyATMmult_cover(1,1,IT),
789     O a_FWbySublimMult(1,1,IT),
790     I bi, bj, myTime, myIter, myThid )
791     ENDDO
792    
793     #ifdef ALLOW_AUTODIFF_TAMC
794     CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
795     CADJ STORE ticeOutMult = comlev1_bibj, key = iicekey, byte = isbyte
796     # ifdef SEAICE_CAP_SUBLIM
797     CADJ STORE latentHeatFluxMaxMult
798     CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
799     # endif
800     CADJ STORE a_QbyATMmult_cover =
801     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
802     CADJ STORE a_QSWbyATMmult_cover =
803     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
804     CADJ STORE a_FWbySublimMult =
805     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
806     #endif /* ALLOW_AUTODIFF_TAMC */
807    
808     DO IT=1,SEAICE_multDim
809     DO J=1,sNy
810     DO I=1,sNx
811     C update TICE & TICES
812     TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
813     & + ticeOutMult(I,J,IT)*recip_multDim
814     TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)
815     C average over categories
816     a_QbyATM_cover (I,J) = a_QbyATM_cover(I,J)
817     & + a_QbyATMmult_cover(I,J,IT)*recip_multDim
818     a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
819     & + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim
820     a_FWbySublim (I,J) = a_FWbySublim(I,J)
821     & + a_FWbySublimMult(I,J,IT)*recip_multDim
822     ENDDO
823     ENDDO
824     ENDDO
825    
826     #ifdef SEAICE_CAP_SUBLIM
827     # ifdef ALLOW_DIAGNOSTICS
828     DO J=1,sNy
829     DO I=1,sNx
830     c The actual latent heat flux realized by SOLVE4TEMP
831     DIAGarrayA(I,J) = a_FWbySublim(I,J) * lhSublim
832     ENDDO
833     ENDDO
834     cif The actual vs. maximum latent heat flux
835     IF ( useDiagnostics ) THEN
836     CALL DIAGNOSTICS_FILL(DIAGarrayA,
837     & 'SIactLHF',0,1,3,bi,bj,myThid)
838     CALL DIAGNOSTICS_FILL(latentHeatFluxMax,
839     & 'SImaxLHF',0,1,3,bi,bj,myThid)
840     ENDIF
841     # endif /* ALLOW_DIAGNOSTICS */
842     #endif /* SEAICE_CAP_SUBLIM */
843    
844     #ifdef ALLOW_AUTODIFF_TAMC
845     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
846     CADJ STORE a_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
847     CADJ STORE a_QSWbyATM_cover= comlev1_bibj, key = iicekey, byte = isbyte
848     CADJ STORE a_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
849     CADJ STORE a_QSWbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
850     CADJ STORE a_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
851     #endif /* ALLOW_AUTODIFF_TAMC */
852    
853     C switch heat fluxes from W/m2 to 'effective' ice meters
854     DO J=1,sNy
855     DO I=1,sNx
856     a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)
857     & * convertQ2HI * AREApreTH(I,J)
858     a_QSWbyATM_cover(I,J) = a_QSWbyATM_cover(I,J)
859     & * convertQ2HI * AREApreTH(I,J)
860     a_QbyATM_open(I,J) = a_QbyATM_open(I,J)
861     & * convertQ2HI * ( ONE - AREApreTH(I,J) )
862     a_QSWbyATM_open(I,J) = a_QSWbyATM_open(I,J)
863     & * convertQ2HI * ( ONE - AREApreTH(I,J) )
864     C and initialize r_QbyATM_cover/r_QbyATM_open
865     r_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
866     r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
867     C Convert fresh water flux by sublimation to 'effective' ice meters.
868     C Negative sublimation is resublimation and will be added as snow.
869     #ifdef SEAICE_DISABLE_SUBLIM
870     cgf just for those who may need to omit this term to reproduce old results
871     a_FWbySublim(I,J) = ZERO
872     #endif /* SEAICE_CAP_SUBLIM */
873     a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce
874     & * a_FWbySublim(I,J)*AREApreTH(I,J)
875     r_FWbySublim(I,J)=a_FWbySublim(I,J)
876     ENDDO
877     ENDDO
878    
879     #ifdef ALLOW_AUTODIFF_TAMC
880     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
881     CADJ STORE a_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
882     CADJ STORE a_QSWbyATM_cover= comlev1_bibj, key = iicekey, byte = isbyte
883     CADJ STORE a_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
884     CADJ STORE a_QSWbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
885     CADJ STORE a_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
886     CADJ STORE r_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
887     CADJ STORE r_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
888     CADJ STORE r_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
889     #endif /* ALLOW_AUTODIFF_TAMC */
890    
891     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
892     Cgf no additional dependency through ice cover
893     IF ( SEAICEadjMODE.GE.3 ) THEN
894     DO J=1,sNy
895     DO I=1,sNx
896     a_QbyATM_cover(I,J) = 0. _d 0
897     r_QbyATM_cover(I,J) = 0. _d 0
898     a_QSWbyATM_cover(I,J) = 0. _d 0
899     ENDDO
900     ENDDO
901     ENDIF
902     #endif
903    
904     C determine available heat due to the ice pack tying the
905     C underlying surface water temperature to freezing point
906     C ======================================================
907    
908     #ifdef ALLOW_AUTODIFF_TAMC
909     CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
910     CADJ & key = iicekey, byte = isbyte
911     CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
912     CADJ & key = iicekey, byte = isbyte
913     #endif
914    
915     DO J=1,sNy
916     DO I=1,sNx
917     c FREEZING TEMP. OF SEA WATER (deg C)
918     tempFrz = SEAICE_tempFrz0 +
919     & SEAICE_dTempFrz_dS *salt(I,J,kSurface,bi,bj)
920     c efficiency of turbulent fluxes : dependency to sign of THETA-TBC
921     IF ( theta(I,J,kSurface,bi,bj) .GE. tempFrz ) THEN
922     tmpscal1 = SEAICE_mcPheePiston
923     ELSE
924     tmpscal1 =SEAICE_frazilFrac*drF(kSurface)/SEAICE_deltaTtherm
925     ENDIF
926     c efficiency of turbulent fluxes : dependency to AREA (McPhee cases)
927     IF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.
928     & (.NOT.SEAICE_mcPheeStepFunc) ) THEN
929     MixedLayerTurbulenceFactor = ONE -
930     & SEAICE_mcPheeTaper * AREApreTH(I,J)
931     ELSEIF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.
932     & (SEAICE_mcPheeStepFunc) ) THEN
933     MixedLayerTurbulenceFactor = ONE - SEAICE_mcPheeTaper
934     ELSE
935     MixedLayerTurbulenceFactor = ONE
936     ENDIF
937     c maximum turbulent flux, in ice meters
938     tmpscal2= - (HeatCapacity_Cp*rhoConst * recip_QI)
939     & * (theta(I,J,kSurface,bi,bj)-tempFrz)
940     & * SEAICE_deltaTtherm * maskC(i,j,kSurface,bi,bj)
941     c available turbulent flux
942     a_QbyOCN(i,j) =
943     & tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor
944     r_QbyOCN(i,j) = a_QbyOCN(i,j)
945     ENDDO
946     ENDDO
947    
948     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
949     CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
950     #endif
951    
952    
953     C ===================================================================
954     C =========PART 3: determine effective thicknesses increments========
955     C ===================================================================
956    
957     C compute snow/ice tendency due to sublimation
958     C ============================================
959    
960     #ifdef ALLOW_AUTODIFF_TAMC
961     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
962     CADJ STORE r_FWbySublim = comlev1_bibj,key=iicekey,byte=isbyte
963     #endif /* ALLOW_AUTODIFF_TAMC */
964     DO J=1,sNy
965     DO I=1,sNx
966     C First sublimate/deposite snow
967     tmpscal2 =
968     & MAX(MIN(r_FWbySublim(I,J),HSNOW(I,J,bi,bj)*SNOW2ICE),ZERO)
969     d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
970     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) - tmpscal2*ICE2SNOW
971     r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2
972     ENDDO
973     ENDDO
974     #ifdef ALLOW_AUTODIFF_TAMC
975     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
976     CADJ STORE r_FWbySublim = comlev1_bibj,key=iicekey,byte=isbyte
977     #endif /* ALLOW_AUTODIFF_TAMC */
978     DO J=1,sNy
979     DO I=1,sNx
980     C If anything is left, sublimate ice
981     tmpscal2 =
982     & MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
983     d_HEFFbySublim(I,J) = - tmpscal2
984     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) - tmpscal2
985     r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2
986     ENDDO
987     ENDDO
988     DO J=1,sNy
989     DO I=1,sNx
990     C If anything is left, it will be evaporated from the ocean rather than sublimated.
991     C Since a_QbyATM_cover was computed for sublimation, not simple evapation, we need to
992     C remove the fusion part for the residual (that happens to be precisely r_FWbySublim).
993     a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)-r_FWbySublim(I,J)
994     r_QbyATM_cover(I,J) = r_QbyATM_cover(I,J)-r_FWbySublim(I,J)
995     ENDDO
996     ENDDO
997    
998     C compute ice thickness tendency due to ice-ocean interaction
999     C ===========================================================
1000    
1001     #ifdef ALLOW_AUTODIFF_TAMC
1002     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1003     CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1004     #endif /* ALLOW_AUTODIFF_TAMC */
1005    
1006     DO J=1,sNy
1007     DO I=1,sNx
1008     d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj))
1009     r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
1010     HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J)
1011     #ifdef ALLOW_SITRACER
1012     SItrHEFF(I,J,bi,bj,2)=HEFF(I,J,bi,bj)
1013     #endif
1014     ENDDO
1015     ENDDO
1016    
1017     C compute snow melt tendency due to snow-atmosphere interaction
1018     C ==================================================================
1019    
1020     #ifdef ALLOW_AUTODIFF_TAMC
1021     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1022     CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1023     #endif /* ALLOW_AUTODIFF_TAMC */
1024    
1025     DO J=1,sNy
1026     DO I=1,sNx
1027     C Convert to standard units (meters of ice) rather than to meters
1028     C of snow. This appears to be more robust.
1029     tmpscal1=MAX(r_QbyATM_cover(I,J),-HSNOW(I,J,bi,bj)*SNOW2ICE)
1030     tmpscal2=MIN(tmpscal1,0. _d 0)
1031     #ifdef SEAICE_MODIFY_GROWTH_ADJ
1032     Cgf no additional dependency through snow
1033     IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1034     #endif
1035     d_HSNWbyATMonSNW(I,J)= tmpscal2*ICE2SNOW
1036     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal2*ICE2SNOW
1037     r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1038     ENDDO
1039     ENDDO
1040    
1041     C compute ice thickness tendency due to the atmosphere
1042     C ====================================================
1043    
1044     #ifdef ALLOW_AUTODIFF_TAMC
1045     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1046     CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1047     #endif /* ALLOW_AUTODIFF_TAMC */
1048    
1049     Cgf note: this block is not actually tested by lab_sea
1050     Cgf where all experiments start in January. So even though
1051     Cgf the v1.81=>v1.82 revision would change results in
1052     Cgf warming conditions, the lab_sea results were not changed.
1053    
1054     DO J=1,sNy
1055     DO I=1,sNx
1056    
1057     #ifdef SEAICE_GROWTH_LEGACY
1058     tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))
1059     #else
1060     tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+
1061     c Limit ice growth by potential melt by ocean
1062     & AREApreTH(I,J) * r_QbyOCN(I,J))
1063     #endif /* SEAICE_GROWTH_LEGACY */
1064    
1065     d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
1066     d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
1067     r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
1068     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2
1069    
1070     #ifdef ALLOW_SITRACER
1071     SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
1072     #endif
1073     ENDDO
1074     ENDDO
1075    
1076     C attribute precip to fresh water or snow stock,
1077     C depending on atmospheric conditions.
1078     C =================================================
1079     #ifdef ALLOW_ATM_TEMP
1080     #ifdef ALLOW_AUTODIFF_TAMC
1081     CADJ STORE a_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1082     CADJ STORE PRECIP(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1083     CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
1084     #endif /* ALLOW_AUTODIFF_TAMC */
1085     DO J=1,sNy
1086     DO I=1,sNx
1087     C possible alternatives to the a_QbyATM_cover criterium
1088     c IF (TICE(I,J,bi,bj) .LT. TMIX) THEN
1089     c IF (atemp(I,J,bi,bj) .LT. celsius2K) THEN
1090     IF ( a_QbyATM_cover(I,J).GE. 0. _d 0 ) THEN
1091     C add precip as snow
1092     d_HFRWbyRAIN(I,J)=0. _d 0
1093     d_HSNWbyRAIN(I,J)=convertPRECIP2HI*ICE2SNOW*
1094     & PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1095     ELSE
1096     C add precip to the fresh water bucket
1097     d_HFRWbyRAIN(I,J)=-convertPRECIP2HI*
1098     & PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1099     d_HSNWbyRAIN(I,J)=0. _d 0
1100     ENDIF
1101     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
1102     ENDDO
1103     ENDDO
1104     Cgf note: this does not affect air-sea heat flux,
1105     Cgf since the implied air heat gain to turn
1106     Cgf rain to snow is not a surface process.
1107     #endif /* ALLOW_ATM_TEMP */
1108    
1109     C compute snow melt due to heat available from ocean.
1110     C =================================================================
1111    
1112     Cgf do we need to keep this comment and cpp bracket?
1113     Cph( very sensitive bit here by JZ
1114     #ifndef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING
1115     #ifdef ALLOW_AUTODIFF_TAMC
1116     CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1117     CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1118     #endif /* ALLOW_AUTODIFF_TAMC */
1119     DO J=1,sNy
1120     DO I=1,sNx
1121     tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))
1122     tmpscal2=MIN(tmpscal1,0. _d 0)
1123     #ifdef SEAICE_MODIFY_GROWTH_ADJ
1124     Cgf no additional dependency through snow
1125     if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1126     #endif
1127     d_HSNWbyOCNonSNW(I,J) = tmpscal2
1128     r_QbyOCN(I,J)=r_QbyOCN(I,J)
1129     & -d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
1130     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_HSNWbyOCNonSNW(I,J)
1131     ENDDO
1132     ENDDO
1133     #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
1134     Cph)
1135    
1136     C gain of new ice over open water
1137     C ===============================
1138     #ifdef ALLOW_AUTODIFF_TAMC
1139     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1140     CADJ STORE r_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1141     CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1142     CADJ STORE a_QSWbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1143     CADJ STORE a_QSWbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1144     #endif /* ALLOW_AUTODIFF_TAMC */
1145    
1146     DO J=1,sNy
1147     DO I=1,sNx
1148     c Initial ice growth is triggered by open water
1149     c heat flux overcoming potential melt by ocean
1150     tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *
1151     & (1.0 _d 0 - AREApreTH(I,J))
1152     c Penetrative shortwave flux beyond first layer
1153     c that is therefore not available to ice growth/melt
1154     tmpscal2=SWFracB * a_QSWbyATM_open(I,J)
1155     C impose -HEFF as the maxmum melting if SEAICE_doOpenWaterMelt
1156     C or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)
1157     tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,
1158     & -HEFF(I,J,bi,bj)*facOpenMelt)*HEFFM(I,J,bi,bj)
1159     d_HEFFbyATMonOCN_open(I,J)=tmpscal3
1160     d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
1161     r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
1162     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
1163     ENDDO
1164     ENDDO
1165    
1166     #ifdef ALLOW_SITRACER
1167     DO J=1,sNy
1168     DO I=1,sNx
1169     c needs to be here to allow use also with LEGACY branch
1170     SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
1171     ENDDO
1172     ENDDO
1173     #endif /* ALLOW_SITRACER */
1174    
1175     C convert snow to ice if submerged.
1176     C =================================
1177    
1178     #ifndef SEAICE_GROWTH_LEGACY
1179     C note: in legacy, this process is done at the end
1180     #ifdef ALLOW_AUTODIFF_TAMC
1181     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1182     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1183     #endif /* ALLOW_AUTODIFF_TAMC */
1184     IF ( SEAICEuseFlooding ) THEN
1185     DO J=1,sNy
1186     DO I=1,sNx
1187     tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1188     & +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
1189     tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
1190     d_HEFFbyFLOODING(I,J)=tmpscal1
1191     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
1192     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1193     & d_HEFFbyFLOODING(I,J)*ICE2SNOW
1194     ENDDO
1195     ENDDO
1196     ENDIF
1197     #endif /* SEAICE_GROWTH_LEGACY */
1198    
1199     C ===================================================================
1200     C ==========PART 4: determine ice cover fraction increments=========-
1201     C ===================================================================
1202    
1203     #ifdef ALLOW_AUTODIFF_TAMC
1204     CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
1205     CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=iicekey,byte=isbyte
1206     CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=iicekey,byte=isbyte
1207     CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
1208     CADJ STORE recip_heffActual = comlev1_bibj,key=iicekey,byte=isbyte
1209     CADJ STORE d_hsnwbyatmonsnw = comlev1_bibj,key=iicekey,byte=isbyte
1210     cph(
1211     cphCADJ STORE d_AREAbyATM = comlev1_bibj,key=iicekey,byte=isbyte
1212     cphCADJ STORE d_AREAbyICE = comlev1_bibj,key=iicekey,byte=isbyte
1213     cphCADJ STORE d_AREAbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1214     cph)
1215     CADJ STORE a_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1216     CADJ STORE heffActual = comlev1_bibj,key=iicekey,byte=isbyte
1217     CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
1218     CADJ STORE HEFF(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1219     CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1220     CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1221     #endif /* ALLOW_AUTODIFF_TAMC */
1222    
1223     DO J=1,sNy
1224     DO I=1,sNx
1225    
1226     IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
1227     recip_HO=1. _d 0 / HO_south
1228     ELSE
1229     recip_HO=1. _d 0 / HO
1230     ENDIF
1231     #ifdef SEAICE_GROWTH_LEGACY
1232     tmpscal0=HEFF(I,J,bi,bj) - d_HEFFbyATMonOCN(I,J)
1233     recip_HH = AREApreTH(I,J) /(tmpscal0+.00001 _d 0)
1234     #else
1235     recip_HH = recip_heffActual(I,J)
1236     #endif
1237    
1238     C gain of ice over open water : computed from
1239     C (SEAICE_areaGainFormula.EQ.1) from growth by ATM
1240     C (SEAICE_areaGainFormula.EQ.2) from predicted growth by ATM
1241     IF (SEAICE_areaGainFormula.EQ.1) THEN
1242     tmpscal4 = MAX(ZERO,d_HEFFbyATMonOCN_open(I,J))
1243     ELSE
1244     tmpscal4=MAX(ZERO,a_QbyATM_open(I,J))
1245     ENDIF
1246    
1247     C loss of ice cover by melting : computed from
1248     C (SEAICE_areaLossFormula.EQ.1) from all but only melt conributions by ATM and OCN
1249     C (SEAICE_areaLossFormula.EQ.2) from net melt-growth>0 by ATM and OCN
1250     C (SEAICE_areaLossFormula.EQ.3) from predicted melt by ATM
1251     IF (SEAICE_areaLossFormula.EQ.1) THEN
1252     tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN_cover(I,J) )
1253     & + MIN( 0. _d 0 , d_HEFFbyATMonOCN_open(I,J) )
1254     & + MIN( 0. _d 0 , d_HEFFbyOCNonICE(I,J) )
1255     ELSEIF (SEAICE_areaLossFormula.EQ.2) THEN
1256     tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN_cover(I,J)
1257     & + d_HEFFbyATMonOCN_open(I,J) + d_HEFFbyOCNonICE(I,J) )
1258     ELSE
1259     C compute heff after ice melt by ocn:
1260     tmpscal0=HEFF(I,J,bi,bj) - d_HEFFbyATMonOCN(I,J)
1261     C compute available heat left after snow melt by atm:
1262     tmpscal1= a_QbyATM_open(I,J)+a_QbyATM_cover(I,J)
1263     & - d_HSNWbyATMonSNW(I,J)*SNOW2ICE
1264     C could not melt more than all the ice
1265     tmpscal2 = MAX(-tmpscal0,tmpscal1)
1266     tmpscal3 = MIN(ZERO,tmpscal2)
1267     ENDIF
1268    
1269     C apply tendency
1270     IF ( (HEFF(i,j,bi,bj).GT.0. _d 0).OR.
1271     & (HSNOW(i,j,bi,bj).GT.0. _d 0) ) THEN
1272     AREA(I,J,bi,bj)=MAX(0. _d 0,
1273     & MIN( SEAICE_area_max, AREA(I,J,bi,bj)
1274     & + recip_HO*tmpscal4+HALF*recip_HH*tmpscal3 ))
1275     ELSE
1276     AREA(I,J,bi,bj)=0. _d 0
1277     ENDIF
1278     #ifdef ALLOW_SITRACER
1279     SItrAREA(I,J,bi,bj,3)=AREA(I,J,bi,bj)
1280     #endif /* ALLOW_SITRACER */
1281     #ifdef ALLOW_DIAGNOSTICS
1282     d_AREAbyATM(I,J)=
1283     & recip_HO*MAX(ZERO,d_HEFFbyATMonOCN_open(I,J))
1284     & +HALF*recip_HH*MIN(0. _d 0,d_HEFFbyATMonOCN_open(I,J))
1285     d_AREAbyICE(I,J)=
1286     & HALF*recip_HH*MIN(0. _d 0,d_HEFFbyATMonOCN_cover(I,J))
1287     d_AREAbyOCN(I,J)=
1288     & HALF*recip_HH*MIN( 0. _d 0,d_HEFFbyOCNonICE(I,J) )
1289     #endif /* ALLOW_DIAGNOSTICS */
1290     ENDDO
1291     ENDDO
1292    
1293     #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1294     Cgf 'bulk' linearization of area=f(HEFF)
1295     IF ( SEAICEadjMODE.GE.1 ) THEN
1296     DO J=1,sNy
1297     DO I=1,sNx
1298     C AREA(I,J,bi,bj) = 0.1 _d 0 * HEFF(I,J,bi,bj)
1299     AREA(I,J,bi,bj) = AREApreTH(I,J) + 0.1 _d 0 *
1300     & ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )
1301     ENDDO
1302     ENDDO
1303     ENDIF
1304     #endif
1305    
1306     C ===================================================================
1307     C =============PART 5: determine ice salinity increments=============
1308     C ===================================================================
1309    
1310     #ifndef SEAICE_VARIABLE_SALINITY
1311     # if (defined ALLOW_AUTODIFF_TAMC && defined ALLOW_SALT_PLUME)
1312     CADJ STORE d_HEFFbyNEG = comlev1_bibj,key=iicekey,byte=isbyte
1313     CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
1314     CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
1315     CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=iicekey,byte=isbyte
1316     CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=iicekey,byte=isbyte
1317     CADJ STORE d_HEFFbyFLOODING = comlev1_bibj,key=iicekey,byte=isbyte
1318     CADJ STORE d_HEFFbySublim = comlev1_bibj,key=iicekey,byte=isbyte
1319     CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
1320     CADJ & key = iicekey, byte = isbyte
1321     # endif /* ALLOW_AUTODIFF_TAMC and ALLOW_SALT_PLUME */
1322     DO J=1,sNy
1323     DO I=1,sNx
1324     tmpscal1 = d_HEFFbyNEG(I,J) + d_HEFFbyOCNonICE(I,J) +
1325     & d_HEFFbyATMonOCN(I,J) + d_HEFFbyFLOODING(I,J)
1326     & + d_HEFFbySublim(I,J)
1327     #ifdef SEAICE_ALLOW_AREA_RELAXATION
1328     + d_HEFFbyRLX(I,J)
1329     #endif
1330     tmpscal2 = tmpscal1 * SEAICE_salt0 * HEFFM(I,J,bi,bj)
1331     & * recip_deltaTtherm * SEAICE_rhoIce
1332     saltFlux(I,J,bi,bj) = tmpscal2
1333     #ifdef ALLOW_SALT_PLUME
1334     tmpscal3 = tmpscal1*salt(I,J,kSurface,bi,bj)*HEFFM(I,J,bi,bj)
1335     & * recip_deltaTtherm * SEAICE_rhoIce
1336     saltPlumeFlux(I,J,bi,bj) = MAX( tmpscal3-tmpscal2 , 0. _d 0)
1337     & *SPsalFRAC
1338     #endif /* ALLOW_SALT_PLUME */
1339     ENDDO
1340     ENDDO
1341     #endif /* ndef SEAICE_VARIABLE_SALINITY */
1342    
1343     #ifdef SEAICE_VARIABLE_SALINITY
1344    
1345     #ifdef SEAICE_GROWTH_LEGACY
1346     # ifdef ALLOW_AUTODIFF_TAMC
1347     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1348     # endif /* ALLOW_AUTODIFF_TAMC */
1349     DO J=1,sNy
1350     DO I=1,sNx
1351     C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean
1352     IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN
1353     saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
1354     & HSALT(I,J,bi,bj) * recip_deltaTtherm
1355     HSALT(I,J,bi,bj) = 0.0 _d 0
1356     ENDIF
1357     ENDDO
1358     ENDDO
1359     #endif /* SEAICE_GROWTH_LEGACY */
1360    
1361     #ifdef ALLOW_AUTODIFF_TAMC
1362     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1363     #endif /* ALLOW_AUTODIFF_TAMC */
1364    
1365     DO J=1,sNy
1366     DO I=1,sNx
1367     C sum up the terms that affect the salt content of the ice pack
1368     tmpscal1=d_HEFFbyOCNonICE(I,J)+d_HEFFbyATMonOCN(I,J)
1369    
1370     C recompute HEFF before thermodynamic updates (which is not AREApreTH in legacy code)
1371     tmpscal2=HEFF(I,J,bi,bj)-tmpscal1-d_HEFFbyFLOODING(I,J)
1372     C tmpscal1 > 0 : m of sea ice that is created
1373     IF ( tmpscal1 .GE. 0.0 ) THEN
1374     saltFlux(I,J,bi,bj) =
1375     & HEFFM(I,J,bi,bj)*recip_deltaTtherm
1376     & *SEAICE_saltFrac*salt(I,J,kSurface,bi,bj)
1377     & *tmpscal1*SEAICE_rhoIce
1378     #ifdef ALLOW_SALT_PLUME
1379     C saltPlumeFlux is defined only during freezing:
1380     saltPlumeFlux(I,J,bi,bj)=
1381     & HEFFM(I,J,bi,bj)*recip_deltaTtherm
1382     & *(ONE-SEAICE_saltFrac)*salt(I,J,kSurface,bi,bj)
1383     & *tmpscal1*SEAICE_rhoIce
1384     & *SPsalFRAC
1385     C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
1386     IF ( .NOT. SaltPlumeSouthernOcean ) THEN
1387     IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )
1388     & saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
1389     ENDIF
1390     #endif /* ALLOW_SALT_PLUME */
1391    
1392     C tmpscal1 < 0 : m of sea ice that is melted
1393     ELSE
1394     saltFlux(I,J,bi,bj) =
1395     & HEFFM(I,J,bi,bj)*recip_deltaTtherm
1396     & *HSALT(I,J,bi,bj)
1397     & *tmpscal1/tmpscal2
1398     #ifdef ALLOW_SALT_PLUME
1399     saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
1400     #endif /* ALLOW_SALT_PLUME */
1401     ENDIF
1402     C update HSALT based on surface saltFlux
1403     HSALT(I,J,bi,bj) = HSALT(I,J,bi,bj) +
1404     & saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
1405     saltFlux(I,J,bi,bj) =
1406     & saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
1407     #ifdef SEAICE_GROWTH_LEGACY
1408     C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean
1409     IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN
1410     saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -
1411     & HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) * recip_deltaTtherm
1412     HSALT(I,J,bi,bj) = 0.0 _d 0
1413     #ifdef ALLOW_SALT_PLUME
1414     saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
1415     #endif /* ALLOW_SALT_PLUME */
1416     ENDIF
1417     #endif /* SEAICE_GROWTH_LEGACY */
1418     ENDDO
1419     ENDDO
1420     #endif /* SEAICE_VARIABLE_SALINITY */
1421    
1422    
1423     C =======================================================================
1424     C ==LEGACY PART 6 (LEGACY) treat pathological cases, then do flooding ===
1425     C =======================================================================
1426    
1427     #ifdef SEAICE_GROWTH_LEGACY
1428    
1429     C treat values of ice cover fraction oustide
1430     C the [0 1] range, and other such issues.
1431     C ===========================================
1432    
1433     Cgf note: this part cannot be heat and water conserving
1434    
1435     #ifdef ALLOW_AUTODIFF_TAMC
1436     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
1437     CADJ & key = iicekey, byte = isbyte
1438     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
1439     CADJ & key = iicekey, byte = isbyte
1440     #endif /* ALLOW_AUTODIFF_TAMC */
1441     DO J=1,sNy
1442     DO I=1,sNx
1443     C NOW SET AREA(I,J,bi,bj)=0 WHERE THERE IS NO ICE
1444     CML replaced "/.0001 _d 0" by "*1. _d 4", 1e-4 is probably
1445     CML meant to be something like a minimum thickness
1446     AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),HEFF(I,J,bi,bj)*1. _d 4)
1447     ENDDO
1448     ENDDO
1449    
1450     #ifdef ALLOW_AUTODIFF_TAMC
1451     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
1452     CADJ & key = iicekey, byte = isbyte
1453     #endif /* ALLOW_AUTODIFF_TAMC */
1454     DO J=1,sNy
1455     DO I=1,sNx
1456     C NOW TRUNCATE AREA
1457     AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj))
1458     ENDDO
1459     ENDDO
1460    
1461     #ifdef ALLOW_AUTODIFF_TAMC
1462     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
1463     CADJ & key = iicekey, byte = isbyte
1464     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
1465     CADJ & key = iicekey, byte = isbyte
1466     #endif /* ALLOW_AUTODIFF_TAMC */
1467     DO J=1,sNy
1468     DO I=1,sNx
1469     AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj))
1470     HSNOW(I,J,bi,bj) = MAX(ZERO,HSNOW(I,J,bi,bj))
1471     AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1472     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1473     #ifdef SEAICE_CAP_HEFF
1474     C This is not energy conserving, but at least it conserves fresh water
1475     tmpscal0 = -MAX(HEFF(I,J,bi,bj)-MAX_HEFF,0. _d 0)
1476     d_HEFFbyNeg(I,J) = d_HEFFbyNeg(I,J) + tmpscal0
1477     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal0
1478     #endif /* SEAICE_CAP_HEFF */
1479     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1480     ENDDO
1481     ENDDO
1482    
1483     C convert snow to ice if submerged.
1484     C =================================
1485    
1486     IF ( SEAICEuseFlooding ) THEN
1487     DO J=1,sNy
1488     DO I=1,sNx
1489     tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1490     & +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
1491     tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
1492     d_HEFFbyFLOODING(I,J)=tmpscal1
1493     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
1494     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1495     & d_HEFFbyFLOODING(I,J)*ICE2SNOW
1496     ENDDO
1497     ENDDO
1498     ENDIF
1499    
1500     #endif /* SEAICE_GROWTH_LEGACY */
1501    
1502     #ifdef ALLOW_SITRACER
1503     DO J=1,sNy
1504     DO I=1,sNx
1505     c needs to be here to allow use also with LEGACY branch
1506     SItrHEFF(I,J,bi,bj,5)=HEFF(I,J,bi,bj)
1507     ENDDO
1508     ENDDO
1509     #endif /* ALLOW_SITRACER */
1510    
1511     C ===================================================================
1512     C ==============PART 7: determine ocean model forcing================
1513     C ===================================================================
1514    
1515     C compute net heat flux leaving/entering the ocean,
1516     C accounting for the part used in melt/freeze processes
1517     C =====================================================
1518    
1519     #ifdef ALLOW_AUTODIFF_TAMC
1520     CADJ STORE d_hsnwbyneg = comlev1_bibj,key=iicekey,byte=isbyte
1521     CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=iicekey,byte=isbyte
1522     #endif /* ALLOW_AUTODIFF_TAMC */
1523    
1524     DO J=1,sNy
1525     DO I=1,sNx
1526     QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + r_QbyATM_open(I,J)
1527     #ifndef SEAICE_GROWTH_LEGACY
1528     C in principle a_QSWbyATM_cover should always be included here, however
1529     C for backward compatibility it is left out of the LEGACY branch
1530     & + a_QSWbyATM_cover(I,J)
1531     #endif /* SEAICE_GROWTH_LEGACY */
1532     & - ( d_HEFFbyOCNonICE(I,J) +
1533     & d_HSNWbyOCNonSNW(I,J)*SNOW2ICE +
1534     & d_HEFFbyNEG(I,J) +
1535     #ifdef SEAICE_ALLOW_AREA_RELAXATION
1536     & d_HEFFbyRLX(I,J) +
1537     #endif
1538     & d_HSNWbyNEG(I,J)*SNOW2ICE )
1539     & * maskC(I,J,kSurface,bi,bj)
1540     QSW(I,J,bi,bj) = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)
1541     ENDDO
1542     ENDDO
1543    
1544     C switch heat fluxes from 'effective' ice meters to W/m2
1545     C ======================================================
1546    
1547     DO J=1,sNy
1548     DO I=1,sNx
1549     QNET(I,J,bi,bj) = QNET(I,J,bi,bj)*convertHI2Q
1550     QSW(I,J,bi,bj) = QSW(I,J,bi,bj)*convertHI2Q
1551     ENDDO
1552     ENDDO
1553    
1554     #ifndef SEAICE_DISABLE_HEATCONSFIX
1555     C treat advective heat flux by ocean to ice water exchange (at 0decC)
1556     C ===================================================================
1557     # ifdef ALLOW_AUTODIFF_TAMC
1558     CADJ STORE d_HEFFbyNEG = comlev1_bibj,key=iicekey,byte=isbyte
1559     CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
1560     CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
1561     CADJ STORE d_HSNWbyNEG = comlev1_bibj,key=iicekey,byte=isbyte
1562     CADJ STORE d_HSNWbyOCNonSNW = comlev1_bibj,key=iicekey,byte=isbyte
1563     CADJ STORE d_HSNWbyATMonSNW = comlev1_bibj,key=iicekey,byte=isbyte
1564     CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
1565     CADJ & key = iicekey, byte = isbyte
1566     # endif /* ALLOW_AUTODIFF_TAMC */
1567     IF ( SEAICEheatConsFix ) THEN
1568     c Unlike for evap and precip, the temperature of gained/lost
1569     c ocean liquid water due to melt/freeze of solid water cannot be chosen
1570     c to be e.g. the ocean SST. It must be done at 0degC. The fix below anticipates
1571     c on external_forcing_surf.F and applies the correction to QNET.
1572     IF ((convertFW2Salt.EQ.-1.).OR.(temp_EvPrRn.EQ.UNSET_RL)) THEN
1573     c I leave alone the exotic case when onvertFW2Salt.NE.-1 and temp_EvPrRn.NE.UNSET_RL and
1574     c the small error of the synchronous time stepping case (see external_forcing_surf.F).
1575     DO J=1,sNy
1576     DO I=1,sNx
1577     #ifdef ALLOW_DIAGNOSTICS
1578     c store unaltered QNET for diagnostic purposes
1579     DIAGarrayA(I,J)=QNET(I,J,bi,bj)
1580     #endif
1581     c compute the ocean water going to ice/snow, in precip units
1582     tmpscal3=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
1583     & ( d_HSNWbyATMonSNW(I,J)*SNOW2ICE
1584     & + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
1585     & + d_HEFFbyOCNonICE(I,J) + d_HEFFbyATMonOCN(I,J)
1586     & + d_HEFFbyNEG(I,J) + d_HSNWbyNEG(I,J)*SNOW2ICE )
1587     & * convertHI2PRECIP
1588     c factor in the heat content that external_forcing_surf.F
1589     c will associate with EMPMR, and remove it from QNET, so that
1590     c melt/freez water is in effect consistently gained/lost at 0degC
1591     IF (temp_EvPrRn.NE.UNSET_RL) THEN
1592     QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*
1593     & HeatCapacity_Cp * temp_EvPrRn
1594     ELSE
1595     QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*
1596     & HeatCapacity_Cp * theta(I,J,kSurface,bi,bj)
1597     ENDIF
1598     #ifdef ALLOW_DIAGNOSTICS
1599     c back out the eventual TFLUX adjustement and fill diag
1600     DIAGarrayA(I,J)=QNET(I,J,bi,bj)-DIAGarrayA(I,J)
1601     #endif
1602     ENDDO
1603     ENDDO
1604     ENDIF
1605     #ifdef ALLOW_DIAGNOSTICS
1606     CALL DIAGNOSTICS_FILL(DIAGarrayA,
1607     & 'SIaaflux',0,1,3,bi,bj,myThid)
1608     #endif
1609     ENDIF
1610     #endif /* ndef SEAICE_DISABLE_HEATCONSFIX */
1611    
1612     C compute net fresh water flux leaving/entering
1613     C the ocean, accounting for fresh/salt water stocks.
1614     C ==================================================
1615    
1616     #ifdef ALLOW_ATM_TEMP
1617     DO J=1,sNy
1618     DO I=1,sNx
1619     tmpscal1= d_HSNWbyATMonSNW(I,J)*SNOW2ICE
1620     & +d_HFRWbyRAIN(I,J)
1621     & +d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
1622     & +d_HEFFbyOCNonICE(I,J)
1623     & +d_HEFFbyATMonOCN(I,J)
1624     & +d_HEFFbyNEG(I,J)
1625     #ifdef SEAICE_ALLOW_AREA_RELAXATION
1626     & +d_HEFFbyRLX(I,J)
1627     #endif
1628     & +d_HSNWbyNEG(I,J)*SNOW2ICE
1629     C If r_FWbySublim>0, then it is evaporated from ocean.
1630     & +r_FWbySublim(I,J)
1631     EmPmR(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
1632     & ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
1633     & * ( ONE - AREApreTH(I,J) )
1634     #ifdef ALLOW_RUNOFF
1635     & - RUNOFF(I,J,bi,bj)
1636     #endif /* ALLOW_RUNOFF */
1637     & + tmpscal1*convertHI2PRECIP
1638     & )*rhoConstFresh
1639     ENDDO
1640     ENDDO
1641    
1642     #ifdef ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION
1643     C--
1644     DO J=1,sNy
1645     DO I=1,sNx
1646     frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
1647     & PRECIP(I,J,bi,bj)
1648     & - EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
1649     # ifdef ALLOW_RUNOFF
1650     & + RUNOFF(I,J,bi,bj)
1651     # endif /* ALLOW_RUNOFF */
1652     & )*rhoConstFresh
1653     # ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
1654     & - a_FWbySublim(I,J)*AREApreTH(I,J)
1655     # endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
1656     ENDDO
1657     ENDDO
1658     C--
1659     #else /* ndef ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION */
1660     C--
1661     # ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
1662     DO J=1,sNy
1663     DO I=1,sNx
1664     frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
1665     & PRECIP(I,J,bi,bj)
1666     & - EVAP(I,J,bi,bj)
1667     & *( ONE - AREApreTH(I,J) )
1668     # ifdef ALLOW_RUNOFF
1669     & + RUNOFF(I,J,bi,bj)
1670     # endif /* ALLOW_RUNOFF */
1671     & )*rhoConstFresh
1672     & - a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
1673     ENDDO
1674     ENDDO
1675     # endif
1676     C--
1677     #endif /* ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION */
1678    
1679     #endif /* ALLOW_ATM_TEMP */
1680    
1681     #ifdef SEAICE_DEBUG
1682     CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
1683     CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
1684     CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
1685     #endif /* SEAICE_DEBUG */
1686    
1687     C Sea Ice Load on the sea surface.
1688     C =================================
1689    
1690     #ifdef ALLOW_AUTODIFF_TAMC
1691     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1692     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1693     #endif /* ALLOW_AUTODIFF_TAMC */
1694    
1695     IF ( useRealFreshWaterFlux ) THEN
1696     DO J=1,sNy
1697     DO I=1,sNx
1698     #ifdef SEAICE_CAP_ICELOAD
1699     tmpscal1 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
1700     & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1701     tmpscal2 = MIN(tmpscal1,heffTooHeavy*rhoConst)
1702     #else
1703     tmpscal2 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
1704     & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1705     #endif
1706     sIceLoad(i,j,bi,bj) = tmpscal2
1707     ENDDO
1708     ENDDO
1709     ENDIF
1710    
1711     C ===================================================================
1712     C ======================PART 8: diagnostics==========================
1713     C ===================================================================
1714    
1715     #ifdef ALLOW_DIAGNOSTICS
1716     IF ( useDiagnostics ) THEN
1717     tmpscal1=1. _d 0 * recip_deltaTtherm
1718     CALL DIAGNOSTICS_SCALE_FILL(a_QbyATM_cover,
1719     & tmpscal1,1,'SIaQbATC',0,1,3,bi,bj,myThid)
1720     CALL DIAGNOSTICS_SCALE_FILL(a_QbyATM_open,
1721     & tmpscal1,1,'SIaQbATO',0,1,3,bi,bj,myThid)
1722     CALL DIAGNOSTICS_SCALE_FILL(a_QbyOCN,
1723     & tmpscal1,1,'SIaQbOCN',0,1,3,bi,bj,myThid)
1724     CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyOCNonICE,
1725     & tmpscal1,1,'SIdHbOCN',0,1,3,bi,bj,myThid)
1726     CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyATMonOCN_cover,
1727     & tmpscal1,1,'SIdHbATC',0,1,3,bi,bj,myThid)
1728     CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyATMonOCN_open,
1729     & tmpscal1,1,'SIdHbATO',0,1,3,bi,bj,myThid)
1730     CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyFLOODING,
1731     & tmpscal1,1,'SIdHbFLO',0,1,3,bi,bj,myThid)
1732     CALL DIAGNOSTICS_SCALE_FILL(d_HSNWbyOCNonSNW,
1733     & tmpscal1,1,'SIdSbOCN',0,1,3,bi,bj,myThid)
1734     CALL DIAGNOSTICS_SCALE_FILL(d_HSNWbyATMonSNW,
1735     & tmpscal1,1,'SIdSbATC',0,1,3,bi,bj,myThid)
1736     CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyATM,
1737     & tmpscal1,1,'SIdAbATO',0,1,3,bi,bj,myThid)
1738     CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyICE,
1739     & tmpscal1,1,'SIdAbATC',0,1,3,bi,bj,myThid)
1740     CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyOCN,
1741     & tmpscal1,1,'SIdAbOCN',0,1,3,bi,bj,myThid)
1742     CALL DIAGNOSTICS_SCALE_FILL(r_QbyATM_open,
1743     & convertHI2Q,1, 'SIqneto ',0,1,3,bi,bj,myThid)
1744     CALL DIAGNOSTICS_SCALE_FILL(r_QbyATM_cover,
1745     & convertHI2Q,1, 'SIqneti ',0,1,3,bi,bj,myThid)
1746     C three that actually need intermediate storage
1747     DO J=1,sNy
1748     DO I=1,sNx
1749     DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
1750     & * d_HSNWbyRAIN(I,J)*SEAICE_rhoSnow*recip_deltaTtherm
1751     DIAGarrayB(I,J) = AREA(I,J,bi,bj)-AREApreTH(I,J)
1752     ENDDO
1753     ENDDO
1754     CALL DIAGNOSTICS_FILL(DIAGarrayA,
1755     & 'SIsnPrcp',0,1,3,bi,bj,myThid)
1756     CALL DIAGNOSTICS_SCALE_FILL(DIAGarrayB,
1757     & tmpscal1,1,'SIdA ',0,1,3,bi,bj,myThid)
1758     #ifdef ALLOW_ATM_TEMP
1759     DO J=1,sNy
1760     DO I=1,sNx
1761     CML If I consider the atmosphere above the ice, the surface flux
1762     CML which is relevant for the air temperature dT/dt Eq
1763     CML accounts for sensible and radiation (with different treatment
1764     CML according to wave-length) fluxes but not for "latent heat flux",
1765     CML since it does not contribute to heating the air.
1766     CML So this diagnostic is only good for heat budget calculations within
1767     CML the ice-ocean system.
1768     DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
1769     #ifndef SEAICE_GROWTH_LEGACY
1770     & a_QSWbyATM_cover(I,J) +
1771     #endif /* SEAICE_GROWTH_LEGACY */
1772     & a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
1773     C
1774     DIAGarrayB(I,J) = maskC(I,J,kSurface,bi,bj) *
1775     & a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
1776     C
1777     DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)*(
1778     & PRECIP(I,J,bi,bj)
1779     & - EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
1780     #ifdef ALLOW_RUNOFF
1781     & + RUNOFF(I,J,bi,bj)
1782     #endif /* ALLOW_RUNOFF */
1783     & )*rhoConstFresh
1784     & - a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
1785     ENDDO
1786     ENDDO
1787     CALL DIAGNOSTICS_FILL(DIAGarrayA,
1788     & 'SIatmQnt',0,1,3,bi,bj,myThid)
1789     CALL DIAGNOSTICS_FILL(DIAGarrayB,
1790     & 'SIfwSubl',0,1,3,bi,bj,myThid)
1791     CALL DIAGNOSTICS_FILL(DIAGarrayC,
1792     & 'SIatmFW ',0,1,3,bi,bj,myThid)
1793     C
1794     DO J=1,sNy
1795     DO I=1,sNx
1796     C the actual Freshwater flux of sublimated ice, >0 decreases ice
1797     DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
1798     & * (a_FWbySublim(I,J)-r_FWbySublim(I,J))
1799     & * SEAICE_rhoIce * recip_deltaTtherm
1800     c the residual Freshwater flux of sublimated ice
1801     DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)
1802     & * r_FWbySublim(I,J)
1803     & * SEAICE_rhoIce * recip_deltaTtherm
1804     C the latent heat flux
1805     tmpscal1= EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
1806     & + r_FWbySublim(I,J)*convertHI2PRECIP
1807     tmpscal2= ( a_FWbySublim(I,J)-r_FWbySublim(I,J) )
1808     & * convertHI2PRECIP
1809     tmpscal3= SEAICE_lhEvap+SEAICE_lhFusion
1810     DIAGarrayB(I,J) = -maskC(I,J,kSurface,bi,bj)*rhoConstFresh
1811     & * ( tmpscal1*SEAICE_lhEvap + tmpscal2*tmpscal3 )
1812     ENDDO
1813     ENDDO
1814     CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIacSubl',0,1,3,bi,bj,myThid)
1815     CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIrsSubl',0,1,3,bi,bj,myThid)
1816     CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIhl ',0,1,3,bi,bj,myThid)
1817    
1818     DO J=1,sNy
1819     DO I=1,sNx
1820     c compute ice/snow water going to atm, in precip units
1821     tmpscal1 = rhoConstFresh*maskC(I,J,kSurface,bi,bj)
1822     & * convertHI2PRECIP * ( - d_HSNWbyRAIN(I,J)*SNOW2ICE
1823     & + a_FWbySublim(I,J) - r_FWbySublim(I,J) )
1824     c compute ocean water going to atm, in precip units
1825     tmpscal2=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
1826     & ( ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
1827     & * ( ONE - AREApreTH(I,J) )
1828     #ifdef ALLOW_RUNOFF
1829     & - RUNOFF(I,J,bi,bj)
1830     #endif /* ALLOW_RUNOFF */
1831     & + ( d_HFRWbyRAIN(I,J) + r_FWbySublim(I,J) )
1832     & *convertHI2PRECIP )
1833     c factor in the advected specific energy (referenced to 0 for 0deC liquid water)
1834     tmpscal1= - tmpscal1*
1835     & ( -SEAICE_lhFusion + HeatCapacity_Cp * ZERO )
1836     IF (temp_EvPrRn.NE.UNSET_RL) THEN
1837     tmpscal2= - tmpscal2*
1838     & ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
1839     ELSE
1840     tmpscal2= - tmpscal2*
1841     & ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )
1842     ENDIF
1843     c add to SIatmQnt, leading to SItflux, which is analogous to TFLUX
1844     DIAGarrayA(I,J)=maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
1845     #ifndef SEAICE_GROWTH_LEGACY
1846     & a_QSWbyATM_cover(I,J) +
1847     #endif
1848     & a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
1849     & -tmpscal1-tmpscal2
1850     ENDDO
1851     ENDDO
1852     CALL DIAGNOSTICS_FILL(DIAGarrayA,
1853     & 'SItflux ',0,1,3,bi,bj,myThid)
1854     #endif /* ALLOW_ATM_TEMP */
1855    
1856     ENDIF
1857     #endif /* ALLOW_DIAGNOSTICS */
1858    
1859     C close bi,bj loops
1860     ENDDO
1861     ENDDO
1862    
1863     RETURN
1864     END
  ViewVC Help
Powered by ViewVC 1.1.22