/[MITgcm]/MITgcm_contrib/atnguyen/verification/lab_sea/code_ad/seaice_growth.F
ViewVC logotype

Annotation of /MITgcm_contrib/atnguyen/verification/lab_sea/code_ad/seaice_growth.F

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

Revision 1.1 - (hide annotations) (download)
Mon Jun 1 22:56:23 2015 UTC (10 years, 2 months ago) by atn
Branch: MAIN
CVS Tags: HEAD 
latest working lab_sea seaice adjoint using IFenty code
1 atn 1.1 C SEAICE GROWTH REBORN v02 2015/01/15
2    
3     #include "SEAICE_OPTIONS.h"
4     #ifdef ALLOW_EXF
5     # include "EXF_OPTIONS.h"
6     #endif
7     #ifdef ALLOW_SALT_PLUME
8     # include "SALT_PLUME_OPTIONS.h"
9     #endif
10     #ifdef ALLOW_AUTODIFF
11     # include "AUTODIFF_OPTIONS.h"
12     #endif
13    
14     CBOP
15     C !ROUTINE: SEAICE_GROWTH
16     C !INTERFACE:
17     SUBROUTINE SEAICE_GROWTH( myTime, myIter, myThid )
18     C !DESCRIPTION: \bv
19     C *==========================================================*
20     C | SUBROUTINE seaice_growth
21     C | o rebirth of seaice_growth_if
22     C *==========================================================*
23     C \ev
24    
25     C !USES:
26     IMPLICIT NONE
27     C === Global variables ===
28     #include "SIZE.h"
29     #include "EEPARAMS.h"
30     #include "PARAMS.h"
31     #include "DYNVARS.h"
32     #include "GRID.h"
33     #include "FFIELDS.h"
34     #include "SEAICE_SIZE.h"
35     #include "SEAICE_PARAMS.h"
36     #include "SEAICE.h"
37     #include "SEAICE_TRACER.h"
38     #ifdef ALLOW_EXF
39     # include "EXF_PARAM.h"
40     # include "EXF_FIELDS.h"
41     #endif
42     #ifdef ALLOW_SALT_PLUME
43     # include "SALT_PLUME.h"
44     #endif
45     #ifdef ALLOW_AUTODIFF_TAMC
46     # include "tamc.h"
47     #endif
48    
49     C !INPUT/OUTPUT PARAMETERS:
50     C === Routine arguments ===
51     C myTime :: Simulation time
52     C myIter :: Simulation timestep number
53     C myThid :: Thread no. that called this routine.
54     _RL myTime
55     INTEGER myIter, myThid
56     CEOP
57    
58     #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
59     C !FUNCTIONS:
60     #ifdef ALLOW_DIAGNOSTICS
61     LOGICAL DIAGNOSTICS_IS_ON
62     EXTERNAL DIAGNOSTICS_IS_ON
63     #endif
64    
65     C !LOCAL VARIABLES:
66     C === Local variables ===
67     C
68     C unit/sign convention:
69    
70    
71     C HEFF is effective Hice thickness (m3/m2)
72     C HSNOW is Heffective snow thickness (m3/m2)
73     C HSALT is Heffective salt content (g/m2)
74     C AREA is the seaice cover fraction (0<=AREA<=1)
75    
76    
77    
78     C i,j,bi,bj :: Loop counters
79     INTEGER i, j, bi, bj
80     C number of surface interface layer
81     INTEGER kSurface
82     C IT :: ice thickness category index (MULTICATEGORIES and ITD code)
83     INTEGER IT
84     C msgBuf :: Informational/error message buffer
85     #ifdef ALLOW_BALANCE_FLUXES
86     CHARACTER*(MAX_LEN_MBUF) msgBuf
87     #elif (defined (SEAICE_DEBUG))
88     CHARACTER*(MAX_LEN_MBUF) msgBuf
89     CHARACTER*12 msgBufForm
90     #endif
91     C constants
92     _RL tempFrz, ICE2SNOW, SNOW2ICE, surf_theta
93     _RL QI, QS, recip_QI
94     _RL RHOSW
95     _RL lhSublim
96    
97     C conversion factors to go from Q (W/m2) to HEFF (ice meters)
98     _RL convertQ2HI, convertHI2Q
99     C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
100     _RL convertPRECIP2HI, convertHI2PRECIP
101    
102    
103    
104     C wind speed square
105     _RL SPEED_SQ
106    
107     C Regularization values squared
108     _RL area_reg_sq, hice_reg_sq
109    
110     C pathological cases thresholds
111     _RL heffTooHeavy
112    
113     C Helper variables: reciprocal of some constants
114     _RL recip_multDim
115     _RL recip_deltaTtherm
116     _RL recip_rhoIce
117    
118     C local value (=1/HO or 1/HO_south)
119     _RL recip_HO
120    
121     C local value (=1/ice thickness)
122     _RL recip_HH
123    
124     #ifndef SEAICE_ITD
125     C facilitate multi-category snow implementation
126     _RL pFac, pFacSnow
127     #endif
128    
129     C temporary variables available for the various computations
130     _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
131    
132     #ifdef ALLOW_SITRACER
133     INTEGER iTr
134     #ifdef ALLOW_DIAGNOSTICS
135     CHARACTER*8 diagName
136     #endif
137    
138    
139     #endif /* ALLOW_SITRACER */
140     #ifdef ALLOW_AUTODIFF_TAMC
141     INTEGER ilockey
142     #endif
143    
144     C== local arrays ==
145     C-- TmixLoc :: ocean surface/mixed-layer temperature (in K)
146     _RL TmixLoc (1:sNx,1:sNy)
147    
148     #ifndef SEAICE_ITD
149     C actual ice thickness (with upper and lower limit)
150     _RL hiceActual (1:sNx,1:sNy)
151     C actual snow thickness
152     _RL hsnowActual (1:sNx,1:sNy)
153     #endif
154     C actual ice thickness (with lower limit only) Reciprocal
155     _RL recip_hiceActual (1:sNx,1:sNy)
156    
157     C AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
158     _RL AREApreTH (1:sNx,1:sNy)
159     _RL HEFFpreTH (1:sNx,1:sNy)
160     _RL HSNWpreTH (1:sNx,1:sNy)
161    
162    
163     C wind speed
164     _RL UG (1:sNx,1:sNy)
165    
166     C temporary variables available for the various computations
167     _RL tmparr1 (1:sNx,1:sNy)
168    
169     _RL ticeInMult (1:sNx,1:sNy,nITD)
170     _RL ticeOutMult (1:sNx,1:sNy,nITD)
171     _RL hiceActualMult (1:sNx,1:sNy,nITD)
172     _RL hsnowActualMult (1:sNx,1:sNy,nITD)
173    
174     _RL F_io_net_mult (1:sNx,1:sNy,nITD)
175     _RL F_ia_net_mult (1:sNx,1:sNy,nITD)
176     _RL F_ia_mult (1:sNx,1:sNy,nITD)
177     _RL QSWI_mult (1:sNx,1:sNy,nITD)
178     _RL FWsublim_mult (1:sNx,1:sNy,nITD)
179    
180     #ifdef ALLOW_DIAGNOSTICS
181     C Helper variables for diagnostics
182     _RL DIAGarrayA (1:sNx,1:sNy)
183     _RL DIAGarrayB (1:sNx,1:sNy)
184     _RL DIAGarrayC (1:sNx,1:sNy)
185     _RL DIAGarrayD (1:sNx,1:sNy)
186     #endif /* ALLOW_DIAGNOSTICS */
187    
188    
189     _RL QSWO_BELOW_FIRST_LAYER (1:sNx,1:sNy)
190     _RL QSWO_IN_FIRST_LAYER (1:sNx,1:sNy)
191     _RL QSWO (1:sNx,1:sNy)
192     _RL QSWI (1:sNx,1:sNy)
193    
194     C Sea ice growth rates (m/s)
195     _RL ActualNewTotalVolumeChange (1:sNx,1:sNy)
196     _RL ActualNewTotalSnowMelt (1:sNx,1:sNy)
197    
198     _RL NetExistingIceGrowthRate (1:sNx,1:sNy)
199     _RL IceGrowthRateUnderExistingIce (1:sNx,1:sNy)
200     _RL IceGrowthRateFromSurface (1:sNx,1:sNy)
201     _RL IceGrowthRateOpenWater (1:sNx,1:sNy)
202     _RL IceGrowthRateMixedLayer (1:sNx,1:sNy)
203    
204     _RL EnergyInNewTotalIceVolume (1:sNx,1:sNy)
205     _RL NetEnergyFluxOutOfOcean (1:sNx,1:sNy)
206     c
207     C The energy taken out of the ocean which is not converted
208     C to sea ice (Joules)
209     _RL ResidualEnergyOutOfOcean (1:sNx,1:sNy)
210    
211     C Snow accumulation rate over ice (m/s)
212     _RL SnowAccRateOverIce (1:sNx,1:sNy)
213    
214     C Total snow accumulation over ice (m)
215     _RL SnowAccOverIce (1:sNx,1:sNy)
216    
217     C The precipitation rate over the ice which goes immediately into the ocean
218     _RL PrecipRateOverIceSurfaceToSea (1:sNx,1:sNy)
219    
220     C The potential snow melt rate if all snow surface heat flux convergences
221     C goes to melting snow (m/s)
222     _RL PotSnowMeltRateFromSurf (1:sNx,1:sNy)
223    
224     C The potential thickness of snow which could be melted by snow surface
225     C heat flux convergence (m)
226     _RL PotSnowMeltFromSurf (1:sNx,1:sNy)
227    
228     C The actual snow melt rate due to snow surface heat flux convergence
229     _RL SnowMeltRateFromSurface (1:sNx,1:sNy)
230    
231     C The actual surface heat flux convergence used to melt snow (W/m^2)
232     _RL SurfHeatFluxConvergToSnowMelt (1:sNx,1:sNy)
233    
234     C The actual thickness of snow to be melted by snow surface
235     C heat flux convergence (m)
236     _RL SnowMeltFromSurface (1:sNx,1:sNy)
237    
238     C The freshwater contribution to the ocean from melting snow (m)
239     _RL FreshwaterContribFromSnowMelt (1:sNx,1:sNy)
240    
241     C The freshwater contribution to (from) the ocean from melting (growing) ice (m)
242     _RL FreshwaterContribFromIce (1:sNx,1:sNy)
243    
244     C S_a : d(AREA)/dt
245     _RL S_a (1:sNx,1:sNy)
246    
247     C S_h : d(HEFF)/dt
248     _RL S_h (1:sNx,1:sNy)
249    
250     C S_hsnow : d(HSNOW)/dt
251     _RL S_hsnow (1:sNx,1:sNy)
252    
253    
254     C F_ia - sea ice/snow surface heat flux with atmosphere (W/m^2)
255     C F_ia > 0, heat loss to atmosphere
256     C F_ia < 0, atmospheric heat flux convergence (ice/snow surface melt)
257     _RL F_ia (1:sNx,1:sNy)
258    
259     C F_ia_net - the net heat flux divergence at the sea ice/snow surface
260     C including sea ice conductive fluxes and atmospheric fluxes (W/m^2)
261     C F_ia_net = 0, sea ice/snow surface energy balance condition met
262     C upward conductive fluxes balance surface heat loss
263     C F_ia_net < 0, net heat flux convergence at ice/snow surface
264     C zero conductive fluxes and net atmospheric convergence
265     _RL F_ia_net (1:sNx,1:sNy)
266    
267     C F_ia_net - the net heat flux divergence at the sea ice/snow surface
268     C before snow is melted with any convergence (W/m^2)
269     C F_ia_net < 0, some snow, if present, will melt
270     _RL F_ia_net_before_snow (1:sNx,1:sNy)
271    
272     C F_io_net - the net upward conductive heat flux through the ice+snow system
273     C realized at the sea ice/snow surface
274     C F_io_net > 0, heat conducting upward from ice base --> basal thickening
275     C F_io_net = 0, no upward heat conduction
276     C ice/snow surface temperature > SEAICE_freeze)
277     _RL F_io_net (1:sNx,1:sNy)
278    
279     C F_ao - heat flux from atmosphere to ocean (W/m^2)
280     C F_ao > 0
281     _RL F_ao (1:sNx,1:sNy)
282    
283     C F_mi - heat flux from ocean to the ice (W/m^2)
284     _RL F_mi (1:sNx,1:sNy)
285    
286     _RL FWsublim (1:sNx,1:sNy)
287    
288     C S_a_from_IGROW : d(AREA)/dt [from ice growth rate from open water fluxes]
289     _RL S_a_IGROW (1:sNx,1:sNy)
290    
291     C S_a_from_IGROW : d(AREA)/dt [from ice growth rate from ocean-ice fluxes]
292     _RL S_a_IGRML (1:sNx,1:sNy)
293    
294     C S_a_from_IGROW : d(AREA)/dt [from ice growth rate from ice-atm fluxes]
295     _RL S_a_IGRNE (1:sNx,1:sNy)
296    
297     C Factor by which we increase the upper ocean friction velocity (u*) when
298    
299     C ice is absent in a grid cell (dimensionless)
300     _RL MLTF (1:sNx,1:sNy)
301    
302     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
303    
304     C ===================================================================
305     C =================PART 0: constants and initializations=============
306     C ===================================================================
307    
308     IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
309     kSurface = Nr
310     ELSE
311     kSurface = 1
312     ENDIF
313    
314     C avoid unnecessary divisions in loops
315     recip_multDim = SEAICE_multDim
316     recip_multDim = ONE / recip_multDim
317     C above/below: double/single precision calculation of recip_multDim
318     c recip_multDim = 1./float(SEAICE_multDim)
319     recip_deltaTtherm = ONE / SEAICE_deltaTtherm
320     recip_rhoIce = ONE / SEAICE_rhoIce
321    
322     C Cutoff for iceload
323     heffTooHeavy=drF(kSurface) / 5. _d 0
324    
325     C RATIO OF SEA ICE DENSITY to SNOW DENSITY
326     ICE2SNOW = SEAICE_rhoIce/SEAICE_rhoSnow
327     SNOW2ICE = ONE / ICE2SNOW
328    
329     C Heat Flux * QI = [W] [m^3/kg] [kg/J] = [m^3/s] or [m/s] per unit area
330     QI = ONE/(SEAICE_rhoIce*SEAICE_lhFusion)
331     C Heat Flux * QS = [W] [m^3/kg] [kg/J] = [m^3/s] or [m/s] per unit area
332     QS = ONE/(SEAICE_rhoSnow*SEAICE_lhFusion)
333    
334     c A reference seawater density (kg m^-3)
335     RHOSW = 1026. _d 0
336     C ICE LATENT HEAT CONSTANT
337     lhSublim = SEAICE_lhEvap + SEAICE_lhFusion
338    
339     C regularization constants
340     area_reg_sq = SEAICE_area_reg * SEAICE_area_reg
341     hice_reg_sq = SEAICE_hice_reg * SEAICE_hice_reg
342    
343    
344     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
345     DO bj=myByLo(myThid),myByHi(myThid)
346     DO bi=myBxLo(myThid),myBxHi(myThid)
347    
348     #ifdef ALLOW_AUTODIFF_TAMC
349     act1 = bi - myBxLo(myThid)
350     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
351     act2 = bj - myByLo(myThid)
352     max2 = myByHi(myThid) - myByLo(myThid) + 1
353     act3 = myThid - 1
354     max3 = nTx*nTy
355     act4 = ikey_dynamics - 1
356     iicekey = (act1 + 1) + act2*max1
357     & + act3*max1*max2
358     & + act4*max1*max2*max3
359     #endif /* ALLOW_AUTODIFF_TAMC */
360    
361     cph-test(
362     #ifdef ALLOW_AUTODIFF_TAMC
363     CADJ STORE area = comlev1_bibj, key = iicekey, byte = isbyte
364     CADJ STORE heff = comlev1_bibj, key = iicekey, byte = isbyte
365     CADJ STORE hsnow = comlev1_bibj, key = iicekey, byte = isbyte
366     CADJ STORE qnet = comlev1_bibj, key = iicekey, byte = isbyte
367     CADJ STORE qsw = comlev1_bibj, key = iicekey, byte = isbyte
368     CADJ STORE tices = comlev1_bibj, key = iicekey, byte = isbyte
369     #endif
370     cph-test)
371    
372     C array initializations
373     C =====================
374    
375     DO J=1,sNy
376     DO I=1,sNx
377    
378     C NEW VARIABLE NAMES
379     NetExistingIceGrowthRate(I,J) = 0.0 _d 0
380     IceGrowthRateOpenWater(I,J) = 0.0 _d 0
381     IceGrowthRateFromSurface(I,J) = 0.0 _d 0
382     IceGrowthRateMixedLayer(I,J) = 0.0 _d 0
383    
384     EnergyInNewTotalIceVolume(I,J) = 0.0 _d 0
385     NetEnergyFluxOutOfOcean(I,J) = 0.0 _d 0
386     ResidualEnergyOutOfOcean(I,J) = 0.0 _d 0
387    
388     PrecipRateOverIceSurfaceToSea(I,J) = 0.0 _d 0
389    
390     DO IT=1,SEAICE_multDim
391     ticeInMult(I,J,IT) = 0.0 _d 0
392     ticeOutMult(I,J,IT) = 0.0 _d 0
393    
394     F_io_net_mult(I,J,IT) = 0.0 _d 0
395     F_ia_net_mult(I,J,IT) = 0.0 _d 0
396     F_ia_mult(I,J,IT) = 0.0 _d 0
397    
398     QSWI_mult(I,J,IT) = 0.0 _d 0
399     FWsublim_mult(I,J,IT) = 0.0 _d 0
400     ENDDO
401    
402     F_io_net(I,J) = 0.0 _d 0
403     F_ia_net(I,J) = 0.0 _d 0
404     F_ia(I,J) = 0.0 _d 0
405    
406     QSWI(I,J) = 0.0 _d 0
407     FWsublim(I,J) = 0.0 _d 0
408    
409     QSWO_BELOW_FIRST_LAYER (I,J) = 0.0 _d 0
410     QSWO_IN_FIRST_LAYER (I,J) = 0.0 _d 0
411     QSWO(I,J) = 0.0 _d 0
412    
413     ActualNewTotalVolumeChange(I,J) = 0.0 _d 0
414     ActualNewTotalSnowMelt(I,J) = 0.0 _d 0
415    
416     SnowAccOverIce(I,J) = 0.0 _d 0
417     SnowAccRateOverIce(I,J) = 0.0 _d 0
418    
419     PotSnowMeltRateFromSurf(I,J) = 0.0 _d 0
420     PotSnowMeltFromSurf(I,J) = 0.0 _d 0
421     SnowMeltRateFromSurface(I,J) = 0.0 _d 0
422     SurfHeatFluxConvergToSnowMelt(I,J) = 0.0 _d 0
423     SnowMeltFromSurface(I,J) = 0.0 _d 0
424    
425     FreshwaterContribFromSnowMelt(I,J) = 0.0 _d 0
426     FreshwaterContribFromIce(I,J) = 0.0 _d 0
427    
428     S_a(I,J) = 0.0 _d 0
429     S_a_IGROW(I,J) = 0.0 _d 0
430     S_a_IGRML(I,J) = 0.0 _d 0
431     S_a_IGRNE(I,J) = 0.0 _d 0
432    
433     S_h(I,J) = 0.0 _d 0
434     S_hsnow(I,J) = 0.0 _d 0
435    
436     MLTF(I,J) = 0.0 _d 0
437     ENDDO
438     ENDDO
439    
440     #if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
441     DO J=1-oLy,sNy+oLy
442     DO I=1-oLx,sNx+oLx
443     frWtrAtm(I,J,bi,bj) = 0.0 _d 0
444     ENDDO
445     ENDDO
446     #endif
447    
448     C =====================================================================
449     C PART 1: Store ice and snow state on onset + regularize actual
450     C ice area and ice thickness
451     C =====================================================================
452    
453     DO J=1,sNy
454     DO I=1,sNx
455    
456     HEFF(I,J,bi,bj) = MAX(ZERO, HEFF(I,J,bi,bj))
457     AREA(I,J,bi,bj) = MAX(ZERO, AREA(I,J,bi,bj))
458    
459     c Cap the area if it exceedes area_max (may also have been
460     c done in do_ridging)
461     AREA(I,J,bi,bj) = MIN(AREA(I,J,bi,bj), SEAICE_area_max)
462    
463     IF (HEFF(I,J,bi,bj) .LE. ZERO) then
464     AREA(I,J, bi,bj) = ZERO
465     HSNOW(I,J, bi,bj) = ZERO
466     ELSEIF (AREA(I,J,bi,bj) .LE. ZERO) then
467     HEFF(I,J,bi,bj) = ZERO
468     HSNOW(I,J,bi,bj) = ZERO
469     ENDIF
470    
471     HEFFpreTH(I,J) = HEFF(I,J,bi,bj)
472     c HSNWpreTH(I,J) = HSNOW(I,J,bi,bj)
473     AREApreTH(I,J) = AREA(I,J,bi,bj)
474    
475     #ifdef ALLOW_DIAGNOSTICS
476     DIAGarrayB(I,J) = AREA(I,J,bi,bj)
477     DIAGarrayC(I,J) = HEFF(I,J,bi,bj)
478     DIAGarrayD(I,J) = HSNOW(I,J,bi,bj)
479     #endif
480     ENDDO
481     ENDDO
482    
483     #ifdef ALLOW_DIAGNOSTICS
484     IF ( useDiagnostics ) THEN
485     CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIareaGA',0,1,3,bi,bj,myThid)
486     CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIheffGA',0,1,3,bi,bj,myThid)
487     CALL DIAGNOSTICS_FILL(DIAGarrayD,'SIhsnoGA',0,1,3,bi,bj,myThid)
488     ENDIF
489     #endif /* ALLOW_DIAGNOSTICS */
490    
491     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
492     #ifdef ALLOW_AUTODIFF_TAMC
493     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
494     CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
495     CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
496     #endif /* ALLOW_AUTODIFF_TAMC */
497     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
498    
499     C COMPUTE ACTUAL ICE/SNOW THICKNESS; USE MIN/MAX VALUES
500     C TO REGULARIZE SEAICE_SOLVE4TEMP/d_AREA COMPUTATIONS
501    
502     DO J=1,sNy
503     DO I=1,sNx
504     IF (HEFFpreTH(I,J) .GT. ZERO) THEN
505     c regularize AREA with SEAICE_area_reg
506     tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)
507    
508     c hiceActual calculated with the regularized AREA
509     tmpscal2 = HEFFpreTH(I,J) / tmpscal1
510    
511     c regularize hiceActual with SEAICE_hice_reg (add lower bound)
512     c hiceActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
513     hiceActual(I,J) = MAX(0.05 _d 0, tmpscal2)
514    
515     c hsnowActual calculated with the regularized AREA (no lower bound)
516     c hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1
517     c actually I do not think we need to regularize this.
518     c hsnowActual(I,J) = HSNWpreTH(I,J) / AREApreTH(I,J)
519    
520     c regularize the inverse of hiceActual by hice_reg
521     recip_hiceActual(I,J) = AREApreTH(I,J) /
522     & sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)
523    
524     c Do not regularize when HEFFpreTH = 0
525     ELSE
526     hiceActual (I,J) = ZERO
527     hsnowActual(I,J) = ZERO
528     recip_hiceActual(I,J) = ZERO
529     ENDIF
530    
531     ENDDO
532     ENDDO
533    
534    
535     C =============================================================================
536     C Part 2: Precipitation as snow or rain over ice
537     C =============================================================================
538    
539     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
540     #ifdef ALLOW_AUTODIFF_TAMC
541     CADJ STORE SnowAccRateOverIce = comlev1_bibj,key=iicekey,byte=isbyte
542     CADJ STORE SnowAccOverIce = comlev1_bibj,key=iicekey,byte=isbyte
543     CADJ STORE PrecipRateOverIceSurfaceToSea =
544     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
545     CADJ STORE PRECIP(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
546     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
547     #endif /* ALLOW_AUTODIFF_TAMC */
548     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
549    
550     DO J=1,sNy
551     DO I=1,sNx
552     c if we have ice and the temperature of the ice is below the freezing point
553     c then the precip falls and accumulates as snow
554     #ifndef FROZEN
555     IF (( AREApreTH(I,J) .GT. ZERO) .AND.
556     & ( TICES(I,J,1,bi,bj) .LT. celsius2k) ) THEN
557     #else
558     IF ( AREApreTH(I,J) .GT. ZERO) THEN
559     #endif
560    
561     c use either prescribed snowfall or PRECIP rate
562     IF ( snowPrecipFile .NE. ' ' ) THEN
563     c rate of actual snow accumulation in m/s over ice
564     c y [m/s] \approx 1.0 [kg/m^3] / 0.9 [m^3/kg] * x [m/s]
565     SnowAccRateOverIce(I,J) = rhoConstFresh/SEAICE_rhoSnow *
566     & snowPrecip(i,j,bi,bj)
567     ELSE
568     SnowAccRateOverIce(I,J) = rhoConstFresh/SEAICE_rhoSnow *
569     & PRECIP(i,j,bi,bj)
570     ENDIF
571     PrecipRateOverIceSurfaceToSea(I,J) = ZERO
572    
573     ELSE
574     c The snow/ice surface is not frozen (wet) so the precipitation
575     c remains wet and runs into the ocean
576     SnowAccRateOverIce(I,J) = ZERO
577     PrecipRateOverIceSurfaceToSea(I,J) = PRECIP(i,j,bi,bj)
578     ENDIF
579    
580     c SnowAccOverIce is the change of mean snow thickness, i.e. HSNOW
581     c over one time step.
582     SnowAccOverIce(I,J) =
583     & SnowAccRateOverIce(I,J) * SEAICE_deltaTtherm * AREApreTH(I,J)
584     c I,J
585     ENDDO
586     ENDDO
587    
588     c acumulate snow on the surface
589     c before we do any surface energy balance calculations
590     DO J=1,sNy
591     DO I=1,sNx
592     IF ( AREApreTH(I,J) .GT. ZERO) THEN
593    
594     c update mean and actual snow thickness.
595     c to this point there are no any thermodynamic calculations,
596     c only potentially accumulated snow based on precip and the
597     c ice surface temp from the previous time step
598    
599     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + SnowAccOverIce(I,J)
600     HSNWpreTH(I,J) = HSNOW(I,J,bi,bj)
601     hsnowActual(I,J) = HSNWpreTH(I,J) / AREApreTH(I,J)
602     ENDIF
603     ENDDO
604     ENDDO
605    
606     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
607     #ifdef ALLOW_AUTODIFF_TAMC
608     CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
609     #endif /* ALLOW_AUTODIFF_TAMC */
610     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
611    
612     C =============================================================================
613     C FIND WIND SPEED
614     C =============================================================================
615    
616     DO j=1,sNy
617     DO i=1,sNx
618     C ocean surface/mixed layer temperature
619     TmixLoc(i,j) = theta(i,j,kSurface,bi,bj) + celsius2K
620     C wind speed from exf
621     UG(I,J) = MAX(SEAICE_EPS, wspeed(I,J,bi,bj))
622     ENDDO
623     ENDDO
624    
625    
626     C =============================================================================
627     c Retrieve the air-sea heat and shortwave radiative fluxes
628     C =============================================================================
629    
630     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
631     #ifdef ALLOW_AUTODIFF_TAMC
632     CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
633     CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
634     cCADJ STORE UG = comlev1_bibj, key = iicekey,byte=isbyte
635     cCADJ STORE TmixLoc = comlev1_bibj, key = iicekey,byte=isbyte
636     #endif /* ALLOW_AUTODIFF_TAMC */
637     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
638    
639     CALL SEAICE_BUDGET_OCEAN(
640     I UG,
641     I TmixLoc,
642     O F_ao, QSWO,
643     I bi, bj, myTime, myIter, myThid )
644    
645    
646     C =============================================================================
647     C Calc air-sea fluxes in the uppermost grid cell
648     C =============================================================================
649    
650     c-- Not all of the sw radiation is absorbed in the uppermost ocean grid cell layer.
651     c Only that fraction which converges in the uppermost ocean grid cell is used to
652     c melt ice.
653     c SWFRACB - the fraction of incoming sw radiation absorbed in the
654     c uppermost ocean grid cell (calculated in seaice_init_vari.F)
655     DO J=1,sNy
656     DO I=1,sNx
657    
658     c The contribution of shortwave heating is
659     c not included without #define SHORTWAVE_HEATING
660     #ifdef SHORTWAVE_HEATING
661     QSWO_BELOW_FIRST_LAYER(I,J)= QSWO(I,J)*SWFRACB
662     QSWO_IN_FIRST_LAYER(I,J) = QSWO(I,J)*(ONE - SWFRACB)
663     #else
664     QSWO_BELOW_FIRST_LAYER(I,J)= ZERO
665     QSWO_IN_FIRST_LAYER(I,J) = ZERO
666     #endif
667     IceGrowthRateOpenWater(I,J) = QI *
668     & (F_ao(I,J) - QSWO(I,J) + QSWO_IN_FIRST_LAYER(I,J))
669    
670     ENDDO
671     ENDDO
672    
673    
674    
675     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
676     #ifdef ALLOW_AUTODIFF_TAMC
677     CADJ STORE hsnowActual = comlev1_bibj, key = iicekey, byte = isbyte
678     CADJ STORE hiceActual = comlev1_bibj, key = iicekey, byte = isbyte
679     CADJ STORE UG = comlev1_bibj, key = iicekey, byte = isbyte
680     CADJ STORE tices(:,:,:,bi,bj)
681     CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
682     CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
683     CADJ & key = iicekey, byte = isbyte
684     #endif /* ALLOW_AUTODIFF_TAMC */
685     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
686    
687     C =============================================================================
688     c calculate heat fluxes within ice (conduction), F_io, and across the
689     c ice/atmosphere interface, F_ia
690     C =============================================================================
691    
692     C-- Start loop over multi-categories
693     DO IT=1,SEAICE_multDim
694    
695     DO J=1,sNy
696     DO I=1,sNx
697     c record prior ice surface temperatures
698     ticeInMult(I,J,IT) = TICES(I,J,IT,bi,bj)
699     ticeOutMult(I,J,IT) = TICES(I,J,IT,bi,bj)
700     TICES(I,J,IT,bi,bj) = ZERO
701     ENDDO
702     ENDDO
703    
704     c set relative thickness of ice categories
705     pFac = (2.0 _d 0*IT - 1.0 _d 0)*recip_multDim
706     pFacSnow = 1. _d 0
707    
708     c find actual snow and ice thickness within categories categories
709     IF ( SEAICE_useMultDimSnow ) pFacSnow=pFac
710    
711     DO J=1,sNy
712     DO I=1,sNx
713     hiceActualMult(I,J,IT) = hiceActual(I,J) *pFac
714     hsnowActualMult(I,J,IT) = hsnowActual(I,J)*pFacSnow
715     ENDDO
716     ENDDO
717    
718     ENDDO /* multDim */
719    
720     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
721     #ifdef ALLOW_AUTODIFF_TAMC
722     CADJ STORE hiceActualMult = comlev1_bibj, key = iicekey, byte = isbyte
723     CADJ STORE hsnowActualMult= comlev1_bibj, key = iicekey, byte = isbyte
724     CADJ STORE ticeOutMult = comlev1_bibj, key = iicekey, byte = isbyte
725     CADJ STORE F_io_net_mult =
726     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
727     CADJ STORE F_ia_net_mult =
728     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
729     CADJ STORE F_ia_mult =
730     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
731     CADJ STORE QSWI_mult =
732     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
733     CADJ STORE FWsublim_mult =
734     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
735     #endif /* ALLOW_AUTODIFF_TAMC */
736     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
737    
738     c =======================================================================
739     c find calculate heat fluxes within ice (conduction) and across the
740     c ice/atmosphere interface for each thickness category
741     c =======================================================================
742    
743     DO IT=1,SEAICE_multDim
744     CALL SEAICE_SOLVE4TEMP(
745     I UG, hiceActualMult(1,1,IT), hsnowActualMult(1,1,IT),
746     U ticeInMult(1,1,IT), ticeOutMult(1,1,IT),
747     O F_io_net_mult(1,1,IT),
748     O F_ia_net_mult(1,1,IT),
749     O F_ia_mult(1,1,IT),
750     O QSWI_mult(1,1,IT),
751     O FWsublim_mult(1,1,IT),
752     I bi, bj, myTime, myIter, myThid )
753     ENDDO
754    
755     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
756     #ifdef ALLOW_AUTODIFF_TAMC
757     CADJ STORE hiceActualMult = comlev1_bibj, key = iicekey, byte = isbyte
758     CADJ STORE hsnowActualMult= comlev1_bibj, key = iicekey, byte = isbyte
759     CADJ STORE ticeOutMult = comlev1_bibj, key = iicekey, byte = isbyte
760     CADJ STORE F_io_net_mult =
761     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
762     CADJ STORE F_ia_net_mult =
763     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
764     CADJ STORE F_ia_mult =
765     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
766     CADJ STORE QSWI_mult =
767     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
768     CADJ STORE FWsublim_mult =
769     CADJ & comlev1_bibj, key = iicekey, byte = isbyte
770     #endif /* ALLOW_AUTODIFF_TAMC */
771     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
772    
773    
774     c =======================================================================
775     c record the ice surface temperature in each category
776     c and find the average of fluxes across each category
777    
778     DO IT=1,SEAICE_multDim
779     DO J=1,sNy
780     DO I=1,sNx
781    
782     C update TICES
783     TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)
784    
785     F_io_net(I,J) = F_io_net(I,J) +
786     & F_io_net_mult(I,J,IT)*recip_multDim
787    
788     F_ia_net(I,J) = F_ia_net(I,J) +
789     & F_ia_net_mult(I,J,IT)*recip_multDim
790    
791     F_ia(I,J) = F_ia(I,J) +
792     & F_ia_mult(I,J,IT)*recip_multDim
793    
794     QSWI(I,J) = QSWI(I,J) + QSWI_mult(I,J,IT)*recip_multDim
795    
796     FWsublim(I,J) = FWsublim(I,J) +
797     & FWsublim_mult(I,J,IT)*recip_multDim
798    
799     ENDDO
800     ENDDO
801     ENDDO
802    
803     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
804     #ifdef ALLOW_AUTODIFF_TAMC
805     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
806     CADJ STORE F_io_net = comlev1_bibj, key = iicekey, byte = isbyte
807     CADJ STORE F_ia_net = comlev1_bibj, key = iicekey, byte = isbyte
808     CADJ STORE F_ia = comlev1_bibj, key = iicekey, byte = isbyte
809     CADJ STORE QSWI = comlev1_bibj, key = iicekey, byte = isbyte
810     CADJ STORE FWSublim = comlev1_bibj, key = iicekey, byte = isbyte
811     #endif /* ALLOW_AUTODIFF_TAMC */
812     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
813    
814     DO J=1,sNy
815     DO I=1,sNx
816     c If there is heat flux convergence at the snow surface,
817     c use that energy to melt snow before melting ice. It is
818     c possible that some snow will remain after melting,
819     c which will drive F_ia_net to zero, or that all of the
820     c snow will be melted, leaving a nonzero F_ia_net to melt
821     c some ice.
822     F_ia_net_before_snow(I,J) = F_ia_net(I,J)
823    
824     c Only continue if there is snow and ice in the cell
825     IF (AREApreTH(I,J) .LE. ZERO) THEN
826     IceGrowthRateUnderExistingIce(I,J) = 0. _d 0
827     IceGrowthRateFromSurface(I,J) = 0. _d 0
828     NetExistingIceGrowthRate(I,J) = 0. _d 0
829     ELSE
830     c The growth rate (m/s) beneath existing ice is given by the upward
831     c ocean-ice conductive flux, F_io_net, and QI.
832     IceGrowthRateUnderExistingIce(I,J) = F_io_net(I,J)*QI
833    
834     c The rate at which snow is melted (m/s) because of surface
835     c heat flux convergence. Note, during snow melt, F_ia_net must
836     c be negative (implying convergence) to make PSMRFW is positive
837     PotSnowMeltRateFromSurf(I,J) = - F_ia_net(I,J)*QS
838    
839     c This is the depth of snow (m) that would be melted in one dt
840     PotSnowMeltFromSurf(I,J) =
841     & PotSnowMeltRateFromSurf(I,J) * SEAICE_deltaTtherm
842    
843     c If we can melt MORE than is actually there, then the melt
844     c rate is reduced so that only that which is there
845     c is melted during the time step. In this case, not all of the
846     c heat flux convergence at the surface is used to melt snow.
847     c Any remaining energy will melt ice.
848    
849     c SurfHeatFluxConvergToSnowMelt is the part of the total heat
850     c flux convergence which melts snow.
851    
852     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
853     CCHECK: HSNOW ACTUAL SHOULDNT BE REGULARIZED FOR THIS
854     IF (PotSnowMeltFromSurf(I,J) .GE. hsnowActual(I,J)) THEN
855    
856     c Snow melt and melt rate [m] (actual snow thickness)
857     SnowMeltFromSurface(I,J) = hsnowActual(I,J)
858    
859     SnowMeltRateFromSurface(I,J) =
860     & SnowMeltFromSurface(I,J) / SEAICE_deltaTtherm
861    
862     SurfHeatFluxConvergToSnowMelt(I,J) =
863     & - hsnowActual(I,J)/QS/SEAICE_deltaTtherm
864     ELSE
865     c In this case there will be snow remaining after melting.
866     c All of the surface heat convergence will be redirected to
867     c this effort.
868     SnowMeltFromSurface(I,J) = PotSnowMeltFromSurf(I,J)
869    
870     SnowMeltRateFromSurface(I,J) =PotSnowMeltRateFromSurf(I,J)
871    
872     SurfHeatFluxConvergToSnowMelt(I,J) = F_ia_net(I,J)
873    
874     ENDIF
875    
876     c Reduce the heat flux convergence available to melt surface
877     c ice by the amount used to melt snow
878     F_ia_net(I,J) = F_ia_net(I,J) -
879     & SurfHeatFluxConvergToSnowMelt(I,J)
880    
881     IceGrowthRateFromSurface(I,J) = F_ia_net(I,J) * QI
882    
883     c The total growth rate (m/s) of the existing ice - the rate of
884     c new ice accretion at the base less the rate due to surface melt
885     NetExistingIceGrowthRate(I,J) =
886     & IceGrowthRateUnderExistingIce(I,J) +
887     & IceGrowthRateFromSurface(I,J)
888     ENDIF
889    
890     ENDDO
891     ENDDO
892    
893     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
894     #ifdef ALLOW_AUTODIFF_TAMC
895     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
896     CADJ STORE SnowMeltRateFromSurface = comlev1_bibj,
897     CADJ & key = iicekey, byte = isbyte
898     CADJ STORE IceGrowthRateUnderExistingIce = comlev1_bibj,
899     CADJ & key = iicekey, byte = isbyte
900     CADJ STORE IceGrowthRateFromSurface = comlev1_bibj,
901     CADJ & key = iicekey, byte = isbyte
902     CADJ STORE NetExistingIceGrowthRate = comlev1_bibj,
903     CADJ & key = iicekey, byte = isbyte
904     #endif /* ALLOW_AUTODIFF_TAMC */
905     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
906    
907    
908     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
909     #ifdef ALLOW_AUTODIFF_TAMC
910     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
911     CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
912     CADJ & key = iicekey, byte = isbyte
913     CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
914     CADJ & key = iicekey, byte = isbyte
915     #endif
916     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
917    
918     C Calcuate the heat fluxes from the ocean to the sea ice
919     DO J=1,sNy
920     DO I=1,sNx
921     c
922     c Bound the ocean temperature to be at or above the freezing point.
923     tempFrz = SEAICE_tempFrz0 +
924     & SEAICE_dTempFrz_dS * salt(I,J,kSurface,bi,bj)
925    
926     surf_theta = max(theta(I,J,kSurface,bi,bj), tempFrz)
927    
928     c inflection point
929     tmpscal0 = 0.4 _d 0
930    
931     c steepness
932     tmpscal1 = 7.0 _d 0
933    
934     MLTF(I,J) = ONE + (12.5 - ONE) *
935     & (ONE + EXP( (AREApreTH(I,J) - tmpscal0)
936     & * tmpscal1/tmpscal0))**(-ONE)
937    
938     c IF (AREApreTH(I,J) .GT. ZERO) THEN
939     c
940     c If ice is present, MixedLayerTurbulenceFactor = 1.0, else 12.50
941     c MLTF = ONE
942     c ELSE
943     c MLTF = 12.5 _d 0
944     c ENDIF
945    
946     F_mi(I,J) = -STANTON_NUMBER * USTAR_BASE * RHOSW *
947     & HeatCapacity_Cp * (surf_theta - tempFrz) * MLTF(I,J)
948    
949     IceGrowthRateMixedLayer (I,J) = F_mi(I,J) * QI
950    
951     ENDDO
952     ENDDO
953    
954    
955     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
956     #ifdef ALLOW_AUTODIFF_TAMC
957     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
958     CADJ STORE F_mi = comlev1_bibj, key = iicekey, byte = isbyte
959     CADJ STORE IceGrowthRateMixedLayer =
960     CADJ & comlev1_bibj,key = iicekey, byte = isbyte
961     #endif /* ALLOW_AUTODIFF_TAMC */
962     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
963    
964     C CALCULATE THICKNESS DERIVATIVES of ice (dhdt) and snow (dhsdt)
965     DO J=1,sNy
966     DO I=1,sNx
967    
968     S_h(I,J) =
969     & NetExistingIceGrowthRate(I,J) * AREApreTH(I,J)
970     & + IceGrowthRateOpenWater(I,J) * (ONE - AREApreTH(I,J))
971     & + IceGrowthRateMixedLayer(I,J)
972    
973     c Both the accumulation and melt rates are in terms
974     c of actual snow thickness. As with ice, multiplying
975     c with area converts to mean snow thickness.
976     c S_hsnow(I,J) = AREApreTH(I,J) *
977     c & ( SnowAccRateOverIce(I,J) - SnowMeltRateFromSurface(I,J))
978    
979     c the only snow tendency term is the surface melting term since
980     c the surface accumulation is taken care of in step 2
981     S_hsnow(I,J) = AREApreTH(I,J) *
982     & ( - SnowMeltRateFromSurface(I,J))
983    
984     ENDDO
985     ENDDO
986    
987     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
988     #ifdef ALLOW_AUTODIFF_TAMC
989     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
990     CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
991     CADJ STORE S_a = comlev1_bibj, key = iicekey, byte = isbyte
992     CADJ STORE S_h = comlev1_bibj, key = iicekey, byte = isbyte
993     CADJ STORE S_hsnow = comlev1_bibj, key = iicekey, byte = isbyte
994     #endif /* ALLOW_AUTODIFF_TAMC */
995     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
996    
997     c Caculate dA/dt (S_a)
998     DO J=1,sNy
999     DO I=1,sNx
1000    
1001     S_a(I,J) = 0. _d 0
1002    
1003     c Caculate the ice area growth rate from the open water fluxes.
1004     c First, determine whether the open water growth rate is positive or
1005     c negative. If positive, make sure that ice is present or that the
1006     c net ice thickness growth rate is positive before extending ice cover
1007    
1008     c this is the geometric term: Area/(2*Heff),
1009     c with Area/Heff regularized as Area/(Heff^2 + epsilon^2)
1010     c Area/(Heff^2 + ep^2) = recip_hiceActual
1011     c epsilon = SEAICE_hice_reg
1012    
1013     tmpscal0 = recip_hiceActual(I,J) / 2.0 _d 0
1014    
1015     C -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
1016     C /* IceGrowthRateOpenWater */
1017    
1018     S_a_IGROW(I,J) = ZERO
1019    
1020     c Expand ice cover if the open water growth rate is positive
1021     IF ( IceGrowthRateOpenWater(I,J) .GT. ZERO) THEN
1022     IF ((AREApreTH(I,J) .GT. ZERO) .OR.
1023     & (S_h(I,J) .GT. ZERO)) THEN
1024     c Determine which hemisphere for hemisphere-dependent
1025     c "lead closing variable", HO
1026     IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
1027     S_a_IGROW(I,J) = (ONE - AREApreTH(I,J)) *
1028     & IceGrowthRateOpenWater(I,J)/HO_south
1029     ELSE
1030     S_a_IGROW(I,J) = (ONE - AREApreTH(I,J)) *
1031     & IceGrowthRateOpenWater(I,J)/HO
1032     ENDIF
1033     ENDIF
1034    
1035     C -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
1036     c Contract ice cover if the open water growth rate is negative
1037     ELSE
1038     S_a_IGROW(I,J) = tmpscal0 *
1039     & IceGrowthRateOpenWater(I,J) * (ONE - AREApreTH(I,J))
1040     ENDIF
1041    
1042     S_a(I,J) = S_a(I,J) + S_a_IGROW(I,J)
1043    
1044    
1045     C -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
1046     c /* IceGrowthRateMixedLayer */
1047     C Contract ice if the IceGrowthRateMixedLayer is negative
1048     S_a_IGRML(I,J) = ZERO
1049    
1050     IF ( IceGrowthRateMixedLayer(I,J) .LE. ZERO) THEN
1051     S_a_IGRML(I,J) = tmpscal0 * IceGrowthRateMixedLayer(I,J)
1052     ENDIF
1053    
1054     S_a(I,J) = S_a(I,J) + S_a_IGRML(I,J)
1055    
1056     c
1057     C -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
1058     C Contract ice if the NetExistingIceGrowthRate is negative
1059     C /* NetExistingIceGrowthRate */
1060     S_a_IGRNE(I,J) = ZERO
1061     IF ( (NetExistingIceGrowthRate(I,J) .LE. ZERO) .AND.
1062     & (HEFFpreTH(I,J) .GT. ZERO) ) THEN
1063    
1064     S_a_IGRNE(I,J) =
1065     & tmpscal0 * NetExistingIceGrowthRate(I,J) * AREApreTH(I,J)
1066    
1067     ENDIF
1068    
1069     S_a(I,J) = S_a(I,J) + S_a_IGRNE(I,J)
1070    
1071     ENDDO /* I,J */
1072     ENDDO /* I,J */
1073    
1074    
1075     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1076     #ifdef ALLOW_AUTODIFF_TAMC
1077     CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
1078     CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
1079     CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
1080     CADJ STORE S_a = comlev1_bibj, key = iicekey, byte = isbyte
1081     CADJ STORE S_h = comlev1_bibj, key = iicekey, byte = isbyte
1082     CADJ STORE S_hsnow = comlev1_bibj, key = iicekey, byte = isbyte
1083     #endif /* ALLOW_AUTODIFF_TAMC */
1084     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1085    
1086     C Update the area, heff, and hsnow
1087     DO J=1,sNy
1088     DO I=1,sNx
1089     HEFF(I,J,bi,bj) = HEFFpreTH(I,J) +
1090     & SEAICE_deltaTtherm * S_h(I,J)
1091    
1092     AREA(I,J,bi,bj) = AREApreTH(I,J) +
1093     & SEAICE_deltaTtherm * S_a(I,J)
1094    
1095     HSNOW(I,J,bi,bj) = HSNWpreTH(I,J) +
1096     & SEAICE_deltaTtherm * S_hsnow(I,J)
1097     ENDDO
1098     ENDDO
1099    
1100     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1101     #ifdef ALLOW_AUTODIFF_TAMC
1102     CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
1103     CADJ STORE HEFFpreTH = comlev1_bibj,key=iicekey,byte=isbyte
1104     CADJ STORE HSNWpreTH = comlev1_bibj,key=iicekey,byte=isbyte
1105     CADJ STORE AREA (:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1106     CADJ STORE HEFF (:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1107     CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1108     #endif /* ALLOW_AUTODIFF_TAMC */
1109     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1110    
1111    
1112     DO J=1,sNy
1113     DO I=1,sNx
1114     c Bound area, heff, and hsnow
1115     AREA(I,J,bi,bj) = MIN(AREA(I,J,bi,bj), SEAICE_area_max)
1116     AREA(I,J,bi,bj) = MAX(ZERO, AREA(I,J,bi,bj))
1117     HEFF(I,J,bi,bj) = MAX(ZERO, HEFF(I,J,bi,bj))
1118     HSNOW(I,J,bi,bj) = MAX(ZERO, HSNOW(I,J,bi,bj))
1119    
1120     c Sanity checks
1121     IF (HEFF(I,J,bi,bj) .LE. ZERO .OR.
1122     & AREA(I,J,bi,bj) .LE. ZERO) THEN
1123    
1124     AREA(I,J,bi,bj) = 0. _d 0
1125     HEFF(I,J,bi,bj) = 0. _d 0
1126     hiceActual(I,J) = 0. _d 0
1127     HSNOW(I,J,bi,bj) = 0. _d 0
1128     hsnowActual(I,J) = 0. _d 0
1129    
1130     ELSE
1131     c recalcuate the actual ice and snow thicknesses
1132     hiceActual(I,J) = HEFF(I,J,bi,bj)/AREA(I,J,bi,bj)
1133     hsnowActual(I,J) = HSNOW(I,J,bi,bj)/AREA(I,J,bi,bj)
1134     ENDIF
1135    
1136     ENDDO
1137     ENDDO
1138    
1139    
1140     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1141     #ifdef ALLOW_AUTODIFF_TAMC
1142     CADJ STORE AREA (:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1143     CADJ STORE HEFF (:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1144     CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1145     CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
1146     CADJ & key = iicekey, byte = isbyte
1147     #endif /* ALLOW_AUTODIFF_TAMC */
1148     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1149    
1150     DO J=1,sNy
1151     DO I=1,sNx
1152    
1153     c THE EFFECTIVE SHORTWAVE HEATING RATE
1154     #ifdef SHORTWAVE_HEATING
1155     QSW(I,J,bi,bj) =
1156     & QSWI(I,J) * ( AREApreTH(I,J)) +
1157     & QSWO(I,J) * (ONE - AREApreTH(I,J))
1158     #else
1159     QSW(I,J,bi,bj) = 0. _d 0
1160     #endif
1161    
1162     c The actual ice volume change over the time step
1163     ActualNewTotalVolumeChange(I,J) =
1164     & HEFF(I,J,bi,bj) - HEFFpreTH(I,J)
1165    
1166     c The net average snow thickness melt that is actually realized. e.g.
1167     c hsnow_orig = 0.25 m (e.g. 1 m of ice over a cell 1/4 covered in snow)
1168     c hsnow_new = 0.20 m
1169     c snow accum = 0.05 m
1170     c melt = 0.25 + 0.05 - 0.2 = 0.1 m
1171    
1172     c snow has been accumulated in HSNWpreTH, so HSNOW can only
1173     c be .LE. HSNWpreTH at this point
1174     ActualNewTotalSnowMelt(I,J) =
1175     & HSNWpreTH(I,J) -
1176     & HSNOW(I,J,bi,bj)
1177     c & SnowAccOverIce(I,J) -
1178    
1179     c The energy required to melt or form the new ice volume
1180     EnergyInNewTotalIceVolume(I,J) =
1181     & ActualNewTotalVolumeChange(I,J)/QI
1182    
1183     c This is the net energy flux out of the ice+ocean system
1184     c Remember -----
1185     c F_ia_net : Under ice/snow surface freezing conditions,
1186     c vertical conductive heat flux convergence (F_c < 0) balances
1187     c heat flux divergence to atmosphere (F_ia > 0)
1188     c Otherwise, F_ia_net = F_ia (pos)
1189     c
1190     c F_io_net : Under ice/snow surface freezing conditions, F_c < 0.
1191     c Under ice surface melting conditions, F_c = 0 (no energy flux
1192     c from the ice to ocean)
1193     c
1194     c So if we are freezing, F_io_net = the conductive flux and there
1195     c is energy balance at ice surface, F_ia_net =0. If we are melting,
1196     c there is a convergence of energy into the ice from above
1197     NetEnergyFluxOutOfOcean(I,J) = SEAICE_deltaTtherm *
1198     & ( AREApreTH(I,J) *
1199     & (F_ia_net(I,J) + F_io_net(I,J) + QSWI(I,J))
1200     & + ( ONE - AREApreTH(I,J)) * F_ao(I,J))
1201    
1202     c THE QUANTITY OF HEAT WHICH IS THE RESIDUAL TO THE QUANTITY OF
1203     c ML temperature. If the net energy flux is exactly balanced by the
1204     c latent energy of fusion in the new ice created then we will not
1205     c change the ML temperature at all.
1206    
1207     ResidualEnergyOutOfOcean(I,J) =
1208     & NetEnergyFluxOutOfOcean(I,J) -
1209     & EnergyInNewTotalIceVolume(I,J)
1210    
1211     C NOW FORMULATE QNET
1212     C THIS QNET DETERMINES THE TEMPERATURE CHANGE
1213     C QNET IS A DEPTH AVERAGED HEAT FLUX FOR THE OCEAN COLUMN
1214    
1215     QNET(I,J,bi,bj) =
1216     & ResidualEnergyOutOfOcean(I,J) / SEAICE_deltaTtherm
1217    
1218    
1219     c Like snow melt, if there is melting, this quantity is positive.
1220     c The change of freshwater content is per unit area over the entire
1221     c cell, not just over the ice covered bits. This term is only used
1222     c to calculate freshwater fluxes for the purpose of changing the
1223     c salinity of the liquid cell. In the case of non-zero ice salinity,
1224     c the amount of freshwater is reduced by the ratio of ice salinity
1225     c to water cell salinity.
1226     IF (salt(I,J,kSurface,bi,bj) .GE. SEAICE_salt0 .AND.
1227     & salt(I,J,kSurface,bi,bj) .GT. 0. _d 0) THEN
1228    
1229     FreshwaterContribFromIce(I,J) =
1230     & - ActualNewTotalVolumeChange(I,J) *
1231     & SEAICE_rhoICE/rhoConstFresh *
1232     & (ONE - SEAICE_salt0/salt(I,J,kSurface,bi,bj))
1233    
1234     ELSE
1235     C If the liquid cell has a lower salinity than the specified
1236     c salinity of sea ice then assume the sea ice is completely fresh
1237     FreshwaterContribFromIce(I,J) =
1238     & -ActualNewTotalVolumeChange(I,J) *
1239     & SEAICE_rhoIce/rhoConstFresh
1240     ENDIF
1241    
1242    
1243     c The freshwater contribution from snow comes only in the form of melt
1244     c unlike ice, which takes freshwater upon growth and yields freshwater
1245     c upon melt. This is why the the actual new average snow melt was determined.
1246     c In m/m^2 over the entire cell.
1247     FreshwaterContribFromSnowMelt(I,J) =
1248     & ActualNewTotalSnowMelt(I,J)*SEAICE_rhoSnow/rhoConstFresh
1249    
1250     c This seems to be in m/s, original time level 2 for area
1251     c Only the precip and evap need to be area weighted. The runoff
1252     c and freshwater contribs from ice and snow melt are already mean
1253     c weighted
1254     EmPmR(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
1255     & ( EVAP(I,J,bi,bj) - PRECIP(I,J,bi,bj) )
1256     & * ( ONE - AREApreTH(I,J) )
1257     & - PrecipRateOverIceSurfaceToSea(I,J)*AREApreTH(I,J)
1258     #ifdef ALLOW_RUNOFF
1259     & - RUNOFF(I,J,bi,bj)
1260     #endif
1261     & - (FreshwaterContribFromIce(I,J) +
1262     & FreshwaterContribFromSnowMelt(I,J)) /
1263     & SEAICE_deltaTtherm ) * rhoConstFresh
1264    
1265     ENDDO
1266     ENDDO
1267    
1268    
1269     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1270     #ifdef ALLOW_AUTODIFF_TAMC
1271     CADJ STORE hiceActual = comlev1_bibj,key=iicekey,byte=isbyte
1272     CADJ STORE hsnowActual = comlev1_bibj,key=iicekey,byte=isbyte
1273     CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
1274     CADJ STORE HEFF(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1275     CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1276     CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1277     CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
1278     CADJ & key = iicekey, byte = isbyte
1279     #endif /* ALLOW_AUTODIFF_TAMC */
1280     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1281    
1282     #ifdef ALLOW_SEAICE_FLOODING
1283     IF(SEAICEuseFlooding) THEN
1284    
1285     DO J = 1,sNy
1286     DO I = 1,sNx
1287     tmpscal0 = FL_C2*( hsnowActual(I,J) -
1288     & hiceActual(I,J)*FL_C3 )
1289    
1290     IF (tmpscal0 .GT. ZERO) THEN
1291     tmpscal1 = FL_C4*tmpscal0
1292    
1293     hiceActual(I,J) = hiceActual(I,J)
1294     & + tmpscal1
1295    
1296     hsnowActual(I,J) = hsnowActual(I,J)
1297     & - tmpscal0
1298    
1299     HEFF(I,J,bi,bj)= hiceActual(I,J) *
1300     & AREA(I,J,bi,bj)
1301    
1302     HSNOW(I,J,bi,bj) = hsnowActual(I,J)*
1303     & AREA(I,J,bi,bj)
1304    
1305     ENDIF /* flooding */
1306     ENDDO
1307     ENDDO
1308     c SEAICEuseFlooding
1309     ENDIF
1310     c ALLOW_SEAICE_FLOODING
1311     #endif
1312    
1313     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1314     #ifdef ALLOW_AUTODIFF_TAMC
1315     CADJ STORE hiceActual = comlev1_bibj,key=iicekey,byte=isbyte
1316     CADJ STORE hsnowActual = comlev1_bibj,key=iicekey,byte=isbyte
1317     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1318     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1319     #endif /* ALLOW_AUTODIFF_TAMC */
1320     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
1321    
1322     C Sea Ice Load on the sea surface.
1323     C =================================
1324     IF ( useRealFreshWaterFlux ) THEN
1325     DO J=1,sNy
1326     DO I=1,sNx
1327     #ifdef SEAICE_CAP_ICELOAD
1328     tmpscal1 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
1329     & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1330     tmpscal2 = MIN(tmpscal1,heffTooHeavy*rhoConst)
1331     #else
1332     tmpscal2 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
1333     & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1334     #endif
1335     sIceLoad(i,j,bi,bj) = tmpscal2
1336     ENDDO
1337     ENDDO
1338     ENDIF
1339    
1340    
1341     #ifdef SEAICE_DEBUG
1342     DO j=1,sNy
1343     DO i=1,sNx
1344    
1345     IF ( (i .EQ. SEAICE_debugPointI) .and.
1346     & (j .EQ. SEAICE_debugPointJ) ) THEN
1347    
1348     print *,'ifice: myTime,myIter:',myTime,myIter
1349    
1350     print '(A,2i4,2(1x,1P3E15.4))',
1351     & 'ifice i j -------------- ',i,j
1352    
1353     print '(A,2i4,2(1x,1P3E15.4))',
1354     & 'ifice i j F(mi ao), rHIA ',
1355     & i,j, F_mi(i,j), F_ao(i,j),
1356     & recip_hiceActual(i,j)
1357    
1358     print '(A,2i4,2(1x,1P3E15.4))',
1359     & 'ifice i j Fi(a,ant2/1 ont)',
1360     & i,j, F_ia(i,j),
1361     & F_ia_net_before_snow(i,j),
1362     & F_ia_net(i,j),
1363     & F_io_net(i,j)
1364    
1365     print '(A,2i4,2(1x,1P3E15.4))',
1366     & 'ifice i j AREA2/1 HEFF2/1 ',i,j,
1367     & AREApreTH(I,J),
1368     & AREA(i,j,bi,bj),
1369     & HEFFpreTH(I,J),
1370     & HEFF(i,j,bi,bj)
1371    
1372     print '(A,2i4,2(1x,1P3E15.4))',
1373     & 'ifice i j HSNOW2/1 TMX ',i,j,
1374     & HSNWpreTH(I,J),
1375     & HSNOW(I,J,bi,bj),
1376     & theta(I,J,kSurface,bi,bj)
1377    
1378     print '(A,2i4,2(1x,1P3E15.4))',
1379     & 'ifice i j TI ATP LWD ',i,j,
1380     & TICES(i,j,1, bi,bj) - celsius2k,
1381     & ATEMP(i,j,bi,bj) - celsius2k,
1382     & LWDOWN(i,j,bi,bj)
1383    
1384     print '(A,2i4,2(1x,1P3E15.4))',
1385     & 'ifice i j S_a(tot,OW,ML,NE',i,j,
1386     & S_a(i,j),
1387     & S_a_IGROW(I,J),
1388     & S_a_IGRML(I,J),
1389     & S_a_IGRNE(I,J)
1390    
1391     print '(A,2i4,2(1x,1P3E15.4))',
1392     & 'ifice i j S_a S_h S_hsnow ',i,j,
1393     & S_a(i,j),
1394     & S_h(i,j),
1395     & S_hsnow(i,j)
1396    
1397     print '(A,2i4,2(1x,1P3E15.4))',
1398     & 'ifice i j IGR(ML OW ICE) ',i,j,
1399     & IceGrowthRateMixedLayer(i,j),
1400     & IceGrowthRateOpenWater(i,j),
1401     & NetExistingIceGrowthRate(i,j)
1402    
1403    
1404     print '(A,2i4,2(1x,1P3E15.4))',
1405     & 'ifice i j THETA/TFRZ/SALT ',i,j,
1406     & theta(I,J,kSurface,bi,bj),
1407     & tmpFrz,
1408     & salt(I,J,kSurface,bi,bj)
1409    
1410     print '(A,2i4,2(1x,1P3E15.4))',
1411     & 'ifice i j IVC(A ENIN) ',i,j,
1412     & ActualNewTotalVolumeChange(i,j),
1413     & EnergyInNewTotalIceVolume(i,j)
1414     c & ExpectedIceVolumeChange(i,j),
1415    
1416     print '(A,2i4,2(1x,1P3E15.4))',
1417     & 'ifice i j EF(NOS RE) QNET ',i,j,
1418     & NetEnergyFluxOutOfOcean(i,j),
1419     & ResidualEnergyOutOfOcean(i,j),
1420     & QNET(I,J,bi,bj)
1421    
1422     print '(A,2i4,3(1x,1P3E15.4))',
1423     & 'ifice i j QSW QSWO QSWI ',i,j,
1424     & QSW(i,j,bi,bj),
1425     & QSWO(i,j),
1426     & QSWI(i,j)
1427    
1428     c print '(A,2i4,3(1x,1P3E15.4))',
1429     c & 'ifice i j SW(BML IML SW) ',i,j,
1430     c & QSW_absorb_below_first_layer(i,j),
1431     c & QSW_absorb_in_first_layer(i,j),
1432     c & SWFRACB
1433    
1434     c print '(A,2i4,3(1x,1P3E15.4))',
1435     c & 'ifice i j ptc(to, qsw, oa)',i,j,
1436     c & PredTempChange(i,j),
1437     c & PredTempChangeFromQSW (i,j),
1438     c & PredTempChangeFromOA_MQNET(i,j)
1439    
1440    
1441     c print '(A,2i4,3(1x,1P3E15.4))',
1442     c & 'ifice i j ptc(fion,ian,ia)',i,j,
1443     c & PredTempChangeFromF_IO_NET(i,j),
1444     c & PredTempChangeFromF_IA_NET(i,j),
1445     c & PredTempChangeFromFIA(i,j)
1446    
1447     c print '(A,2i4,3(1x,1P3E15.4))',
1448     c & 'ifice i j ptc(niv) ',i,j,
1449     c & PredTempChangeFromNewIceVol(i,j)
1450    
1451    
1452     print '(A,2i4,3(1x,1P3E15.4))',
1453     & 'ifice i j EmPmR EVP PRE RU',i,j,
1454     & EmPmR(I,J,bi,bj),
1455     & EVAP(I,J,bi,bj),
1456     & PRECIP(I,J,bi,bj),
1457     & RUNOFF(I,J,bi,bj)
1458    
1459     print '(A,2i4,3(1x,1P3E15.4))',
1460     & 'ifice i j PRROIS,SAOI(R .)',i,j,
1461     & PrecipRateOverIceSurfaceToSea(I,J),
1462     & SnowAccRateOverIce(I,J),
1463     & SnowAccOverIce(I,J)
1464    
1465     print '(A,2i4,4(1x,1P3E15.4))',
1466     & 'ifice i j SM(PM PMR . .R) ',i,j,
1467     & PotSnowMeltFromSurf(I,J),
1468     & PotSnowMeltRateFromSurf(I,J),
1469     & SnowMeltFromSurface(I,J),
1470     & SnowMeltRateFromSurface(I,J)
1471    
1472     print '(A,2i4,4(1x,1P3E15.4))',
1473     & 'ifice i j TotSnwMlt ',i,j,
1474     & ActualNewTotalSnowMelt(I,J)
1475     c & ExpectedSnowVolumeChange(I,J)
1476    
1477     print '(A,2i4,4(1x,1P3E15.4))',
1478     & 'ifice i j fw(CFICE, CFSM) ',i,j,
1479     & FreshwaterContribFromIce(I,J),
1480     & FreshwaterContribFromSnowMelt(I,J)
1481    
1482     print '(A,2i4,2(1x,1P3E15.4))',
1483     & 'ifice i j -------------- ',i,j
1484    
1485     ENDIF
1486    
1487     ENDDO
1488     ENDDO
1489     #endif /* SEAICE_DEBUG */
1490    
1491    
1492     C close bi,bj loops
1493     ENDDO
1494     ENDDO
1495    
1496     #else /* ALLOW_EXF and ALLOW_ATM_TEMP */
1497     STOP 'SEAICE_GROWTH not compiled without EXF and ALLOW_ATM_TEMP'
1498     #endif /* ALLOW_EXF and ALLOW_ATM_TEMP */
1499    
1500     RETURN
1501     END
  ViewVC Help
Powered by ViewVC 1.1.22