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

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

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

Revision 1.4 - (show annotations) (download)
Fri Sep 28 19:01:17 2012 UTC (12 years, 10 months ago) by torge
Branch: MAIN
Changes since 1.3: +44 -27 lines 
added missing lines related to a_QbyATM_open and a_QSWbyATM_open
to SEAICE_ITD part of code;
also removing some write and print statements used for debugging
1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_growth.F,v 1.162 2012/03/15 03:07:31 jmc Exp $
2 C $Name: $
3
4 #include "SEAICE_OPTIONS.h"
5 #ifdef ALLOW_EXF
6 # include "EXF_OPTIONS.h"
7 #endif
8
9 CBOP
10 C !ROUTINE: SEAICE_GROWTH
11 C !INTERFACE:
12 SUBROUTINE SEAICE_GROWTH( myTime, myIter, myThid )
13 C !DESCRIPTION: \bv
14 C *==========================================================*
15 C | SUBROUTINE seaice_growth
16 C | o Updata ice thickness and snow depth
17 C *==========================================================*
18 C \ev
19
20 C !USES:
21 IMPLICIT NONE
22 C === Global variables ===
23 #include "SIZE.h"
24 #include "EEPARAMS.h"
25 #include "PARAMS.h"
26 #include "DYNVARS.h"
27 #include "GRID.h"
28 #include "FFIELDS.h"
29 #include "SEAICE_SIZE.h"
30 #include "SEAICE_PARAMS.h"
31 #include "SEAICE.h"
32 #include "SEAICE_TRACER.h"
33 #ifdef ALLOW_EXF
34 # include "EXF_PARAM.h"
35 # include "EXF_FIELDS.h"
36 #endif
37 #ifdef ALLOW_SALT_PLUME
38 # include "SALT_PLUME.h"
39 #endif
40 #ifdef ALLOW_AUTODIFF_TAMC
41 # include "tamc.h"
42 #endif
43
44 C !INPUT/OUTPUT PARAMETERS:
45 C === Routine arguments ===
46 C myTime :: Simulation time
47 C myIter :: Simulation timestep number
48 C myThid :: Thread no. that called this routine.
49 _RL myTime
50 INTEGER myIter, myThid
51 CEOP
52
53 C !FUNCTIONS:
54 #ifdef ALLOW_DIAGNOSTICS
55 LOGICAL DIAGNOSTICS_IS_ON
56 EXTERNAL DIAGNOSTICS_IS_ON
57 #endif
58
59 C !LOCAL VARIABLES:
60 C === Local variables ===
61 c ToM<<< debug seaice_growth
62 C msgBuf :: Informational/error message buffer
63 CHARACTER*(MAX_LEN_MBUF) msgBuf
64 c ToM>>>
65 C
66 C unit/sign convention:
67 C Within the thermodynamic computation all stocks, except HSNOW,
68 C are in 'effective ice meters' units, and >0 implies more ice.
69 C This holds for stocks due to ocean and atmosphere heat,
70 C at the outset of 'PART 2: determine heat fluxes/stocks'
71 C and until 'PART 7: determine ocean model forcing'
72 C This strategy minimizes the need for multiplications/divisions
73 C by ice fraction, heat capacity, etc. The only conversions that
74 C occurs are for the HSNOW (in effective snow meters) and
75 C PRECIP (fresh water m/s).
76 C
77 C HEFF is effective Hice thickness (m3/m2)
78 C HSNOW is Heffective snow thickness (m3/m2)
79 C HSALT is Heffective salt content (g/m2)
80 C AREA is the seaice cover fraction (0<=AREA<=1)
81 C Q denotes heat stocks -- converted to ice stocks (m3/m2) early on
82 C
83 C For all other stocks/increments, such as d_HEFFbyATMonOCN
84 C or a_QbyATM_cover, the naming convention is as follows:
85 C The prefix 'a_' means available, the prefix 'd_' means delta
86 C (i.e. increment), and the prefix 'r_' means residual.
87 C The suffix '_cover' denotes a value for the ice covered fraction
88 C of the grid cell, whereas '_open' is for the open water fraction.
89 C The main part of the name states what ice/snow stock is concerned
90 C (e.g. QbyATM or HEFF), and how it is affected (e.g. d_HEFFbyATMonOCN
91 C is the increment of HEFF due to the ATMosphere extracting heat from the
92 C OCeaN surface, or providing heat to the OCeaN surface).
93
94 C i,j,bi,bj :: Loop counters
95 INTEGER i, j, bi, bj
96 C number of surface interface layer
97 INTEGER kSurface
98 C constants
99 _RL tempFrz, ICE2SNOW, SNOW2ICE
100 _RL QI, QS, recip_QI
101
102 C-- TmixLoc :: ocean surface/mixed-layer temperature (in K)
103 _RL TmixLoc (1:sNx,1:sNy)
104
105 C a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
106 C the atmosphere and the ocean surface - for ice covered water
107 C a_QbyATM_open :: same but for open water
108 C r_QbyATM_cover :: residual of a_QbyATM_cover after freezing/melting processes
109 C r_QbyATM_open :: same but for open water
110 _RL a_QbyATM_cover (1:sNx,1:sNy)
111 _RL a_QbyATM_open (1:sNx,1:sNy)
112 _RL r_QbyATM_cover (1:sNx,1:sNy)
113 _RL r_QbyATM_open (1:sNx,1:sNy)
114 C a_QSWbyATM_open - short wave heat flux over ocean in W/m^2
115 C a_QSWbyATM_cover - short wave heat flux under ice in W/m^2
116 _RL a_QSWbyATM_open (1:sNx,1:sNy)
117 _RL a_QSWbyATM_cover (1:sNx,1:sNy)
118 C a_QbyOCN :: available heat (in in W/m^2) due to the
119 C interaction of the ice pack and the ocean surface
120 C r_QbyOCN :: residual of a_QbyOCN after freezing/melting
121 C processes have been accounted for
122 _RL a_QbyOCN (1:sNx,1:sNy)
123 _RL r_QbyOCN (1:sNx,1:sNy)
124
125 C conversion factors to go from Q (W/m2) to HEFF (ice meters)
126 _RL convertQ2HI, convertHI2Q
127 C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
128 _RL convertPRECIP2HI, convertHI2PRECIP
129
130 #ifdef ALLOW_DIAGNOSTICS
131 C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
132 _RL d_AREAbyATM (1:sNx,1:sNy)
133 _RL d_AREAbyOCN (1:sNx,1:sNy)
134 _RL d_AREAbyICE (1:sNx,1:sNy)
135 #endif
136
137 #ifdef SEAICE_ALLOW_AREA_RELAXATION
138 C ICE/SNOW stocks tendency associated with relaxation towards observation
139 _RL d_AREAbyRLX (1:sNx,1:sNy)
140 c The change of mean ice thickness due to relaxation
141 _RL d_HEFFbyRLX (1:sNx,1:sNy)
142 #endif
143
144 #ifdef SEAICE_ITD
145 c The change of mean ice area due to out-of-bounds values following
146 c sea ice dynamics
147 _RL d_AREAbyNEG (1:sNx,1:sNy)
148 #endif
149 c The change of mean ice thickness due to out-of-bounds values following
150 c sea ice dynamics
151 _RL d_HEFFbyNEG (1:sNx,1:sNy)
152
153 c The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
154 _RL d_HEFFbyOCNonICE (1:sNx,1:sNy)
155
156 c The sum of mean ice thickness increments due to atmospheric fluxes over the open water
157 c fraction and ice-covered fractions of the grid cell
158 _RL d_HEFFbyATMonOCN (1:sNx,1:sNy)
159 c The change of mean ice thickness due to flooding by snow
160 _RL d_HEFFbyFLOODING (1:sNx,1:sNy)
161
162 c The mean ice thickness increments due to atmospheric fluxes over the open water
163 c fraction and ice-covered fractions of the grid cell, respectively
164 _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)
165 _RL d_HEFFbyATMonOCN_cover(1:sNx,1:sNy)
166
167 _RL d_HSNWbyNEG (1:sNx,1:sNy)
168 _RL d_HSNWbyATMonSNW (1:sNx,1:sNy)
169 _RL d_HSNWbyOCNonSNW (1:sNx,1:sNy)
170 _RL d_HSNWbyRAIN (1:sNx,1:sNy)
171
172 _RL d_HFRWbyRAIN (1:sNx,1:sNy)
173 C
174 C a_FWbySublim :: fresh water flux implied by latent heat of
175 C sublimation to atmosphere, same sign convention
176 C as EVAP (positive upward)
177 _RL a_FWbySublim (1:sNx,1:sNy)
178 _RL r_FWbySublim (1:sNx,1:sNy)
179 _RL d_HEFFbySublim (1:sNx,1:sNy)
180 _RL d_HSNWbySublim (1:sNx,1:sNy)
181
182 #if (defined(SEAICE_CAP_SUBLIM) || defined(SEAICE_ITD))
183 C The latent heat flux which will sublimate all snow and ice
184 C over one time step
185 _RL latentHeatFluxMax (1:sNx,1:sNy)
186 _RL latentHeatFluxMaxMult (1:sNx,1:sNy,MULTDIM)
187 #endif
188
189 C actual ice thickness (with upper and lower limit)
190 _RL heffActual (1:sNx,1:sNy)
191 C actual snow thickness
192 _RL hsnowActual (1:sNx,1:sNy)
193 C actual ice thickness (with lower limit only) Reciprocal
194 _RL recip_heffActual (1:sNx,1:sNy)
195 C local value (=1/HO or 1/HO_south)
196 _RL recip_HO
197 C local value (=1/ice thickness)
198 _RL recip_HH
199
200 C AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
201 _RL AREApreTH (1:sNx,1:sNy)
202 _RL HEFFpreTH (1:sNx,1:sNy)
203 _RL HSNWpreTH (1:sNx,1:sNy)
204 #ifdef SEAICE_ITD
205 _RL AREAITDpreTH (1:sNx,1:sNy,1:nITD)
206 _RL HEFFITDpreTH (1:sNx,1:sNy,1:nITD)
207 _RL HSNWITDpreTH (1:sNx,1:sNy,1:nITD)
208 _RL areaFracFactor (1:sNx,1:sNy,1:nITD)
209 _RL heffFracFactor (1:sNx,1:sNy,1:nITD)
210 #endif
211
212 C wind speed
213 _RL UG (1:sNx,1:sNy)
214 #ifdef ALLOW_ATM_WIND
215 _RL SPEED_SQ
216 #endif
217
218 C Regularization values squared
219 _RL area_reg_sq, hice_reg_sq
220
221 C pathological cases thresholds
222 _RL heffTooHeavy
223
224 _RL lhSublim
225
226 C temporary variables available for the various computations
227 _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
228 _RL tmparr1 (1:sNx,1:sNy)
229
230 #ifdef SEAICE_VARIABLE_SALINITY
231 _RL saltFluxAdjust (1:sNx,1:sNy)
232 #endif
233
234 INTEGER ilockey
235 INTEGER it
236 #ifdef SEAICE_ITD
237 INTEGER K
238 #endif
239 _RL pFac
240 _RL ticeInMult (1:sNx,1:sNy,MULTDIM)
241 _RL ticeOutMult (1:sNx,1:sNy,MULTDIM)
242 _RL heffActualMult (1:sNx,1:sNy,MULTDIM)
243 #ifdef SEAICE_ITD
244 _RL hsnowActualMult (1:sNx,1:sNy,MULTDIM)
245 _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
246 #endif
247 _RL a_QbyATMmult_cover (1:sNx,1:sNy,MULTDIM)
248 _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
249 _RL a_FWbySublimMult (1:sNx,1:sNy,MULTDIM)
250 #ifdef SEAICE_ITD
251 _RL r_QbyATMmult_cover (1:sNx,1:sNy,MULTDIM)
252 _RL r_FWbySublimMult (1:sNx,1:sNy,MULTDIM)
253 #endif
254 C Helper variables: reciprocal of some constants
255 _RL recip_multDim
256 _RL recip_deltaTtherm
257 _RL recip_rhoIce
258
259 C Factor by which we increase the upper ocean friction velocity (u*) when
260 C ice is absent in a grid cell (dimensionless)
261 _RL MixedLayerTurbulenceFactor
262
263 #ifdef ALLOW_SITRACER
264 INTEGER iTr
265 CHARACTER*8 diagName
266 #endif
267 #ifdef ALLOW_DIAGNOSTICS
268 c Helper variables for diagnostics
269 _RL DIAGarrayA (1:sNx,1:sNy)
270 _RL DIAGarrayB (1:sNx,1:sNy)
271 _RL DIAGarrayC (1:sNx,1:sNy)
272 _RL DIAGarrayD (1:sNx,1:sNy)
273 #endif
274
275
276 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
277
278 C ===================================================================
279 C =================PART 0: constants and initializations=============
280 C ===================================================================
281
282 IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
283 kSurface = Nr
284 ELSE
285 kSurface = 1
286 ENDIF
287
288 C avoid unnecessary divisions in loops
289 #ifdef SEAICE_ITD
290 CToM SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
291 #endif
292 recip_multDim = SEAICE_multDim
293 recip_multDim = ONE / recip_multDim
294 C above/below: double/single precision calculation of recip_multDim
295 c recip_multDim = 1./float(SEAICE_multDim)
296 recip_deltaTtherm = ONE / SEAICE_deltaTtherm
297 recip_rhoIce = ONE / SEAICE_rhoIce
298
299 C Cutoff for iceload
300 heffTooHeavy=drF(kSurface) / 5. _d 0
301
302 C RATIO OF SEA ICE DENSITY to SNOW DENSITY
303 ICE2SNOW = SEAICE_rhoIce/SEAICE_rhoSnow
304 SNOW2ICE = ONE / ICE2SNOW
305
306 C HEAT OF FUSION OF ICE (J/m^3)
307 QI = SEAICE_rhoIce*SEAICE_lhFusion
308 recip_QI = ONE / QI
309 C HEAT OF FUSION OF SNOW (J/m^3)
310 QS = SEAICE_rhoSnow*SEAICE_lhFusion
311
312 C ICE LATENT HEAT CONSTANT
313 lhSublim = SEAICE_lhEvap + SEAICE_lhFusion
314
315 C regularization constants
316 area_reg_sq = SEAICE_area_reg * SEAICE_area_reg
317 hice_reg_sq = SEAICE_hice_reg * SEAICE_hice_reg
318
319 C conversion factors to go from Q (W/m2) to HEFF (ice meters)
320 convertQ2HI=SEAICE_deltaTtherm/QI
321 convertHI2Q = ONE/convertQ2HI
322 C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
323 convertPRECIP2HI=SEAICE_deltaTtherm*rhoConstFresh/SEAICE_rhoIce
324 convertHI2PRECIP = ONE/convertPRECIP2HI
325
326 DO bj=myByLo(myThid),myByHi(myThid)
327 DO bi=myBxLo(myThid),myBxHi(myThid)
328
329 #ifdef ALLOW_AUTODIFF_TAMC
330 act1 = bi - myBxLo(myThid)
331 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
332 act2 = bj - myByLo(myThid)
333 max2 = myByHi(myThid) - myByLo(myThid) + 1
334 act3 = myThid - 1
335 max3 = nTx*nTy
336 act4 = ikey_dynamics - 1
337 iicekey = (act1 + 1) + act2*max1
338 & + act3*max1*max2
339 & + act4*max1*max2*max3
340 #endif /* ALLOW_AUTODIFF_TAMC */
341
342
343 C array initializations
344 C =====================
345
346 DO J=1,sNy
347 DO I=1,sNx
348 a_QbyATM_cover (I,J) = 0.0 _d 0
349 a_QbyATM_open(I,J) = 0.0 _d 0
350 r_QbyATM_cover (I,J) = 0.0 _d 0
351 r_QbyATM_open (I,J) = 0.0 _d 0
352
353 a_QSWbyATM_open (I,J) = 0.0 _d 0
354 a_QSWbyATM_cover (I,J) = 0.0 _d 0
355
356 a_QbyOCN (I,J) = 0.0 _d 0
357 r_QbyOCN (I,J) = 0.0 _d 0
358
359 #ifdef ALLOW_DIAGNOSTICS
360 d_AREAbyATM(I,J) = 0.0 _d 0
361 d_AREAbyICE(I,J) = 0.0 _d 0
362 d_AREAbyOCN(I,J) = 0.0 _d 0
363 #endif
364
365 #ifdef SEAICE_ALLOW_AREA_RELAXATION
366 d_AREAbyRLX(I,J) = 0.0 _d 0
367 d_HEFFbyRLX(I,J) = 0.0 _d 0
368 #endif
369
370 #ifdef SEAICE_ITD
371 d_AREAbyNEG(I,J) = 0.0 _d 0
372 #endif
373 d_HEFFbyNEG(I,J) = 0.0 _d 0
374 d_HEFFbyOCNonICE(I,J) = 0.0 _d 0
375 d_HEFFbyATMonOCN(I,J) = 0.0 _d 0
376 d_HEFFbyFLOODING(I,J) = 0.0 _d 0
377
378 d_HEFFbyATMonOCN_open(I,J) = 0.0 _d 0
379 d_HEFFbyATMonOCN_cover(I,J) = 0.0 _d 0
380
381 d_HSNWbyNEG(I,J) = 0.0 _d 0
382 d_HSNWbyATMonSNW(I,J) = 0.0 _d 0
383 d_HSNWbyOCNonSNW(I,J) = 0.0 _d 0
384 d_HSNWbyRAIN(I,J) = 0.0 _d 0
385 a_FWbySublim(I,J) = 0.0 _d 0
386 r_FWbySublim(I,J) = 0.0 _d 0
387 d_HEFFbySublim(I,J) = 0.0 _d 0
388 d_HSNWbySublim(I,J) = 0.0 _d 0
389 #ifdef SEAICE_CAP_SUBLIM
390 latentHeatFluxMax(I,J) = 0.0 _d 0
391 #endif
392 c
393 d_HFRWbyRAIN(I,J) = 0.0 _d 0
394
395 tmparr1(I,J) = 0.0 _d 0
396
397 #ifdef SEAICE_VARIABLE_SALINITY
398 saltFluxAdjust(I,J) = 0.0 _d 0
399 #endif
400 DO IT=1,SEAICE_multDim
401 ticeInMult(I,J,IT) = 0.0 _d 0
402 ticeOutMult(I,J,IT) = 0.0 _d 0
403 a_QbyATMmult_cover(I,J,IT) = 0.0 _d 0
404 a_QSWbyATMmult_cover(I,J,IT) = 0.0 _d 0
405 a_FWbySublimMult(I,J,IT) = 0.0 _d 0
406 #ifdef SEAICE_ITD
407 r_QbyATMmult_cover (I,J,IT) = 0.0 _d 0
408 r_FWbySublimMult(I,J,IT) = 0.0 _d 0
409 #endif
410 #if (defined(SEAICE_CAP_SUBLIM) || defined(SEAICE_ITD))
411 latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
412 #endif
413 ENDDO
414 ENDDO
415 ENDDO
416 #if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
417 DO J=1-oLy,sNy+oLy
418 DO I=1-oLx,sNx+oLx
419 frWtrAtm(I,J,bi,bj) = 0.0 _d 0
420 ENDDO
421 ENDDO
422 #endif
423
424
425 C =====================================================================
426 C ===========PART 1: treat pathological cases (post advdiff)===========
427 C =====================================================================
428
429 #ifdef SEAICE_GROWTH_LEGACY
430
431 DO J=1,sNy
432 DO I=1,sNx
433 HEFFpreTH(I,J)=HEFFNM1(I,J,bi,bj)
434 HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
435 AREApreTH(I,J)=AREANM1(I,J,bi,bj)
436 d_HEFFbyNEG(I,J)=0. _d 0
437 d_HSNWbyNEG(I,J)=0. _d 0
438 #ifdef ALLOW_DIAGNOSTICS
439 DIAGarrayA(I,J) = AREANM1(I,J,bi,bj)
440 DIAGarrayB(I,J) = AREANM1(I,J,bi,bj)
441 DIAGarrayC(I,J) = HEFFNM1(I,J,bi,bj)
442 DIAGarrayD(I,J) = HSNOW(I,J,bi,bj)
443 #endif
444 ENDDO
445 ENDDO
446 #ifdef SEAICE_ITD
447 DO K=1,nITD
448 DO J=1,sNy
449 DO I=1,sNx
450 HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)
451 HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)
452 AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)
453 ENDDO
454 ENDDO
455 ENDDO
456 #endif
457
458 #else /* SEAICE_GROWTH_LEGACY */
459
460 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
461 Cgf no dependency through pathological cases treatment
462 IF ( SEAICEadjMODE.EQ.0 ) THEN
463 #endif
464
465 #ifdef SEAICE_ALLOW_AREA_RELAXATION
466 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
467 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
468 C 0) relax sea ice concentration towards observation
469 IF (SEAICE_tauAreaObsRelax .GT. 0.) THEN
470 DO J=1,sNy
471 DO I=1,sNx
472 C d_AREAbyRLX(i,j) = 0. _d 0
473 C d_HEFFbyRLX(i,j) = 0. _d 0
474 IF ( obsSIce(I,J,bi,bj).GT.AREA(I,J,bi,bj)) THEN
475 d_AREAbyRLX(i,j) =
476 & SEAICE_deltaTtherm/SEAICE_tauAreaObsRelax
477 & * (obsSIce(I,J,bi,bj) - AREA(I,J,bi,bj))
478 ENDIF
479 IF ( obsSIce(I,J,bi,bj).GT.0. _d 0 .AND.
480 & AREA(I,J,bi,bj).EQ.0. _d 0) THEN
481 C d_HEFFbyRLX(i,j) = 1. _d 1 * siEps * d_AREAbyRLX(i,j)
482 d_HEFFbyRLX(i,j) = 1. _d 1 * siEps
483 ENDIF
484 #ifdef SEAICE_ITD
485 AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
486 & + d_AREAbyRLX(i,j)
487 HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
488 & + d_HEFFbyRLX(i,j)
489 #endif
490 AREA(I,J,bi,bj) = AREA(I,J,bi,bj) + d_AREAbyRLX(i,j)
491 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + d_HEFFbyRLX(i,j)
492 ENDDO
493 ENDDO
494 ENDIF
495 #endif /* SEAICE_ALLOW_AREA_RELAXATION */
496
497 C 1) treat the case of negative values:
498
499 #ifdef ALLOW_AUTODIFF_TAMC
500 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
501 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
502 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
503 #endif /* ALLOW_AUTODIFF_TAMC */
504 DO J=1,sNy
505 DO I=1,sNx
506 #ifdef SEAICE_ITD
507 DO K=1,nITD
508 tmpscal1=0. _d 0
509 tmpscal2=0. _d 0
510 tmpscal3=0. _d 0
511 tmpscal2=MAX(-HEFFITD(I,J,K,bi,bj),0. _d 0)
512 HEFFITD(I,J,K,bi,bj)=HEFFITD(I,J,K,bi,bj)+tmpscal2
513 d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
514 tmpscal3=MAX(-HSNOWITD(I,J,K,bi,bj),0. _d 0)
515 HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3
516 d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
517 tmpscal1=MAX(-AREAITD(I,J,K,bi,bj),0. _d 0)
518 AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal1
519 d_AREAbyNEG(I,J)=d_AREAbyNEG(I,J)+tmpscal1
520 ENDDO
521 CToM AREA, HEFF, and HSNOW will be updated at end of PART 1
522 C by calling SEAICE_ITD_SUM
523 #else
524 d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
525 d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
526 AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
527 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
528 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
529 #endif
530 ENDDO
531 ENDDO
532
533 C 1.25) treat the case of very thin ice:
534
535 #ifdef ALLOW_AUTODIFF_TAMC
536 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
537 #endif /* ALLOW_AUTODIFF_TAMC */
538 DO J=1,sNy
539 DO I=1,sNx
540 #ifdef SEAICE_ITD
541 DO K=1,nITD
542 #endif
543 tmpscal2=0. _d 0
544 tmpscal3=0. _d 0
545 #ifdef SEAICE_ITD
546 IF (HEFFITD(I,J,K,bi,bj).LE.siEps) THEN
547 tmpscal2=-HEFFITD(I,J,K,bi,bj)
548 tmpscal3=-HSNOWITD(I,J,K,bi,bj)
549 TICES(I,J,K,bi,bj)=celsius2K
550 CToM TICE will be updated at end of Part 1 together with AREA and HEFF
551 ENDIF
552 HEFFITD(I,J,K,bi,bj) =HEFFITD(I,J,K,bi,bj) +tmpscal2
553 HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3
554 #else
555 IF (HEFF(I,J,bi,bj).LE.siEps) THEN
556 tmpscal2=-HEFF(I,J,bi,bj)
557 tmpscal3=-HSNOW(I,J,bi,bj)
558 TICE(I,J,bi,bj)=celsius2K
559 DO IT=1,SEAICE_multDim
560 TICES(I,J,IT,bi,bj)=celsius2K
561 ENDDO
562 ENDIF
563 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal2
564 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
565 #endif
566 d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
567 d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
568 #ifdef SEAICE_ITD
569 ENDDO
570 #endif
571 ENDDO
572 ENDDO
573
574 C 1.5) treat the case of area but no ice/snow:
575
576 #ifdef ALLOW_AUTODIFF_TAMC
577 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
578 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
579 #endif /* ALLOW_AUTODIFF_TAMC */
580 DO J=1,sNy
581 DO I=1,sNx
582 #ifdef SEAICE_ITD
583 DO K=1,nITD
584 IF ((HEFFITD(i,j,k,bi,bj).EQ.0. _d 0).AND.
585 & (HSNOWITD(i,j,k,bi,bj).EQ.0. _d 0))
586 & AREAITD(I,J,K,bi,bj)=0. _d 0
587 ENDDO
588 #else
589 IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
590 & (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
591 #endif
592 ENDDO
593 ENDDO
594
595 C 2) treat the case of very small area:
596
597 #ifndef DISABLE_AREA_FLOOR
598 #ifdef ALLOW_AUTODIFF_TAMC
599 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
600 #endif /* ALLOW_AUTODIFF_TAMC */
601 DO J=1,sNy
602 DO I=1,sNx
603 #ifdef SEAICE_ITD
604 DO K=1,nITD
605 IF ((HEFFITD(i,j,k,bi,bj).GT.0).OR.
606 & (HSNOWITD(i,j,k,bi,bj).GT.0)) THEN
607 CToM SEAICE_area_floor*nITD cannot be allowed to exceed 1
608 C hence use SEAICE_area_floor devided by nITD
609 C (or install a warning in e.g. seaice_readparms.F)
610 AREAITD(I,J,K,bi,bj)=
611 & MAX(AREAITD(I,J,K,bi,bj),SEAICE_area_floor/float(nITD))
612 ENDIF
613 ENDDO
614 #else
615 IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN
616 AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor)
617 ENDIF
618 #endif
619 ENDDO
620 ENDDO
621 #endif /* DISABLE_AREA_FLOOR */
622
623 C 2.5) treat case of excessive ice cover, e.g., due to ridging:
624
625 CToM for SEAICE_ITD this case is treated in SEAICE_ITD_REDIST,
626 C which is called at end of PART 1 below
627 #ifndef SEAICE_ITD
628 #ifdef ALLOW_AUTODIFF_TAMC
629 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
630 #endif /* ALLOW_AUTODIFF_TAMC */
631 DO J=1,sNy
632 DO I=1,sNx
633 #ifdef ALLOW_DIAGNOSTICS
634 DIAGarrayA(I,J) = AREA(I,J,bi,bj)
635 #endif
636 #ifdef ALLOW_SITRACER
637 SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
638 #endif
639 AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)
640 ENDDO
641 ENDDO
642 #endif /* SEAICE_ITD */
643
644 #ifdef SEAICE_ITD
645 CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
646 C first, update AREA and HEFF:
647 CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
648 C
649 DO J=1,sNy
650 DO I=1,sNx
651 C TICES was changed above (item 1.25), now update TICE as ice volume
652 C weighted average of TICES
653 tmpscal1 = 0. _d 0
654 tmpscal2 = 0. _d 0
655 DO K=1,nITD
656 tmpscal1=tmpscal1 + TICES(I,J,K,bi,bj)*HEFFITD(I,J,K,bi,bj)
657 tmpscal2=tmpscal2 + HEFFITD(I,J,K,bi,bj)
658 ENDDO
659 TICE(I,J,bi,bj)=tmpscal1/tmpscal2
660 C lines of item 2.5 that were omitted:
661 C in 2.5 these lines are executed before "ridging" is applied to AREA
662 C hence we execute them here before SEAICE_ITD_REDIST is called
663 C although this means that AREA has not been completely regularized
664 #ifdef ALLOW_DIAGNOSTICS
665 DIAGarrayA(I,J) = AREA(I,J,bi,bj)
666 #endif
667 #ifdef ALLOW_SITRACER
668 SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
669 #endif
670 ENDDO
671 ENDDO
672
673 CToM finally make sure that all categories meet their thickness limits
674 C which includes ridging as in item 2.5
675 C and update AREA, HEFF, and HSNOW
676 CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
677 CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
678
679 #endif
680 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
681 C ENDIF SEAICEadjMODE.EQ.0
682 ENDIF
683 #endif
684
685 C 3) store regularized values of heff, hsnow, area at the onset of thermo.
686 DO J=1,sNy
687 DO I=1,sNx
688 HEFFpreTH(I,J)=HEFF(I,J,bi,bj)
689 HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
690 AREApreTH(I,J)=AREA(I,J,bi,bj)
691 #ifdef ALLOW_DIAGNOSTICS
692 DIAGarrayB(I,J) = AREA(I,J,bi,bj)
693 DIAGarrayC(I,J) = HEFF(I,J,bi,bj)
694 DIAGarrayD(I,J) = HSNOW(I,J,bi,bj)
695 #endif
696 #ifdef ALLOW_SITRACER
697 SItrHEFF(I,J,bi,bj,1)=HEFF(I,J,bi,bj)
698 SItrAREA(I,J,bi,bj,2)=AREA(I,J,bi,bj)
699 #endif
700 ENDDO
701 ENDDO
702 #ifdef SEAICE_ITD
703 DO K=1,nITD
704 DO J=1,sNy
705 DO I=1,sNx
706 HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)
707 HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)
708 AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)
709
710 C memorize areal and volume fraction of each ITD category
711 IF (AREA(I,J,bi,bj).GT.0) THEN
712 areaFracFactor(I,J,K)=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj)
713 ELSE
714 areaFracFactor(I,J,K)=ZERO
715 ENDIF
716 IF (HEFF(I,J,bi,bj).GT.0) THEN
717 heffFracFactor(I,J,K)=HEFFITD(I,J,K,bi,bj)/HEFF(I,J,bi,bj)
718 ELSE
719 heffFracFactor(I,J,K)=ZERO
720 ENDIF
721 ENDDO
722 ENDDO
723 ENDDO
724 C prepare SItrHEFF to be computed as cumulative sum
725 DO K=2,5
726 DO J=1,sNy
727 DO I=1,sNx
728 SItrHEFF(I,J,bi,bj,K)=ZERO
729 ENDDO
730 ENDDO
731 ENDDO
732 C prepare SItrAREA to be computed as cumulative sum
733 DO J=1,sNy
734 DO I=1,sNx
735 SItrAREA(I,J,bi,bj,3)=ZERO
736 ENDDO
737 ENDDO
738 #endif
739
740 C 4) treat sea ice salinity pathological cases
741 #ifdef SEAICE_VARIABLE_SALINITY
742 #ifdef ALLOW_AUTODIFF_TAMC
743 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
744 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
745 #endif /* ALLOW_AUTODIFF_TAMC */
746 DO J=1,sNy
747 DO I=1,sNx
748 IF ( (HSALT(I,J,bi,bj) .LT. 0.0).OR.
749 & (HEFF(I,J,bi,bj) .EQ. 0.0) ) THEN
750 saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
751 & HSALT(I,J,bi,bj) * recip_deltaTtherm
752 HSALT(I,J,bi,bj) = 0.0 _d 0
753 ENDIF
754 ENDDO
755 ENDDO
756 #endif /* SEAICE_VARIABLE_SALINITY */
757
758 #endif /* SEAICE_GROWTH_LEGACY */
759
760 #ifdef ALLOW_DIAGNOSTICS
761 IF ( useDiagnostics ) THEN
762 CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIareaPR',0,1,3,bi,bj,myThid)
763 CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIareaPT',0,1,3,bi,bj,myThid)
764 CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIheffPT',0,1,3,bi,bj,myThid)
765 CALL DIAGNOSTICS_FILL(DIAGarrayD,'SIhsnoPT',0,1,3,bi,bj,myThid)
766 #ifdef ALLOW_SITRACER
767 DO iTr = 1, SItrNumInUse
768 WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTr,'PT'
769 IF (SItrMate(iTr).EQ.'HEFF') THEN
770 CALL DIAGNOSTICS_FRACT_FILL(
771 I SItracer(1-OLx,1-OLy,bi,bj,iTr),HEFF(1-OLx,1-OLy,bi,bj),
772 I ONE, 1, diagName,0,1,2,bi,bj,myThid )
773 ELSE
774 CALL DIAGNOSTICS_FRACT_FILL(
775 I SItracer(1-OLx,1-OLy,bi,bj,iTr),AREA(1-OLx,1-OLy,bi,bj),
776 I ONE, 1, diagName,0,1,2,bi,bj,myThid )
777 ENDIF
778 ENDDO
779 #endif /* ALLOW_SITRACER */
780 ENDIF
781 #endif /* ALLOW_DIAGNOSTICS */
782
783 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
784 Cgf no additional dependency of air-sea fluxes to ice
785 IF ( SEAICEadjMODE.GE.1 ) THEN
786 DO J=1,sNy
787 DO I=1,sNx
788 HEFFpreTH(I,J) = 0. _d 0
789 HSNWpreTH(I,J) = 0. _d 0
790 AREApreTH(I,J) = 0. _d 0
791 ENDDO
792 ENDDO
793 #ifdef SEAICE_ITD
794 DO K=1,nITD
795 DO J=1,sNy
796 DO I=1,sNx
797 HEFFITDpreTH(I,J,K) = 0. _d 0
798 HSNWITDpreTH(I,J,K) = 0. _d 0
799 AREAITDpreTH(I,J,K) = 0. _d 0
800 ENDDO
801 ENDDO
802 ENDDO
803 #endif
804 ENDIF
805 #endif
806
807 #if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
808 DO J=1,sNy
809 DO I=1,sNx
810 AREAforAtmFW(I,J,bi,bj) = AREApreTH(I,J)
811 ENDDO
812 ENDDO
813 #endif
814
815 C 4) COMPUTE ACTUAL ICE/SNOW THICKNESS; USE MIN/MAX VALUES
816 C TO REGULARIZE SEAICE_SOLVE4TEMP/d_AREA COMPUTATIONS
817
818 #ifdef ALLOW_AUTODIFF_TAMC
819 CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
820 CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
821 CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
822 #endif /* ALLOW_AUTODIFF_TAMC */
823 #ifdef SEAICE_ITD
824 DO K=1,nITD
825 #endif
826 DO J=1,sNy
827 DO I=1,sNx
828
829 #ifdef SEAICE_ITD
830 IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN
831 #ifdef SEAICE_GROWTH_LEGACY
832 tmpscal1 = MAX(SEAICE_area_reg/float(nITD),
833 & AREAITDpreTH(I,J,K))
834 hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K)/tmpscal1
835 tmpscal2 = HEFFITDpreTH(I,J,K)/tmpscal1
836 heffActualMult(I,J,K) = MAX(tmpscal2,SEAICE_hice_reg)
837 #else /* SEAICE_GROWTH_LEGACY */
838 cif regularize AREA with SEAICE_area_reg
839 tmpscal1 = SQRT(AREAITDpreTH(I,J,K) * AREAITDpreTH(I,J,K)
840 & + area_reg_sq)
841 cif heffActual calculated with the regularized AREA
842 tmpscal2 = HEFFITDpreTH(I,J,K) / tmpscal1
843 cif regularize heffActual with SEAICE_hice_reg (add lower bound)
844 heffActualMult(I,J,K) = SQRT(tmpscal2 * tmpscal2
845 & + hice_reg_sq)
846 cif hsnowActual calculated with the regularized AREA
847 hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K) / tmpscal1
848 #endif /* SEAICE_GROWTH_LEGACY */
849 cif regularize the inverse of heffActual by hice_reg
850 recip_heffActualMult(I,J,K) = AREAITDpreTH(I,J,K) /
851 & sqrt(HEFFITDpreTH(I,J,K) * HEFFITDpreTH(I,J,K)
852 & + hice_reg_sq)
853 cif Do not regularize when HEFFpreTH = 0
854 ELSE
855 heffActualMult(I,J,K) = ZERO
856 hsnowActualMult(I,J,K) = ZERO
857 recip_heffActualMult(I,J,K) = ZERO
858 ENDIF
859 #else /* SEAICE_ITD */
860 IF (HEFFpreTH(I,J) .GT. ZERO) THEN
861 #ifdef SEAICE_GROWTH_LEGACY
862 tmpscal1 = MAX(SEAICE_area_reg,AREApreTH(I,J))
863 hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1
864 tmpscal2 = HEFFpreTH(I,J)/tmpscal1
865 heffActual(I,J) = MAX(tmpscal2,SEAICE_hice_reg)
866 #else /* SEAICE_GROWTH_LEGACY */
867 cif regularize AREA with SEAICE_area_reg
868 tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)
869 cif heffActual calculated with the regularized AREA
870 tmpscal2 = HEFFpreTH(I,J) / tmpscal1
871 cif regularize heffActual with SEAICE_hice_reg (add lower bound)
872 heffActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
873 cif hsnowActual calculated with the regularized AREA
874 hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1
875 #endif /* SEAICE_GROWTH_LEGACY */
876 cif regularize the inverse of heffActual by hice_reg
877 recip_heffActual(I,J) = AREApreTH(I,J) /
878 & sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)
879 cif Do not regularize when HEFFpreTH = 0
880 ELSE
881 heffActual(I,J) = ZERO
882 hsnowActual(I,J) = ZERO
883 recip_heffActual(I,J) = ZERO
884 ENDIF
885 #endif /* SEAICE_ITD */
886
887 ENDDO
888 ENDDO
889 #ifdef SEAICE_ITD
890 ENDDO
891 #endif
892
893 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
894 CALL ZERO_ADJ_1D( sNx*sNy, heffActual, myThid)
895 CALL ZERO_ADJ_1D( sNx*sNy, hsnowActual, myThid)
896 CALL ZERO_ADJ_1D( sNx*sNy, recip_heffActual, myThid)
897 #endif
898
899 #ifdef SEAICE_CAP_SUBLIM
900 C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
901 C AND SNOW THICKNESS
902 #ifdef SEAICE_ITD
903 DO K=1,nITD
904 #endif
905 DO J=1,sNy
906 DO I=1,sNx
907 c The latent heat flux over the sea ice which
908 c will sublimate all of the snow and ice over one time
909 c step (W/m^2)
910 #ifdef SEAICE_ITD
911 IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN
912 latentHeatFluxMaxMult(I,J,K) = lhSublim*recip_deltaTtherm *
913 & (HEFFITDpreTH(I,J,K)*SEAICE_rhoIce +
914 & HSNWITDpreTH(I,J,K)*SEAICE_rhoSnow)/AREAITDpreTH(I,J,K)
915 ELSE
916 latentHeatFluxMaxMult(I,J,K) = ZERO
917 ENDIF
918 #else
919 IF (HEFFpreTH(I,J) .GT. ZERO) THEN
920 latentHeatFluxMax(I,J) = lhSublim * recip_deltaTtherm *
921 & (HEFFpreTH(I,J) * SEAICE_rhoIce +
922 & HSNWpreTH(I,J) * SEAICE_rhoSnow)/AREApreTH(I,J)
923 ELSE
924 latentHeatFluxMax(I,J) = ZERO
925 ENDIF
926 #endif
927 ENDDO
928 ENDDO
929 #ifdef SEAICE_ITD
930 ENDDO
931 #endif
932 #endif /* SEAICE_CAP_SUBLIM */
933
934 C ===================================================================
935 C ================PART 2: determine heat fluxes/stocks===============
936 C ===================================================================
937
938 C determine available heat due to the atmosphere -- for open water
939 C ================================================================
940
941 DO j=1,sNy
942 DO i=1,sNx
943 C ocean surface/mixed layer temperature
944 TmixLoc(i,j) = theta(i,j,kSurface,bi,bj)+celsius2K
945 C wind speed from exf
946 UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))
947 ENDDO
948 ENDDO
949
950 #ifdef ALLOW_AUTODIFF_TAMC
951 CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
952 CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
953 cCADJ STORE UG = comlev1_bibj, key = iicekey,byte=isbyte
954 cCADJ STORE TmixLoc = comlev1_bibj, key = iicekey,byte=isbyte
955 #endif /* ALLOW_AUTODIFF_TAMC */
956
957 CALL SEAICE_BUDGET_OCEAN(
958 I UG,
959 I TmixLoc,
960 O a_QbyATM_open, a_QSWbyATM_open,
961 I bi, bj, myTime, myIter, myThid )
962
963 C determine available heat due to the atmosphere -- for ice covered water
964 C =======================================================================
965
966 #ifdef ALLOW_ATM_WIND
967 IF (useRelativeWind) THEN
968 C Compute relative wind speed over sea ice.
969 DO J=1,sNy
970 DO I=1,sNx
971 SPEED_SQ =
972 & (uWind(I,J,bi,bj)
973 & +0.5 _d 0*(uVel(i,j,kSurface,bi,bj)
974 & +uVel(i+1,j,kSurface,bi,bj))
975 & -0.5 _d 0*(uice(i,j,bi,bj)+uice(i+1,j,bi,bj)))**2
976 & +(vWind(I,J,bi,bj)
977 & +0.5 _d 0*(vVel(i,j,kSurface,bi,bj)
978 & +vVel(i,j+1,kSurface,bi,bj))
979 & -0.5 _d 0*(vice(i,j,bi,bj)+vice(i,j+1,bi,bj)))**2
980 IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
981 UG(I,J)=SEAICE_EPS
982 ELSE
983 UG(I,J)=SQRT(SPEED_SQ)
984 ENDIF
985 ENDDO
986 ENDDO
987 ENDIF
988 #endif /* ALLOW_ATM_WIND */
989
990 #ifdef ALLOW_AUTODIFF_TAMC
991 CADJ STORE tice(:,:,bi,bj)
992 CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
993 CADJ STORE hsnowActual = comlev1_bibj, key = iicekey, byte = isbyte
994 CADJ STORE heffActual = comlev1_bibj, key = iicekey, byte = isbyte
995 CADJ STORE UG = comlev1_bibj, key = iicekey, byte = isbyte
996 CADJ STORE tices(:,:,:,bi,bj)
997 CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
998 CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
999 CADJ & key = iicekey, byte = isbyte
1000 #endif /* ALLOW_AUTODIFF_TAMC */
1001
1002 C-- Start loop over multi-categories
1003 #ifdef SEAICE_ITD
1004 CToM SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
1005 #endif
1006 DO IT=1,SEAICE_multDim
1007 c homogeneous distribution between 0 and 2 x heffActual
1008 #ifndef SEAICE_ITD
1009 pFac = (2.0 _d 0*real(IT)-1.0 _d 0)*recip_multDim
1010 #endif
1011 DO J=1,sNy
1012 DO I=1,sNx
1013 #ifndef SEAICE_ITD
1014 CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
1015 C (instead of heffActual and latentHeatFluxMax)
1016 heffActualMult(I,J,IT)= heffActual(I,J)*pFac
1017 #ifdef SEAICE_CAP_SUBLIM
1018 latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac
1019 #endif
1020 #endif
1021 ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1022 ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
1023 TICE(I,J,bi,bj) = ZERO
1024 TICES(I,J,IT,bi,bj) = ZERO
1025 ENDDO
1026 ENDDO
1027 ENDDO
1028
1029 #ifdef ALLOW_AUTODIFF_TAMC
1030 CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
1031 CADJ STORE ticeInMult = comlev1_bibj, key = iicekey, byte = isbyte
1032 # ifdef SEAICE_CAP_SUBLIM
1033 CADJ STORE latentHeatFluxMaxMult
1034 CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
1035 # endif
1036 CADJ STORE a_QbyATMmult_cover =
1037 CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1038 CADJ STORE a_QSWbyATMmult_cover =
1039 CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1040 CADJ STORE a_FWbySublimMult =
1041 CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1042 #endif /* ALLOW_AUTODIFF_TAMC */
1043
1044 DO IT=1,SEAICE_multDim
1045 CALL SEAICE_SOLVE4TEMP(
1046 #ifdef SEAICE_ITD
1047 I UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),
1048 #else
1049 I UG, heffActualMult(1,1,IT), hsnowActual,
1050 #endif
1051 #ifdef SEAICE_CAP_SUBLIM
1052 I latentHeatFluxMaxMult(1,1,IT),
1053 #endif
1054 U ticeInMult(1,1,IT), ticeOutMult(1,1,IT),
1055 O a_QbyATMmult_cover(1,1,IT), a_QSWbyATMmult_cover(1,1,IT),
1056 O a_FWbySublimMult(1,1,IT),
1057 I bi, bj, myTime, myIter, myThid )
1058 ENDDO
1059
1060 #ifdef ALLOW_AUTODIFF_TAMC
1061 CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
1062 CADJ STORE ticeOutMult = comlev1_bibj, key = iicekey, byte = isbyte
1063 # ifdef SEAICE_CAP_SUBLIM
1064 CADJ STORE latentHeatFluxMaxMult
1065 CADJ & = comlev1_bibj, key = iicekey, byte = isbyte
1066 # endif
1067 CADJ STORE a_QbyATMmult_cover =
1068 CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1069 CADJ STORE a_QSWbyATMmult_cover =
1070 CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1071 CADJ STORE a_FWbySublimMult =
1072 CADJ & comlev1_bibj, key = iicekey, byte = isbyte
1073 #endif /* ALLOW_AUTODIFF_TAMC */
1074
1075 DO IT=1,SEAICE_multDim
1076 DO J=1,sNy
1077 DO I=1,sNx
1078 C update TICE & TICES
1079 #ifdef SEAICE_ITD
1080 C calculate area weighted mean
1081 C (although the ice's temperature relates to its energy content
1082 C and hence should be averaged weighted by ice volume [heffFracFactor],
1083 C the temperature here is a result of the fluxes through the ice surface
1084 C computed individually for each single category in SEAICE_SOLVE4TEMP
1085 C and hence is averaged area weighted [areaFracFactor])
1086 TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1087 & + ticeOutMult(I,J,IT)*areaFracFactor(I,J,K)
1088 #else
1089 TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
1090 & + ticeOutMult(I,J,IT)*recip_multDim
1091 #endif
1092 TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)
1093 C average over categories
1094 #ifdef SEAICE_ITD
1095 C calculate area weighted mean
1096 C (fluxes are per unit (ice surface) area and are thus area weighted)
1097 a_QbyATM_cover (I,J) = a_QbyATM_cover(I,J)
1098 & + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)
1099 a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1100 & + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)
1101 a_FWbySublim (I,J) = a_FWbySublim(I,J)
1102 & + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,K)
1103 #else
1104 a_QbyATM_cover (I,J) = a_QbyATM_cover(I,J)
1105 & + a_QbyATMmult_cover(I,J,IT)*recip_multDim
1106 a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
1107 & + a_QSWbyATMmult_cover(I,J,IT)*recip_multDim
1108 a_FWbySublim (I,J) = a_FWbySublim(I,J)
1109 & + a_FWbySublimMult(I,J,IT)*recip_multDim
1110 #endif
1111 ENDDO
1112 ENDDO
1113 ENDDO
1114
1115 #ifdef SEAICE_CAP_SUBLIM
1116 # ifdef ALLOW_DIAGNOSTICS
1117 DO J=1,sNy
1118 DO I=1,sNx
1119 c The actual latent heat flux realized by SOLVE4TEMP
1120 DIAGarrayA(I,J) = a_FWbySublim(I,J) * lhSublim
1121 ENDDO
1122 ENDDO
1123 cif The actual vs. maximum latent heat flux
1124 IF ( useDiagnostics ) THEN
1125 CALL DIAGNOSTICS_FILL(DIAGarrayA,
1126 & 'SIactLHF',0,1,3,bi,bj,myThid)
1127 CALL DIAGNOSTICS_FILL(latentHeatFluxMax,
1128 & 'SImaxLHF',0,1,3,bi,bj,myThid)
1129 ENDIF
1130 # endif /* ALLOW_DIAGNOSTICS */
1131 #endif /* SEAICE_CAP_SUBLIM */
1132
1133 #ifdef ALLOW_AUTODIFF_TAMC
1134 CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
1135 CADJ STORE a_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
1136 CADJ STORE a_QSWbyATM_cover= comlev1_bibj, key = iicekey, byte = isbyte
1137 CADJ STORE a_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
1138 CADJ STORE a_QSWbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
1139 CADJ STORE a_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
1140 #endif /* ALLOW_AUTODIFF_TAMC */
1141
1142 C switch heat fluxes from W/m2 to 'effective' ice meters
1143 #ifdef SEAICE_ITD
1144 DO K=1,nITD
1145 DO J=1,sNy
1146 DO I=1,sNx
1147 a_QbyATMmult_cover(I,J,K) = a_QbyATMmult_cover(I,J,K)
1148 & * convertQ2HI * AREAITDpreTH(I,J,K)
1149 a_QSWbyATMmult_cover(I,J,K) = a_QSWbyATMmult_cover(I,J,K)
1150 & * convertQ2HI * AREAITDpreTH(I,J,K)
1151 C and initialize r_QbyATM_cover
1152 r_QbyATMmult_cover(I,J,K)=a_QbyATMmult_cover(I,J,K)
1153 C Convert fresh water flux by sublimation to 'effective' ice meters.
1154 C Negative sublimation is resublimation and will be added as snow.
1155 #ifdef SEAICE_DISABLE_SUBLIM
1156 a_FWbySublimMult(I,J,K) = ZERO
1157 #endif
1158 a_FWbySublimMult(I,J,K) = SEAICE_deltaTtherm*recip_rhoIce
1159 & * a_FWbySublimMult(I,J,K)*AREAITDpreTH(I,J,K)
1160 r_FWbySublimMult(I,J,K)=a_FWbySublimMult(I,J,K)
1161 ENDDO
1162 ENDDO
1163 ENDDO
1164 DO J=1,sNy
1165 DO I=1,sNx
1166 a_QbyATM_open(I,J) = a_QbyATM_open(I,J)
1167 & * convertQ2HI * ( ONE - AREApreTH(I,J) )
1168 a_QSWbyATM_open(I,J) = a_QSWbyATM_open(I,J)
1169 & * convertQ2HI * ( ONE - AREApreTH(I,J) )
1170 C and initialize r_QbyATM_open
1171 r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
1172 ENDDO
1173 ENDDO
1174 #else /* SEAICE_ITD */
1175 DO J=1,sNy
1176 DO I=1,sNx
1177 a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)
1178 & * convertQ2HI * AREApreTH(I,J)
1179 a_QSWbyATM_cover(I,J) = a_QSWbyATM_cover(I,J)
1180 & * convertQ2HI * AREApreTH(I,J)
1181 a_QbyATM_open(I,J) = a_QbyATM_open(I,J)
1182 & * convertQ2HI * ( ONE - AREApreTH(I,J) )
1183 a_QSWbyATM_open(I,J) = a_QSWbyATM_open(I,J)
1184 & * convertQ2HI * ( ONE - AREApreTH(I,J) )
1185 C and initialize r_QbyATM_cover/r_QbyATM_open
1186 r_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
1187 r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
1188 C Convert fresh water flux by sublimation to 'effective' ice meters.
1189 C Negative sublimation is resublimation and will be added as snow.
1190 #ifdef SEAICE_DISABLE_SUBLIM
1191 cgf just for those who may need to omit this term to reproduce old results
1192 a_FWbySublim(I,J) = ZERO
1193 #endif
1194 a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce
1195 & * a_FWbySublim(I,J)*AREApreTH(I,J)
1196 r_FWbySublim(I,J)=a_FWbySublim(I,J)
1197 ENDDO
1198 ENDDO
1199 #endif /* SEAICE_ITD */
1200
1201 #ifdef ALLOW_AUTODIFF_TAMC
1202 CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
1203 CADJ STORE a_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
1204 CADJ STORE a_QSWbyATM_cover= comlev1_bibj, key = iicekey, byte = isbyte
1205 CADJ STORE a_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
1206 CADJ STORE a_QSWbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
1207 CADJ STORE a_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
1208 CADJ STORE r_QbyATM_cover = comlev1_bibj, key = iicekey, byte = isbyte
1209 CADJ STORE r_QbyATM_open = comlev1_bibj, key = iicekey, byte = isbyte
1210 CADJ STORE r_FWbySublim = comlev1_bibj, key = iicekey, byte = isbyte
1211 #endif /* ALLOW_AUTODIFF_TAMC */
1212
1213 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1214 Cgf no additional dependency through ice cover
1215 IF ( SEAICEadjMODE.GE.3 ) THEN
1216 #ifdef SEAICE_ITD
1217 DO K=1,nITD
1218 DO J=1,sNy
1219 DO I=1,sNx
1220 a_QbyATMmult_cover(I,J,K) = 0. _d 0
1221 r_QbyATMmult_cover(I,J,K) = 0. _d 0
1222 a_QSWbyATMmult_cover(I,J,K) = 0. _d 0
1223 ENDDO
1224 ENDDO
1225 ENDDO
1226 #else
1227 DO J=1,sNy
1228 DO I=1,sNx
1229 a_QbyATM_cover(I,J) = 0. _d 0
1230 r_QbyATM_cover(I,J) = 0. _d 0
1231 a_QSWbyATM_cover(I,J) = 0. _d 0
1232 ENDDO
1233 ENDDO
1234 #endif
1235 ENDIF
1236 #endif
1237
1238 C determine available heat due to the ice pack tying the
1239 C underlying surface water temperature to freezing point
1240 C ======================================================
1241
1242 #ifdef ALLOW_AUTODIFF_TAMC
1243 CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
1244 CADJ & key = iicekey, byte = isbyte
1245 CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
1246 CADJ & key = iicekey, byte = isbyte
1247 #endif
1248
1249 DO J=1,sNy
1250 DO I=1,sNx
1251 c FREEZING TEMP. OF SEA WATER (deg C)
1252 tempFrz = SEAICE_tempFrz0 +
1253 & SEAICE_dTempFrz_dS *salt(I,J,kSurface,bi,bj)
1254 c efficiency of turbulent fluxes : dependency to sign of THETA-TBC
1255 IF ( theta(I,J,kSurface,bi,bj) .GE. tempFrz ) THEN
1256 tmpscal1 = SEAICE_mcPheePiston
1257 ELSE
1258 tmpscal1 =SEAICE_frazilFrac*drF(kSurface)/SEAICE_deltaTtherm
1259 ENDIF
1260 c efficiency of turbulent fluxes : dependency to AREA (McPhee cases)
1261 IF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.
1262 & (.NOT.SEAICE_mcPheeStepFunc) ) THEN
1263 MixedLayerTurbulenceFactor = ONE -
1264 & SEAICE_mcPheeTaper * AREApreTH(I,J)
1265 ELSEIF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.
1266 & (SEAICE_mcPheeStepFunc) ) THEN
1267 MixedLayerTurbulenceFactor = ONE - SEAICE_mcPheeTaper
1268 ELSE
1269 MixedLayerTurbulenceFactor = ONE
1270 ENDIF
1271 c maximum turbulent flux, in ice meters
1272 tmpscal2= - (HeatCapacity_Cp*rhoConst * recip_QI)
1273 & * (theta(I,J,kSurface,bi,bj)-tempFrz)
1274 & * SEAICE_deltaTtherm * maskC(i,j,kSurface,bi,bj)
1275 c available turbulent flux
1276 a_QbyOCN(i,j) =
1277 & tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor
1278 r_QbyOCN(i,j) = a_QbyOCN(i,j)
1279 ENDDO
1280 ENDDO
1281
1282 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1283 CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
1284 #endif
1285
1286
1287 C ===================================================================
1288 C =========PART 3: determine effective thicknesses increments========
1289 C ===================================================================
1290
1291 C compute snow/ice tendency due to sublimation
1292 C ============================================
1293
1294 #ifdef ALLOW_AUTODIFF_TAMC
1295 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1296 CADJ STORE r_FWbySublim = comlev1_bibj,key=iicekey,byte=isbyte
1297 #endif /* ALLOW_AUTODIFF_TAMC */
1298 #ifdef SEAICE_ITD
1299 DO K=1,nITD
1300 #endif
1301 DO J=1,sNy
1302 DO I=1,sNx
1303 C First sublimate/deposite snow
1304 tmpscal2 =
1305 #ifdef SEAICE_ITD
1306 & MAX(MIN(r_FWbySublimMult(I,J,K),HSNOWITD(I,J,K,bi,bj)
1307 & *SNOW2ICE),ZERO)
1308 C accumulate change over ITD categories
1309 d_HSNWbySublim(I,J) = d_HSNWbySublim(I,J) - tmpscal2
1310 & *ICE2SNOW
1311 HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal2
1312 & *ICE2SNOW
1313 r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
1314 C keep total up to date, too
1315 r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2
1316 #else
1317 & MAX(MIN(r_FWbySublim(I,J),HSNOW(I,J,bi,bj)*SNOW2ICE),ZERO)
1318 d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
1319 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) - tmpscal2*ICE2SNOW
1320 r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2
1321 #endif
1322 ENDDO
1323 ENDDO
1324 #ifdef ALLOW_AUTODIFF_TAMC
1325 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1326 CADJ STORE r_FWbySublim = comlev1_bibj,key=iicekey,byte=isbyte
1327 #endif /* ALLOW_AUTODIFF_TAMC */
1328 DO J=1,sNy
1329 DO I=1,sNx
1330 C If anything is left, sublimate ice
1331 tmpscal2 =
1332 #ifdef SEAICE_ITD
1333 & MAX(MIN(r_FWbySublimMult(I,J,K),HEFFITD(I,J,K,bi,bj)),ZERO)
1334 C accumulate change over ITD categories
1335 d_HSNWbySublim(I,J) = d_HSNWbySublim(I,J) - tmpscal2
1336 HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) - tmpscal2
1337 r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
1338 C keep total up to date, too
1339 r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2
1340 #else
1341 & MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
1342 d_HEFFbySublim(I,J) = - tmpscal2
1343 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) - tmpscal2
1344 r_FWbySublim(I,J) = r_FWbySublim(I,J) - tmpscal2
1345 #endif
1346 ENDDO
1347 ENDDO
1348 DO J=1,sNy
1349 DO I=1,sNx
1350 C If anything is left, it will be evaporated from the ocean rather than sublimated.
1351 C Since a_QbyATM_cover was computed for sublimation, not simple evaporation, we need to
1352 C remove the fusion part for the residual (that happens to be precisely r_FWbySublim).
1353 #ifdef SEAICE_ITD
1354 a_QbyATMmult_cover(I,J,K) = a_QbyATMmult_cover(I,J,K)
1355 & - r_FWbySublimMult(I,J,K)
1356 r_QbyATMmult_cover(I,J,K) = r_QbyATMmult_cover(I,J,K)
1357 & - r_FWbySublimMult(I,J,K)
1358 ENDDO
1359 ENDDO
1360 C end K loop
1361 ENDDO
1362 C then update totals
1363 DO J=1,sNy
1364 DO I=1,sNx
1365 #endif
1366 a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)-r_FWbySublim(I,J)
1367 r_QbyATM_cover(I,J) = r_QbyATM_cover(I,J)-r_FWbySublim(I,J)
1368 ENDDO
1369 ENDDO
1370 c ToM<<< debug seaice_growth
1371 WRITE(msgBuf,'(A,7F6.2)')
1372 #ifdef SEAICE_ITD
1373 & ' SEAICE_GROWTH: Heff increments 1, HEFFITD = ',
1374 & HEFFITD(20,20,:,bi,bj)
1375 #else
1376 & ' SEAICE_GROWTH: Heff increments 1, HEFF = ',
1377 & HEFF(20,20,bi,bj)
1378 #endif
1379 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1380 & SQUEEZE_RIGHT , myThid)
1381 c ToM>>>
1382
1383 C compute ice thickness tendency due to ice-ocean interaction
1384 C ===========================================================
1385
1386 #ifdef ALLOW_AUTODIFF_TAMC
1387 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1388 CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1389 #endif /* ALLOW_AUTODIFF_TAMC */
1390
1391 #ifdef SEAICE_ITD
1392 DO K=1,nITD
1393 DO J=1,sNy
1394 DO I=1,sNx
1395 C ice growth/melt due to ocean heat is equally distributed under the ice
1396 C and hence weighted by fractional area of each thickness category
1397 tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,K),
1398 & -HEFFITD(I,J,K,bi,bj))
1399 d_HEFFbyOCNonICE(I,J)= d_HEFFbyOCNonICE(I,J) + tmpscal1
1400 r_QbyOCN(I,J) = r_QbyOCN(I,J) - tmpscal1
1401 HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal1
1402 #ifdef ALLOW_SITRACER
1403 SItrHEFF(I,J,bi,bj,2) = SItrHEFF(I,J,bi,bj,2)
1404 & + HEFFITD(I,J,K,bi,bj)
1405 #endif
1406 ENDDO
1407 ENDDO
1408 ENDDO
1409 #else /* SEAICE_ITD */
1410 DO J=1,sNy
1411 DO I=1,sNx
1412 d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj))
1413 r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
1414 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J)
1415 #ifdef ALLOW_SITRACER
1416 SItrHEFF(I,J,bi,bj,2)=HEFF(I,J,bi,bj)
1417 #endif
1418 ENDDO
1419 ENDDO
1420 #endif /* SEAICE_ITD */
1421 c ToM<<< debug seaice_growth
1422 WRITE(msgBuf,'(A,7F6.2)')
1423 #ifdef SEAICE_ITD
1424 & ' SEAICE_GROWTH: Heff increments 2, HEFFITD = ',
1425 & HEFFITD(20,20,:,bi,bj)
1426 #else
1427 & ' SEAICE_GROWTH: Heff increments 2, HEFF = ',
1428 & HEFF(20,20,bi,bj)
1429 #endif
1430 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1431 & SQUEEZE_RIGHT , myThid)
1432 c ToM>>>
1433
1434 C compute snow melt tendency due to snow-atmosphere interaction
1435 C ==================================================================
1436
1437 #ifdef ALLOW_AUTODIFF_TAMC
1438 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1439 CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1440 #endif /* ALLOW_AUTODIFF_TAMC */
1441
1442 #ifdef SEAICE_ITD
1443 DO K=1,nITD
1444 DO J=1,sNy
1445 DO I=1,sNx
1446 C Convert to standard units (meters of ice) rather than to meters
1447 C of snow. This appears to be more robust.
1448 tmpscal1=MAX(r_QbyATMmult_cover(I,J,K),-HSNOWITD(I,J,K,bi,bj)
1449 & *SNOW2ICE)
1450 tmpscal2=MIN(tmpscal1,0. _d 0)
1451 #ifdef SEAICE_MODIFY_GROWTH_ADJ
1452 Cgf no additional dependency through snow
1453 IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1454 #endif
1455 d_HSNWbyATMonSNW(I,J)=d_HSNWbyATMonSNW(I,J)+tmpscal2*ICE2SNOW
1456 HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal2*ICE2SNOW
1457 r_QbyATMmult_cover(I,J,K)=r_QbyATMmult_cover(I,J,K) - tmpscal2
1458 C keep the total up to date, too
1459 r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1460 ENDDO
1461 ENDDO
1462 ENDDO
1463 #else /* SEAICE_ITD */
1464 DO J=1,sNy
1465 DO I=1,sNx
1466 C Convert to standard units (meters of ice) rather than to meters
1467 C of snow. This appears to be more robust.
1468 tmpscal1=MAX(r_QbyATM_cover(I,J),-HSNOW(I,J,bi,bj)*SNOW2ICE)
1469 tmpscal2=MIN(tmpscal1,0. _d 0)
1470 #ifdef SEAICE_MODIFY_GROWTH_ADJ
1471 Cgf no additional dependency through snow
1472 IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1473 #endif
1474 d_HSNWbyATMonSNW(I,J)= tmpscal2*ICE2SNOW
1475 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal2*ICE2SNOW
1476 r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
1477 ENDDO
1478 ENDDO
1479 #endif /* SEAICE_ITD */
1480 c ToM<<< debug seaice_growth
1481 WRITE(msgBuf,'(A,7F6.2)')
1482 #ifdef SEAICE_ITD
1483 & ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',
1484 & HEFFITD(20,20,:,bi,bj)
1485 #else
1486 & ' SEAICE_GROWTH: Heff increments 3, HEFF = ',
1487 & HEFF(20,20,bi,bj)
1488 #endif
1489 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1490 & SQUEEZE_RIGHT , myThid)
1491 c ToM>>>
1492
1493 C compute ice thickness tendency due to the atmosphere
1494 C ====================================================
1495
1496 #ifdef ALLOW_AUTODIFF_TAMC
1497 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1498 CADJ STORE r_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1499 #endif /* ALLOW_AUTODIFF_TAMC */
1500
1501 Cgf note: this block is not actually tested by lab_sea
1502 Cgf where all experiments start in January. So even though
1503 Cgf the v1.81=>v1.82 revision would change results in
1504 Cgf warming conditions, the lab_sea results were not changed.
1505
1506 #ifdef SEAICE_ITD
1507 DO K=1,nITD
1508 DO J=1,sNy
1509 DO I=1,sNx
1510 #ifdef SEAICE_GROWTH_LEGACY
1511 tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K))
1512 #else
1513 tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)
1514 c Limit ice growth by potential melt by ocean
1515 & + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J))
1516 #endif /* SEAICE_GROWTH_LEGACY */
1517 d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
1518 & + tmpscal2
1519 d_HEFFbyATMonOCN(I,J) = d_HEFFbyATMonOCN(I,J)
1520 & + tmpscal2
1521 r_QbyATM_cover(I,J) = r_QbyATM_cover(I,J)
1522 & - tmpscal2
1523 HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal2
1524
1525 #ifdef ALLOW_SITRACER
1526 SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,3)
1527 & + HEFFITD(I,J,K,bi,bj)
1528 #endif
1529 ENDDO
1530 ENDDO
1531 ENDDO
1532 #else /* SEAICE_ITD */
1533 DO J=1,sNy
1534 DO I=1,sNx
1535
1536 #ifdef SEAICE_GROWTH_LEGACY
1537 tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))
1538 #else
1539 tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+
1540 c Limit ice growth by potential melt by ocean
1541 & AREApreTH(I,J) * r_QbyOCN(I,J))
1542 #endif /* SEAICE_GROWTH_LEGACY */
1543
1544 d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
1545 d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
1546 r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
1547 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2
1548
1549 #ifdef ALLOW_SITRACER
1550 SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
1551 #endif
1552 ENDDO
1553 ENDDO
1554 #endif /* SEAICE_ITD */
1555 c ToM<<< debug seaice_growth
1556 WRITE(msgBuf,'(A,7F6.2)')
1557 #ifdef SEAICE_ITD
1558 & ' SEAICE_GROWTH: Heff increments 4, HEFFITD = ',
1559 & HEFFITD(20,20,:,bi,bj)
1560 #else
1561 & ' SEAICE_GROWTH: Heff increments 4, HEFF = ',
1562 & HEFF(20,20,bi,bj)
1563 #endif
1564 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1565 & SQUEEZE_RIGHT , myThid)
1566 c ToM>>>
1567
1568 C attribute precip to fresh water or snow stock,
1569 C depending on atmospheric conditions.
1570 C =================================================
1571 #ifdef ALLOW_ATM_TEMP
1572 #ifdef ALLOW_AUTODIFF_TAMC
1573 CADJ STORE a_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1574 CADJ STORE PRECIP(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1575 CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
1576 #endif /* ALLOW_AUTODIFF_TAMC */
1577 DO J=1,sNy
1578 DO I=1,sNx
1579 C possible alternatives to the a_QbyATM_cover criterium
1580 c IF (TICE(I,J,bi,bj) .LT. TMIX) THEN
1581 c IF (atemp(I,J,bi,bj) .LT. celsius2K) THEN
1582 IF ( a_QbyATM_cover(I,J).GE. 0. _d 0 ) THEN
1583 C add precip as snow
1584 d_HFRWbyRAIN(I,J)=0. _d 0
1585 d_HSNWbyRAIN(I,J)=convertPRECIP2HI*ICE2SNOW*
1586 & PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1587 ELSE
1588 C add precip to the fresh water bucket
1589 d_HFRWbyRAIN(I,J)=-convertPRECIP2HI*
1590 & PRECIP(I,J,bi,bj)*AREApreTH(I,J)
1591 d_HSNWbyRAIN(I,J)=0. _d 0
1592 ENDIF
1593 ENDDO
1594 ENDDO
1595 #ifdef SEAICE_ITD
1596 DO K=1,nITD
1597 DO J=1,sNy
1598 DO I=1,sNx
1599 HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj)
1600 & + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,K)
1601 ENDDO
1602 ENDDO
1603 ENDDO
1604 #else
1605 DO J=1,sNy
1606 DO I=1,sNx
1607 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
1608 ENDDO
1609 ENDDO
1610 #endif
1611 Cgf note: this does not affect air-sea heat flux,
1612 Cgf since the implied air heat gain to turn
1613 Cgf rain to snow is not a surface process.
1614 #endif /* ALLOW_ATM_TEMP */
1615 c ToM<<< debug seaice_growth
1616 WRITE(msgBuf,'(A,7F6.2)')
1617 #ifdef SEAICE_ITD
1618 & ' SEAICE_GROWTH: Heff increments 5, HEFFITD = ',
1619 & HEFFITD(20,20,:,bi,bj)
1620 #else
1621 & ' SEAICE_GROWTH: Heff increments 5, HEFF = ',
1622 & HEFF(20,20,bi,bj)
1623 #endif
1624 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1625 & SQUEEZE_RIGHT , myThid)
1626 c ToM>>>
1627
1628 C compute snow melt due to heat available from ocean.
1629 C =================================================================
1630
1631 Cgf do we need to keep this comment and cpp bracket?
1632 Cph( very sensitive bit here by JZ
1633 #ifndef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING
1634 #ifdef ALLOW_AUTODIFF_TAMC
1635 CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1636 CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1637 #endif /* ALLOW_AUTODIFF_TAMC */
1638
1639 #ifdef SEAICE_ITD
1640 DO K=1,nITD
1641 DO J=1,sNy
1642 DO I=1,sNx
1643 tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,K),
1644 & -HSNOWITD(I,J,K,bi,bj))
1645 tmpscal2=MIN(tmpscal1,0. _d 0)
1646 #ifdef SEAICE_MODIFY_GROWTH_ADJ
1647 Cgf no additional dependency through snow
1648 if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1649 #endif
1650 d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
1651 r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE
1652 HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) + tmpscal2
1653 ENDDO
1654 ENDDO
1655 ENDDO
1656 #else /* SEAICE_ITD */
1657 DO J=1,sNy
1658 DO I=1,sNx
1659 tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))
1660 tmpscal2=MIN(tmpscal1,0. _d 0)
1661 #ifdef SEAICE_MODIFY_GROWTH_ADJ
1662 Cgf no additional dependency through snow
1663 if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
1664 #endif
1665 d_HSNWbyOCNonSNW(I,J) = tmpscal2
1666 r_QbyOCN(I,J)=r_QbyOCN(I,J)
1667 & -d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
1668 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_HSNWbyOCNonSNW(I,J)
1669 ENDDO
1670 ENDDO
1671 #endif /* SEAICE_ITD */
1672 #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
1673 Cph)
1674 c ToM<<< debug seaice_growth
1675 WRITE(msgBuf,'(A,7F6.2)')
1676 #ifdef SEAICE_ITD
1677 & ' SEAICE_GROWTH: Heff increments 6, HEFFITD = ',
1678 & HEFFITD(20,20,:,bi,bj)
1679 #else
1680 & ' SEAICE_GROWTH: Heff increments 6, HEFF = ',
1681 & HEFF(20,20,bi,bj)
1682 #endif
1683 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1684 & SQUEEZE_RIGHT , myThid)
1685 c ToM>>>
1686
1687 C gain of new ice over open water
1688 C ===============================
1689 #ifdef ALLOW_AUTODIFF_TAMC
1690 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1691 CADJ STORE r_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1692 CADJ STORE r_QbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1693 CADJ STORE a_QSWbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
1694 CADJ STORE a_QSWbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1695 #endif /* ALLOW_AUTODIFF_TAMC */
1696
1697 DO J=1,sNy
1698 DO I=1,sNx
1699 c Initial ice growth is triggered by open water
1700 c heat flux overcoming potential melt by ocean
1701 tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *
1702 & (1.0 _d 0 - AREApreTH(I,J))
1703 c Penetrative shortwave flux beyond first layer
1704 c that is therefore not available to ice growth/melt
1705 tmpscal2=SWFracB * a_QSWbyATM_open(I,J)
1706 C impose -HEFF as the maxmum melting if SEAICE_doOpenWaterMelt
1707 C or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)
1708 tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,
1709 & -HEFF(I,J,bi,bj)*facOpenMelt)*HEFFM(I,J,bi,bj)
1710 d_HEFFbyATMonOCN_open(I,J)=tmpscal3
1711 d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
1712 r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
1713 #ifdef SEAICE_ITD
1714 C open water area fraction
1715 tmpscal0 = ONE-AREApreTH(I,J)
1716 C determine thickness of new ice
1717 C considering the entire open water area to refreeze
1718 tmpscal1 = tmpscal3/tmpscal0
1719 C then add new ice volume to appropriate thickness category
1720 DO K=1,nITD
1721 IF (tmpscal1.LT.Hlimit(K)) THEN
1722 HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal3
1723 tmpscal3=ZERO
1724 C not sure if AREAITD should be changed here since AREA is incremented
1725 C in PART 4 below in non-itd code
1726 C in this scenario all open water is covered by new ice instantaneously,
1727 C i.e. no delayed lead closing is concidered such as is associated with
1728 C Hibler's h_0 parameter
1729 AREAITD(I,J,K,bi,bj) = AREAITD(I,J,K,bi,bj)
1730 & + tmpscal0
1731 tmpscal0=ZERO
1732 ENDIF
1733 ENDDO
1734 #else
1735 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
1736 #endif
1737 ENDDO
1738 ENDDO
1739
1740 #ifdef ALLOW_SITRACER
1741 #ifdef SEAICE_ITD
1742 DO K=1,nITD
1743 DO J=1,sNy
1744 DO I=1,sNx
1745 c needs to be here to allow use also with LEGACY branch
1746 SItrHEFF(I,J,bi,bj,4) = SItrHEFF(I,J,bi,bj,4)
1747 & + HEFFITD(I,J,K,bi,bj)
1748 ENDDO
1749 ENDDO
1750 ENDDO
1751 #else
1752 DO J=1,sNy
1753 DO I=1,sNx
1754 c needs to be here to allow use also with LEGACY branch
1755 SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
1756 ENDDO
1757 ENDDO
1758 #endif
1759 #endif /* ALLOW_SITRACER */
1760 c ToM<<< debug seaice_growth
1761 WRITE(msgBuf,'(A,7F6.2)')
1762 #ifdef SEAICE_ITD
1763 & ' SEAICE_GROWTH: Heff increments 7, HEFFITD = ',
1764 & HEFFITD(20,20,:,bi,bj)
1765 #else
1766 & ' SEAICE_GROWTH: Heff increments 7, HEFF = ',
1767 & HEFF(20,20,bi,bj)
1768 #endif
1769 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1770 & SQUEEZE_RIGHT , myThid)
1771 c ToM>>>
1772
1773 C convert snow to ice if submerged.
1774 C =================================
1775
1776 #ifndef SEAICE_GROWTH_LEGACY
1777 C note: in legacy, this process is done at the end
1778 #ifdef ALLOW_AUTODIFF_TAMC
1779 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1780 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1781 #endif /* ALLOW_AUTODIFF_TAMC */
1782 IF ( SEAICEuseFlooding ) THEN
1783 #ifdef SEAICE_ITD
1784 DO K=1,nITD
1785 DO J=1,sNy
1786 DO I=1,sNx
1787 tmpscal0 = (HSNOWITD(I,J,K,bi,bj)*SEAICE_rhoSnow
1788 & +HEFFITD(I,J,K,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
1789 tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,K,bi,bj))
1790 d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1
1791 HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal1
1792 HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal1
1793 & * ICE2SNOW
1794 ENDDO
1795 ENDDO
1796 ENDDO
1797 #else
1798 DO J=1,sNy
1799 DO I=1,sNx
1800 tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1801 & +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
1802 tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
1803 d_HEFFbyFLOODING(I,J)=tmpscal1
1804 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
1805 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
1806 & d_HEFFbyFLOODING(I,J)*ICE2SNOW
1807 ENDDO
1808 ENDDO
1809 #endif
1810 ENDIF
1811 #endif /* SEAICE_GROWTH_LEGACY */
1812 c ToM<<< debug seaice_growth
1813 WRITE(msgBuf,'(A,7F6.2)')
1814 #ifdef SEAICE_ITD
1815 & ' SEAICE_GROWTH: Heff increments 8, HEFFITD = ',
1816 & HEFFITD(20,20,:,bi,bj)
1817 #else
1818 & ' SEAICE_GROWTH: Heff increments 8, HEFF = ',
1819 & HEFF(20,20,bi,bj)
1820 #endif
1821 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
1822 & SQUEEZE_RIGHT , myThid)
1823 c ToM>>>
1824
1825 C ===================================================================
1826 C ==========PART 4: determine ice cover fraction increments=========-
1827 C ===================================================================
1828
1829 #ifdef ALLOW_AUTODIFF_TAMC
1830 CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
1831 CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=iicekey,byte=isbyte
1832 CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=iicekey,byte=isbyte
1833 CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
1834 CADJ STORE recip_heffActual = comlev1_bibj,key=iicekey,byte=isbyte
1835 CADJ STORE d_hsnwbyatmonsnw = comlev1_bibj,key=iicekey,byte=isbyte
1836 cph(
1837 cphCADJ STORE d_AREAbyATM = comlev1_bibj,key=iicekey,byte=isbyte
1838 cphCADJ STORE d_AREAbyICE = comlev1_bibj,key=iicekey,byte=isbyte
1839 cphCADJ STORE d_AREAbyOCN = comlev1_bibj,key=iicekey,byte=isbyte
1840 cph)
1841 CADJ STORE a_QbyATM_open = comlev1_bibj,key=iicekey,byte=isbyte
1842 CADJ STORE heffActual = comlev1_bibj,key=iicekey,byte=isbyte
1843 CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
1844 CADJ STORE HEFF(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1845 CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1846 CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
1847 #endif /* ALLOW_AUTODIFF_TAMC */
1848
1849 #ifdef SEAICE_ITD
1850 C replaces Hibler '79 scheme and lead closing parameter
1851 C because ITD accounts explicitly for lead openings and
1852 C different melt rates due to varying ice thickness
1853 C
1854 C only consider ice area loss due to total ice thickness loss
1855 C ice area gain due to freezing of open water as handled above
1856 C under "gain of new ice over open water"
1857 C
1858 C does not account for lateral melt of ice floes
1859 C
1860 C AREAITD is incremented in section "gain of new ice over open water" above
1861 C
1862 DO K=1,nITD
1863 DO J=1,sNy
1864 DO I=1,sNx
1865 IF (HEFFITD(I,J,K,bi,bj).LE.ZERO) THEN
1866 AREAITD(I,J,K,bi,bj)=ZERO
1867 ENDIF
1868 #ifdef ALLOW_SITRACER
1869 SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)
1870 & + AREAITD(I,J,K,bi,bj)
1871 #endif /* ALLOW_SITRACER */
1872 ENDDO
1873 ENDDO
1874 ENDDO
1875 #else /* SEAICE_ITD */
1876 DO J=1,sNy
1877 DO I=1,sNx
1878
1879 IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
1880 recip_HO=1. _d 0 / HO_south
1881 ELSE
1882 recip_HO=1. _d 0 / HO
1883 ENDIF
1884 #ifdef SEAICE_GROWTH_LEGACY
1885 tmpscal0=HEFF(I,J,bi,bj) - d_HEFFbyATMonOCN(I,J)
1886 recip_HH = AREApreTH(I,J) /(tmpscal0+.00001 _d 0)
1887 #else
1888 recip_HH = recip_heffActual(I,J)
1889 #endif
1890
1891 C gain of ice over open water : computed from
1892 C (SEAICE_areaGainFormula.EQ.1) from growth by ATM
1893 C (SEAICE_areaGainFormula.EQ.2) from predicted growth by ATM
1894 IF (SEAICE_areaGainFormula.EQ.1) THEN
1895 tmpscal4 = MAX(ZERO,d_HEFFbyATMonOCN_open(I,J))
1896 ELSE
1897 tmpscal4=MAX(ZERO,a_QbyATM_open(I,J))
1898 ENDIF
1899
1900 C loss of ice cover by melting : computed from
1901 C (SEAICE_areaLossFormula.EQ.1) from all but only melt conributions by ATM and OCN
1902 C (SEAICE_areaLossFormula.EQ.2) from net melt-growth>0 by ATM and OCN
1903 C (SEAICE_areaLossFormula.EQ.3) from predicted melt by ATM
1904 IF (SEAICE_areaLossFormula.EQ.1) THEN
1905 tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN_cover(I,J) )
1906 & + MIN( 0. _d 0 , d_HEFFbyATMonOCN_open(I,J) )
1907 & + MIN( 0. _d 0 , d_HEFFbyOCNonICE(I,J) )
1908 ELSEIF (SEAICE_areaLossFormula.EQ.2) THEN
1909 tmpscal3 = MIN( 0. _d 0 , d_HEFFbyATMonOCN_cover(I,J)
1910 & + d_HEFFbyATMonOCN_open(I,J) + d_HEFFbyOCNonICE(I,J) )
1911 ELSE
1912 C compute heff after ice melt by ocn:
1913 tmpscal0=HEFF(I,J,bi,bj) - d_HEFFbyATMonOCN(I,J)
1914 C compute available heat left after snow melt by atm:
1915 tmpscal1= a_QbyATM_open(I,J)+a_QbyATM_cover(I,J)
1916 & - d_HSNWbyATMonSNW(I,J)*SNOW2ICE
1917 C could not melt more than all the ice
1918 tmpscal2 = MAX(-tmpscal0,tmpscal1)
1919 tmpscal3 = MIN(ZERO,tmpscal2)
1920 ENDIF
1921
1922 C apply tendency
1923 IF ( (HEFF(i,j,bi,bj).GT.0. _d 0).OR.
1924 & (HSNOW(i,j,bi,bj).GT.0. _d 0) ) THEN
1925 AREA(I,J,bi,bj)=MAX(0. _d 0,
1926 & MIN( SEAICE_area_max, AREA(I,J,bi,bj)
1927 & + recip_HO*tmpscal4+HALF*recip_HH*tmpscal3 ))
1928 ELSE
1929 AREA(I,J,bi,bj)=0. _d 0
1930 ENDIF
1931 #ifdef ALLOW_SITRACER
1932 SItrAREA(I,J,bi,bj,3)=AREA(I,J,bi,bj)
1933 #endif /* ALLOW_SITRACER */
1934 #ifdef ALLOW_DIAGNOSTICS
1935 d_AREAbyATM(I,J)=
1936 & recip_HO*MAX(ZERO,d_HEFFbyATMonOCN_open(I,J))
1937 & +HALF*recip_HH*MIN(0. _d 0,d_HEFFbyATMonOCN_open(I,J))
1938 d_AREAbyICE(I,J)=
1939 & HALF*recip_HH*MIN(0. _d 0,d_HEFFbyATMonOCN_cover(I,J))
1940 d_AREAbyOCN(I,J)=
1941 & HALF*recip_HH*MIN( 0. _d 0,d_HEFFbyOCNonICE(I,J) )
1942 #endif /* ALLOW_DIAGNOSTICS */
1943 ENDDO
1944 ENDDO
1945 #endif /* SEAICE_ITD */
1946
1947 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
1948 Cgf 'bulk' linearization of area=f(HEFF)
1949 IF ( SEAICEadjMODE.GE.1 ) THEN
1950 #ifdef SEAICE_ITD
1951 DO K=1,nITD
1952 DO J=1,sNy
1953 DO I=1,sNx
1954 AREAITD(I,J,K,bi,bj) = AREAITDpreTH(I,J,K) + 0.1 _d 0 *
1955 & ( HEFFITD(I,J,K,bi,bj) - HEFFITDpreTH(I,J,K) )
1956 ENDDO
1957 ENDDO
1958 ENDDO
1959 #else
1960 DO J=1,sNy
1961 DO I=1,sNx
1962 C AREA(I,J,bi,bj) = 0.1 _d 0 * HEFF(I,J,bi,bj)
1963 AREA(I,J,bi,bj) = AREApreTH(I,J) + 0.1 _d 0 *
1964 & ( HEFF(I,J,bi,bj) - HEFFpreTH(I,J) )
1965 ENDDO
1966 ENDDO
1967 #endif
1968 ENDIF
1969 #endif
1970 #ifdef SEAICE_ITD
1971 C check categories for consistency with limits after growth/melt
1972 CALL SEAICE_ITD_REDIST(bi, bj, myTime,myIter,myThid)
1973 C finally update total AREA, HEFF, HSNOW
1974 CALL SEAICE_ITD_SUM(bi, bj, myTime,myIter,myThid)
1975 #endif
1976
1977 C ===================================================================
1978 C =============PART 5: determine ice salinity increments=============
1979 C ===================================================================
1980
1981 #ifndef SEAICE_VARIABLE_SALINITY
1982 # if (defined ALLOW_AUTODIFF_TAMC && defined ALLOW_SALT_PLUME)
1983 CADJ STORE d_HEFFbyNEG = comlev1_bibj,key=iicekey,byte=isbyte
1984 CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
1985 CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
1986 CADJ STORE d_HEFFbyATMonOCN_open = comlev1_bibj,key=iicekey,byte=isbyte
1987 CADJ STORE d_HEFFbyATMonOCN_cover = comlev1_bibj,key=iicekey,byte=isbyte
1988 CADJ STORE d_HEFFbyFLOODING = comlev1_bibj,key=iicekey,byte=isbyte
1989 CADJ STORE d_HEFFbySublim = comlev1_bibj,key=iicekey,byte=isbyte
1990 CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
1991 CADJ & key = iicekey, byte = isbyte
1992 # endif /* ALLOW_AUTODIFF_TAMC and ALLOW_SALT_PLUME */
1993 DO J=1,sNy
1994 DO I=1,sNx
1995 tmpscal1 = d_HEFFbyNEG(I,J) + d_HEFFbyOCNonICE(I,J) +
1996 & d_HEFFbyATMonOCN(I,J) + d_HEFFbyFLOODING(I,J)
1997 & + d_HEFFbySublim(I,J)
1998 #ifdef SEAICE_ALLOW_AREA_RELAXATION
1999 + d_HEFFbyRLX(I,J)
2000 #endif
2001 tmpscal2 = tmpscal1 * SEAICE_salt0 * HEFFM(I,J,bi,bj)
2002 & * recip_deltaTtherm * SEAICE_rhoIce
2003 saltFlux(I,J,bi,bj) = tmpscal2
2004 #ifdef ALLOW_SALT_PLUME
2005 tmpscal3 = tmpscal1*salt(I,J,kSurface,bi,bj)*HEFFM(I,J,bi,bj)
2006 & * recip_deltaTtherm * SEAICE_rhoIce
2007 saltPlumeFlux(I,J,bi,bj) = MAX( tmpscal3-tmpscal2 , 0. _d 0)
2008 & *SPsalFRAC
2009 #endif /* ALLOW_SALT_PLUME */
2010 ENDDO
2011 ENDDO
2012 #endif /* ndef SEAICE_VARIABLE_SALINITY */
2013
2014 #ifdef SEAICE_VARIABLE_SALINITY
2015
2016 #ifdef SEAICE_GROWTH_LEGACY
2017 # ifdef ALLOW_AUTODIFF_TAMC
2018 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
2019 # endif /* ALLOW_AUTODIFF_TAMC */
2020 DO J=1,sNy
2021 DO I=1,sNx
2022 C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean
2023 IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN
2024 saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
2025 & HSALT(I,J,bi,bj) * recip_deltaTtherm
2026 HSALT(I,J,bi,bj) = 0.0 _d 0
2027 ENDIF
2028 ENDDO
2029 ENDDO
2030 #endif /* SEAICE_GROWTH_LEGACY */
2031
2032 #ifdef ALLOW_AUTODIFF_TAMC
2033 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
2034 #endif /* ALLOW_AUTODIFF_TAMC */
2035
2036 DO J=1,sNy
2037 DO I=1,sNx
2038 C sum up the terms that affect the salt content of the ice pack
2039 tmpscal1=d_HEFFbyOCNonICE(I,J)+d_HEFFbyATMonOCN(I,J)
2040
2041 C recompute HEFF before thermodynamic updates (which is not AREApreTH in legacy code)
2042 tmpscal2=HEFF(I,J,bi,bj)-tmpscal1-d_HEFFbyFLOODING(I,J)
2043 C tmpscal1 > 0 : m of sea ice that is created
2044 IF ( tmpscal1 .GE. 0.0 ) THEN
2045 saltFlux(I,J,bi,bj) =
2046 & HEFFM(I,J,bi,bj)*recip_deltaTtherm
2047 & *SEAICE_saltFrac*salt(I,J,kSurface,bi,bj)
2048 & *tmpscal1*SEAICE_rhoIce
2049 #ifdef ALLOW_SALT_PLUME
2050 C saltPlumeFlux is defined only during freezing:
2051 saltPlumeFlux(I,J,bi,bj)=
2052 & HEFFM(I,J,bi,bj)*recip_deltaTtherm
2053 & *(ONE-SEAICE_saltFrac)*salt(I,J,kSurface,bi,bj)
2054 & *tmpscal1*SEAICE_rhoIce
2055 & *SPsalFRAC
2056 C if SaltPlumeSouthernOcean=.FALSE. turn off salt plume in Southern Ocean
2057 IF ( .NOT. SaltPlumeSouthernOcean ) THEN
2058 IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 )
2059 & saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
2060 ENDIF
2061 #endif /* ALLOW_SALT_PLUME */
2062
2063 C tmpscal1 < 0 : m of sea ice that is melted
2064 ELSE
2065 saltFlux(I,J,bi,bj) =
2066 & HEFFM(I,J,bi,bj)*recip_deltaTtherm
2067 & *HSALT(I,J,bi,bj)
2068 & *tmpscal1/tmpscal2
2069 #ifdef ALLOW_SALT_PLUME
2070 saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
2071 #endif /* ALLOW_SALT_PLUME */
2072 ENDIF
2073 C update HSALT based on surface saltFlux
2074 HSALT(I,J,bi,bj) = HSALT(I,J,bi,bj) +
2075 & saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
2076 saltFlux(I,J,bi,bj) =
2077 & saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
2078 #ifdef SEAICE_GROWTH_LEGACY
2079 C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean
2080 IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN
2081 saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -
2082 & HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) * recip_deltaTtherm
2083 HSALT(I,J,bi,bj) = 0.0 _d 0
2084 #ifdef ALLOW_SALT_PLUME
2085 saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
2086 #endif /* ALLOW_SALT_PLUME */
2087 ENDIF
2088 #endif /* SEAICE_GROWTH_LEGACY */
2089 ENDDO
2090 ENDDO
2091 #endif /* SEAICE_VARIABLE_SALINITY */
2092
2093
2094 C =======================================================================
2095 C ==LEGACY PART 6 (LEGACY) treat pathological cases, then do flooding ===
2096 C =======================================================================
2097
2098 #ifdef SEAICE_GROWTH_LEGACY
2099
2100 C treat values of ice cover fraction oustide
2101 C the [0 1] range, and other such issues.
2102 C ===========================================
2103
2104 Cgf note: this part cannot be heat and water conserving
2105
2106 #ifdef ALLOW_AUTODIFF_TAMC
2107 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
2108 CADJ & key = iicekey, byte = isbyte
2109 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
2110 CADJ & key = iicekey, byte = isbyte
2111 #endif /* ALLOW_AUTODIFF_TAMC */
2112 DO J=1,sNy
2113 DO I=1,sNx
2114 C NOW SET AREA(I,J,bi,bj)=0 WHERE THERE IS NO ICE
2115 CML replaced "/.0001 _d 0" by "*1. _d 4", 1e-4 is probably
2116 CML meant to be something like a minimum thickness
2117 AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),HEFF(I,J,bi,bj)*1. _d 4)
2118 ENDDO
2119 ENDDO
2120
2121 #ifdef ALLOW_AUTODIFF_TAMC
2122 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
2123 CADJ & key = iicekey, byte = isbyte
2124 #endif /* ALLOW_AUTODIFF_TAMC */
2125 DO J=1,sNy
2126 DO I=1,sNx
2127 C NOW TRUNCATE AREA
2128 AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj))
2129 ENDDO
2130 ENDDO
2131
2132 #ifdef ALLOW_AUTODIFF_TAMC
2133 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
2134 CADJ & key = iicekey, byte = isbyte
2135 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
2136 CADJ & key = iicekey, byte = isbyte
2137 #endif /* ALLOW_AUTODIFF_TAMC */
2138 DO J=1,sNy
2139 DO I=1,sNx
2140 AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj))
2141 HSNOW(I,J,bi,bj) = MAX(ZERO,HSNOW(I,J,bi,bj))
2142 AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
2143 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)
2144 #ifdef SEAICE_CAP_HEFF
2145 C This is not energy conserving, but at least it conserves fresh water
2146 tmpscal0 = -MAX(HEFF(I,J,bi,bj)-MAX_HEFF,0. _d 0)
2147 d_HEFFbyNeg(I,J) = d_HEFFbyNeg(I,J) + tmpscal0
2148 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal0
2149 #endif /* SEAICE_CAP_HEFF */
2150 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
2151 ENDDO
2152 ENDDO
2153
2154 C convert snow to ice if submerged.
2155 C =================================
2156
2157 IF ( SEAICEuseFlooding ) THEN
2158 DO J=1,sNy
2159 DO I=1,sNx
2160 tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
2161 & +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
2162 tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
2163 d_HEFFbyFLOODING(I,J)=tmpscal1
2164 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
2165 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
2166 & d_HEFFbyFLOODING(I,J)*ICE2SNOW
2167 ENDDO
2168 ENDDO
2169 ENDIF
2170
2171 #endif /* SEAICE_GROWTH_LEGACY */
2172
2173 #ifdef ALLOW_SITRACER
2174 DO J=1,sNy
2175 DO I=1,sNx
2176 c needs to be here to allow use also with LEGACY branch
2177 SItrHEFF(I,J,bi,bj,5)=HEFF(I,J,bi,bj)
2178 ENDDO
2179 ENDDO
2180 #endif /* ALLOW_SITRACER */
2181
2182 C ===================================================================
2183 C ==============PART 7: determine ocean model forcing================
2184 C ===================================================================
2185
2186 C compute net heat flux leaving/entering the ocean,
2187 C accounting for the part used in melt/freeze processes
2188 C =====================================================
2189
2190 #ifdef ALLOW_AUTODIFF_TAMC
2191 CADJ STORE d_hsnwbyneg = comlev1_bibj,key=iicekey,byte=isbyte
2192 CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=iicekey,byte=isbyte
2193 #endif /* ALLOW_AUTODIFF_TAMC */
2194
2195 DO J=1,sNy
2196 DO I=1,sNx
2197 QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + r_QbyATM_open(I,J)
2198 #ifndef SEAICE_GROWTH_LEGACY
2199 C in principle a_QSWbyATM_cover should always be included here, however
2200 C for backward compatibility it is left out of the LEGACY branch
2201 & + a_QSWbyATM_cover(I,J)
2202 #endif /* SEAICE_GROWTH_LEGACY */
2203 & - ( d_HEFFbyOCNonICE(I,J) +
2204 & d_HSNWbyOCNonSNW(I,J)*SNOW2ICE +
2205 & d_HEFFbyNEG(I,J) +
2206 #ifdef SEAICE_ALLOW_AREA_RELAXATION
2207 & d_HEFFbyRLX(I,J) +
2208 #endif
2209 & d_HSNWbyNEG(I,J)*SNOW2ICE )
2210 & * maskC(I,J,kSurface,bi,bj)
2211 QSW(I,J,bi,bj) = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)
2212 ENDDO
2213 ENDDO
2214
2215 C switch heat fluxes from 'effective' ice meters to W/m2
2216 C ======================================================
2217
2218 DO J=1,sNy
2219 DO I=1,sNx
2220 QNET(I,J,bi,bj) = QNET(I,J,bi,bj)*convertHI2Q
2221 QSW(I,J,bi,bj) = QSW(I,J,bi,bj)*convertHI2Q
2222 ENDDO
2223 ENDDO
2224
2225 #ifndef SEAICE_DISABLE_HEATCONSFIX
2226 C treat advective heat flux by ocean to ice water exchange (at 0decC)
2227 C ===================================================================
2228 # ifdef ALLOW_AUTODIFF_TAMC
2229 CADJ STORE d_HEFFbyNEG = comlev1_bibj,key=iicekey,byte=isbyte
2230 CADJ STORE d_HEFFbyOCNonICE = comlev1_bibj,key=iicekey,byte=isbyte
2231 CADJ STORE d_HEFFbyATMonOCN = comlev1_bibj,key=iicekey,byte=isbyte
2232 CADJ STORE d_HSNWbyNEG = comlev1_bibj,key=iicekey,byte=isbyte
2233 CADJ STORE d_HSNWbyOCNonSNW = comlev1_bibj,key=iicekey,byte=isbyte
2234 CADJ STORE d_HSNWbyATMonSNW = comlev1_bibj,key=iicekey,byte=isbyte
2235 CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
2236 CADJ & key = iicekey, byte = isbyte
2237 # endif /* ALLOW_AUTODIFF_TAMC */
2238 IF ( SEAICEheatConsFix ) THEN
2239 c Unlike for evap and precip, the temperature of gained/lost
2240 c ocean liquid water due to melt/freeze of solid water cannot be chosen
2241 c to be e.g. the ocean SST. It must be done at 0degC. The fix below anticipates
2242 c on external_forcing_surf.F and applies the correction to QNET.
2243 IF ((convertFW2Salt.EQ.-1.).OR.(temp_EvPrRn.EQ.UNSET_RL)) THEN
2244 c I leave alone the exotic case when onvertFW2Salt.NE.-1 and temp_EvPrRn.NE.UNSET_RL and
2245 c the small error of the synchronous time stepping case (see external_forcing_surf.F).
2246 DO J=1,sNy
2247 DO I=1,sNx
2248 #ifdef ALLOW_DIAGNOSTICS
2249 c store unaltered QNET for diagnostic purposes
2250 DIAGarrayA(I,J)=QNET(I,J,bi,bj)
2251 #endif
2252 c compute the ocean water going to ice/snow, in precip units
2253 tmpscal3=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
2254 & ( d_HSNWbyATMonSNW(I,J)*SNOW2ICE
2255 & + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
2256 & + d_HEFFbyOCNonICE(I,J) + d_HEFFbyATMonOCN(I,J)
2257 & + d_HEFFbyNEG(I,J) + d_HSNWbyNEG(I,J)*SNOW2ICE )
2258 & * convertHI2PRECIP
2259 c factor in the heat content that external_forcing_surf.F
2260 c will associate with EMPMR, and remove it from QNET, so that
2261 c melt/freez water is in effect consistently gained/lost at 0degC
2262 IF (temp_EvPrRn.NE.UNSET_RL) THEN
2263 QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*
2264 & HeatCapacity_Cp * temp_EvPrRn
2265 ELSE
2266 QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*
2267 & HeatCapacity_Cp * theta(I,J,kSurface,bi,bj)
2268 ENDIF
2269 #ifdef ALLOW_DIAGNOSTICS
2270 c back out the eventual TFLUX adjustement and fill diag
2271 DIAGarrayA(I,J)=QNET(I,J,bi,bj)-DIAGarrayA(I,J)
2272 #endif
2273 ENDDO
2274 ENDDO
2275 ENDIF
2276 #ifdef ALLOW_DIAGNOSTICS
2277 CALL DIAGNOSTICS_FILL(DIAGarrayA,
2278 & 'SIaaflux',0,1,3,bi,bj,myThid)
2279 #endif
2280 ENDIF
2281 #endif /* ndef SEAICE_DISABLE_HEATCONSFIX */
2282
2283 C compute net fresh water flux leaving/entering
2284 C the ocean, accounting for fresh/salt water stocks.
2285 C ==================================================
2286
2287 #ifdef ALLOW_ATM_TEMP
2288 DO J=1,sNy
2289 DO I=1,sNx
2290 tmpscal1= d_HSNWbyATMonSNW(I,J)*SNOW2ICE
2291 & +d_HFRWbyRAIN(I,J)
2292 & +d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
2293 & +d_HEFFbyOCNonICE(I,J)
2294 & +d_HEFFbyATMonOCN(I,J)
2295 & +d_HEFFbyNEG(I,J)
2296 #ifdef SEAICE_ALLOW_AREA_RELAXATION
2297 & +d_HEFFbyRLX(I,J)
2298 #endif
2299 & +d_HSNWbyNEG(I,J)*SNOW2ICE
2300 C If r_FWbySublim>0, then it is evaporated from ocean.
2301 & +r_FWbySublim(I,J)
2302 EmPmR(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
2303 & ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
2304 & * ( ONE - AREApreTH(I,J) )
2305 #ifdef ALLOW_RUNOFF
2306 & - RUNOFF(I,J,bi,bj)
2307 #endif /* ALLOW_RUNOFF */
2308 & + tmpscal1*convertHI2PRECIP
2309 & )*rhoConstFresh
2310 ENDDO
2311 ENDDO
2312
2313 #ifdef ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION
2314 C--
2315 DO J=1,sNy
2316 DO I=1,sNx
2317 frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
2318 & PRECIP(I,J,bi,bj)
2319 & - EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
2320 # ifdef ALLOW_RUNOFF
2321 & + RUNOFF(I,J,bi,bj)
2322 # endif /* ALLOW_RUNOFF */
2323 & )*rhoConstFresh
2324 # ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
2325 & - a_FWbySublim(I,J)*AREApreTH(I,J)
2326 # endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
2327 ENDDO
2328 ENDDO
2329 C--
2330 #else /* ndef ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION */
2331 C--
2332 # ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
2333 DO J=1,sNy
2334 DO I=1,sNx
2335 frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
2336 & PRECIP(I,J,bi,bj)
2337 & - EVAP(I,J,bi,bj)
2338 & *( ONE - AREApreTH(I,J) )
2339 # ifdef ALLOW_RUNOFF
2340 & + RUNOFF(I,J,bi,bj)
2341 # endif /* ALLOW_RUNOFF */
2342 & )*rhoConstFresh
2343 & - a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
2344 ENDDO
2345 ENDDO
2346 # endif
2347 C--
2348 #endif /* ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION */
2349
2350 #endif /* ALLOW_ATM_TEMP */
2351
2352 #ifdef SEAICE_DEBUG
2353 CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
2354 CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
2355 CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
2356 #endif /* SEAICE_DEBUG */
2357
2358 C Sea Ice Load on the sea surface.
2359 C =================================
2360
2361 #ifdef ALLOW_AUTODIFF_TAMC
2362 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
2363 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
2364 #endif /* ALLOW_AUTODIFF_TAMC */
2365
2366 IF ( useRealFreshWaterFlux ) THEN
2367 DO J=1,sNy
2368 DO I=1,sNx
2369 #ifdef SEAICE_CAP_ICELOAD
2370 tmpscal1 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
2371 & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
2372 tmpscal2 = MIN(tmpscal1,heffTooHeavy*rhoConst)
2373 #else
2374 tmpscal2 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
2375 & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
2376 #endif
2377 sIceLoad(i,j,bi,bj) = tmpscal2
2378 ENDDO
2379 ENDDO
2380 ENDIF
2381
2382 C ===================================================================
2383 C ======================PART 8: diagnostics==========================
2384 C ===================================================================
2385
2386 #ifdef ALLOW_DIAGNOSTICS
2387 IF ( useDiagnostics ) THEN
2388 tmpscal1=1. _d 0 * recip_deltaTtherm
2389 CALL DIAGNOSTICS_SCALE_FILL(a_QbyATM_cover,
2390 & tmpscal1,1,'SIaQbATC',0,1,3,bi,bj,myThid)
2391 CALL DIAGNOSTICS_SCALE_FILL(a_QbyATM_open,
2392 & tmpscal1,1,'SIaQbATO',0,1,3,bi,bj,myThid)
2393 CALL DIAGNOSTICS_SCALE_FILL(a_QbyOCN,
2394 & tmpscal1,1,'SIaQbOCN',0,1,3,bi,bj,myThid)
2395 CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyOCNonICE,
2396 & tmpscal1,1,'SIdHbOCN',0,1,3,bi,bj,myThid)
2397 CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyATMonOCN_cover,
2398 & tmpscal1,1,'SIdHbATC',0,1,3,bi,bj,myThid)
2399 CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyATMonOCN_open,
2400 & tmpscal1,1,'SIdHbATO',0,1,3,bi,bj,myThid)
2401 CALL DIAGNOSTICS_SCALE_FILL(d_HEFFbyFLOODING,
2402 & tmpscal1,1,'SIdHbFLO',0,1,3,bi,bj,myThid)
2403 CALL DIAGNOSTICS_SCALE_FILL(d_HSNWbyOCNonSNW,
2404 & tmpscal1,1,'SIdSbOCN',0,1,3,bi,bj,myThid)
2405 CALL DIAGNOSTICS_SCALE_FILL(d_HSNWbyATMonSNW,
2406 & tmpscal1,1,'SIdSbATC',0,1,3,bi,bj,myThid)
2407 CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyATM,
2408 & tmpscal1,1,'SIdAbATO',0,1,3,bi,bj,myThid)
2409 CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyICE,
2410 & tmpscal1,1,'SIdAbATC',0,1,3,bi,bj,myThid)
2411 CALL DIAGNOSTICS_SCALE_FILL(d_AREAbyOCN,
2412 & tmpscal1,1,'SIdAbOCN',0,1,3,bi,bj,myThid)
2413 CALL DIAGNOSTICS_SCALE_FILL(r_QbyATM_open,
2414 & convertHI2Q,1, 'SIqneto ',0,1,3,bi,bj,myThid)
2415 CALL DIAGNOSTICS_SCALE_FILL(r_QbyATM_cover,
2416 & convertHI2Q,1, 'SIqneti ',0,1,3,bi,bj,myThid)
2417 C three that actually need intermediate storage
2418 DO J=1,sNy
2419 DO I=1,sNx
2420 DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
2421 & * d_HSNWbyRAIN(I,J)*SEAICE_rhoSnow*recip_deltaTtherm
2422 DIAGarrayB(I,J) = AREA(I,J,bi,bj)-AREApreTH(I,J)
2423 ENDDO
2424 ENDDO
2425 CALL DIAGNOSTICS_FILL(DIAGarrayA,
2426 & 'SIsnPrcp',0,1,3,bi,bj,myThid)
2427 CALL DIAGNOSTICS_SCALE_FILL(DIAGarrayB,
2428 & tmpscal1,1,'SIdA ',0,1,3,bi,bj,myThid)
2429 #ifdef ALLOW_ATM_TEMP
2430 DO J=1,sNy
2431 DO I=1,sNx
2432 CML If I consider the atmosphere above the ice, the surface flux
2433 CML which is relevant for the air temperature dT/dt Eq
2434 CML accounts for sensible and radiation (with different treatment
2435 CML according to wave-length) fluxes but not for "latent heat flux",
2436 CML since it does not contribute to heating the air.
2437 CML So this diagnostic is only good for heat budget calculations within
2438 CML the ice-ocean system.
2439 DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
2440 #ifndef SEAICE_GROWTH_LEGACY
2441 & a_QSWbyATM_cover(I,J) +
2442 #endif /* SEAICE_GROWTH_LEGACY */
2443 & a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
2444 C
2445 DIAGarrayB(I,J) = maskC(I,J,kSurface,bi,bj) *
2446 & a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
2447 C
2448 DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)*(
2449 & PRECIP(I,J,bi,bj)
2450 & - EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
2451 #ifdef ALLOW_RUNOFF
2452 & + RUNOFF(I,J,bi,bj)
2453 #endif /* ALLOW_RUNOFF */
2454 & )*rhoConstFresh
2455 & - a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
2456 ENDDO
2457 ENDDO
2458 CALL DIAGNOSTICS_FILL(DIAGarrayA,
2459 & 'SIatmQnt',0,1,3,bi,bj,myThid)
2460 CALL DIAGNOSTICS_FILL(DIAGarrayB,
2461 & 'SIfwSubl',0,1,3,bi,bj,myThid)
2462 CALL DIAGNOSTICS_FILL(DIAGarrayC,
2463 & 'SIatmFW ',0,1,3,bi,bj,myThid)
2464 C
2465 DO J=1,sNy
2466 DO I=1,sNx
2467 C the actual Freshwater flux of sublimated ice, >0 decreases ice
2468 DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
2469 & * (a_FWbySublim(I,J)-r_FWbySublim(I,J))
2470 & * SEAICE_rhoIce * recip_deltaTtherm
2471 c the residual Freshwater flux of sublimated ice
2472 DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)
2473 & * r_FWbySublim(I,J)
2474 & * SEAICE_rhoIce * recip_deltaTtherm
2475 C the latent heat flux
2476 tmpscal1= EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
2477 & + r_FWbySublim(I,J)*convertHI2PRECIP
2478 tmpscal2= ( a_FWbySublim(I,J)-r_FWbySublim(I,J) )
2479 & * convertHI2PRECIP
2480 tmpscal3= SEAICE_lhEvap+SEAICE_lhFusion
2481 DIAGarrayB(I,J) = -maskC(I,J,kSurface,bi,bj)*rhoConstFresh
2482 & * ( tmpscal1*SEAICE_lhEvap + tmpscal2*tmpscal3 )
2483 ENDDO
2484 ENDDO
2485 CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIacSubl',0,1,3,bi,bj,myThid)
2486 CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIrsSubl',0,1,3,bi,bj,myThid)
2487 CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIhl ',0,1,3,bi,bj,myThid)
2488
2489 DO J=1,sNy
2490 DO I=1,sNx
2491 c compute ice/snow water going to atm, in precip units
2492 tmpscal1 = rhoConstFresh*maskC(I,J,kSurface,bi,bj)
2493 & * convertHI2PRECIP * ( - d_HSNWbyRAIN(I,J)*SNOW2ICE
2494 & + a_FWbySublim(I,J) - r_FWbySublim(I,J) )
2495 c compute ocean water going to atm, in precip units
2496 tmpscal2=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
2497 & ( ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
2498 & * ( ONE - AREApreTH(I,J) )
2499 #ifdef ALLOW_RUNOFF
2500 & - RUNOFF(I,J,bi,bj)
2501 #endif /* ALLOW_RUNOFF */
2502 & + ( d_HFRWbyRAIN(I,J) + r_FWbySublim(I,J) )
2503 & *convertHI2PRECIP )
2504 c factor in the advected specific energy (referenced to 0 for 0deC liquid water)
2505 tmpscal1= - tmpscal1*
2506 & ( -SEAICE_lhFusion + HeatCapacity_Cp * ZERO )
2507 IF (temp_EvPrRn.NE.UNSET_RL) THEN
2508 tmpscal2= - tmpscal2*
2509 & ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
2510 ELSE
2511 tmpscal2= - tmpscal2*
2512 & ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )
2513 ENDIF
2514 c add to SIatmQnt, leading to SItflux, which is analogous to TFLUX
2515 DIAGarrayA(I,J)=maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
2516 #ifndef SEAICE_GROWTH_LEGACY
2517 & a_QSWbyATM_cover(I,J) +
2518 #endif
2519 & a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
2520 & -tmpscal1-tmpscal2
2521 ENDDO
2522 ENDDO
2523 CALL DIAGNOSTICS_FILL(DIAGarrayA,
2524 & 'SItflux ',0,1,3,bi,bj,myThid)
2525 #endif /* ALLOW_ATM_TEMP */
2526
2527 ENDIF
2528 #endif /* ALLOW_DIAGNOSTICS */
2529
2530 C close bi,bj loops
2531 ENDDO
2532 ENDDO
2533
2534 RETURN
2535 END
  ViewVC Help
Powered by ViewVC 1.1.22