/[MITgcm]/MITgcm/pkg/seaice/seaice_growth.F
ViewVC logotype

Annotation of /MITgcm/pkg/seaice/seaice_growth.F

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


Revision 1.98 - (hide annotations) (download)
Thu Nov 11 16:37:39 2010 UTC (13 years, 7 months ago) by mlosch
Branch: MAIN
Changes since 1.97: +38 -45 lines
a few cosmetic changes that reduce the number of lines of code

1 mlosch 1.98 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_growth.F,v 1.97 2010/11/05 08:10:57 mlosch Exp $
2 jmc 1.91 C $Name: $
3    
4 dimitri 1.69 #include "SEAICE_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: SEAICE_GROWTH
8     C !INTERFACE:
9     SUBROUTINE SEAICE_GROWTH( myTime, myIter, myThid )
10     C !DESCRIPTION: \bv
11     C *==========================================================*
12     C | SUBROUTINE seaice_growth
13     C | o Updata ice thickness and snow depth
14     C *==========================================================*
15     C \ev
16    
17     C !USES:
18     IMPLICIT NONE
19     C === Global variables ===
20     #include "SIZE.h"
21     #include "EEPARAMS.h"
22     #include "PARAMS.h"
23     #include "DYNVARS.h"
24     #include "GRID.h"
25     #include "FFIELDS.h"
26     #include "SEAICE_PARAMS.h"
27     #include "SEAICE.h"
28     #ifdef ALLOW_EXF
29     # include "EXF_OPTIONS.h"
30     # include "EXF_FIELDS.h"
31     # include "EXF_PARAM.h"
32     #endif
33     #ifdef ALLOW_SALT_PLUME
34     # include "SALT_PLUME.h"
35     #endif
36     #ifdef ALLOW_AUTODIFF_TAMC
37     # include "tamc.h"
38     #endif
39    
40     C !INPUT/OUTPUT PARAMETERS:
41     C === Routine arguments ===
42     C myTime :: Simulation time
43     C myIter :: Simulation timestep number
44     C myThid :: Thread no. that called this routine.
45     _RL myTime
46     INTEGER myIter, myThid
47    
48     C !LOCAL VARIABLES:
49     C === Local variables ===
50 gforget 1.89 C
51 jmc 1.91 C unit/sign convention:
52     C Within the thermodynamic computation all stocks, except HSNOW,
53 gforget 1.89 C are in 'effective ice meters' units, and >0 implies more ice.
54 jmc 1.91 C This holds for stocks due to ocean and atmosphere heat,
55     C at the outset of 'PART 2: determine heat fluxes/stocks'
56 gforget 1.89 C and until 'PART 7: determine ocean model forcing'
57     C This strategy minimizes the need for multiplications/divisions
58 jmc 1.91 C by ice fraction, heat capacity, etc. The only conversions that
59     C occurs are for the HSNOW (in effective snow meters) and
60     C PRECIP (fresh water m/s).
61 gforget 1.89 c
62     c HEFF is effective Hice thickness (m3/m2)
63     c HSNOW is Heffective snow thickness (m3/m2)
64     c HSALT is Heffective salt content (g/m2)
65     c AREA is the seaice cover fraction (0<=AREA<=1)
66     c Q denotes heat stocks -- converted to ice stocks (m3/m2) early on
67     c
68 jmc 1.91 c For all other stocks/increments, such as d_HEFFbyATMonOCN
69 gforget 1.89 c or a_QbyATM_cover, the naming convention is as follows:
70 jmc 1.91 c The prefix 'a_' means available, the prefix 'd_' means delta
71     c (i.e. increment), and the prefix 'r_' means residual.
72 gforget 1.89 c The suffix '_cover' denotes a value for the ice covered fraction
73 mlosch 1.97 c of the grid cell, whereas '_open' is for the open water fraction.
74 jmc 1.91 c The main part of the name states what ice/snow stock is concerned
75 gforget 1.89 c (e.g. QbyATM or HEFF), and how it is affected (e.g. d_HEFFbyATMonOCN
76 jmc 1.91 c is the increment of HEFF due to the ATMosphere extracting heat from the
77 gforget 1.89 c OCeaN surface, or providing heat to the OCeaN surface).
78    
79     CEOP
80 dimitri 1.69 C i,j,bi,bj :: Loop counters
81     INTEGER i, j, bi, bj
82     C number of surface interface layer
83     INTEGER kSurface
84     C constants
85 gforget 1.80 _RL TBC, ICE2SNOW
86 gforget 1.85 _RL QI, QS, Lfusion
87 dimitri 1.69
88 jmc 1.91 C a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
89 gforget 1.75 C the atmosphere and the ocean surface - for ice covered water
90 gforget 1.89 C a_QbyATM_open :: same but for open water
91 gforget 1.84 C r_QbyATM_cover :: residual of a_QbyATM_cover after freezing/melting processes
92 gforget 1.89 C r_QbyATM_open :: same but for open water
93 gforget 1.76 _RL a_QbyATM_cover (1:sNx,1:sNy)
94     _RL a_QbyATM_open (1:sNx,1:sNy)
95     _RL r_QbyATM_cover (1:sNx,1:sNy)
96 gforget 1.89 _RL r_QbyATM_open (1:sNx,1:sNy)
97 gforget 1.75 C a_QSWbyATM_open - short wave heat flux over ocean in W/m^2
98     C a_QSWbyATM_cover - short wave heat flux under ice in W/m^2
99 gforget 1.76 _RL a_QSWbyATM_open (1:sNx,1:sNy)
100     _RL a_QSWbyATM_cover (1:sNx,1:sNy)
101 jmc 1.91 C a_QbyOCN :: available heat (in in W/m^2) due to the
102 gforget 1.76 C interaction of the ice pack and the ocean surface
103 jmc 1.91 C r_QbyOCN :: residual of a_QbyOCN after freezing/melting
104 gforget 1.75 C processes have been accounted for
105 gforget 1.84 _RL a_QbyOCN (1:sNx,1:sNy)
106     _RL r_QbyOCN (1:sNx,1:sNy)
107 gforget 1.76
108     c conversion factors to go from Q (W/m2) to HEFF (ice meters)
109     _RL convertQ2HI, convertHI2Q
110     c conversion factors to go from precip (m/s) unit to HEFF (ice meters)
111     _RL convertPRECIP2HI, convertHI2PRECIP
112 gforget 1.75
113     c ICE/SNOW stocks tendencies associated with the various melt/freeze processes
114 gforget 1.76 _RL d_AREAbyATM (1:sNx,1:sNy)
115 gforget 1.71 c
116 gforget 1.87 _RL d_HEFFbyNEG (1:sNx,1:sNy)
117 gforget 1.84 _RL d_HEFFbyOCNonICE (1:sNx,1:sNy)
118 gforget 1.75 _RL d_HEFFbyATMonOCN (1:sNx,1:sNy)
119 gforget 1.84 _RL d_HEFFbyFLOODING (1:sNx,1:sNy)
120 gforget 1.73 c
121 gforget 1.90 _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)
122     c
123 gforget 1.87 _RL d_HSNWbyNEG (1:sNx,1:sNy)
124 gforget 1.84 _RL d_HSNWbyATMonSNW (1:sNx,1:sNy)
125 gforget 1.75 _RL d_HSNWbyOCNonSNW (1:sNx,1:sNy)
126 gforget 1.84 _RL d_HSNWbyRAIN (1:sNx,1:sNy)
127     c
128     _RL d_HFRWbyRAIN (1:sNx,1:sNy)
129 gforget 1.71
130 dimitri 1.69 C actual ice thickness with upper and lower limit
131 gforget 1.83 _RL heffActual (1:sNx,1:sNy)
132 dimitri 1.69 C actual snow thickness
133 gforget 1.83 _RL hsnowActual (1:sNx,1:sNy)
134 gforget 1.87
135     C AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
136     _RL AREApreTH (1:sNx,1:sNy)
137     _RL HEFFpreTH (1:sNx,1:sNy)
138     _RL HSNWpreTH (1:sNx,1:sNy)
139    
140 dimitri 1.69 C wind speed
141 gforget 1.76 _RL UG (1:sNx,1:sNy)
142 dimitri 1.69 _RL SPEED_SQ
143 gforget 1.96
144     C pathological cases thresholds
145 gforget 1.92 _RL heffTooThin, heffTooHeavy
146 dimitri 1.69
147 gforget 1.75 c temporary variables available for the various computations
148 gforget 1.88 _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
149 gforget 1.76 _RL tmparr1 (1:sNx,1:sNy)
150 gforget 1.75
151    
152     #ifdef ALLOW_SEAICE_FLOODING
153     _RL hDraft
154     #endif /* ALLOW_SEAICE_FLOODING */
155 jmc 1.91
156 gforget 1.75 #ifdef SEAICE_SALINITY
157 gforget 1.83 _RL saltFluxAdjust (1:sNx,1:sNy)
158 gforget 1.75 #endif
159    
160 mlosch 1.98 INTEGER nDim
161 dimitri 1.69 #ifdef SEAICE_MULTICATEGORY
162 mlosch 1.98 INTEGER ilockey
163     PARAMETER ( nDim = MULTDIM )
164     #else
165     PARAMETER ( nDim = 1 )
166     #endif /* SEAICE_MULTICATEGORY */
167 dimitri 1.69 INTEGER it
168 mlosch 1.98 _RL pFac
169 gforget 1.83 _RL heffActualP (1:sNx,1:sNy)
170     _RL a_QbyATMmult_cover (1:sNx,1:sNy)
171     _RL a_QSWbyATMmult_cover(1:sNx,1:sNy)
172 dimitri 1.69
173     #ifdef ALLOW_DIAGNOSTICS
174     _RL DIAGarray (1:sNx,1:sNy)
175     LOGICAL DIAGNOSTICS_IS_ON
176     EXTERNAL DIAGNOSTICS_IS_ON
177     #endif
178    
179 gforget 1.88 c ===================================================================
180     c =================PART 0: constants and initializations=============
181     c ===================================================================
182    
183 dimitri 1.69 IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
184     kSurface = Nr
185     ELSE
186     kSurface = 1
187     ENDIF
188    
189 gforget 1.92 c Cutoff for very thin ice
190     heffTooThin=1. _d -5
191     c Cutoff for iceload
192     heffTooHeavy=dRf(kSurface) / 5. _d 0
193    
194 dimitri 1.69 C FREEZING TEMP. OF SEA WATER (deg C)
195     TBC = SEAICE_freeze
196 gforget 1.80
197 dimitri 1.69 C RATIO OF SEA ICE DENSITY to SNOW DENSITY
198 gforget 1.80 ICE2SNOW = SEAICE_rhoIce/SEAICE_rhoSnow
199 gforget 1.85
200 gforget 1.89 C FRESH ICE/SNOW LATENT HEAT OF FUSION (J/kg)
201 gforget 1.85 Lfusion = 3.34 _d +05
202 dimitri 1.69 C HEAT OF FUSION OF ICE (J/m^3)
203 gforget 1.85 QI = SEAICE_rhoIce*Lfusion
204 dimitri 1.69 C HEAT OF FUSION OF SNOW (J/m^3)
205 gforget 1.85 QS = SEAICE_rhoSnow*Lfusion
206     C
207     C note that QI/QS=ICE2SNOW -- except it wasnt in old code
208    
209    
210 dimitri 1.69
211 gforget 1.76 c conversion factors to go from Q (W/m2) to HEFF (ice meters)
212 gforget 1.84 convertQ2HI=SEAICE_deltaTtherm/QI
213 gforget 1.76 convertHI2Q=1/convertQ2HI
214     c conversion factors to go from precip (m/s) unit to HEFF (ice meters)
215     convertPRECIP2HI=SEAICE_deltaTtherm*rhoConstFresh/SEAICE_rhoIce
216     convertHI2PRECIP=1./convertPRECIP2HI
217    
218 dimitri 1.69 DO bj=myByLo(myThid),myByHi(myThid)
219     DO bi=myBxLo(myThid),myBxHi(myThid)
220     c
221     #ifdef ALLOW_AUTODIFF_TAMC
222     act1 = bi - myBxLo(myThid)
223     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
224     act2 = bj - myByLo(myThid)
225     max2 = myByHi(myThid) - myByLo(myThid) + 1
226     act3 = myThid - 1
227     max3 = nTx*nTy
228     act4 = ikey_dynamics - 1
229     iicekey = (act1 + 1) + act2*max1
230     & + act3*max1*max2
231     & + act4*max1*max2*max3
232     #endif /* ALLOW_AUTODIFF_TAMC */
233 gforget 1.75
234 dimitri 1.69 #ifdef ALLOW_AUTODIFF_TAMC
235     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
236     CADJ & key = iicekey, byte = isbyte
237     CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj,
238     CADJ & key = iicekey, byte = isbyte
239     CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj,
240     CADJ & key = iicekey, byte = isbyte
241     #endif /* ALLOW_AUTODIFF_TAMC */
242 gforget 1.75
243    
244     C array initializations
245     C =====================
246    
247 dimitri 1.69 DO J=1,sNy
248     DO I=1,sNx
249 gforget 1.90 a_QbyATM_cover (I,J) = 0.0 _d 0
250     a_QbyATM_open(I,J) = 0.0 _d 0
251     r_QbyATM_cover (I,J) = 0.0 _d 0
252     r_QbyATM_open (I,J) = 0.0 _d 0
253 gforget 1.75 c
254 gforget 1.71 a_QSWbyATM_open (I,J) = 0.0 _d 0
255 gforget 1.90 a_QSWbyATM_cover (I,J) = 0.0 _d 0
256 gforget 1.75 c
257 gforget 1.90 a_QbyOCN (I,J) = 0.0 _d 0
258     r_QbyOCN (I,J) = 0.0 _d 0
259 gforget 1.75 c
260 gforget 1.90 d_AREAbyATM(I,J) = 0.0 _d 0
261 gforget 1.75 c
262 gforget 1.84 d_HEFFbyOCNonICE(I,J) = 0.0 _d 0
263 gforget 1.75 d_HEFFbyATMonOCN(I,J) = 0.0 _d 0
264 gforget 1.84 d_HEFFbyFLOODING(I,J) = 0.0 _d 0
265 gforget 1.75 c
266 gforget 1.90 d_HEFFbyATMonOCN_open(I,J) = 0.0 _d 0
267     c
268 gforget 1.84 d_HSNWbyATMonSNW(I,J) = 0.0 _d 0
269 gforget 1.75 d_HSNWbyOCNonSNW(I,J) = 0.0 _d 0
270 gforget 1.90 d_HSNWbyRAIN(I,J) = 0.0 _d 0
271 gforget 1.84 c
272 gforget 1.90 d_HFRWbyRAIN(I,J) = 0.0 _d 0
273 gforget 1.75 c
274 gforget 1.90 tmparr1(I,J) = 0.0 _d 0
275 gforget 1.75 c
276 dimitri 1.69 #ifdef SEAICE_SALINITY
277 gforget 1.90 saltFluxAdjust(I,J) = 0.0 _d 0
278 dimitri 1.69 #endif
279 gforget 1.90 a_QbyATMmult_cover(I,J) = 0.0 _d 0
280     a_QSWbyATMmult_cover(I,J) = 0.0 _d 0
281 dimitri 1.69 ENDDO
282     ENDDO
283 gforget 1.84 #ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
284 dimitri 1.69 DO J=1-oLy,sNy+oLy
285     DO I=1-oLx,sNx+oLx
286 gforget 1.90 frWtrAtm(I,J,bi,bj) = 0.0 _d 0
287 dimitri 1.69 ENDDO
288     ENDDO
289 gforget 1.84 #endif
290 dimitri 1.69
291 gforget 1.87
292 gforget 1.88
293     c =====================================================================
294     c ===========PART 1: treat pathological cases (post advdiff)===========
295     c =====================================================================
296    
297    
298 gforget 1.87
299     #ifndef SEAICE_GROWTH_LEGACY
300    
301     c 1) treat the case of negative values:
302     c
303     #ifdef ALLOW_AUTODIFF_TAMC
304     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
305     CADJ & key = iicekey, byte = isbyte
306     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
307     CADJ & key = iicekey, byte = isbyte
308     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
309     CADJ & key = iicekey, byte = isbyte
310     #endif /* ALLOW_AUTODIFF_TAMC */
311     DO J=1,sNy
312     DO I=1,sNx
313 gforget 1.89 d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
314 gforget 1.87 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
315 gforget 1.89 d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
316 gforget 1.87 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
317     AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
318     ENDDO
319     ENDDO
320    
321 gforget 1.92 c 1.25) treat the case of very thin ice:
322     c
323     #ifdef ALLOW_AUTODIFF_TAMC
324     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
325     CADJ & key = iicekey, byte = isbyte
326     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
327     CADJ & key = iicekey, byte = isbyte
328     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
329     CADJ & key = iicekey, byte = isbyte
330     #endif /* ALLOW_AUTODIFF_TAMC */
331     DO J=1,sNy
332     DO I=1,sNx
333     if (HEFF(I,J,bi,bj).LE.heffTooThin) then
334     tmpscal2=-HEFF(I,J,bi,bj)
335     tmpscal3=-HSNOW(I,J,bi,bj)
336 gforget 1.96 TICE(I,J,bi,bj)=celsius2K
337 gforget 1.92 else
338     tmpscal2=0. _d 0
339     tmpscal3=0. _d 0
340     endif
341     HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal2
342     d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
343     HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
344     d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
345     ENDDO
346     ENDDO
347    
348 gforget 1.90 c 1.5) treat the case of area but no ice/snow:
349     c
350     #ifdef ALLOW_AUTODIFF_TAMC
351     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
352     CADJ & key = iicekey, byte = isbyte
353     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
354     CADJ & key = iicekey, byte = isbyte
355     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
356     CADJ & key = iicekey, byte = isbyte
357     #endif /* ALLOW_AUTODIFF_TAMC */
358     DO J=1,sNy
359     DO I=1,sNx
360     IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
361     & (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
362     ENDDO
363     ENDDO
364    
365 gforget 1.87 c 2) treat the case of very small area:
366     c
367     #ifdef ALLOW_AUTODIFF_TAMC
368     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
369     CADJ & key = iicekey, byte = isbyte
370     #endif /* ALLOW_AUTODIFF_TAMC */
371     DO J=1,sNy
372     DO I=1,sNx
373 jmc 1.91 IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0))
374 gforget 1.87 & AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),areaMin)
375     ENDDO
376     ENDDO
377    
378 gforget 1.89 c 2.5) treat case of excessive ice cover:
379     c
380     #ifdef ALLOW_AUTODIFF_TAMC
381     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
382     CADJ & key = iicekey, byte = isbyte
383     #endif /* ALLOW_AUTODIFF_TAMC */
384     DO J=1,sNy
385     DO I=1,sNx
386     AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),areaMax)
387     ENDDO
388     ENDDO
389    
390 gforget 1.87 c 3) store regularized values of heff, hsnow, area at the onset of thermo.
391     DO J=1,sNy
392     DO I=1,sNx
393     HEFFpreTH(I,J)=HEFF(I,J,bi,bj)
394     HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
395     AREApreTH(I,J)=AREA(I,J,bi,bj)
396     ENDDO
397     ENDDO
398    
399     c 4) treat sea ice salinity pathological cases
400     #ifdef SEAICE_SALINITY
401     #ifdef ALLOW_AUTODIFF_TAMC
402     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
403     CADJ & key = iicekey, byte = isbyte
404     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
405     CADJ & key = iicekey, byte = isbyte
406     #endif /* ALLOW_AUTODIFF_TAMC */
407     DO J=1,sNy
408     DO I=1,sNx
409     IF ( (HSALT(I,J,bi,bj) .LT. 0.0).OR.
410     & (HEFF(I,J,bi,bj) .EQ. 0.0) ) THEN
411     saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
412     & HSALT(I,J,bi,bj) / SEAICE_deltaTtherm
413     HSALT(I,J,bi,bj) = 0.0 _d 0
414     ENDIF
415     ENDDO
416     ENDDO
417     #endif /* SEAICE_SALINITY */
418    
419     c 5) treat sea ice age pathological cases
420     c ...
421 jmc 1.91
422 gforget 1.87 #else
423     #ifdef ALLOW_AUTODIFF_TAMC
424     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
425     CADJ & key = iicekey, byte = isbyte
426     #endif /* ALLOW_AUTODIFF_TAMC */
427 dimitri 1.69 DO J=1,sNy
428     DO I=1,sNx
429 gforget 1.87 HEFFpreTH(I,J)=HEFFNM1(I,J,bi,bj)
430     HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
431     AREApreTH(I,J)=AREANM1(I,J,bi,bj)
432     d_HEFFbyNEG(I,J)=0. _d 0
433     d_HSNWbyNEG(I,J)=0. _d 0
434 dimitri 1.69 ENDDO
435     ENDDO
436 gforget 1.87 #endif /* SEAICE_GROWTH_LEGACY */
437 dimitri 1.69
438    
439 gforget 1.89 c 4) COMPUTE ACTUAL ICE/SNOW THICKNESS; USE MIN/MAX VALUES
440     c TO REGULARIZE SEAICE_SOLVE4TEMP/d_AREA COMPUTATIONS
441 gforget 1.87 c
442 dimitri 1.69 #ifdef ALLOW_AUTODIFF_TAMC
443 gforget 1.87 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
444     CADJ & key = iicekey, byte = isbyte
445 dimitri 1.69 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
446     CADJ & key = iicekey, byte = isbyte
447     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
448     CADJ & key = iicekey, byte = isbyte
449     #endif /* ALLOW_AUTODIFF_TAMC */
450     DO J=1,sNy
451     DO I=1,sNx
452 gforget 1.87 tmpscal1 = MAX(areaMin,AREApreTH(I,J))
453     hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1
454     tmpscal2 = HEFFpreTH(I,J)/tmpscal1
455     heffActual(I,J) = MAX(tmpscal2,hiceMin)
456 gforget 1.89 cgf do we need to keep this comment?
457     C
458 dimitri 1.69 C Capping the actual ice thickness effectively enforces a
459     C minimum of heat flux through the ice and helps getting rid of
460     C very thick ice.
461     cdm actually, this does exactly the opposite, i.e., ice is thicker
462 gforget 1.83 cdm when heffActual is capped, so I am commenting out
463     cdm heffActual(I,J) = MIN(heffActual(I,J),9.0 _d +00)
464 dimitri 1.69 ENDDO
465     ENDDO
466    
467 gforget 1.95 #ifdef ALLOW_AUTODIFF_TAMC
468 gforget 1.96 #ifdef SEAICE_SIMPLIFY_GROWTH_ADJ
469 gforget 1.95 CALL ZERO_ADJ_1D( sNx*sNy, heffActual, myThid)
470 gforget 1.96 CALL ZERO_ADJ_1D( sNx*sNy, hsnowActual, myThid)
471 gforget 1.95 #endif
472     #endif
473 gforget 1.87
474 gforget 1.88 c ===================================================================
475 gforget 1.89 c ================PART 2: determine heat fluxes/stocks===============
476 gforget 1.88 c ===================================================================
477    
478    
479    
480 gforget 1.75 C determine available heat due to the atmosphere -- for open water
481     C ================================================================
482    
483     C ocean surface/mixed layer temperature
484 dimitri 1.69 DO J=1,sNy
485     DO I=1,sNx
486 gforget 1.83 TMIX(I,J,bi,bj)=theta(I,J,kSurface,bi,bj)+celsius2K
487 dimitri 1.69 ENDDO
488     ENDDO
489    
490 gforget 1.75 C wind speed from exf
491 dimitri 1.69 DO J=1,sNy
492     DO I=1,sNx
493     UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))
494     ENDDO
495     ENDDO
496    
497 gforget 1.75 CALL SEAICE_BUDGET_OCEAN(
498     I UG,
499     U TMIX,
500     O a_QbyATM_open, a_QSWbyATM_open,
501     I bi, bj, myTime, myIter, myThid )
502    
503    
504     C determine available heat due to the atmosphere -- for ice covered water
505     C =======================================================================
506 dimitri 1.69
507 gforget 1.89 #ifdef ALLOW_ATM_WIND
508 dimitri 1.69 IF (useRelativeWind) THEN
509     C Compute relative wind speed over sea ice.
510     DO J=1,sNy
511     DO I=1,sNx
512     SPEED_SQ =
513     & (uWind(I,J,bi,bj)
514     & +0.5 _d 0*(uVel(i,j,kSurface,bi,bj)
515     & +uVel(i+1,j,kSurface,bi,bj))
516     & -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj)))**2
517     & +(vWind(I,J,bi,bj)
518     & +0.5 _d 0*(vVel(i,j,kSurface,bi,bj)
519     & +vVel(i,j+1,kSurface,bi,bj))
520     & -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj)))**2
521     IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
522     UG(I,J)=SEAICE_EPS
523     ELSE
524     UG(I,J)=SQRT(SPEED_SQ)
525     ENDIF
526     ENDDO
527     ENDDO
528     ENDIF
529 gforget 1.89 #endif
530 gforget 1.75
531 mlosch 1.98 #ifdef ALLOW_AUTODIFF_TAMC
532     CADJ STORE tice = comlev1, key = ikey_dynamics, byte = isbyte
533     # ifdef SEAICE_MULTICATEGORY
534     CADJ STORE tices = comlev1, key = ikey_dynamics, byte = isbyte
535     # endif /* SEAICE_MULTICATEGORY */
536     #endif /* ALLOW_AUTODIFF_TAMC */
537    
538     C-- Start loop over muli-categories, if SEAICE_MULTICATEGORY is undefined
539     C nDim = 1, and there is only one loop iteration
540     DO IT=1,nDim
541     #ifdef ALLOW_AUTODIFF_TAMC
542 dimitri 1.69 #ifdef SEAICE_MULTICATEGORY
543 mlosch 1.98 C Why do we need this store directive when we have just stored
544     C TICES before the loop?
545     ilockey = (iicekey-1)*nDim + IT
546 dimitri 1.69 CADJ STORE tices(:,:,it,bi,bj) = comlev1_multdim,
547     CADJ & key = ilockey, byte = isbyte
548 mlosch 1.98 #endif /* SEAICE_MULTICATEGORY */
549 dimitri 1.69 #endif /* ALLOW_AUTODIFF_TAMC */
550 mlosch 1.98 pFac = (2.0 _d 0*real(IT)-1.0 _d 0)/nDim
551 dimitri 1.69 DO J=1,sNy
552     DO I=1,sNx
553 mlosch 1.98 heffActualP(I,J)= heffActual(I,J)*pFac
554     #ifdef SEAICE_MULTICATEGORY
555 dimitri 1.69 TICE(I,J,bi,bj)=TICES(I,J,IT,bi,bj)
556 mlosch 1.98 #endif /* SEAICE_MULTICATEGORY */
557 dimitri 1.69 ENDDO
558     ENDDO
559 jmc 1.70 CALL SEAICE_SOLVE4TEMP(
560 gforget 1.83 I UG, heffActualP, hsnowActual,
561 dimitri 1.69 U TICE,
562 gforget 1.71 O a_QbyATMmult_cover, a_QSWbyATMmult_cover,
563 dimitri 1.69 I bi, bj, myTime, myIter, myThid )
564     DO J=1,sNy
565     DO I=1,sNx
566 gforget 1.75 C average over categories
567 mlosch 1.98 a_QbyATM_cover (I,J) = a_QbyATM_cover(I,J)
568     & + a_QbyATMmult_cover(I,J)/nDim
569     a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
570     & + a_QSWbyATMmult_cover(I,J)/nDim
571     #ifdef SEAICE_MULTICATEGORY
572 dimitri 1.69 TICES(I,J,IT,bi,bj) = TICE(I,J,bi,bj)
573 mlosch 1.98 #endif /* SEAICE_MULTICATEGORY */
574 dimitri 1.69 ENDDO
575     ENDDO
576     ENDDO
577     C-- End loop over multi-categories
578    
579     #ifdef ALLOW_DIAGNOSTICS
580     IF ( useDiagnostics ) THEN
581     IF ( DIAGNOSTICS_IS_ON('SIatmQnt',myThid) ) THEN
582     DO J=1,sNy
583     DO I=1,sNx
584     DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj) * (
585 gforget 1.87 & a_QbyATM_cover(I,J) * AREApreTH(I,J) +
586     & a_QbyATM_open(I,J) * ( ONE - AREApreTH(I,J) ) )
587 dimitri 1.69 ENDDO
588     ENDDO
589     CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmQnt',0,1,3,bi,bj,myThid)
590     ENDIF
591     ENDIF
592     #endif
593    
594 gforget 1.83 c switch heat fluxes from W/m2 to 'effective' ice meters
595     DO J=1,sNy
596     DO I=1,sNx
597 jmc 1.91 a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)
598 gforget 1.87 & * convertQ2HI * AREApreTH(I,J)
599 jmc 1.91 a_QSWbyATM_cover(I,J) = a_QSWbyATM_cover(I,J)
600 gforget 1.87 & * convertQ2HI * AREApreTH(I,J)
601 jmc 1.91 a_QbyATM_open(I,J) = a_QbyATM_open(I,J)
602 gforget 1.87 & * convertQ2HI * ( ONE - AREApreTH(I,J) )
603 jmc 1.91 a_QSWbyATM_open(I,J) = a_QSWbyATM_open(I,J)
604 gforget 1.87 & * convertQ2HI * ( ONE - AREApreTH(I,J) )
605 gforget 1.89 c and initialize r_QbyATM_cover/r_QbyATM_open
606 gforget 1.84 r_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
607 gforget 1.89 r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
608 gforget 1.83 ENDDO
609     ENDDO
610 gforget 1.75
611    
612 jmc 1.91 C determine available heat due to the ice pack tying the
613 gforget 1.75 C underlying surface water temperature to freezing point
614     C ======================================================
615    
616 dimitri 1.69 #ifdef ALLOW_AUTODIFF_TAMC
617     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj,
618     CADJ & key = iicekey, byte = isbyte
619     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
620     CADJ & key = iicekey, byte = isbyte
621     #endif
622 gforget 1.72
623 dimitri 1.69 DO J=1,sNy
624     DO I=1,sNx
625     IF ( .NOT. inAdMode ) THEN
626     #ifdef SEAICE_VARIABLE_FREEZING_POINT
627     TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0
628     #endif /* SEAICE_VARIABLE_FREEZING_POINT */
629     IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN
630 gforget 1.84 a_QbyOCN(i,j) = -SEAICE_availHeatFrac
631 dimitri 1.69 & * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
632 gforget 1.77 & * hFacC(i,j,kSurface,bi,bj) *
633 gforget 1.83 & (HeatCapacity_Cp*rhoConst/QI)
634 dimitri 1.69 ELSE
635 gforget 1.84 a_QbyOCN(i,j) = -SEAICE_availHeatFracFrz
636 dimitri 1.69 & * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
637 gforget 1.77 & * hFacC(i,j,kSurface,bi,bj) *
638 gforget 1.83 & (HeatCapacity_Cp*rhoConst/QI)
639 dimitri 1.69 ENDIF
640     ELSE
641 gforget 1.84 a_QbyOCN(i,j) = 0.
642 dimitri 1.69 ENDIF
643 gforget 1.77 cgf heat and water conservation: ok -- since rid of 72.0764 factor
644 gforget 1.72 ENDDO
645     ENDDO
646    
647 gforget 1.96 #ifdef ALLOW_AUTODIFF_TAMC
648     #ifdef SEAICE_SIMPLIFY_GROWTH_ADJ
649     CALL ZERO_ADJ_1D( sNx*sNy, a_QbyOCN, myThid)
650     #endif
651     #endif
652 gforget 1.88
653     c ===================================================================
654     c =========PART 3: determine effective thicknesses increments========
655     c ===================================================================
656    
657    
658    
659 gforget 1.75 C compute ice thickness tendency due to ice-ocean interaction
660     C ===========================================================
661    
662 gforget 1.72 DO J=1,sNy
663     DO I=1,sNx
664 gforget 1.85 d_HEFFbyOCNonICE(I,J)=MAX(a_QbyOCN(i,j), -HEFF(I,J,bi,bj))
665 gforget 1.84 r_QbyOCN(I,J)=a_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
666     HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J)
667 dimitri 1.69 ENDDO
668     ENDDO
669    
670 gforget 1.96 #ifdef SEAICE_GROWTH_LEGACY
671 dimitri 1.69 #ifdef ALLOW_DIAGNOSTICS
672     IF ( useDiagnostics ) THEN
673     IF ( DIAGNOSTICS_IS_ON('SIyneg ',myThid) ) THEN
674 gforget 1.84 CALL DIAGNOSTICS_FILL(d_HEFFbyOCNonICE,
675 gforget 1.71 & 'SIyneg ',0,1,1,bi,bj,myThid)
676 dimitri 1.69 ENDIF
677     ENDIF
678     #endif
679 gforget 1.96 #endif
680 gforget 1.75
681 gforget 1.89 C compute snow melt tendency due to snow-atmosphere interaction
682 gforget 1.75 C ==================================================================
683    
684 dimitri 1.69 #ifdef ALLOW_AUTODIFF_TAMC
685     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
686     CADJ & key = iicekey, byte = isbyte
687     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
688     CADJ & key = iicekey, byte = isbyte
689 gforget 1.71 CADJ STORE a_QbyATM_cover(:,:) = comlev1_bibj,
690 dimitri 1.69 CADJ & key = iicekey, byte = isbyte
691     #endif /* ALLOW_AUTODIFF_TAMC */
692 gforget 1.75
693 dimitri 1.69 DO J=1,sNy
694     DO I=1,sNx
695 gforget 1.89 IF ( a_QbyATM_cover(I,J).LT. 0. _d 0 ) THEN
696 gforget 1.85 tmpscal1=
697     & MAX(r_QbyATM_cover(I,J)*ICE2SNOW, -HSNOW(I,J,bi,bj))
698 jmc 1.91 ELSE
699 gforget 1.89 tmpscal1=0. _d 0
700 gforget 1.73 ENDIF
701 gforget 1.89 d_HSNWbyATMonSNW(I,J)= tmpscal1
702 jmc 1.91 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal1
703 gforget 1.89 r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal1/ICE2SNOW
704 dimitri 1.69 ENDDO
705     ENDDO
706    
707 gforget 1.89 C compute ice thickness tendency due to the atmosphere
708     C ====================================================
709 gforget 1.75
710 dimitri 1.69 #ifdef ALLOW_AUTODIFF_TAMC
711 gforget 1.89 CADJ STORE AREApreTH(:,:) = comlev1_bibj,
712 dimitri 1.69 CADJ & key = iicekey, byte = isbyte
713     #endif
714 gforget 1.75
715 gforget 1.82 cgf note: this block is not actually tested by lab_sea
716     cgf where all experiments start in January. So even though
717     cgf the v1.81=>v1.82 revision would change results in
718     cgf warming conditions, the lab_sea results were not changed.
719    
720 dimitri 1.69 DO J=1,sNy
721     DO I=1,sNx
722 gforget 1.84 tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))
723 gforget 1.89 d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
724     r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
725     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2
726 gforget 1.74 ENDDO
727     ENDDO
728 dimitri 1.69
729 gforget 1.75
730 jmc 1.91 C attribute precip to fresh water or snow stock,
731 gforget 1.75 C depending on atmospheric conditions.
732     C =================================================
733 dimitri 1.69 #ifdef ALLOW_ATM_TEMP
734 gforget 1.74 DO J=1,sNy
735     DO I=1,sNx
736 gforget 1.89 c possible alternatives to the a_QbyATM_cover criterium
737     c IF (TICE(I,J,bi,bj) .LT. TMIX) THEN
738     c IF (atemp(I,J,bi,bj) .LT. celsius2K) THEN
739     IF ( a_QbyATM_cover(I,J).GE. 0. _d 0 ) THEN
740 gforget 1.73 C add precip as snow
741 gforget 1.84 d_HFRWbyRAIN(I,J)=0. _d 0
742     d_HSNWbyRAIN(I,J)=convertPRECIP2HI*ICE2SNOW*
743 gforget 1.87 & PRECIP(I,J,bi,bj)*AREApreTH(I,J)
744 jmc 1.91 ELSE
745 gforget 1.73 c add precip to the fresh water bucket
746 gforget 1.84 d_HFRWbyRAIN(I,J)=-convertPRECIP2HI*
747 gforget 1.87 & PRECIP(I,J,bi,bj)*AREApreTH(I,J)
748 gforget 1.84 d_HSNWbyRAIN(I,J)=0. _d 0
749 dimitri 1.69 ENDIF
750 jmc 1.91 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
751 dimitri 1.69 ENDDO
752     ENDDO
753 gforget 1.89 cgf note: this does not affect air-sea heat flux,
754 jmc 1.91 cgf since the implied air heat gain to turn
755 gforget 1.89 cgf rain to snow is not a surface process.
756 gforget 1.74 #endif /* ALLOW_ATM_TEMP */
757 dimitri 1.69
758 gforget 1.75
759 gforget 1.89 C compute snow melt due to heat available from ocean.
760 gforget 1.75 C =================================================================
761 dimitri 1.69
762 gforget 1.93 #ifdef ALLOW_AUTODIFF_TAMC
763     CADJ STORE r_QbyOCN(:,:) = comlev1_bibj,
764     CADJ & key = iicekey, byte = isbyte
765     CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,
766     CADJ & key = iicekey, byte = isbyte
767     #endif
768 gforget 1.89 cgf do we need to keep this comment and cpp bracket?
769     c
770 dimitri 1.69 cph( very sensitive bit here by JZ
771     #ifndef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING
772     DO J=1,sNy
773     DO I=1,sNx
774 gforget 1.84 IF( r_QbyOCN(i,j) .LT. ZERO .AND.
775 dimitri 1.69 & HSNOW(I,J,bi,bj) .GT. ZERO ) THEN
776 gforget 1.85 tmpscal2= MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))
777 gforget 1.72 ELSE
778 gforget 1.81 tmpscal2= 0. _d 0
779 dimitri 1.69 ENDIF
780 gforget 1.85 d_HSNWbyOCNonSNW(I,J) = tmpscal2
781 gforget 1.84 r_QbyOCN(I,J)=r_QbyOCN(I,J)
782     & -d_HSNWbyOCNonSNW(I,J)/ICE2SNOW
783 gforget 1.81 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_HSNWbyOCNonSNW(I,J)
784 dimitri 1.69 ENDDO
785     ENDDO
786 gforget 1.75 #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
787 dimitri 1.69 cph)
788    
789 gforget 1.90
790     C gain of new ice over open water
791     C ===============================
792     #ifndef SEAICE_GROWTH_LEGACY
793     #ifdef SEAICE_DO_OPEN_WATER_GROWTH
794     #ifdef ALLOW_AUTODIFF_TAMC
795     CADJ STORE AREApreTH(:,:) = comlev1_bibj,
796     CADJ & key = iicekey, byte = isbyte
797     #endif
798     DO J=1,sNy
799     DO I=1,sNx
800     if ( (r_QbyATM_open(I,J).GT.0. _d 0).AND.
801     & (AREApreTH(I,J).GT.0. _d 0) ) then
802     tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j)
803 jmc 1.91 c at this point r_QbyOCN(i,j)<=0 and represents the heat
804 gforget 1.90 c that is still needed to get to the first layer to freezing point
805     tmpscal2=SWFRACB*(a_QSWbyATM_cover(I,J)
806     & +a_QSWbyATM_open(I,J))
807 jmc 1.91 c SWFRACB*tmpscal2<=0 is the heat (out of qnet) that is not
808 gforget 1.90 c going to the first layer, which favors its freezing
809     tmpscal3=MAX(0. _d 0, tmpscal1-tmpscal2)
810     else
811     tmpscal3=0. _d 0
812     endif
813     d_HEFFbyATMonOCN_open(I,J)=tmpscal3
814     c The distinct d_HEFFbyATMonOCN_open array is only needed for d_AREA computation.
815 jmc 1.91 c For the rest it is treated as another contribution to d_HEFFbyATMonOCN.
816 gforget 1.90 d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
817     r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
818     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
819     ENDDO
820     ENDDO
821     #endif /* SEAICE_DO_OPEN_WATER_GROWTH */
822     #endif /* SEAICE_GROWTH_LEGACY */
823    
824    
825 gforget 1.87 C convert snow to ice if submerged.
826     C =================================
827    
828 gforget 1.89 #ifndef SEAICE_GROWTH_LEGACY
829     c note: in legacy, this process is done at the end
830 gforget 1.87 #ifdef ALLOW_SEAICE_FLOODING
831     IF ( SEAICEuseFlooding ) THEN
832     DO J=1,sNy
833     DO I=1,sNx
834     hDraft = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
835     & +HEFF(I,J,bi,bj)*SEAICE_rhoIce)/rhoConst
836     tmpscal1 = MAX( 0. _d 0, hDraft - HEFF(I,J,bi,bj))
837     d_HEFFbyFLOODING(I,J)=tmpscal1
838     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
839     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
840 jmc 1.91 & d_HEFFbyFLOODING(I,J)*ICE2SNOW
841 gforget 1.87 ENDDO
842     ENDDO
843     ENDIF
844     #endif /* ALLOW_SEAICE_FLOODING */
845     #endif /* SEAICE_GROWTH_LEGACY */
846    
847    
848 gforget 1.75
849 gforget 1.88 c ===================================================================
850     c ==========PART 4: determine ice cover fraction increments=========-
851     c ===================================================================
852    
853    
854 gforget 1.75
855 dimitri 1.69 #ifdef ALLOW_AUTODIFF_TAMC
856 gforget 1.88 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
857     CADJ & key = iicekey, byte = isbyte
858 gforget 1.90 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
859     CADJ & key = iicekey, byte = isbyte
860 gforget 1.88 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
861     CADJ & key = iicekey, byte = isbyte
862     CADJ STORE r_QbyATM_cover(:,:) = comlev1_bibj,
863     CADJ & key = iicekey, byte = isbyte
864 gforget 1.89 CADJ STORE r_QbyATM_open(:,:) = comlev1_bibj,
865     CADJ & key = iicekey, byte = isbyte
866 gforget 1.88 CADJ STORE a_QbyATM_cover(:,:) = comlev1_bibj,
867     CADJ & key = iicekey, byte = isbyte
868     CADJ STORE a_QbyATM_open(:,:) = comlev1_bibj,
869     CADJ & key = iicekey, byte = isbyte
870     CADJ STORE a_QSWbyATM_cover(:,:) = comlev1_bibj,
871     CADJ & key = iicekey, byte = isbyte
872     CADJ STORE a_QSWbyATM_open(:,:) = comlev1_bibj,
873     CADJ & key = iicekey, byte = isbyte
874 dimitri 1.69 #endif /* ALLOW_AUTODIFF_TAMC */
875    
876     DO J=1,sNy
877     DO I=1,sNx
878 gforget 1.89 #ifndef SEAICE_GROWTH_LEGACY
879 gforget 1.90 c compute ice melt due to ATM (and OCN) heat stocks
880 gforget 1.89 c
881 jmc 1.91 # ifdef SEAICE_OCN_MELT_ACT_ON_AREA
882 gforget 1.90 c ice cover reduction by joint OCN+ATM melt
883 jmc 1.91 tmpscal3 = MIN( 0. _d 0 ,
884     & d_HEFFbyATMonOCN(I,J)+d_HEFFbyOCNonICE(I,J) )
885 gforget 1.90 # else
886     c ice cover reduction by ATM melt only -- as in legacy code
887 jmc 1.91 tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN(I,J) )
888     # endif
889 gforget 1.89 C gain of new ice over open water
890     c
891 gforget 1.90 # ifdef SEAICE_DO_OPEN_WATER_GROWTH
892     c the one effectively used to increment HEFF
893     tmpscal4 = d_HEFFbyATMonOCN_open(I,J)
894     # else
895     c the virtual one -- as in legcy code
896     tmpscal4 = MAX(ZERO,a_QbyATM_open(I,J))
897     # endif
898 gforget 1.89 #else
899 gforget 1.88 c compute heff after ice melt by ocn:
900     tmpscal0=HEFF(I,J,bi,bj)
901     & - d_HEFFbyATMonOCN(I,J) - d_HEFFbyFLOODING(I,J)
902     c compute available heat left after snow melt by atm:
903     tmpscal1= a_QbyATM_open(I,J)+a_QbyATM_cover(I,J)
904     & - d_HSNWbyATMonSNW(I,J)/ICE2SNOW
905     c (cannot melt more than all the ice)
906     tmpscal2 = MAX(-tmpscal0,tmpscal1)
907     tmpscal3 = MIN(ZERO,tmpscal2)
908     #ifdef ALLOW_DIAGNOSTICS
909     DIAGarray(I,J) = tmpscal2
910     #endif
911 gforget 1.89 C gain of new ice over open water
912 gforget 1.88 tmpscal4 = MAX(ZERO,a_QbyATM_open(I,J))
913 gforget 1.89 #endif /* SEAICE_GROWTH_LEGACY */
914 gforget 1.88 c compute cover fraction tendency
915     IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
916     d_AREAbyATM(I,J)=tmpscal4/HO_south
917     ELSE
918     d_AREAbyATM(I,J)=tmpscal4/HO
919 mlosch 1.98 ENDIF
920     d_AREAbyATM(I,J)=d_AREAbyATM(I,J)
921 gforget 1.88 #ifndef SEAICE_GROWTH_LEGACY
922 mlosch 1.98 & +HALF*tmpscal3/heffActual(I,J)
923 gforget 1.88 #else
924 mlosch 1.98 & +HALF*tmpscal3*AREApreTH(I,J)
925     & /(tmpscal0+.00001 _d 0)
926 dimitri 1.69 #endif
927 gforget 1.88 c apply tendency
928 gforget 1.90 IF ( (HEFF(i,j,bi,bj).GT.0. _d 0).OR.
929     & (HSNOW(i,j,bi,bj).GT.0. _d 0) ) THEN
930 jmc 1.91 AREA(I,J,bi,bj)=max(0. _d 0 , min( 1. _d 0,
931 gforget 1.88 & AREA(I,J,bi,bj)+d_AREAbyATM(I,J) ) )
932 gforget 1.90 ELSE
933     AREA(I,J,bi,bj)=0. _d 0
934     ENDIF
935 dimitri 1.69 ENDDO
936     ENDDO
937    
938 gforget 1.96 #ifdef SEAICE_GROWTH_LEGACY
939 dimitri 1.69 #ifdef ALLOW_DIAGNOSTICS
940     IF ( useDiagnostics ) THEN
941 gforget 1.88 IF ( DIAGNOSTICS_IS_ON('SIfice ',myThid) ) THEN
942     CALL DIAGNOSTICS_FILL(DIAGarray,'SIfice ',0,1,3,bi,bj,myThid)
943 dimitri 1.69 ENDIF
944     ENDIF
945     #endif
946 gforget 1.96 #endif
947 gforget 1.75
948 dimitri 1.69
949 gforget 1.88 c ===================================================================
950     c =============PART 5: determine ice salinity increments=============
951     c ===================================================================
952    
953    
954     #ifdef ALLOW_ATM_TEMP
955 dimitri 1.69 #ifdef SEAICE_SALINITY
956    
957 gforget 1.87 #ifdef SEAICE_GROWTH_LEGACY
958 gforget 1.88 # ifdef ALLOW_AUTODIFF_TAMC
959     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
960     CADJ & key = iicekey, byte = isbyte
961     # endif /* ALLOW_AUTODIFF_TAMC */
962 dimitri 1.69 DO J=1,sNy
963     DO I=1,sNx
964     C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean
965     IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN
966     saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
967     & HSALT(I,J,bi,bj) / SEAICE_deltaTtherm
968     HSALT(I,J,bi,bj) = 0.0 _d 0
969     ENDIF
970     ENDDO
971     ENDDO
972 gforget 1.87 #endif /* SEAICE_GROWTH_LEGACY */
973 dimitri 1.69
974     #ifdef ALLOW_AUTODIFF_TAMC
975     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
976     CADJ & key = iicekey, byte = isbyte
977     #endif /* ALLOW_AUTODIFF_TAMC */
978    
979     DO J=1,sNy
980     DO I=1,sNx
981 gforget 1.87 c sum up the terms that affect the salt content of the ice pack
982 gforget 1.84 tmpscal1=d_HEFFbyOCNonICE(I,J)+d_HEFFbyATMonOCN(I,J)
983 gforget 1.87 c recompute HEFF before thermodyncamic updates (which is not AREApreTH in legacy code)
984     tmpscal2=HEFF(I,J,bi,bj)-tmpscal1-d_HEFFbyFLOODING(I,J)
985 gforget 1.84 C tmpscal1 > 0 : m of sea ice that is created
986     IF ( tmpscal1 .GE. 0.0 ) THEN
987 dimitri 1.69 saltFlux(I,J,bi,bj) =
988 gforget 1.87 & HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm
989     & *SEAICE_salinity*salt(I,j,kSurface,bi,bj)
990     & *tmpscal1*ICE2WATR*rhoConstFresh
991 dimitri 1.69 #ifdef ALLOW_SALT_PLUME
992     C saltPlumeFlux is defined only during freezing:
993     saltPlumeFlux(I,J,bi,bj)=
994 gforget 1.87 & HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm
995     & *(1-SEAICE_salinity)*salt(I,j,kSurface,bi,bj)
996     & *tmpscal1*ICE2WATR*rhoConstFresh
997 dimitri 1.69 C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
998     IF ( .NOT. SaltPlumeSouthernOcean ) THEN
999     IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )
1000     & saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
1001     ENDIF
1002    
1003     #endif /* ALLOW_SALT_PLUME */
1004 gforget 1.84 C tmpscal1 < 0 : m of sea ice that is melted
1005 dimitri 1.69 ELSE
1006     saltFlux(I,J,bi,bj) =
1007 gforget 1.87 & HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm
1008     & *HSALT(I,J,bi,bj)
1009     & *tmpscal1/tmpscal2
1010 dimitri 1.69 #ifdef ALLOW_SALT_PLUME
1011     saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
1012     #endif /* ALLOW_SALT_PLUME */
1013     ENDIF
1014     C update HSALT based on surface saltFlux
1015     HSALT(I,J,bi,bj) = HSALT(I,J,bi,bj) +
1016     & saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
1017     saltFlux(I,J,bi,bj) =
1018     & saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
1019 gforget 1.87 #ifdef SEAICE_GROWTH_LEGACY
1020 dimitri 1.69 C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean
1021     IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN
1022     saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -
1023     & HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) /
1024     & SEAICE_deltaTtherm
1025     HSALT(I,J,bi,bj) = 0.0 _d 0
1026     #ifdef ALLOW_SALT_PLUME
1027     saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
1028     #endif /* ALLOW_SALT_PLUME */
1029     ENDIF
1030 gforget 1.87 #endif /* SEAICE_GROWTH_LEGACY */
1031 dimitri 1.69 ENDDO
1032     ENDDO
1033     #endif /* SEAICE_SALINITY */
1034     #endif /* ALLOW_ATM_TEMP */
1035    
1036 gforget 1.75
1037    
1038 gforget 1.88 c =======================================================================
1039     c =====LEGACY PART 5.5: treat pathological cases, then do flooding ======
1040     c =======================================================================
1041 dimitri 1.69
1042    
1043 gforget 1.72
1044 gforget 1.87 #ifdef SEAICE_GROWTH_LEGACY
1045    
1046 jmc 1.91 C treat values of ice cover fraction oustide
1047 gforget 1.75 C the [0 1] range, and other such issues.
1048     C ===========================================
1049    
1050 gforget 1.89 cgf note: this part cannot be heat and water conserving
1051 gforget 1.76
1052 dimitri 1.69 #ifdef ALLOW_AUTODIFF_TAMC
1053     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
1054     CADJ & key = iicekey, byte = isbyte
1055     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
1056     CADJ & key = iicekey, byte = isbyte
1057     #endif /* ALLOW_AUTODIFF_TAMC */
1058     DO J=1,sNy
1059     DO I=1,sNx
1060     C NOW SET AREA(I,J,bi,bj)=0 WHERE NO ICE IS
1061     AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj)
1062     & ,HEFF(I,J,bi,bj)/.0001 _d 0)
1063     ENDDO
1064     ENDDO
1065 gforget 1.75
1066 dimitri 1.69 #ifdef ALLOW_AUTODIFF_TAMC
1067     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
1068     CADJ & key = iicekey, byte = isbyte
1069     #endif /* ALLOW_AUTODIFF_TAMC */
1070     DO J=1,sNy
1071     DO I=1,sNx
1072     C NOW TRUNCATE AREA
1073     AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj))
1074     ENDDO
1075     ENDDO
1076 gforget 1.75
1077 dimitri 1.69 #ifdef ALLOW_AUTODIFF_TAMC
1078     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
1079     CADJ & key = iicekey, byte = isbyte
1080     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
1081     CADJ & key = iicekey, byte = isbyte
1082     #endif /* ALLOW_AUTODIFF_TAMC */
1083     DO J=1,sNy
1084     DO I=1,sNx
1085     AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj))
1086     HSNOW(I,J,bi,bj) = MAX(ZERO,HSNOW(I,J,bi,bj))
1087     AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1088     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1089     #ifdef SEAICE_CAP_HEFF
1090 mlosch 1.97 C This is not energy conserving, but at least it conserves fresh water
1091     tmpscal0 = MAX(HEFF(I,J,bi,bj)-MAX_HEFF,0. _d 0)
1092     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) - tmpscal0
1093     EmPmR(I,J,bi,bj) = EmPmR(I,J,bi,bj)
1094     & - tmpscal0*SEAICE_rhoIce/SEAICE_deltaTtherm
1095 dimitri 1.69 #endif /* SEAICE_CAP_HEFF */
1096     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1097     ENDDO
1098     ENDDO
1099    
1100     #ifdef ALLOW_DIAGNOSTICS
1101     IF ( useDiagnostics ) THEN
1102     IF ( DIAGNOSTICS_IS_ON('SIthdgrh',myThid) ) THEN
1103     DO J=1,sNy
1104     DO I=1,sNx
1105 gforget 1.87 tmparr1(I,J) = (HEFF(I,J,bi,bj)-HEFFpreTH(I,J))
1106 dimitri 1.69 & /SEAICE_deltaTtherm
1107     ENDDO
1108     ENDDO
1109 gforget 1.73 CALL DIAGNOSTICS_FILL(tmparr1,'SIthdgrh',0,1,3,bi,bj,myThid)
1110 dimitri 1.69 ENDIF
1111     ENDIF
1112     #endif /* ALLOW_DIAGNOSTICS */
1113    
1114 gforget 1.75
1115     C convert snow to ice if submerged.
1116     C =================================
1117    
1118 dimitri 1.69 #ifdef ALLOW_SEAICE_FLOODING
1119     IF ( SEAICEuseFlooding ) THEN
1120     DO J=1,sNy
1121     DO I=1,sNx
1122     hDraft = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1123 gforget 1.85 & +HEFF(I,J,bi,bj)*SEAICE_rhoIce)/rhoConst
1124 jmc 1.86 tmpscal1 = MAX( 0. _d 0, hDraft - HEFF(I,J,bi,bj))
1125 gforget 1.84 d_HEFFbyFLOODING(I,J)=tmpscal1
1126     HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
1127     HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1128 jmc 1.91 & d_HEFFbyFLOODING(I,J)*ICE2SNOW
1129 dimitri 1.69 ENDDO
1130     ENDDO
1131     #ifdef ALLOW_DIAGNOSTICS
1132     IF ( useDiagnostics ) THEN
1133     IF ( DIAGNOSTICS_IS_ON('SIsnwice',myThid) ) THEN
1134     DO J=1,sNy
1135     DO I=1,sNx
1136 gforget 1.84 tmparr1(I,J) = d_HEFFbyFLOODING(I,J)/SEAICE_deltaTtherm
1137 dimitri 1.69 ENDDO
1138     ENDDO
1139 gforget 1.73 CALL DIAGNOSTICS_FILL(tmparr1,'SIsnwice',0,1,3,bi,bj,myThid)
1140 dimitri 1.69 ENDIF
1141     ENDIF
1142     #endif /* ALLOW_DIAGNOSTICS */
1143     ENDIF
1144     #endif /* ALLOW_SEAICE_FLOODING */
1145    
1146 gforget 1.87 #endif /* SEAICE_GROWTH_LEGACY */
1147 gforget 1.75
1148    
1149 gforget 1.88
1150     c ===================================================================
1151     c ===============PART 6: determine ice age increments================
1152     c ===================================================================
1153 dimitri 1.69
1154 gforget 1.75
1155    
1156 dimitri 1.69 #ifdef SEAICE_AGE
1157     # ifndef SEAICE_AGE_VOL
1158     C Sources and sinks for sea ice age:
1159     C assume that a) freezing: new ice fraction forms with zero age
1160     C b) melting: ice fraction vanishes with current age
1161     DO J=1,sNy
1162     DO I=1,sNx
1163     IF ( AREA(I,J,bi,bj) .GT. 0.15 ) THEN
1164 gforget 1.87 IF ( AREA(i,j,bi,bj) .LT. AREApreTH(i,j) ) THEN
1165 dimitri 1.69 C-- scale effective ice-age to account for ice-age sink associated with melting
1166     IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)
1167 gforget 1.87 & *AREA(i,j,bi,bj)/AREApreTH(i,j)
1168 dimitri 1.69 ENDIF
1169     C-- account for aging:
1170     IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)
1171     & + AREA(i,j,bi,bj) * SEAICE_deltaTtherm
1172     ELSE
1173     IceAge(i,j,bi,bj) = ZERO
1174     ENDIF
1175     ENDDO
1176     ENDDO
1177     # else /* ifdef SEAICE_AGE_VOL */
1178     C Sources and sinks for sea ice age:
1179     C assume that a) freezing: new ice volume forms with zero age
1180     C b) melting: ice volume vanishes with current age
1181     DO J=1,sNy
1182     DO I=1,sNx
1183     C-- compute actual age from effective age:
1184 gforget 1.87 IF (AREApreTH(i,j).GT.0. _d 0) THEN
1185     tmpscal1=IceAge(i,j,bi,bj)/AREApreTH(i,j)
1186 dimitri 1.69 ELSE
1187 gforget 1.87 tmpscal1=0. _d 0
1188 jmc 1.91 ENDIF
1189 gforget 1.87 IF ( (HEFFpreTH(i,j).LT.HEFF(i,j,bi,bj)).AND.
1190 dimitri 1.69 & (AREA(i,j,bi,bj).GT.0.15) ) THEN
1191 gforget 1.87 tmpscal2=tmpscal1*HEFFpreTH(i,j)/
1192 dimitri 1.69 & HEFF(i,j,bi,bj)+SEAICE_deltaTtherm
1193     ELSEIF (AREA(i,j,bi,bj).LE.0.15) THEN
1194 gforget 1.87 tmpscal2=0. _d 0
1195 dimitri 1.69 ELSE
1196 gforget 1.87 tmpscal2=tmpscal1+SEAICE_deltaTtherm
1197 dimitri 1.69 ENDIF
1198     C-- re-scale to effective age:
1199 gforget 1.87 IceAge(i,j,bi,bj) = tmpscal2*AREA(i,j,bi,bj)
1200 dimitri 1.69 ENDDO
1201     ENDDO
1202     # endif /* SEAICE_AGE_VOL */
1203     #endif /* SEAICE_AGE */
1204    
1205 gforget 1.88
1206    
1207     c ===================================================================
1208     c ==============PART 7: determine ocean model forcing================
1209     c ===================================================================
1210    
1211    
1212    
1213 jmc 1.91 C compute net heat flux leaving/entering the ocean,
1214 gforget 1.88 C accounting for the part used in melt/freeze processes
1215     C =====================================================
1216    
1217     DO J=1,sNy
1218     DO I=1,sNx
1219 gforget 1.89 QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + r_QbyATM_open(I,J)
1220 jmc 1.91 & - ( d_HEFFbyOCNonICE(I,J) +
1221 gforget 1.88 & d_HSNWbyOCNonSNW(I,J)/ICE2SNOW +
1222 jmc 1.91 & d_HEFFbyNEG(I,J) +
1223 gforget 1.88 & d_HSNWbyNEG(I,J)/ICE2SNOW )
1224     & * maskC(I,J,kSurface,bi,bj)
1225     QSW(I,J,bi,bj) = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)
1226     ENDDO
1227     ENDDO
1228    
1229     #ifdef ALLOW_DIAGNOSTICS
1230     IF ( useDiagnostics ) THEN
1231     IF ( DIAGNOSTICS_IS_ON('SIqneto ',myThid) ) THEN
1232     DO J=1,sNy
1233     DO I=1,sNx
1234 gforget 1.96 DIAGarray(I,J) = r_QbyATM_open(I,J) * convertHI2Q
1235 gforget 1.88 ENDDO
1236     ENDDO
1237     CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneto ',0,1,3,bi,bj,myThid)
1238     ENDIF
1239     IF ( DIAGNOSTICS_IS_ON('SIqneti ',myThid) ) THEN
1240     DO J=1,sNy
1241     DO I=1,sNx
1242 gforget 1.96 DIAGarray(I,J) = r_QbyATM_cover(I,J) * convertHI2Q
1243 gforget 1.88 ENDDO
1244     ENDDO
1245     CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneti ',0,1,3,bi,bj,myThid)
1246     ENDIF
1247     ENDIF
1248     #endif
1249    
1250 gforget 1.83 c switch heat fluxes from 'effective' ice meters to W/m2
1251 gforget 1.88 C ======================================================
1252 gforget 1.83
1253     DO J=1,sNy
1254     DO I=1,sNx
1255     QNET(I,J,bi,bj) = QNET(I,J,bi,bj)*convertHI2Q
1256     QSW(I,J,bi,bj) = QSW(I,J,bi,bj)*convertHI2Q
1257     ENDDO
1258     ENDDO
1259 jmc 1.91
1260     C compute net fresh water flux leaving/entering
1261 gforget 1.88 C the ocean, accounting for fresh/salt water stocks.
1262     C ==================================================
1263    
1264     #ifdef ALLOW_ATM_TEMP
1265     DO J=1,sNy
1266     DO I=1,sNx
1267     tmpscal1= d_HSNWbyATMonSNW(I,J)/ICE2SNOW
1268     & +d_HFRWbyRAIN(I,J)
1269     & +d_HSNWbyOCNonSNW(I,J)/ICE2SNOW
1270     & +d_HEFFbyOCNonICE(I,J)
1271     & +d_HEFFbyATMonOCN(I,J)
1272 jmc 1.91 & +d_HEFFbyNEG(I,J)
1273 gforget 1.88 & +d_HSNWbyNEG(I,J)/ICE2SNOW
1274     EmPmR(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
1275     & ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
1276     & * ( ONE - AREApreTH(I,J) )
1277     #ifdef ALLOW_RUNOFF
1278     & - RUNOFF(I,J,bi,bj)
1279     #endif
1280     & + tmpscal1*convertHI2PRECIP
1281     & )*rhoConstFresh
1282     ENDDO
1283     ENDDO
1284    
1285     #ifdef ALLOW_DIAGNOSTICS
1286     IF ( useDiagnostics ) THEN
1287     IF ( DIAGNOSTICS_IS_ON('SIatmFW ',myThid) ) THEN
1288     DO J=1,sNy
1289     DO I=1,sNx
1290     DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj)*(
1291     & PRECIP(I,J,bi,bj)
1292     & - EVAP(I,J,bi,bj)
1293     & *( ONE - AREApreTH(I,J) )
1294     & + RUNOFF(I,J,bi,bj)
1295     & )*rhoConstFresh
1296     ENDDO
1297     ENDDO
1298     CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmFW ',0,1,3,bi,bj,myThid)
1299     ENDIF
1300     ENDIF
1301     #endif
1302     #ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
1303     DO J=1,sNy
1304     DO I=1,sNx
1305     frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
1306     & PRECIP(I,J,bi,bj)
1307     & - EVAP(I,J,bi,bj)
1308     & *( ONE - AREApreTH(I,J) )
1309     & + RUNOFF(I,J,bi,bj)
1310     & )*rhoConstFresh
1311     ENDDO
1312     ENDDO
1313     #endif
1314     #endif /* ALLOW_ATM_TEMP */
1315    
1316 gforget 1.96 #ifdef SEAICE_DEBUG
1317     CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
1318     CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
1319     CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
1320     #endif /* SEAICE_DEBUG */
1321 gforget 1.88
1322     C Sea Ice Load on the sea surface.
1323     C =================================
1324    
1325     IF ( useRealFreshWaterFlux ) THEN
1326     DO J=1,sNy
1327     DO I=1,sNx
1328 gforget 1.92 #ifdef SEAICE_CAP_ICELOAD
1329     tmpscal1 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
1330     & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1331 jmc 1.94 tmpscal2 = min(tmpscal1,heffTooHeavy*rhoConst)
1332     #else
1333 gforget 1.92 tmpscal2 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
1334     & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1335     #endif
1336     sIceLoad(i,j,bi,bj) = tmpscal2
1337 gforget 1.88 ENDDO
1338     ENDDO
1339     ENDIF
1340    
1341    
1342     C close bi,bj loops
1343 dimitri 1.69 ENDDO
1344     ENDDO
1345    
1346     RETURN
1347     END

  ViewVC Help
Powered by ViewVC 1.1.22