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

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

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


Revision 1.96 - (hide annotations) (download)
Sun Oct 31 20:07:01 2010 UTC (13 years, 7 months ago) by gforget
Branch: MAIN
Changes since 1.95: +22 -38 lines
- move areaMin, hiceMin, areaMax to common blocks.
- sort out diagnostics in seaice_growth.F.
- practical approximations in seaice_growth.F adjoint.
- reset TICE to celsius2K when no ice is present.

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

  ViewVC Help
Powered by ViewVC 1.1.22