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

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

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

Revision 1.75 - (show annotations) (download)
Fri Oct 1 21:36:45 2010 UTC (13 years, 8 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint62l
Changes since 1.74: +234 -202 lines 
Merging seaice_growth codes -- preliminary steps.

Done with this part!

This part, as intended, has not changed anything to the actual computation
results (lab_sea ad&fwd). It was "simply" a matter untying the many
logical knots and rendering what we have been using understandable.

I will give it a rest for a couple days, and let the testing proceed
unaltered over the week end. Then we will make a tag and move on.

The next part will consist of actually changing the computations, which will
change verification experiments results. We will most likely let go of rigorous
backward compatibility soon after the upcoming tagged version.

Cheers.
1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_growth.F,v 1.72 2010/10/01 01:45:15 gforget Exp $
2 C $Name: $
3
4 #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 CEOP
48
49 C !LOCAL VARIABLES:
50 C === Local variables ===
51 C i,j,bi,bj :: Loop counters
52 INTEGER i, j, bi, bj
53 C number of surface interface layer
54 INTEGER kSurface
55 C constants
56 _RL TBC, SDF, ICE2SNOW
57 _RL QI, recip_QI, QS
58
59 C note: for all heat stocks, the sign convention is positive towards the atmosphere
60 C
61 C a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
62 C the atmosphere and the ocean surface - for ice covered water
63 C a_QbyATM_open :: same but for open water
64 C a_QbyATM :: weighted average (depending on the ice cover fraction)
65 C of this available heat over the grid cell.
66 C r_QbyATM_cover :: residual of a_QbyATM_cover after freezing/melting
67 C processes have been accounted for
68 _RL a_QbyATM_cover (1:sNx,1:sNy)
69 _RL a_QbyATM_open (1:sNx,1:sNy)
70 _RL a_QbyATM (1:sNx,1:sNy)
71 _RL r_QbyATM_cover (1:sNx,1:sNy)
72 C a_QSWbyATM_open - short wave heat flux over ocean in W/m^2
73 C a_QSWbyATM_cover - short wave heat flux under ice in W/m^2
74 _RL a_QSWbyATM_open (1:sNx,1:sNy)
75 _RL a_QSWbyATM_cover (1:sNx,1:sNy)
76 C a_QbyICE :: available heat (in equivalent m of ice) due to
77 C the interaction of the ice pack and the ocean surface
78 C r_QbyICE :: residual of a_QbyICE after freezing/melting
79 C processes have been accounted for
80 _RL a_QbyICE (1:sNx,1:sNy)
81 _RL r_QbyICE (1:sNx,1:sNy)
82
83 c FRWfromSNW :: flag that states whether the atmostpheric conditions are
84 c prone to generate SNOW (FRWfromSNW=1) or fresh water (FRWfromSNW=2)
85 integer FRWfromSNW (1:sNx,1:sNy)
86
87 C Available Heat tendencies associated with melt/freeze processes
88 _RL d_QbyICE (1:sNx,1:sNy)
89 _RL d_QbySNW (1:sNx,1:sNy)
90 _RL d_QbyATMonOCN (1:sNx,1:sNy)
91 _RL d_QbyATMonSNW (1:sNx,1:sNy)
92
93 c ICE/SNOW stocks tendencies associated with the various melt/freeze processes
94 _RL d_AREAbyATM (1:sNx,1:sNy)
95 c
96 _RL d_HEFFbyICEonOCN (1:sNx,1:sNy)
97 _RL d_HEFFbyATMonOCN (1:sNx,1:sNy)
98 _RL d_HEFFfromSNWflood (1:sNx,1:sNy)
99 c
100 _RL d_HSNWfromFRW (1:sNx,1:sNy)
101 _RL d_HSNWbyOCNonSNW (1:sNx,1:sNy)
102 _RL d_SNWintoICEflood (1:sNx,1:sNy)
103 _RL d_HFRWfromSNW (1:sNx,1:sNy)
104
105 C actual ice thickness with upper and lower limit
106 _RL HICE (1:sNx,1:sNy)
107 C actual snow thickness
108 _RL hSnwLoc (1:sNx,1:sNy)
109 C wind speed
110 _RL UG (1:sNx,1:sNy)
111 _RL SPEED_SQ
112 C local copy of AREA
113 _RL areaLoc
114
115 c temporary variables available for the various computations
116 _RL tmpscal1, tmpscal2, tmpscal3, tmpscal4
117 _RL tmparr1 (1:sNx,1:sNy)
118
119 C auxillary variables used for specific processes
120 _RL snowEnergy
121
122 #ifdef ALLOW_SEAICE_FLOODING
123 _RL hDraft
124 #endif /* ALLOW_SEAICE_FLOODING */
125
126 #ifdef SEAICE_SALINITY
127 _RL saltFluxAdjust(1:sNx,1:sNy)
128 #endif
129
130 #ifdef SEAICE_MULTICATEGORY
131 INTEGER it
132 INTEGER ilockey
133 _RL RK
134 _RL HICEP (1:sNx,1:sNy)
135 _RL a_QbyATMmult_cover (1:sNx,1:sNy)
136 _RL a_QSWbyATMmult_cover (1:sNx,1:sNy)
137 #endif
138
139 #ifdef SEAICE_AGE
140 C old_AREA :: hold sea-ice fraction field before any seaice-thermo update
141 _RL old_AREA (1:sNx,1:sNy)
142 # ifdef SEAICE_AGE_VOL
143 C old_HEFF :: hold sea-ice effective thickness field before any seaice-thermo update
144 _RL old_HEFF (1:sNx,1:sNy)
145 _RL age_actual
146 # endif /* SEAICE_AGE_VOL */
147 #endif /* SEAICE_AGE */
148
149 #ifdef ALLOW_DIAGNOSTICS
150 _RL DIAGarray (1:sNx,1:sNy)
151 LOGICAL DIAGNOSTICS_IS_ON
152 EXTERNAL DIAGNOSTICS_IS_ON
153 #endif
154
155 IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
156 kSurface = Nr
157 ELSE
158 kSurface = 1
159 ENDIF
160
161 C FREEZING TEMP. OF SEA WATER (deg C)
162 TBC = SEAICE_freeze
163 C RATIO OF WATER DESITY TO SNOW DENSITY
164 SDF = 1000.0 _d 0/SEAICE_rhoSnow
165 C RATIO OF SEA ICE DENSITY to SNOW DENSITY
166 ICE2SNOW = SDF * ICE2WATR
167 C HEAT OF FUSION OF ICE (J/m^3)
168 QI = 302.0 _d +06
169 recip_QI = 1.0 _d 0 / QI
170 C HEAT OF FUSION OF SNOW (J/m^3)
171 QS = 1.1 _d +08
172
173
174 DO bj=myByLo(myThid),myByHi(myThid)
175 DO bi=myBxLo(myThid),myBxHi(myThid)
176 c
177 #ifdef ALLOW_AUTODIFF_TAMC
178 act1 = bi - myBxLo(myThid)
179 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
180 act2 = bj - myByLo(myThid)
181 max2 = myByHi(myThid) - myByLo(myThid) + 1
182 act3 = myThid - 1
183 max3 = nTx*nTy
184 act4 = ikey_dynamics - 1
185 iicekey = (act1 + 1) + act2*max1
186 & + act3*max1*max2
187 & + act4*max1*max2*max3
188 #endif /* ALLOW_AUTODIFF_TAMC */
189
190 #ifdef ALLOW_AUTODIFF_TAMC
191 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
192 CADJ & key = iicekey, byte = isbyte
193 CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj,
194 CADJ & key = iicekey, byte = isbyte
195 CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj,
196 CADJ & key = iicekey, byte = isbyte
197 #endif /* ALLOW_AUTODIFF_TAMC */
198
199
200 C array initializations
201 C =====================
202
203 DO J=1,sNy
204 DO I=1,sNx
205 a_QbyATM_cover (I,J) = 0.0 _d 0
206 a_QbyATM_open(I,J) = 0.0 _d 0
207 a_QbyATM(I,J) = 0.0 _d 0
208 r_QbyATM_cover (I,J) = 0.0 _d 0
209 c
210 a_QSWbyATM_open (I,J) = 0.0 _d 0
211 a_QSWbyATM_cover (I,J) = 0.0 _d 0
212 c
213 a_QbyICE (I,J) = 0.0 _d 0
214 r_QbyICE (I,J) = 0.0 _d 0
215 c
216 FRWfromSNW (I,J) = 0
217 c
218 d_QbyICE (I,J) = 0.0 _d 0
219 d_QbySNW (I,J) = 0.0 _d 0
220 d_QbyATMonOCN (I,J) = 0.0 _d 0
221 d_QbyATMonSNW (I,J) = 0.0 _d 0
222 c
223 d_AREAbyATM(I,J) = 0.0 _d 0
224 c
225 d_HEFFbyICEonOCN(I,J) = 0.0 _d 0
226 d_HEFFbyATMonOCN(I,J) = 0.0 _d 0
227 d_HEFFfromSNWflood(I,J) = 0.0 _d 0
228 c
229 d_HSNWfromFRW(I,J) = 0.0 _d 0
230 d_HSNWbyOCNonSNW(I,J) = 0.0 _d 0
231 d_SNWintoICEflood(I,J) = 0.0 _d 0
232 d_HFRWfromSNW(I,J) = 0.0 _d 0
233 c
234 tmparr1(I,J) = 0.0 _d 0
235 c
236 #ifdef SEAICE_SALINITY
237 saltFluxAdjust(I,J) = 0.0 _d 0
238 #endif
239 #ifdef SEAICE_MULTICATEGORY
240 a_QbyATMmult_cover(I,J) = 0.0 _d 0
241 a_QSWbyATMmult_cover(I,J) = 0.0 _d 0
242 #endif
243 ENDDO
244 ENDDO
245 DO J=1-oLy,sNy+oLy
246 DO I=1-oLx,sNx+oLx
247 saltWtrIce(I,J,bi,bj) = 0.0 _d 0
248 frWtrIce(I,J,bi,bj) = 0.0 _d 0
249 #ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
250 frWtrAtm(I,J,bi,bj) = 0.0 _d 0
251 #endif
252 ENDDO
253 ENDDO
254
255 #ifdef SEAICE_AGE
256 C store the initial ice fraction over the domain
257 DO J=1,sNy
258 DO I=1,sNx
259 old_AREA(i,j) = AREA(I,J,bi,bj)
260 # ifdef SEAICE_AGE_VOL
261 old_HEFF(i,j) = HEFF(I,J,bi,bj)
262 # endif
263 ENDDO
264 ENDDO
265 #endif /* SEAICE_AGE */
266
267
268 #ifdef ALLOW_AUTODIFF_TAMC
269 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
270 CADJ & key = iicekey, byte = isbyte
271 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
272 CADJ & key = iicekey, byte = isbyte
273 #endif /* ALLOW_AUTODIFF_TAMC */
274 DO J=1,sNy
275 DO I=1,sNx
276 C COMPUTE ACTUAL ICE THICKNESS AND PUT MINIMUM/MAXIMUM
277 C ON ICE THICKNESS FOR BUDGET COMPUTATION
278 C The default of A22 = 0.15 is a common threshold for defining
279 C the ice edge. This ice concentration usually does not occur
280 C due to thermodynamics but due to advection.
281 areaLoc = MAX(A22,AREANm1(I,J,bi,bj))
282 HICE(I,J) = HEFFNm1(I,J,bi,bj)/areaLoc
283 C Do we know what this is for?
284 HICE(I,J) = MAX(HICE(I,J),0.05 _d +00)
285 C Capping the actual ice thickness effectively enforces a
286 C minimum of heat flux through the ice and helps getting rid of
287 C very thick ice.
288 cdm actually, this does exactly the opposite, i.e., ice is thicker
289 cdm when HICE is capped, so I am commenting out
290 cdm HICE(I,J) = MIN(HICE(I,J),9.0 _d +00)
291 hSnwLoc(I,J) = HSNOW(I,J,bi,bj)/areaLoc
292 ENDDO
293 ENDDO
294
295
296 C determine available heat due to the atmosphere -- for open water
297 C ================================================================
298
299 C ocean surface/mixed layer temperature
300 DO J=1,sNy
301 DO I=1,sNx
302 TMIX(I,J,bi,bj)=theta(I,J,kSurface,bi,bj)+273.16 _d +00
303 #ifdef SEAICE_DEBUG
304 TMIX(I,J,bi,bj)=MAX(TMIX(I,J,bi,bj),271.2 _d +00)
305 #endif
306 ENDDO
307 ENDDO
308
309 C wind speed from exf
310 DO J=1,sNy
311 DO I=1,sNx
312 UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))
313 ENDDO
314 ENDDO
315
316 CALL SEAICE_BUDGET_OCEAN(
317 I UG,
318 U TMIX,
319 O a_QbyATM_open, a_QSWbyATM_open,
320 I bi, bj, myTime, myIter, myThid )
321
322
323 C determine available heat due to the atmosphere -- for ice covered water
324 C =======================================================================
325
326 #ifdef ALLOW_AUTODIFF_TAMC
327 CADJ STORE tice = comlev1, key = ikey_dynamics, byte = isbyte
328 # ifdef SEAICE_MULTICATEGORY
329 CADJ STORE tices = comlev1, key = ikey_dynamics, byte = isbyte
330 # endif
331 #endif /* ALLOW_AUTODIFF_TAMC */
332
333 IF (useRelativeWind) THEN
334 C Compute relative wind speed over sea ice.
335 DO J=1,sNy
336 DO I=1,sNx
337 SPEED_SQ =
338 & (uWind(I,J,bi,bj)
339 & +0.5 _d 0*(uVel(i,j,kSurface,bi,bj)
340 & +uVel(i+1,j,kSurface,bi,bj))
341 & -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj)))**2
342 & +(vWind(I,J,bi,bj)
343 & +0.5 _d 0*(vVel(i,j,kSurface,bi,bj)
344 & +vVel(i,j+1,kSurface,bi,bj))
345 & -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj)))**2
346 IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
347 UG(I,J)=SEAICE_EPS
348 ELSE
349 UG(I,J)=SQRT(SPEED_SQ)
350 ENDIF
351 ENDDO
352 ENDDO
353 ENDIF
354
355 #ifdef SEAICE_MULTICATEGORY
356 C-- Start loop over muli-categories
357 DO IT=1,MULTDIM
358 #ifdef ALLOW_AUTODIFF_TAMC
359 ilockey = (iicekey-1)*MULTDIM + IT
360 CADJ STORE tices(:,:,it,bi,bj) = comlev1_multdim,
361 CADJ & key = ilockey, byte = isbyte
362 #endif /* ALLOW_AUTODIFF_TAMC */
363 RK=REAL(IT)
364 DO J=1,sNy
365 DO I=1,sNx
366 HICEP(I,J)=(HICE(I,J)/MULTDIM)*((2.0 _d 0*RK)-1.0 _d 0)
367 TICE(I,J,bi,bj)=TICES(I,J,IT,bi,bj)
368 ENDDO
369 ENDDO
370 CALL SEAICE_SOLVE4TEMP(
371 I UG, HICEP, hSnwLoc,
372 U TICE,
373 O a_QbyATMmult_cover, a_QSWbyATMmult_cover,
374 I bi, bj, myTime, myIter, myThid )
375 DO J=1,sNy
376 DO I=1,sNx
377 C average over categories
378 a_QbyATM_cover (I,J) =
379 & a_QbyATM_cover(I,J) + a_QbyATMmult_cover(I,J)/MULTDIM
380 a_QSWbyATM_cover (I,J) =
381 & a_QSWbyATM_cover(I,J) + a_QSWbyATMmult_cover(I,J)/MULTDIM
382 TICES(I,J,IT,bi,bj) = TICE(I,J,bi,bj)
383 ENDDO
384 ENDDO
385 ENDDO
386 C-- End loop over multi-categories
387 #else /* SEAICE_MULTICATEGORY */
388 CALL SEAICE_SOLVE4TEMP(
389 I UG, HICE, hSnwLoc,
390 U TICE,
391 O a_QbyATM_cover, a_QSWbyATM_cover,
392 I bi, bj, myTime, myIter, myThid )
393 #endif /* SEAICE_MULTICATEGORY */
394
395 #ifdef ALLOW_DIAGNOSTICS
396 IF ( useDiagnostics ) THEN
397 IF ( DIAGNOSTICS_IS_ON('SIatmQnt',myThid) ) THEN
398 DO J=1,sNy
399 DO I=1,sNx
400 DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj) * (
401 & a_QbyATM_cover(I,J) * areaNm1(I,J,bi,bj) +
402 & a_QbyATM_open(I,J) * ( ONE - areaNm1(I,J,bi,bj) ) )
403 ENDDO
404 ENDDO
405 CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmQnt',0,1,3,bi,bj,myThid)
406 ENDIF
407 ENDIF
408 #endif
409
410
411 C determine whether the atmostpheric conditions are prone
412 C to generate SNOW (FRWfromSNW=1) or fresh water (FRWfromSNW=2)
413 C ==============================================================
414
415 DO J=1,sNy
416 DO I=1,sNx
417 IF (a_QbyATM_cover(I,J).LT.ZERO.AND.
418 & AREANm1(I,J,bi,bj).GT.ZERO) THEN
419 FRWfromSNW(I,J)=2
420 ELSE
421 FRWfromSNW(I,J)=1
422 ENDIF
423 ENDDO
424 ENDDO
425
426
427 C determine available heat due to the ice pack tying the
428 C underlying surface water temperature to freezing point
429 C ======================================================
430
431 #ifdef ALLOW_AUTODIFF_TAMC
432 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj,
433 CADJ & key = iicekey, byte = isbyte
434 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
435 CADJ & key = iicekey, byte = isbyte
436 CADJ STORE FRWfromSNW = comlev1_bibj,
437 CADJ & key = iicekey, byte = isbyte
438 #endif
439
440 cgf a_QbyICE: available heat that may be extracted by the ICE acting
441 cgf on the OCN surface/mixed layer to bring it back to freezig point
442 cgf (in m, negative out of the ocean -- different convention than a_QbyATM)
443 cgf
444 cgf Using the same convention (W/m2, positive out of the ocean) as
445 cgf for the other heat stocks breaks the lab_sea test at this point.
446 cgf We should revise these conventions in a later revision though.
447 c
448 DO J=1,sNy
449 DO I=1,sNx
450 IF ( .NOT. inAdMode ) THEN
451 #ifdef SEAICE_VARIABLE_FREEZING_POINT
452 TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0
453 #endif /* SEAICE_VARIABLE_FREEZING_POINT */
454 IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN
455 a_QbyICE(i,j) = SEAICE_availHeatFrac
456 & * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
457 & * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0
458 #ifdef SEAICE_QBYICE_IS_WPERM2
459 c using W/m2 as the Q unit
460 & * ( - QI / SEAICE_deltaTtherm )
461 #endif
462 ELSE
463 a_QbyICE(i,j) = SEAICE_availHeatFracFrz
464 & * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
465 & * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0
466 #ifdef SEAICE_QBYICE_IS_WPERM2
467 & * ( - QI / SEAICE_deltaTtherm )
468 #endif
469 ENDIF
470 ELSE
471 a_QbyICE(i,j) = 0.
472 ENDIF
473 ENDDO
474 ENDDO
475
476
477 C compute ice thickness tendency due to ice-ocean interaction
478 C ===========================================================
479
480 DO J=1,sNy
481 DO I=1,sNx
482 tmpscal1=a_QbyICE(i,j)
483 #ifdef SEAICE_QBYICE_IS_WPERM2
484 & / ( - QI / SEAICE_deltaTtherm )
485 #endif
486 d_HEFFbyICEonOCN(I,J) =
487 & MAX(ZERO, HEFF(I,J,bi,bj)-tmpscal1)- HEFF(I,J,bi,bj)
488 d_QbyICE(I,J)=d_HEFFbyICEonOCN(I,J)
489 #ifdef SEAICE_QBYICE_IS_WPERM2
490 & * ( - QI / SEAICE_deltaTtherm )
491 #endif
492 c apply tendency
493 r_QbyICE(I,J)=a_QbyICE(I,J)+d_QbyICE(I,J)
494 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyICEonOCN(I,J)
495 saltWtrIce(I,J,bi,bj) = saltWtrIce(I,J,bi,bj)
496 & + d_HEFFbyICEonOCN(I,J)
497 ENDDO
498 ENDDO
499
500 #ifdef ALLOW_DIAGNOSTICS
501 IF ( useDiagnostics ) THEN
502 IF ( DIAGNOSTICS_IS_ON('SIyneg ',myThid) ) THEN
503 CALL DIAGNOSTICS_FILL(d_HEFFbyICEonOCN,
504 & 'SIyneg ',0,1,1,bi,bj,myThid)
505 ENDIF
506 ENDIF
507 #endif
508
509
510 C compute snow thickness tendency due to snow-atmosphere interaction
511 C ==================================================================
512
513 #ifdef ALLOW_AUTODIFF_TAMC
514 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
515 CADJ & key = iicekey, byte = isbyte
516 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
517 CADJ & key = iicekey, byte = isbyte
518 CADJ STORE a_QbyATM_cover(:,:) = comlev1_bibj,
519 CADJ & key = iicekey, byte = isbyte
520 #endif /* ALLOW_AUTODIFF_TAMC */
521
522 DO J=1,sNy
523 DO I=1,sNx
524 tmpscal1=-SEAICE_deltaTtherm*
525 & a_QbyATM_cover(I,J)*AREANm1(I,J,bi,bj)
526 IF ( FRWfromSNW(I,J).EQ.2 ) THEN
527 snowEnergy=HSNOW(I,J,bi,bj)*QS
528 IF(tmpscal1.LE.snowEnergy) THEN
529 C not enough heat to melt all snow; use up all of a_QbyATM_cover
530 d_HSNWfromFRW(I,J)=-tmpscal1/QS
531 C SNOW CONVERTED INTO WATER AND THEN INTO equivalent m of ICE melt
532 C The factor 1/ICE2SNOW converts m of snow to m of sea-ice
533 d_HFRWfromSNW(I,J)= - tmpscal1/(QS*ICE2SNOW)
534 d_QbyATMonSNW(I,J) = -a_QbyATM_cover(I,J)
535 ELSE
536 C enough heat to melt snow completely;
537 C compute remaining heat that will melt ice
538 d_QbyATMonSNW(I,J)=-(tmpscal1-snowEnergy)/
539 & SEAICE_deltaTtherm/AREANm1(I,J,bi,bj)-a_QbyATM_cover(I,J)
540 C convert all snow to melt water (fresh water flux)
541 d_HFRWfromSNW(I,J)=-HSNOW(I,J,bi,bj)/ICE2SNOW
542 d_HSNWfromFRW(I,J)=-HSNOW(I,J,bi,bj)
543 END IF
544 ELSE
545 c this process is inactive under cold conditions
546 d_QbyATMonSNW(I,J)=0. _d 0
547 d_HFRWfromSNW(I,J)=0. _d 0
548 d_HSNWfromFRW(I,J)=0. _d 0
549 ENDIF
550 c apply tendency
551 frWtrIce(I,J,bi,bj) = frWtrIce(I,J,bi,bj) +
552 & d_HFRWfromSNW(I,J)
553 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWfromFRW(I,J)
554 a_QbyATM_cover(I,J)= a_QbyATM_cover(I,J) + d_QbyATMonSNW(I,J)
555 ENDDO
556 ENDDO
557
558
559 C compute heat due to the atmosphere that
560 C remain available after melt/freeze processes
561 C ============================================
562
563 #ifdef ALLOW_AUTODIFF_TAMC
564 CADJ STORE a_QbyATM_cover(:,:) = comlev1_bibj,
565 CADJ & key = iicekey, byte = isbyte
566 #endif /* ALLOW_AUTODIFF_TAMC */
567
568 DO J=1,sNy
569 DO I=1,sNx
570 a_QbyATM(I,J)= a_QbyATM_cover(I,J) * AREANm1(I,J,bi,bj)
571 & + a_QbyATM_open(I,J) * (ONE-AREANm1(I,J,bi,bj))
572 ENDDO
573 ENDDO
574
575
576 C compute ice cover fraction tendency from a_QbyATM
577 C =================================================
578
579 #ifdef ALLOW_AUTODIFF_TAMC
580 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
581 CADJ & key = iicekey, byte = isbyte
582 CADJ STORE a_QbyATM_cover(:,:) = comlev1_bibj,
583 CADJ & key = iicekey, byte = isbyte
584 CADJ STORE a_QbyATM(:,:) = comlev1_bibj,
585 CADJ & key = iicekey, byte = isbyte
586 CADJ STORE a_QbyATM_open(:,:) = comlev1_bibj,
587 CADJ & key = iicekey, byte = isbyte
588 CADJ STORE a_QSWbyATM_cover(:,:) = comlev1_bibj,
589 CADJ & key = iicekey, byte = isbyte
590 CADJ STORE a_QSWbyATM_open(:,:) = comlev1_bibj,
591 CADJ & key = iicekey, byte = isbyte
592 #endif /* ALLOW_AUTODIFF_TAMC */
593
594 DO J=1,sNy
595 DO I=1,sNx
596 tmpscal1 = -SEAICE_deltaTtherm*a_QbyATM(I,J)*recip_QI
597 c (cannot melt more than all the ice)
598 tmpscal2 = -ONE*MIN(HEFF(I,J,bi,bj),tmpscal1)
599 tmpscal3 = MIN(ZERO,tmpscal2)
600 #ifdef ALLOW_DIAGNOSTICS
601 DIAGarray(I,J) = tmpscal2
602 #endif
603 C gain of new ice over open water (>0 by definition)
604 tmpscal4 = MAX(ZERO,SEAICE_deltaTtherm*
605 & a_QbyATM_open(I,J)*recip_QI)
606 c compute cover fraction tendency
607 IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
608 d_AREAbyATM(I,J)=
609 & (ONE-AREANm1(I,J,bi,bj))*tmpscal4/HO_south
610 & +HALF*tmpscal3*AREANm1(I,J,bi,bj)
611 & /(HEFF(I,J,bi,bj)+.00001 _d 0)
612 ELSE
613 d_AREAbyATM(I,J)=
614 & (ONE-AREANm1(I,J,bi,bj))*tmpscal4/HO
615 & +HALF*tmpscal3*AREANm1(I,J,bi,bj)
616 & /(HEFF(I,J,bi,bj)+.00001 _d 0)
617 ENDIF
618 c apply tendency
619 AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+d_AREAbyATM(I,J)
620 ENDDO
621 ENDDO
622
623 #ifdef ALLOW_DIAGNOSTICS
624 IF ( useDiagnostics ) THEN
625 IF ( DIAGNOSTICS_IS_ON('SIfice ',myThid) ) THEN
626 CALL DIAGNOSTICS_FILL(DIAGarray,'SIfice ',0,1,3,bi,bj,myThid)
627 ENDIF
628 ENDIF
629 #endif
630
631
632 C compute ice thickness tendency due to the atmosphere.
633 C Freezing ocean water/melting ice also affects
634 C the salt water and heat stocks.
635 C =================================================
636
637 #ifdef ALLOW_AUTODIFF_TAMC
638 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
639 CADJ & key = iicekey, byte = isbyte
640 #endif
641
642 DO J=1,sNy
643 DO I=1,sNx
644 tmpscal1 = SEAICE_deltaTtherm*
645 & a_QbyATM_cover(I,J)*recip_QI*AREANm1(I,J,bi,bj)
646 C (cannot melt more than all the ice)
647 tmpscal2 = MAX(-HEFF(I,J,bi,bj),tmpscal1)
648 d_HEFFbyATMonOCN(I,J)=tmpscal2
649 C compute the r_QbyATM_cover residual as the difference between
650 C the available heat tmpscal1 and the used tmpscal2;
651 d_QbyATMonOCN(I,J)=(tmpscal2 - tmpscal1)
652 & *QI/SEAICE_deltaTtherm-a_QbyATM_cover(I,J)
653 c apply tendency
654 r_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)+d_QbyATMonOCN(I,J)
655 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + d_HEFFbyATMonOCN(I,J)
656 saltWtrIce(I,J,bi,bj) = saltWtrIce(I,J,bi,bj) + tmpscal2
657 ENDDO
658 ENDDO
659
660
661 C attribute precip to fresh water or snow stock,
662 C depending on atmospheric conditions.
663 C =================================================
664 #ifdef ALLOW_ATM_TEMP
665 DO J=1,sNy
666 DO I=1,sNx
667 IF ( FRWfromSNW(I,J).EQ.1 ) THEN
668 C add precip as snow
669 d_HFRWfromSNW(I,J)=0. _d 0
670 d_HSNWfromFRW(I,J)=SEAICE_deltaTtherm*
671 & PRECIP(I,J,bi,bj)*AREANm1(I,J,bi,bj)*SDF
672 ELSE
673 c add precip to the fresh water bucket
674 d_HFRWfromSNW(I,J)=-PRECIP(I,J,bi,bj)*AREANm1(I,J,bi,bj)*
675 & SEAICE_deltaTtherm/ICE2WATR
676 d_HSNWfromFRW(I,J)=0. _d 0
677 ENDIF
678 c apply tendency
679 frWtrIce(I,J,bi,bj) = frWtrIce(I,J,bi,bj) +
680 & d_HFRWfromSNW(I,J)
681 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWfromFRW(I,J)
682 ENDDO
683 ENDDO
684 #endif /* ALLOW_ATM_TEMP */
685
686
687 C compute snow tendency due to heat available from atmosphere
688 C can only reduce HSNOW. melted snow is added to fresh water stock.
689 C =================================================================
690
691 cph( very sensitive bit here by JZ
692 #ifndef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING
693 DO J=1,sNy
694 DO I=1,sNx
695 tmpscal1=r_QbyICE(i,j)
696 #ifdef SEAICE_QBYICE_IS_WPERM2
697 & / ( - QI / SEAICE_deltaTtherm )
698 #endif
699 IF( tmpscal1 .GT. ZERO .AND.
700 & HSNOW(I,J,bi,bj) .GT. ZERO ) THEN
701 d_HSNWbyOCNonSNW(I,J) =
702 & - MIN( HSNOW(I,J,bi,bj)/SDF/ICE2WATR , tmpscal1 )
703 ELSE
704 d_HSNWbyOCNonSNW(I,J) = 0. _d 0
705 ENDIF
706 d_QbySNW(I,J)=d_HSNWbyOCNonSNW(I,J)
707 #ifdef SEAICE_QBYICE_IS_WPERM2
708 & * ( - QI / SEAICE_deltaTtherm )
709 #endif
710 c apply tendency
711 r_QbyICE(I,J)=a_QbyICE(I,J)+d_QbySNW(I,J)
712 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)
713 & +d_HSNWbyOCNonSNW(I,J)*SDF*ICE2WATR
714 frWtrIce(I,J,bi,bj) = frWtrIce(I,J,bi,bj)
715 & +d_HSNWbyOCNonSNW(I,J)
716 ENDDO
717 ENDDO
718 #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
719 cph)
720
721
722 C compute net fresh water flux leaving/entering
723 C the ocean, accounting for fresh/salt water stocks.
724 C ==================================================
725
726 #ifdef ALLOW_AUTODIFF_TAMC
727 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
728 CADJ & key = iicekey, byte = isbyte
729 # ifdef SEAICE_SALINITY
730 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
731 CADJ & key = iicekey, byte = isbyte
732 # endif /* SEAICE_SALINITY */
733 #endif /* ALLOW_AUTODIFF_TAMC */
734
735 #ifdef ALLOW_ATM_TEMP
736 DO J=1,sNy
737 DO I=1,sNx
738 C NOW GET FRESH WATER FLUX
739 EmPmR(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
740 & ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
741 & * ( ONE - AREANm1(I,J,bi,bj) )
742 #ifdef ALLOW_RUNOFF
743 & - RUNOFF(I,J,bi,bj)
744 #endif
745 & + frWtrIce(I,J,bi,bj)*ICE2WATR/SEAICE_deltaTtherm
746 & + saltWtrIce(I,J,bi,bj)*ICE2WATR/SEAICE_deltaTtherm
747 & )*rhoConstFresh
748 ENDDO
749 ENDDO
750
751 #ifdef ALLOW_DIAGNOSTICS
752 IF ( useDiagnostics ) THEN
753 IF ( DIAGNOSTICS_IS_ON('SIatmFW ',myThid) ) THEN
754 DO J=1,sNy
755 DO I=1,sNx
756 DIAGarray(I,J) = maskC(I,J,kSurface,bi,bj)*(
757 & PRECIP(I,J,bi,bj)
758 & - EVAP(I,J,bi,bj)
759 & *( ONE - AREANm1(I,J,bi,bj) )
760 & + RUNOFF(I,J,bi,bj)
761 & )*rhoConstFresh
762 ENDDO
763 ENDDO
764 CALL DIAGNOSTICS_FILL(DIAGarray,'SIatmFW ',0,1,3,bi,bj,myThid)
765 ENDIF
766 ENDIF
767 #endif
768 #ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
769 DO J=1,sNy
770 DO I=1,sNx
771 frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
772 & PRECIP(I,J,bi,bj)
773 & - EVAP(I,J,bi,bj)
774 & *( ONE - AREANm1(I,J,bi,bj) )
775 & + RUNOFF(I,J,bi,bj)
776 & )*rhoConstFresh
777 ENDDO
778 ENDDO
779 #endif
780
781
782 C COMPUTE SURFACE SALT FLUX AND ADJUST ICE SALINITY
783 C ==================================================
784
785 #ifdef SEAICE_SALINITY
786
787 DO J=1,sNy
788 DO I=1,sNx
789 C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean
790 IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN
791 saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
792 & HSALT(I,J,bi,bj) / SEAICE_deltaTtherm
793 HSALT(I,J,bi,bj) = 0.0 _d 0
794 ENDIF
795 ENDDO
796 ENDDO
797
798 #ifdef ALLOW_AUTODIFF_TAMC
799 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
800 CADJ & key = iicekey, byte = isbyte
801 #endif /* ALLOW_AUTODIFF_TAMC */
802
803 DO J=1,sNy
804 DO I=1,sNx
805 C saltWtrIce > 0 : m of sea ice that is created
806 IF ( saltWtrIce(I,J,bi,bj) .GE. 0.0 ) THEN
807 saltFlux(I,J,bi,bj) =
808 & HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*
809 & ICE2WATR*rhoConstFresh*SEAICE_salinity*
810 & salt(I,j,kSurface,bi,bj)/SEAICE_deltaTtherm
811 #ifdef ALLOW_SALT_PLUME
812 C saltPlumeFlux is defined only during freezing:
813 saltPlumeFlux(I,J,bi,bj)=
814 & HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*
815 & ICE2WATR*rhoConstFresh*(1-SEAICE_salinity)*
816 & salt(I,j,kSurface,bi,bj)/SEAICE_deltaTtherm
817 C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
818 IF ( .NOT. SaltPlumeSouthernOcean ) THEN
819 IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )
820 & saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
821 ENDIF
822
823 #endif /* ALLOW_SALT_PLUME */
824 C saltWtrIce < 0 : m of sea ice that is melted
825 ELSE
826 saltFlux(I,J,bi,bj) =
827 & HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*
828 & HSALT(I,J,bi,bj)/
829 & (HEFF(I,J,bi,bj)-saltWtrIce(I,J,bi,bj))/
830 & SEAICE_deltaTtherm
831 #ifdef ALLOW_SALT_PLUME
832 saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
833 #endif /* ALLOW_SALT_PLUME */
834 ENDIF
835 C update HSALT based on surface saltFlux
836 HSALT(I,J,bi,bj) = HSALT(I,J,bi,bj) +
837 & saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
838 saltFlux(I,J,bi,bj) =
839 & saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
840 C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean
841 IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN
842 saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -
843 & HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) /
844 & SEAICE_deltaTtherm
845 HSALT(I,J,bi,bj) = 0.0 _d 0
846 #ifdef ALLOW_SALT_PLUME
847 saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
848 #endif /* ALLOW_SALT_PLUME */
849 ENDIF
850 ENDDO
851 ENDDO
852 #endif /* SEAICE_SALINITY */
853 #endif /* ALLOW_ATM_TEMP */
854
855
856 C compute net heat flux leaving/entering the ocean,
857 C accounting for the part used in melt/freeze processes
858 C =====================================================
859
860 DO J=1,sNy
861 DO I=1,sNx
862 QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) * AREANm1(I,J,bi,bj)
863 & +a_QbyATM_open(I,J) * (ONE-AREANm1(I,J,bi,bj))
864 QSW(I,J,bi,bj) = a_QSWbyATM_cover(I,J) * AREANm1(I,J,bi,bj)
865 & +a_QSWbyATM_open(I,J) * (ONE-AREANm1(I,J,bi,bj))
866 ENDDO
867 ENDDO
868
869 #ifdef ALLOW_DIAGNOSTICS
870 IF ( useDiagnostics ) THEN
871 IF ( DIAGNOSTICS_IS_ON('SIqneto ',myThid) ) THEN
872 DO J=1,sNy
873 DO I=1,sNx
874 DIAGarray(I,J) = a_QbyATM_open(I,J) *
875 & (ONE-areaNm1(I,J,bi,bj))
876 ENDDO
877 ENDDO
878 CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneto ',0,1,3,bi,bj,myThid)
879 ENDIF
880 IF ( DIAGNOSTICS_IS_ON('SIqneti ',myThid) ) THEN
881 DO J=1,sNy
882 DO I=1,sNx
883 DIAGarray(I,J) = r_QbyATM_cover(I,J) * areaNm1(I,J,bi,bj)
884 ENDDO
885 ENDDO
886 CALL DIAGNOSTICS_FILL(DIAGarray,'SIqneti ',0,1,3,bi,bj,myThid)
887 ENDIF
888 ENDIF
889 #endif
890
891
892 C account for contribution due to ocean-ice interaction.
893 C ======================================================
894
895 DO J=1,sNy
896 DO I=1,sNx
897 #ifdef SEAICE_QBYICE_IS_WPERM2
898 tmpscal1=
899 & -maskC(I,J,kSurface,bi,bj)
900 & *(HeatCapacity_Cp*rUnit2mass/QI)
901 & *72.0764 _d 0
902 QNET(I,J,bi,bj) = QNET(I,J,bi,bj)
903 & - ( d_QbyICE(I,J) + d_QbySNW(I,J) )
904 & * tmpscal1
905 #else
906 tmpscal1 =
907 & - ( d_HEFFbyICEonOCN(I,J)+d_HSNWbyOCNonSNW(I,J) )
908 & *recip_dRf(kSurface)*recip_hFacC(i,j,kSurface,bi,bj)
909 & *72.0764 _d 0
910 QNET(I,J,bi,bj) = QNET(I,J,bi,bj)
911 & +tmpscal1
912 & /SEAICE_deltaTtherm*maskC(I,J,kSurface,bi,bj)
913 & *HeatCapacity_Cp*rUnit2mass
914 & *drF(kSurface)*hFacC(i,j,kSurface,bi,bj)
915 #endif
916
917 ENDDO
918 ENDDO
919
920 #ifdef SEAICE_DEBUG
921 CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
922 CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
923 CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
924 #endif /* SEAICE_DEBUG */
925
926
927 C treat values of ice cover fraction oustide
928 C the [0 1] range, and other such issues.
929 C ===========================================
930
931 #ifdef ALLOW_AUTODIFF_TAMC
932 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
933 CADJ & key = iicekey, byte = isbyte
934 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
935 CADJ & key = iicekey, byte = isbyte
936 #endif /* ALLOW_AUTODIFF_TAMC */
937 DO J=1,sNy
938 DO I=1,sNx
939 C NOW SET AREA(I,J,bi,bj)=0 WHERE NO ICE IS
940 AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj)
941 & ,HEFF(I,J,bi,bj)/.0001 _d 0)
942 ENDDO
943 ENDDO
944
945 #ifdef ALLOW_AUTODIFF_TAMC
946 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
947 CADJ & key = iicekey, byte = isbyte
948 #endif /* ALLOW_AUTODIFF_TAMC */
949 DO J=1,sNy
950 DO I=1,sNx
951 C NOW TRUNCATE AREA
952 AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj))
953 ENDDO
954 ENDDO
955
956 #ifdef ALLOW_AUTODIFF_TAMC
957 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
958 CADJ & key = iicekey, byte = isbyte
959 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
960 CADJ & key = iicekey, byte = isbyte
961 #endif /* ALLOW_AUTODIFF_TAMC */
962 DO J=1,sNy
963 DO I=1,sNx
964 AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj))
965 HSNOW(I,J,bi,bj) = MAX(ZERO,HSNOW(I,J,bi,bj))
966 AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
967 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)
968 #ifdef SEAICE_CAP_HEFF
969 HEFF(I,J,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,bi,bj))
970 #endif /* SEAICE_CAP_HEFF */
971 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
972 ENDDO
973 ENDDO
974
975 #ifdef ALLOW_DIAGNOSTICS
976 IF ( useDiagnostics ) THEN
977 IF ( DIAGNOSTICS_IS_ON('SIthdgrh',myThid) ) THEN
978 C use (abuse) tmparr1 to diagnose the total thermodynamic growth rate
979 DO J=1,sNy
980 DO I=1,sNx
981 tmparr1(I,J) = (HEFF(I,J,bi,bj)-HEFFNm1(I,J,bi,bj))
982 & /SEAICE_deltaTtherm
983 ENDDO
984 ENDDO
985 CALL DIAGNOSTICS_FILL(tmparr1,'SIthdgrh',0,1,3,bi,bj,myThid)
986 ENDIF
987 ENDIF
988 #endif /* ALLOW_DIAGNOSTICS */
989
990
991 C convert snow to ice if submerged.
992 C =================================
993
994 #ifdef ALLOW_SEAICE_FLOODING
995 IF ( SEAICEuseFlooding ) THEN
996 DO J=1,sNy
997 DO I=1,sNx
998 hDraft = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
999 & +HEFF(I,J,bi,bj)*SEAICE_rhoIce)/1000. _d 0
1000 tmparr1(I,J) = hDraft - MIN(hDraft,HEFF(I,J,bi,bj))
1001 C
1002 d_HEFFfromSNWflood(I,J)=tmparr1(I,J)
1003 d_SNWintoICEflood(I,J)=MAX(0. _d 0, HSNOW(I,J,bi,bj)
1004 & -tmparr1(I,J)*ICE2SNOW ) - HSNOW(I,J,bi,bj)
1005 c apply tendency
1006 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFfromSNWflood(I,J)
1007 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_SNWintoICEflood(I,J)
1008 ENDDO
1009 ENDDO
1010 #ifdef ALLOW_DIAGNOSTICS
1011 IF ( useDiagnostics ) THEN
1012 IF ( DIAGNOSTICS_IS_ON('SIsnwice',myThid) ) THEN
1013 C turn tmparr1 into a rate
1014 DO J=1,sNy
1015 DO I=1,sNx
1016 tmparr1(I,J) = tmparr1(I,J)/SEAICE_deltaTtherm
1017 ENDDO
1018 ENDDO
1019 CALL DIAGNOSTICS_FILL(tmparr1,'SIsnwice',0,1,3,bi,bj,myThid)
1020 ENDIF
1021 ENDIF
1022 #endif /* ALLOW_DIAGNOSTICS */
1023 ENDIF
1024 #endif /* ALLOW_SEAICE_FLOODING */
1025
1026
1027 C Sea Ice Load on the sea surface.
1028 C =================================
1029
1030 IF ( useRealFreshWaterFlux ) THEN
1031 DO J=1,sNy
1032 DO I=1,sNx
1033 sIceLoad(i,j,bi,bj) = HEFF(I,J,bi,bj)*SEAICE_rhoIce
1034 & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1035 ENDDO
1036 ENDDO
1037 ENDIF
1038
1039
1040 C Sea Ice Age Tracer.
1041 C ===================
1042
1043 #ifdef SEAICE_AGE
1044 # ifndef SEAICE_AGE_VOL
1045 C Sources and sinks for sea ice age:
1046 C assume that a) freezing: new ice fraction forms with zero age
1047 C b) melting: ice fraction vanishes with current age
1048 DO J=1,sNy
1049 DO I=1,sNx
1050 IF ( AREA(I,J,bi,bj) .GT. 0.15 ) THEN
1051 IF ( AREA(i,j,bi,bj) .LT. old_AREA(i,j) ) THEN
1052 C-- scale effective ice-age to account for ice-age sink associated with melting
1053 IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)
1054 & *AREA(i,j,bi,bj)/old_AREA(i,j)
1055 ENDIF
1056 C-- account for aging:
1057 IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj)
1058 & + AREA(i,j,bi,bj) * SEAICE_deltaTtherm
1059 ELSE
1060 IceAge(i,j,bi,bj) = ZERO
1061 ENDIF
1062 ENDDO
1063 ENDDO
1064 # else /* ifdef SEAICE_AGE_VOL */
1065 C Sources and sinks for sea ice age:
1066 C assume that a) freezing: new ice volume forms with zero age
1067 C b) melting: ice volume vanishes with current age
1068 DO J=1,sNy
1069 DO I=1,sNx
1070 C-- compute actual age from effective age:
1071 IF (OLD_AREA(i,j).GT.0. _d 0) THEN
1072 age_actual=IceAge(i,j,bi,bj)/OLD_AREA(i,j)
1073 ELSE
1074 age_actual=0. _d 0
1075 ENDIF
1076 IF ( (OLD_HEFF(i,j).LT.HEFF(i,j,bi,bj)).AND.
1077 & (AREA(i,j,bi,bj).GT.0.15) ) THEN
1078 age_actual=age_actual*OLD_HEFF(i,j)/
1079 & HEFF(i,j,bi,bj)+SEAICE_deltaTtherm
1080 ELSEIF (AREA(i,j,bi,bj).LE.0.15) THEN
1081 age_actual=0. _d 0
1082 ELSE
1083 age_actual=age_actual+SEAICE_deltaTtherm
1084 ENDIF
1085 C-- re-scale to effective age:
1086 IceAge(i,j,bi,bj) = age_actual*AREA(i,j,bi,bj)
1087 ENDDO
1088 ENDDO
1089 # endif /* SEAICE_AGE_VOL */
1090 #endif /* SEAICE_AGE */
1091
1092 ENDDO
1093 ENDDO
1094
1095 RETURN
1096 END
  ViewVC Help
Powered by ViewVC 1.1.22