/[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.55 - (hide annotations) (download)
Wed Mar 18 14:35:37 2009 UTC (15 years, 2 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint61n, checkpoint61l, checkpoint61m
Changes since 1.54: +3 -3 lines 
forgot a few instances of 330., when replace this number by a runtime
parameter SEAICE_rhoSnow. There is still one instance in
cost_ice_test.F, which I am not touching
1 mlosch 1.55 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_growth.F,v 1.54 2009/03/18 14:02:25 mlosch Exp $
2 heimbach 1.2 C $Name: $
3 mlosch 1.1
4     #include "SEAICE_OPTIONS.h"
5    
6     CStartOfInterface
7     SUBROUTINE SEAICE_GROWTH( myTime, myIter, myThid )
8     C /==========================================================\
9     C | SUBROUTINE seaice_growth |
10     C | o Updata ice thickness and snow depth |
11     C |==========================================================|
12     C \==========================================================/
13     IMPLICIT NONE
14    
15     C === Global variables ===
16     #include "SIZE.h"
17     #include "EEPARAMS.h"
18     #include "PARAMS.h"
19     #include "DYNVARS.h"
20     #include "GRID.h"
21     #include "FFIELDS.h"
22     #include "SEAICE_PARAMS.h"
23     #include "SEAICE.h"
24 dimitri 1.38 #ifdef ALLOW_EXF
25     # include "EXF_OPTIONS.h"
26     # include "EXF_FIELDS.h"
27 dimitri 1.40 # include "EXF_PARAM.h"
28 dimitri 1.38 #endif
29 dimitri 1.37 #ifdef ALLOW_SALT_PLUME
30 dimitri 1.38 # include "SALT_PLUME.h"
31     #endif
32 mlosch 1.1 #ifdef ALLOW_AUTODIFF_TAMC
33     # include "tamc.h"
34     #endif
35 dimitri 1.38
36 mlosch 1.1 C === Routine arguments ===
37     C myTime - Simulation time
38     C myIter - Simulation timestep number
39     C myThid - Thread no. that called this routine.
40     _RL myTime
41     INTEGER myIter, myThid
42     CEndOfInterface
43    
44     C === Local variables ===
45     C i,j,bi,bj - Loop counters
46    
47     INTEGER i, j, bi, bj
48 mlosch 1.3 C number of surface interface layer
49     INTEGER kSurface
50 mlosch 1.8 C constants
51 dimitri 1.26 _RL TBC, SDF, ICE2SNOW
52 mlosch 1.8 _RL QI, recip_QI, QS, recip_QS
53     C auxillary variables
54 heimbach 1.45 _RL snowEnergy
55 heimbach 1.32 _RL growthHEFF, growthNeg
56 mlosch 1.1 #ifdef ALLOW_SEAICE_FLOODING
57 mlosch 1.36 _RL hDraft
58 mlosch 1.1 #endif /* ALLOW_SEAICE_FLOODING */
59 heimbach 1.45 _RL availHeat (1:sNx,1:sNy)
60     _RL hEffOld (1:sNx,1:sNy)
61 dimitri 1.39 _RL GAREA (1:sNx,1:sNy)
62     _RL GHEFF (1:sNx,1:sNy)
63 mlosch 1.1 C RESID_HEAT is residual heat above freezing in equivalent m of ice
64 dimitri 1.39 _RL RESID_HEAT (1:sNx,1:sNy)
65 heimbach 1.32 #ifdef SEAICE_SALINITY
66 dimitri 1.39 _RL saltFluxAdjust(1:sNx,1:sNy)
67 heimbach 1.32 #endif
68 mlosch 1.1
69     C FICE - thermodynamic ice growth rate over sea ice in W/m^2
70     C >0 causes ice growth, <0 causes snow and sea ice melt
71     C FHEFF - effective thermodynamic ice growth rate over sea ice in W/m^2
72     C >0 causes ice growth, <0 causes snow and sea ice melt
73     C QNETO - thermodynamic ice growth rate over open water in W/m^2
74     C ( = surface heat flux )
75     C >0 causes ice growth, <0 causes snow and sea ice melt
76     C QNETI - net surface heat flux under ice in W/m^2
77     C QSWO - short wave heat flux over ocean in W/m^2
78     C QSWI - short wave heat flux under ice in W/m^2
79 dimitri 1.39 _RL FHEFF (1:sNx,1:sNy)
80     _RL FICE (1:sNx,1:sNy)
81     _RL QNETO (1:sNx,1:sNy)
82     _RL QNETI (1:sNx,1:sNy)
83     _RL QSWO (1:sNx,1:sNy)
84     _RL QSWI (1:sNx,1:sNy)
85 mlosch 1.1 C
86 dimitri 1.39 _RL HCORR (1:sNx,1:sNy)
87 dimitri 1.6 C actual ice thickness with upper and lower limit
88 dimitri 1.39 _RL HICE (1:sNx,1:sNy)
89 mlosch 1.1 C actual snow thickness
90 dimitri 1.39 _RL hSnwLoc (1:sNx,1:sNy)
91 mlosch 1.1 C wind speed
92 dimitri 1.39 _RL UG (1:sNx,1:sNy)
93 mlosch 1.1 _RL SPEED_SQ
94 mlosch 1.3 C local copy of AREA
95 mlosch 1.7 _RL areaLoc
96 mlosch 1.1
97 mlosch 1.7 #ifdef SEAICE_MULTICATEGORY
98 mlosch 1.1 INTEGER it
99     INTEGER ilockey
100     _RL RK
101 dimitri 1.39 _RL HICEP (1:sNx,1:sNy)
102     _RL FICEP (1:sNx,1:sNy)
103     _RL QSWIP (1:sNx,1:sNy)
104 mlosch 1.1 #endif
105    
106 mlosch 1.36 #ifdef ALLOW_DIAGNOSTICS
107     LOGICAL DIAGNOSTICS_IS_ON
108     EXTERNAL DIAGNOSTICS_IS_ON
109     #endif
110    
111 mlosch 1.1 if ( buoyancyRelation .eq. 'OCEANICP' ) then
112     kSurface = Nr
113     else
114     kSurface = 1
115     endif
116    
117 mlosch 1.3 C FREEZING TEMP. OF SEA WATER (deg C)
118     TBC = SEAICE_freeze
119     C RATIO OF WATER DESITY TO SNOW DENSITY
120 mlosch 1.53 SDF = 1000.0 _d 0/SEAICE_rhoSnow
121 mlosch 1.8 C RATIO OF SEA ICE DENSITY to SNOW DENSITY
122 mlosch 1.10 ICE2SNOW = SDF * ICE2WATR
123 mlosch 1.8 C HEAT OF FUSION OF ICE (m^3/J)
124     QI = 302.0 _d +06
125     recip_QI = 1.0 _d 0 / QI
126 mlosch 1.3 C HEAT OF FUSION OF SNOW (J/m^3)
127     QS = 1.1 _d +08
128 mlosch 1.8 recip_QS = 1.1 _d 0 / QS
129 mlosch 1.1
130     DO bj=myByLo(myThid),myByHi(myThid)
131     DO bi=myBxLo(myThid),myBxHi(myThid)
132     c
133     #ifdef ALLOW_AUTODIFF_TAMC
134     act1 = bi - myBxLo(myThid)
135     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
136     act2 = bj - myByLo(myThid)
137     max2 = myByHi(myThid) - myByLo(myThid) + 1
138     act3 = myThid - 1
139     max3 = nTx*nTy
140     act4 = ikey_dynamics - 1
141     iicekey = (act1 + 1) + act2*max1
142     & + act3*max1*max2
143     & + act4*max1*max2*max3
144     #endif /* ALLOW_AUTODIFF_TAMC */
145     C
146     C initialise a few fields
147     C
148 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
149     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
150     CADJ & key = iicekey, byte = isbyte
151     CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj,
152     CADJ & key = iicekey, byte = isbyte
153     CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj,
154     CADJ & key = iicekey, byte = isbyte
155     #endif /* ALLOW_AUTODIFF_TAMC */
156 mlosch 1.1 DO J=1,sNy
157     DO I=1,sNx
158 mlosch 1.8 FHEFF(I,J) = 0.0 _d 0
159     FICE (I,J) = 0.0 _d 0
160     QNETO(I,J) = 0.0 _d 0
161     QNETI(I,J) = 0.0 _d 0
162     QSWO (I,J) = 0.0 _d 0
163     QSWI (I,J) = 0.0 _d 0
164     HCORR(I,J) = 0.0 _d 0
165     RESID_HEAT(I,J) = 0.0 _d 0
166 heimbach 1.32 #ifdef SEAICE_SALINITY
167 heimbach 1.33 saltFluxAdjust(I,J) = 0.0 _d 0
168 heimbach 1.32 #endif
169 dimitri 1.39 #ifdef SEAICE_MULTICATEGORY
170     FICEP(I,J) = 0.0 _d 0
171     QSWIP(I,J) = 0.0 _d 0
172     #endif
173 mlosch 1.1 ENDDO
174     ENDDO
175 heimbach 1.49 DO J=1-oLy,sNy+oLy
176     DO I=1-oLx,sNx+oLx
177     saltWtrIce(I,J,bi,bj) = 0.0 _d 0
178     frWtrIce(I,J,bi,bj) = 0.0 _d 0
179     frWtrAtm(I,J,bi,bj) = 0.0 _d 0
180     ENDDO
181     ENDDO
182 mlosch 1.1 #ifdef ALLOW_AUTODIFF_TAMC
183 heimbach 1.2 CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
184     CADJ & key = iicekey, byte = isbyte
185     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
186     CADJ & key = iicekey, byte = isbyte
187 mlosch 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
188 heimbach 1.2 DO J=1,sNy
189     DO I=1,sNx
190 dimitri 1.6 C COMPUTE ACTUAL ICE THICKNESS AND PUT MINIMUM/MAXIMUM
191     C ON ICE THICKNESS FOR BUDGET COMPUTATION
192 mlosch 1.8 C The default of A22 = 0.15 is a common threshold for defining
193     C the ice edge. This ice concentration usually does not occur
194     C due to thermodynamics but due to advection.
195 mlosch 1.7 areaLoc = MAX(A22,AREA(I,J,2,bi,bj))
196     HICE(I,J) = HEFF(I,J,2,bi,bj)/areaLoc
197 mlosch 1.8 C Do we know what this is for?
198 dimitri 1.6 HICE(I,J) = MAX(HICE(I,J),0.05 _d +00)
199 mlosch 1.8 C Capping the actual ice thickness effectively enforces a
200     C minimum of heat flux through the ice and helps getting rid of
201     C very thick ice.
202 dimitri 1.20 cdm actually, this does exactly the opposite, i.e., ice is thicker
203     cdm when HICE is capped, so I am commenting out
204     cdm HICE(I,J) = MIN(HICE(I,J),9.0 _d +00)
205 mlosch 1.7 hSnwLoc(I,J) = HSNOW(I,J,bi,bj)/areaLoc
206 heimbach 1.2 ENDDO
207     ENDDO
208 mlosch 1.1
209     C NOW DETERMINE MIXED LAYER TEMPERATURE
210     DO J=1,sNy
211     DO I=1,sNx
212     TMIX(I,J,bi,bj)=theta(I,J,kSurface,bi,bj)+273.16 _d +00
213     #ifdef SEAICE_DEBUG
214     TMIX(I,J,bi,bj)=MAX(TMIX(I,J,bi,bj),271.2 _d +00)
215     #endif
216     ENDDO
217     ENDDO
218    
219     C THERMAL WIND OF ATMOSPHERE
220     DO J=1,sNy
221     DO I=1,sNx
222 mlosch 1.54 C copy the wind speed computed in exf_wind.F to UG
223     UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))
224     CML this is the old code, which does the same
225     CML SPEED_SQ = UWIND(I,J,bi,bj)**2 + VWIND(I,J,bi,bj)**2
226     CML IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
227     CML UG(I,J)=SEAICE_EPS
228     CML ELSE
229     CML UG(I,J)=SQRT(SPEED_SQ)
230     CML ENDIF
231 mlosch 1.1 ENDDO
232     ENDDO
233    
234    
235     #ifdef ALLOW_AUTODIFF_TAMC
236 heimbach 1.52 cphCADJ STORE heff = comlev1, key = ikey_dynamics, byte = isbyte
237     cphCADJ STORE hsnow = comlev1, key = ikey_dynamics, byte = isbyte
238     cphCADJ STORE uwind = comlev1, key = ikey_dynamics, byte = isbyte
239     cphCADJ STORE vwind = comlev1, key = ikey_dynamics, byte = isbyte
240 heimbach 1.2 c
241 heimbach 1.52 CADJ STORE tice = comlev1, key = ikey_dynamics, byte = isbyte
242 mlosch 1.7 # ifdef SEAICE_MULTICATEGORY
243 heimbach 1.52 CADJ STORE tices = comlev1, key = ikey_dynamics, byte = isbyte
244 mlosch 1.1 # endif
245     #endif /* ALLOW_AUTODIFF_TAMC */
246    
247     C NOW DETERMINE GROWTH RATES
248     C FIRST DO OPEN WATER
249     CALL SEAICE_BUDGET_OCEAN(
250     I UG,
251     U TMIX,
252     O QNETO, QSWO,
253 mlosch 1.11 I bi, bj, myThid )
254 mlosch 1.1
255     C NOW DO ICE
256 dimitri 1.40 IF (useRelativeWind) THEN
257     C Compute relative wind speed over sea ice.
258     DO J=1,sNy
259     DO I=1,sNx
260     SPEED_SQ =
261     & (uWind(I,J,bi,bj)
262     & +0.5 _d 0*(uVel(i,j,1,bi,bj)+uVel(i+1,j,1,bi,bj))
263     & -0.5 _d 0*(uice(i,j,1,bi,bj)+uice(i+1,j,1,bi,bj)))**2
264     & +(vWind(I,J,bi,bj)
265     & +0.5 _d 0*(vVel(i,j,1,bi,bj)+vVel(i,j+1,1,bi,bj))
266     & -0.5 _d 0*(vice(i,j,1,bi,bj)+vice(i,j+1,1,bi,bj)))**2
267     IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
268     UG(I,J)=SEAICE_EPS
269     ELSE
270     UG(I,J)=SQRT(SPEED_SQ)
271     ENDIF
272     ENDDO
273     ENDDO
274     ENDIF
275 mlosch 1.7 #ifdef SEAICE_MULTICATEGORY
276     C-- Start loop over muli-categories
277 mlosch 1.1 DO IT=1,MULTDIM
278     #ifdef ALLOW_AUTODIFF_TAMC
279     ilockey = (iicekey-1)*MULTDIM + IT
280     CADJ STORE tices(:,:,it,bi,bj) = comlev1_multdim,
281     CADJ & key = ilockey, byte = isbyte
282     #endif /* ALLOW_AUTODIFF_TAMC */
283 mlosch 1.7 RK=REAL(IT)
284 mlosch 1.1 DO J=1,sNy
285     DO I=1,sNx
286 mlosch 1.7 HICEP(I,J)=(HICE(I,J)/MULTDIM)*((2.0 _d 0*RK)-1.0 _d 0)
287 mlosch 1.1 TICE(I,J,bi,bj)=TICES(I,J,IT,bi,bj)
288     ENDDO
289     ENDDO
290     CALL SEAICE_BUDGET_ICE(
291 mlosch 1.5 I UG, HICEP, hSnwLoc,
292 mlosch 1.1 U TICE,
293 mlosch 1.7 O FICEP, QSWIP,
294 mlosch 1.11 I bi, bj, myThid )
295 mlosch 1.1 DO J=1,sNy
296     DO I=1,sNx
297 mlosch 1.7 C average surface heat fluxes/growth rates
298     FICE (I,J) = FICE(I,J) + FICEP(I,J)/MULTDIM
299     QSWI (I,J) = QSWI(I,J) + QSWIP(I,J)/MULTDIM
300     TICES(I,J,IT,bi,bj) = TICE(I,J,bi,bj)
301 mlosch 1.1 ENDDO
302     ENDDO
303     ENDDO
304 mlosch 1.7 C-- End loop over multi-categories
305     #else /* SEAICE_MULTICATEGORY */
306 mlosch 1.1 CALL SEAICE_BUDGET_ICE(
307     I UG, HICE, hSnwLoc,
308     U TICE,
309     O FICE, QSWI,
310 mlosch 1.11 I bi, bj, myThid )
311 mlosch 1.7 #endif /* SEAICE_MULTICATEGORY */
312 mlosch 1.1
313 mlosch 1.3 #ifdef ALLOW_AUTODIFF_TAMC
314     CADJ STORE theta(:,:,:,bi,bj)= comlev1_bibj,
315     CADJ & key = iicekey, byte = isbyte
316     CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
317     CADJ & key = iicekey, byte = isbyte
318     #endif /* ALLOW_AUTODIFF_TAMC */
319 mlosch 1.8 C
320     C-- compute and apply ice growth due to oceanic heat flux from below
321     C
322 mlosch 1.3 DO J=1,sNy
323     DO I=1,sNx
324     C-- Create or melt sea-ice so that first-level oceanic temperature
325     C is approximately at the freezing point when there is sea-ice.
326 mlosch 1.8 C Initially the units of YNEG/availHeat are m of sea-ice.
327 mlosch 1.3 C The factor dRf(1)/72.0764, used to convert temperature
328     C change in deg K to m of sea-ice, is approximately:
329     C dRf(1) * (sea water heat capacity = 3996 J/kg/K)
330     C * (density of sea-water = 1026 kg/m^3)
331     C / (latent heat of fusion of sea-ice = 334000 J/kg)
332     C / (density of sea-ice = 910 kg/m^3)
333 mlosch 1.8 C Negative YNEG/availHeat leads to ice growth.
334     C Positive YNEG/availHeat leads to ice melting.
335 mlosch 1.3 IF ( .NOT. inAdMode ) THEN
336     #ifdef SEAICE_VARIABLE_FREEZING_POINT
337     TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0
338     #endif /* SEAICE_VARIABLE_FREEZING_POINT */
339 dimitri 1.43 IF ( theta(I,J,kSurface,bi,bj) .GE. TBC ) THEN
340 heimbach 1.45 availHeat(i,j) = SEAICE_availHeatFrac
341 dimitri 1.42 & * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
342     & * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0
343     ELSE
344 heimbach 1.45 availHeat(i,j) = SEAICE_availHeatFracFrz
345 dimitri 1.42 & * (theta(I,J,kSurface,bi,bj)-TBC) * dRf(kSurface)
346     & * hFacC(i,j,kSurface,bi,bj) / 72.0764 _d 0
347     ENDIF
348 mlosch 1.3 ELSE
349 heimbach 1.45 availHeat(i,j) = 0.
350 mlosch 1.3 ENDIF
351 mlosch 1.8 C local copy of old effective ice thickness
352 heimbach 1.45 hEffOld(I,J) = HEFF(I,J,1,bi,bj)
353 mlosch 1.8 C Melt (availHeat>0) or create (availHeat<0) sea ice
354 heimbach 1.45 HEFF(I,J,1,bi,bj) = MAX(ZERO,HEFF(I,J,1,bi,bj)-availHeat(I,J))
355 mlosch 1.8 C
356 heimbach 1.45 YNEG(I,J,bi,bj) = hEffOld(I,J) - HEFF(I,J,1,bi,bj)
357 mlosch 1.8 C
358 heimbach 1.49 saltWtrIce(I,J,bi,bj) = saltWtrIce(I,J,bi,bj)
359     & - YNEG(I,J,bi,bj)
360 heimbach 1.45 RESID_HEAT(I,J) = availHeat(I,J) - YNEG(I,J,bi,bj)
361 mlosch 1.3 C YNEG now contains m of ice melted (>0) or created (<0)
362 dimitri 1.16 C saltWtrIce contains m of ice melted (<0) or created (>0)
363 mlosch 1.3 C RESID_HEAT is residual heat above freezing in equivalent m of ice
364     ENDDO
365     ENDDO
366    
367 mlosch 1.1 cph(
368     #ifdef ALLOW_AUTODIFF_TAMC
369 heimbach 1.52 cphCADJ STORE heff = comlev1, key = ikey_dynamics, byte = isbyte
370     cphCADJ STORE hsnow = comlev1, key = ikey_dynamics, byte = isbyte
371 mlosch 1.1 #endif
372     cph)
373     c
374     #ifdef ALLOW_AUTODIFF_TAMC
375     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
376     CADJ & key = iicekey, byte = isbyte
377     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
378     CADJ & key = iicekey, byte = isbyte
379 heimbach 1.2 CADJ STORE fice(:,:) = comlev1_bibj,
380 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
381     #endif /* ALLOW_AUTODIFF_TAMC */
382     cph)
383 mlosch 1.8 C
384     C-- compute and apply ice growth due to atmospheric fluxes from above
385     C
386 mlosch 1.1 DO J=1,sNy
387     DO I=1,sNx
388 mlosch 1.3 C NOW CALCULATE CORRECTED effective growth in J/m^2 (>0=melt)
389     GHEFF(I,J)=-SEAICE_deltaTtherm*FICE(I,J)*AREA(I,J,2,bi,bj)
390 mlosch 1.1 ENDDO
391     ENDDO
392 heimbach 1.2
393 mlosch 1.1 #ifdef ALLOW_AUTODIFF_TAMC
394 mlosch 1.3 CADJ STORE fice(:,:) = comlev1_bibj,
395     CADJ & key = iicekey, byte = isbyte
396 mlosch 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
397    
398     DO J=1,sNy
399     DO I=1,sNx
400     IF(FICE(I,J).LT.ZERO.AND.AREA(I,J,2,bi,bj).GT.ZERO) THEN
401 mlosch 1.8 C use FICE to melt snow and CALCULATE CORRECTED GROWTH
402     C effective snow thickness in J/m^2
403     snowEnergy=HSNOW(I,J,bi,bj)*QS
404     IF(GHEFF(I,J).LE.snowEnergy) THEN
405     C not enough heat to melt all snow; use up all heat flux FICE
406 mlosch 1.1 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)-GHEFF(I,J)/QS
407 mlosch 1.8 C SNOW CONVERTED INTO WATER AND THEN INTO equivalent m of ICE melt
408     C The factor 1/ICE2SNOW converts m of snow to m of sea-ice
409 heimbach 1.49 frWtrIce(I,J,bi,bj) =
410     & frWtrIce(I,J,bi,bj) - GHEFF(I,J)/(QS*ICE2SNOW)
411 mlosch 1.8 FICE (I,J) = ZERO
412 mlosch 1.1 ELSE
413 mlosch 1.8 C enought heat to melt snow completely;
414     C compute remaining heat flux that will melt ice
415     FICE(I,J)=-(GHEFF(I,J)-snowEnergy)/
416 mlosch 1.1 & SEAICE_deltaTtherm/AREA(I,J,2,bi,bj)
417     C convert all snow to melt water (fresh water flux)
418 heimbach 1.49 frWtrIce(I,J,bi,bj) = frWtrIce(I,J,bi,bj)
419 mlosch 1.8 & -HSNOW(I,J,bi,bj)/ICE2SNOW
420 dimitri 1.37 HSNOW(I,J,bi,bj)=0.0 _d 0
421 mlosch 1.1 END IF
422     END IF
423 heimbach 1.2 ENDDO
424     ENDDO
425 mlosch 1.1
426 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
427 mlosch 1.3 CADJ STORE fice(:,:) = comlev1_bibj,
428     CADJ & key = iicekey, byte = isbyte
429 heimbach 1.2 #endif /* ALLOW_AUTODIFF_TAMC */
430    
431     DO J=1,sNy
432     DO I=1,sNx
433 mlosch 1.8 C now get cell averaged growth rate in W/m^2, >0 causes ice growth
434 mlosch 1.1 FHEFF(I,J)= FICE(I,J) * AREA(I,J,2,bi,bj)
435     & + QNETO(I,J) * (ONE-AREA(I,J,2,bi,bj))
436     ENDDO
437     ENDDO
438     cph(
439     #ifdef ALLOW_AUTODIFF_TAMC
440     CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
441     CADJ & key = iicekey, byte = isbyte
442 mlosch 1.3 CADJ STORE fice(:,:) = comlev1_bibj,
443 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
444 mlosch 1.3 CADJ STORE fheff(:,:) = comlev1_bibj,
445 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
446 mlosch 1.3 CADJ STORE qneto(:,:) = comlev1_bibj,
447 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
448 mlosch 1.3 CADJ STORE qswi(:,:) = comlev1_bibj,
449 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
450 mlosch 1.3 CADJ STORE qswo(:,:) = comlev1_bibj,
451 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
452     #endif /* ALLOW_AUTODIFF_TAMC */
453     cph)
454 mlosch 1.8 C
455     C First update (freeze or melt) ice area
456     C
457 mlosch 1.1 DO J=1,sNy
458     DO I=1,sNx
459 mlosch 1.8 C negative growth in meters of ice (>0 for melting)
460     growthNeg = -SEAICE_deltaTtherm*FHEFF(I,J)*recip_QI
461     C negative growth must not exceed effective ice thickness (=volume)
462     C (that is, cannot melt more than all the ice)
463     growthHEFF = -ONE*MIN(HEFF(I,J,1,bi,bj),growthNeg)
464     C growthHEFF < 0 means melting
465     HCORR(I,J) = MIN(ZERO,growthHEFF)
466     C gain of new effective ice thickness over open water (>0 by definition)
467     GAREA(I,J) = MAX(ZERO,SEAICE_deltaTtherm*QNETO(I,J)*recip_QI)
468     CML removed these loops and moved TAMC store directive up
469     CML ENDDO
470     CML ENDDO
471     CML#ifdef ALLOW_AUTODIFF_TAMC
472     CMLCADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
473     CMLCADJ & key = iicekey, byte = isbyte
474     CML#endif
475     CML DO J=1,sNy
476     CML DO I=1,sNx
477     C Here we finally compute the new AREA
478     AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)+
479     & (ONE-AREA(I,J,2,bi,bj))*GAREA(I,J)/HO
480     & +HALF*HCORR(I,J)*AREA(I,J,2,bi,bj)
481     & /(HEFF(I,J,1,bi,bj)+.00001 _d 0)
482 mlosch 1.1 ENDDO
483     ENDDO
484     #ifdef ALLOW_AUTODIFF_TAMC
485     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
486     CADJ & key = iicekey, byte = isbyte
487     #endif
488 mlosch 1.8 C
489     C now update (freeze or melt) HEFF
490     C
491 mlosch 1.1 DO J=1,sNy
492     DO I=1,sNx
493 mlosch 1.8 C negative growth (>0 for melting) of existing ice in meters
494     growthNeg = -SEAICE_deltaTtherm*
495     & FICE(I,J)*recip_QI*AREA(I,J,2,bi,bj)
496     C negative growth must not exceed effective ice thickness (=volume)
497     C (that is, cannot melt more than all the ice)
498     growthHEFF = -ONE*MIN(HEFF(I,J,1,bi,bj),growthNeg)
499     C growthHEFF < 0 means melting
500     HEFF(I,J,1,bi,bj)= HEFF(I,J,1,bi,bj) + growthHEFF
501     C add effective growth to fresh water of ice
502 heimbach 1.49 saltWtrIce(I,J,bi,bj) = saltWtrIce(I,J,bi,bj) + growthHEFF
503 mlosch 1.8
504     C now calculate QNETI under ice (if any) as the difference between
505     C the available "heat flux" growthNeg and the actual growthHEFF;
506     C keep in mind that growthNeg and growthHEFF have different signs
507     C by construction
508     QNETI(I,J) = (growthHEFF + growthNeg)*QI/SEAICE_deltaTtherm
509 mlosch 1.1
510 mlosch 1.8 C now update other things
511 mlosch 1.1
512 jmc 1.13 #ifdef ALLOW_ATM_TEMP
513 mlosch 1.1 IF(FICE(I,J).GT.ZERO) THEN
514 mlosch 1.8 C freezing, add precip as snow
515 mlosch 1.3 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm*
516 mlosch 1.1 & PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*SDF
517     ELSE
518 mlosch 1.8 C add precip as rain, water converted into equivalent m of
519 mlosch 1.10 C ice by 1/ICE2WATR.
520 mlosch 1.9 C Do not get confused by the sign:
521     C precip > 0 for downward flux of fresh water
522     C frWtrIce > 0 for more ice (corresponds to an upward "fresh water flux"),
523     C so that here the rain is added *as if* it is melted ice (which is not
524     C true, but just a trick; physically the rain just runs as water
525     C through the ice into the ocean)
526 heimbach 1.49 frWtrIce(I,J,bi,bj) = frWtrIce(I,J,bi,bj)
527 mlosch 1.1 & -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*
528 mlosch 1.10 & SEAICE_deltaTtherm/ICE2WATR
529 mlosch 1.1 ENDIF
530 jmc 1.13 #else /* ALLOW_ATM_TEMP */
531     STOP 'ABNORMAL END: S/R THSICE_GROWTH: ATM_TEMP undef'
532     #endif /* ALLOW_ATM_TEMP */
533 mlosch 1.1
534 heimbach 1.21 ENDDO
535     ENDDO
536    
537     #ifdef ALLOW_AUTODIFF_TAMC
538 heimbach 1.45 cphCADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
539     cphCADJ & key = iicekey, byte = isbyte
540 heimbach 1.21 #endif
541 heimbach 1.34
542     cph( very sensitive bit here by JZ
543     #ifndef SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING
544 heimbach 1.21 DO J=1,sNy
545     DO I=1,sNx
546 mlosch 1.8 C Now melt snow if there is residual heat left in surface level
547     C Note that units of YNEG and frWtrIce are m of ice
548     IF( RESID_HEAT(I,J) .GT. ZERO .AND.
549     & HSNOW(I,J,bi,bj) .GT. ZERO ) THEN
550 mlosch 1.10 GHEFF(I,J) = MIN( HSNOW(I,J,bi,bj)/SDF/ICE2WATR,
551 mlosch 1.3 & RESID_HEAT(I,J) )
552 mlosch 1.8 YNEG(I,J,bi,bj) = YNEG(I,J,bi,bj) +GHEFF(I,J)
553 heimbach 1.21 ENDIF
554     ENDDO
555     ENDDO
556    
557     #ifdef ALLOW_AUTODIFF_TAMC
558     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
559     CADJ & key = iicekey, byte = isbyte
560     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
561     CADJ & key = iicekey, byte = isbyte
562     #endif
563     DO J=1,sNy
564     DO I=1,sNx
565     IF( RESID_HEAT(I,J) .GT. ZERO .AND.
566     & HSNOW(I,J,bi,bj) .GT. ZERO ) THEN
567 mlosch 1.10 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-GHEFF(I,J)*SDF*ICE2WATR
568 heimbach 1.49 frWtrIce(I,J,bi,bj) = frWtrIce(I,J,bi,bj)-GHEFF(I,J)
569 mlosch 1.1 ENDIF
570 heimbach 1.21 ENDDO
571     ENDDO
572 heimbach 1.34 #endif
573     cph)
574 heimbach 1.21
575     #ifdef ALLOW_AUTODIFF_TAMC
576     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
577     CADJ & key = iicekey, byte = isbyte
578 heimbach 1.32 # ifdef SEAICE_SALINITY
579     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
580     CADJ & key = iicekey, byte = isbyte
581 dimitri 1.37 # endif /* SEAICE_SALINITY */
582     #endif /* ALLOW_AUTODIFF_TAMC */
583    
584 heimbach 1.47 #ifdef ALLOW_ATM_TEMP
585 heimbach 1.21 DO J=1,sNy
586     DO I=1,sNx
587    
588 mlosch 1.1 C NOW GET FRESH WATER FLUX
589 mlosch 1.3 EmPmR(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
590 mlosch 1.9 & ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
591     & * ( ONE - AREA(I,J,2,bi,bj) )
592 jmc 1.14 #ifdef ALLOW_RUNOFF
593 mlosch 1.9 & - RUNOFF(I,J,bi,bj)
594 dimitri 1.37 #endif /* ALLOW_RUNOFF */
595 heimbach 1.49 & + frWtrIce(I,J,bi,bj)*ICE2WATR/SEAICE_deltaTtherm
596     & + saltWtrIce(I,J,bi,bj)*ICE2WATR/SEAICE_deltaTtherm
597 jmc 1.35 & )*rhoConstFresh
598 jmc 1.15 #ifdef ALLOW_DIAGNOSTICS
599 heimbach 1.49 frWtrAtm(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
600 mlosch 1.48 & PRECIP(I,J,bi,bj)
601 jmc 1.15 & - EVAP(I,J,bi,bj)
602     & *( ONE - AREA(I,J,2,bi,bj) )
603     & + RUNOFF(I,J,bi,bj)
604 mlosch 1.48 & )*rhoConstFresh
605 dimitri 1.37 #endif /* ALLOW_DIAGNOSTICS */
606 dimitri 1.16
607 heimbach 1.47 ENDDO
608     ENDDO
609    
610 dimitri 1.24 C COMPUTE SURFACE SALT FLUX AND ADJUST ICE SALINITY
611 heimbach 1.47
612 dimitri 1.23 #ifdef SEAICE_SALINITY
613 heimbach 1.47
614     DO J=1,sNy
615     DO I=1,sNx
616 dimitri 1.27 C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean
617     IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN
618 heimbach 1.32 saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
619 dimitri 1.27 & HSALT(I,J,bi,bj) / SEAICE_deltaTtherm
620 dimitri 1.37 HSALT(I,J,bi,bj) = 0.0 _d 0
621 dimitri 1.27 ENDIF
622 heimbach 1.32 ENDDO
623     ENDDO
624    
625     #ifdef ALLOW_AUTODIFF_TAMC
626     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
627     CADJ & key = iicekey, byte = isbyte
628 dimitri 1.37 #endif /* ALLOW_AUTODIFF_TAMC */
629    
630 heimbach 1.32 DO J=1,sNy
631     DO I=1,sNx
632 dimitri 1.27 C saltWtrIce > 0 : m of sea ice that is created
633 heimbach 1.49 IF ( saltWtrIce(I,J,bi,bj) .GE. 0.0 ) THEN
634     saltFlux(I,J,bi,bj) =
635     & HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*
636 dimitri 1.28 & ICE2WATR*rhoConstFresh*SEAICE_salinity*
637 dimitri 1.31 & salt(I,j,kSurface,bi,bj)/SEAICE_deltaTtherm
638     #ifdef ALLOW_SALT_PLUME
639     C saltPlumeFlux is defined only during freezing:
640 heimbach 1.49 saltPlumeFlux(I,J,bi,bj)=
641     & HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*
642 dimitri 1.31 & ICE2WATR*rhoConstFresh*(1-SEAICE_salinity)*
643     & salt(I,j,kSurface,bi,bj)/SEAICE_deltaTtherm
644     #endif /* ALLOW_SALT_PLUME */
645 dimitri 1.27 C saltWtrIce < 0 : m of sea ice that is melted
646 dimitri 1.25 ELSE
647 heimbach 1.49 saltFlux(I,J,bi,bj) =
648     & HEFFM(I,J,bi,bj)*saltWtrIce(I,J,bi,bj)*
649     & HSALT(I,J,bi,bj)/
650     & (HEFF(I,J,1,bi,bj)-saltWtrIce(I,J,bi,bj))/
651 dimitri 1.25 & SEAICE_deltaTtherm
652 dimitri 1.37 #ifdef ALLOW_SALT_PLUME
653     saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
654     #endif /* ALLOW_SALT_PLUME */
655 dimitri 1.25 ENDIF
656 dimitri 1.37 C update HSALT based on surface saltFlux
657 dimitri 1.24 HSALT(I,J,bi,bj) = HSALT(I,J,bi,bj) +
658     & saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
659 heimbach 1.32 saltFlux(I,J,bi,bj) =
660     & saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
661 dimitri 1.27 C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean
662 dimitri 1.24 IF ( HEFF(I,J,1,bi,bj) .EQ. 0.0 ) THEN
663     saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -
664     & HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) /
665     & SEAICE_deltaTtherm
666 dimitri 1.37 HSALT(I,J,bi,bj) = 0.0 _d 0
667     #ifdef ALLOW_SALT_PLUME
668     saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
669     #endif /* ALLOW_SALT_PLUME */
670 dimitri 1.24 ENDIF
671 heimbach 1.47 ENDDO
672     ENDDO
673    
674 dimitri 1.27 #endif /* SEAICE_SALINITY */
675 dimitri 1.16
676 jmc 1.13 #else /* ALLOW_ATM_TEMP */
677 heimbach 1.47 STOP 'ABNORMAL END: S/R THSICE_GROWTH: ATM_TEMP undef'
678 jmc 1.13 #endif /* ALLOW_ATM_TEMP */
679 mlosch 1.1
680 heimbach 1.47 DO J=1,sNy
681     DO I=1,sNx
682 mlosch 1.1 C NOW GET TOTAL QNET AND QSW
683 mlosch 1.3 QNET(I,J,bi,bj) = QNETI(I,J) * AREA(I,J,2,bi,bj)
684     & +QNETO(I,J) * (ONE-AREA(I,J,2,bi,bj))
685     QSW(I,J,bi,bj) = QSWI(I,J) * AREA(I,J,2,bi,bj)
686     & +QSWO(I,J) * (ONE-AREA(I,J,2,bi,bj))
687 heimbach 1.21 ENDDO
688     ENDDO
689    
690     #ifdef ALLOW_AUTODIFF_TAMC
691 heimbach 1.45 cphCADJ STORE yneg(:,:,bi,bj) = comlev1_bibj,
692     cphCADJ & key = iicekey, byte = isbyte
693 heimbach 1.21 #endif
694     DO J=1,sNy
695     DO I=1,sNx
696 mlosch 1.1 C Now convert YNEG back to deg K.
697 dimitri 1.19 YNEG(I,J,bi,bj) = YNEG(I,J,bi,bj)*recip_dRf(kSurface) *
698     & recip_hFacC(i,j,kSurface,bi,bj)*72.0764 _d 0
699 heimbach 1.21 ENDDO
700     ENDDO
701 mlosch 1.1
702 heimbach 1.21 #ifdef ALLOW_AUTODIFF_TAMC
703     CADJ STORE yneg(:,:,bi,bj) = comlev1_bibj,
704     CADJ & key = iicekey, byte = isbyte
705     #endif
706     DO J=1,sNy
707     DO I=1,sNx
708 mlosch 1.1 C Add YNEG contribution to QNET
709 mlosch 1.3 QNET(I,J,bi,bj) = QNET(I,J,bi,bj)
710 mlosch 1.1 & +YNEG(I,J,bi,bj)/SEAICE_deltaTtherm
711     & *maskC(I,J,kSurface,bi,bj)
712 jmc 1.22 & *HeatCapacity_Cp*rUnit2mass
713 mlosch 1.1 & *drF(kSurface)*hFacC(i,j,kSurface,bi,bj)
714     ENDDO
715     ENDDO
716 heimbach 1.21
717 mlosch 1.1 #ifdef SEAICE_DEBUG
718     CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
719     CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
720     CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
721     #endif /* SEAICE_DEBUG */
722    
723     crg Added by Ralf Giering: do we need DO_WE_NEED_THIS ?
724     #define DO_WE_NEED_THIS
725     C NOW ZERO OUTSIDE POINTS
726     #ifdef ALLOW_AUTODIFF_TAMC
727     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
728     CADJ & key = iicekey, byte = isbyte
729     CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
730     CADJ & key = iicekey, byte = isbyte
731     #endif /* ALLOW_AUTODIFF_TAMC */
732     DO J=1,sNy
733     DO I=1,sNx
734     C NOW SET AREA(I,J,1,bi,bj)=0 WHERE NO ICE IS
735     AREA(I,J,1,bi,bj)=MIN(AREA(I,J,1,bi,bj)
736     & ,HEFF(I,J,1,bi,bj)/.0001 _d 0)
737     ENDDO
738     ENDDO
739     #ifdef ALLOW_AUTODIFF_TAMC
740     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
741     CADJ & key = iicekey, byte = isbyte
742     #endif /* ALLOW_AUTODIFF_TAMC */
743     DO J=1,sNy
744     DO I=1,sNx
745     C NOW TRUNCATE AREA
746     #ifdef DO_WE_NEED_THIS
747     AREA(I,J,1,bi,bj)=MIN(ONE,AREA(I,J,1,bi,bj))
748     ENDDO
749     ENDDO
750     #ifdef ALLOW_AUTODIFF_TAMC
751     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
752     CADJ & key = iicekey, byte = isbyte
753     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
754     CADJ & key = iicekey, byte = isbyte
755     #endif /* ALLOW_AUTODIFF_TAMC */
756     DO J=1,sNy
757     DO I=1,sNx
758 mlosch 1.3 AREA(I,J,1,bi,bj) = MAX(ZERO,AREA(I,J,1,bi,bj))
759     HSNOW(I,J,bi,bj) = MAX(ZERO,HSNOW(I,J,bi,bj))
760 dimitri 1.20 #endif /* DO_WE_NEED_THIS */
761 mlosch 1.3 AREA(I,J,1,bi,bj) = AREA(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)
762     HEFF(I,J,1,bi,bj) = HEFF(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)
763 dimitri 1.20 #ifdef SEAICE_CAP_HEFF
764     HEFF(I,J,1,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,1,bi,bj))
765     #endif /* SEAICE_CAP_HEFF */
766 mlosch 1.3 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
767 mlosch 1.1 ENDDO
768     ENDDO
769    
770 mlosch 1.36 #ifdef ALLOW_DIAGNOSTICS
771     IF ( useDiagnostics ) THEN
772     IF ( DIAGNOSTICS_IS_ON('SIthdgrh',myThid) ) THEN
773     C use (abuse) GHEFF to diagnose the total thermodynamic growth rate
774     DO J=1,sNy
775     DO I=1,sNx
776     GHEFF(I,J) = (HEFF(I,J,1,bi,bj)-HEFF(I,J,2,bi,bj))
777     & /SEAICE_deltaTtherm
778     ENDDO
779     ENDDO
780 mlosch 1.46 CALL DIAGNOSTICS_FILL(GHEFF,'SIthdgrh',0,1,3,bi,bj,myThid)
781 mlosch 1.36 ENDIF
782     ENDIF
783     #endif /* ALLOW_DIAGNOSTICS */
784    
785 mlosch 1.1 #ifdef ALLOW_SEAICE_FLOODING
786     IF ( SEAICEuseFlooding ) THEN
787     C convert snow to ice if submerged
788     DO J=1,sNy
789     DO I=1,sNx
790 mlosch 1.55 hDraft = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
791 mlosch 1.1 & +HEFF(I,J,1,bi,bj)*SEAICE_rhoIce)/1000. _d 0
792 mlosch 1.36 C here GEFF is the gain of ice due to flooding
793     GHEFF(I,J) = hDraft - MIN(hDraft,HEFF(I,J,1,bi,bj))
794     HEFF(I,J,1,bi,bj) = HEFF(I,J,1,bi,bj) + GHEFF(I,J)
795 mlosch 1.10 HSNOW(I,J,bi,bj) = MAX(0. _d 0,
796 mlosch 1.36 & HSNOW(I,J,bi,bj)-GHEFF(I,J)*ICE2SNOW)
797 mlosch 1.1 ENDDO
798     ENDDO
799 mlosch 1.36 #ifdef ALLOW_DIAGNOSTICS
800     IF ( useDiagnostics ) THEN
801     IF ( DIAGNOSTICS_IS_ON('SIsnwice',myThid) ) THEN
802     C turn GHEFF into a rate
803     DO J=1,sNy
804     DO I=1,sNx
805     GHEFF(I,J) = GHEFF(I,J)/SEAICE_deltaTtherm
806     ENDDO
807     ENDDO
808 mlosch 1.46 CALL DIAGNOSTICS_FILL(GHEFF,'SIsnwice',0,1,3,bi,bj,myThid)
809 mlosch 1.36 ENDIF
810     ENDIF
811     #endif /* ALLOW_DIAGNOSTICS */
812 mlosch 1.1 ENDIF
813     #endif /* ALLOW_SEAICE_FLOODING */
814    
815     IF ( useRealFreshWaterFlux ) THEN
816     DO J=1,sNy
817     DO I=1,sNx
818     sIceLoad(i,j,bi,bj) = HEFF(I,J,1,bi,bj)*SEAICE_rhoIce
819 mlosch 1.55 & + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
820 mlosch 1.1 ENDDO
821     ENDDO
822     ENDIF
823    
824 dimitri 1.50 #ifdef SEAICE_AGE
825     C Sources and sinks for sea ice age
826     DO J=1,sNy
827     DO I=1,sNx
828     IF ( AREA(I,J,1,bi,bj) .GT. 0.15 ) THEN
829 dimitri 1.51 IceAge(i,j,bi,bj) = IceAge(i,j,bi,bj) +
830     & AREA(I,J,1,bi,bj) * SEAICE_deltaTtherm
831 dimitri 1.50 ELSE
832     IceAge(i,j,bi,bj) = ZERO
833     ENDIF
834     ENDDO
835     ENDDO
836     #endif
837    
838 mlosch 1.1 ENDDO
839     ENDDO
840    
841     RETURN
842     END
  ViewVC Help
Powered by ViewVC 1.1.22