/[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.101 - (hide annotations) (download)
Fri Nov 19 16:21:08 2010 UTC (13 years, 6 months ago) by mlosch
Branch: MAIN
Changes since 1.100: +4 -6 lines
  - replace irritating parameters SEAICE_latentWater/Ice and SEAICE_sensHeat
    by something more sensible (parameters that are what their name implies)
  - change some defaults, so that by default exf-parameters are used for
    things like rhoAir, cpAir,latent/sensible heat parameters

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

  ViewVC Help
Powered by ViewVC 1.1.22