/[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.18 - (hide annotations) (download)
Sat Jun 30 21:21:46 2007 UTC (16 years, 10 months ago) by dimitri
Branch: MAIN
Changes since 1.17: +3 -2 lines 
removing bug fix:
  24-Jun-07: bug fix for SEAICE_salinity: salt rejection was being double-counted
recovering previous version of verification/lab_sea/results/output* files
  SEAICE_salinity=0 should not change results if fix is correct
changing sign of saltWtrIce to saltFlux
  it should oppose contribution to EmPmR, which has opposite sign convention
1 dimitri 1.18 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_growth.F,v 1.17 2007/06/24 14:21:36 dimitri 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     #include "SEAICE_FFIELDS.h"
25    
26     #ifdef ALLOW_AUTODIFF_TAMC
27     # include "tamc.h"
28     #endif
29     C === Routine arguments ===
30     C myTime - Simulation time
31     C myIter - Simulation timestep number
32     C myThid - Thread no. that called this routine.
33     _RL myTime
34     INTEGER myIter, myThid
35     CEndOfInterface
36    
37     C === Local variables ===
38     C i,j,bi,bj - Loop counters
39    
40     INTEGER i, j, bi, bj
41 mlosch 1.3 C number of surface interface layer
42     INTEGER kSurface
43 mlosch 1.8 C constants
44 dimitri 1.16 _RL TBC, SDF, ICE2WATR, ICE2SNOW
45 mlosch 1.8 _RL QI, recip_QI, QS, recip_QS
46     C auxillary variables
47     _RL availHeat, hEffOld, snowEnergy, snowAsIce
48     _RL growthHEFF, growthNeg
49 mlosch 1.1 #ifdef ALLOW_SEAICE_FLOODING
50     _RL hDraft, hFlood
51     #endif /* ALLOW_SEAICE_FLOODING */
52     _RL GAREA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
53     _RL GHEFF ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
54     C RESID_HEAT is residual heat above freezing in equivalent m of ice
55     _RL RESID_HEAT ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
56    
57     C FICE - thermodynamic ice growth rate over sea ice in W/m^2
58     C >0 causes ice growth, <0 causes snow and sea ice melt
59     C FHEFF - effective thermodynamic ice growth rate over sea ice in W/m^2
60     C >0 causes ice growth, <0 causes snow and sea ice melt
61     C QNETO - thermodynamic ice growth rate over open water in W/m^2
62     C ( = surface heat flux )
63     C >0 causes ice growth, <0 causes snow and sea ice melt
64     C QNETI - net surface heat flux under ice in W/m^2
65     C QSWO - short wave heat flux over ocean in W/m^2
66     C QSWI - short wave heat flux under ice in W/m^2
67     _RL FHEFF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68     _RL FICE (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69     _RL QNETO (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70     _RL QNETI (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71     _RL QSWO (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72     _RL QSWI (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73     C
74     _RL HCORR (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 dimitri 1.16 C saltWtrIce contains m of salty ice melted (<0) or created (>0)
76     C for time being its salinity is assumed constant SEAICE_salinity
77     _RL saltWtrIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78     C frWtrIce contains m of freshwater ice melted (<0) or created (>0)
79     C that is, ice due to precipitation or snow
80 mlosch 1.8 _RL frWtrIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81 jmc 1.15 C frWtrAtm contains freshwater flux from the atmosphere
82     _RL frWtrAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83 dimitri 1.6 C actual ice thickness with upper and lower limit
84 mlosch 1.1 _RL HICE (1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
85     C actual snow thickness
86     _RL hSnwLoc(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
87     C wind speed
88     _RL UG (1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
89     _RL SPEED_SQ
90 mlosch 1.3 C local copy of AREA
91 mlosch 1.7 _RL areaLoc
92 mlosch 1.1
93 mlosch 1.7 #ifdef SEAICE_MULTICATEGORY
94 mlosch 1.1 INTEGER it
95     INTEGER ilockey
96     _RL RK
97     _RL HICEP(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
98     _RL FICEP(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
99 mlosch 1.7 _RL QSWIP(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
100 mlosch 1.1 #endif
101    
102     if ( buoyancyRelation .eq. 'OCEANICP' ) then
103     kSurface = Nr
104     else
105     kSurface = 1
106     endif
107    
108 mlosch 1.3 C FREEZING TEMP. OF SEA WATER (deg C)
109     TBC = SEAICE_freeze
110     C RATIO OF WATER DESITY TO SNOW DENSITY
111     SDF = 1000.0 _d 0/330.0 _d 0
112     C RATIO OF SEA ICE DESITY TO WATER DENSITY
113 mlosch 1.10 ICE2WATR = 0.920 _d 0
114     C this makes more sense, but changes the results
115     C ICE2WATR = SEAICE_rhoIce * 1. _d -03
116 mlosch 1.8 C RATIO OF SEA ICE DENSITY to SNOW DENSITY
117 mlosch 1.10 ICE2SNOW = SDF * ICE2WATR
118 mlosch 1.8 C HEAT OF FUSION OF ICE (m^3/J)
119     QI = 302.0 _d +06
120     recip_QI = 1.0 _d 0 / QI
121 mlosch 1.3 C HEAT OF FUSION OF SNOW (J/m^3)
122     QS = 1.1 _d +08
123 mlosch 1.8 recip_QS = 1.1 _d 0 / QS
124 mlosch 1.1
125     DO bj=myByLo(myThid),myByHi(myThid)
126     DO bi=myBxLo(myThid),myBxHi(myThid)
127     c
128     #ifdef ALLOW_AUTODIFF_TAMC
129     act1 = bi - myBxLo(myThid)
130     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
131     act2 = bj - myByLo(myThid)
132     max2 = myByHi(myThid) - myByLo(myThid) + 1
133     act3 = myThid - 1
134     max3 = nTx*nTy
135     act4 = ikey_dynamics - 1
136     iicekey = (act1 + 1) + act2*max1
137     & + act3*max1*max2
138     & + act4*max1*max2*max3
139     #endif /* ALLOW_AUTODIFF_TAMC */
140     C
141     C initialise a few fields
142     C
143 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
144     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
145     CADJ & key = iicekey, byte = isbyte
146     CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj,
147     CADJ & key = iicekey, byte = isbyte
148     CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj,
149     CADJ & key = iicekey, byte = isbyte
150     #endif /* ALLOW_AUTODIFF_TAMC */
151 mlosch 1.1 DO J=1,sNy
152     DO I=1,sNx
153 mlosch 1.8 FHEFF(I,J) = 0.0 _d 0
154     FICE (I,J) = 0.0 _d 0
155 mlosch 1.7 #ifdef SEAICE_MULTICATEGORY
156 mlosch 1.8 FICEP(I,J) = 0.0 _d 0
157     QSWIP(I,J) = 0.0 _d 0
158 mlosch 1.1 #endif
159 mlosch 1.8 FHEFF(I,J) = 0.0 _d 0
160     FICE (I,J) = 0.0 _d 0
161     QNETO(I,J) = 0.0 _d 0
162     QNETI(I,J) = 0.0 _d 0
163     QSWO (I,J) = 0.0 _d 0
164     QSWI (I,J) = 0.0 _d 0
165     HCORR(I,J) = 0.0 _d 0
166 dimitri 1.16 saltWtrIce(I,J) = 0.0 _d 0
167 mlosch 1.8 frWtrIce(I,J) = 0.0 _d 0
168 jmc 1.15 frWtrAtm(I,J) = 0.0 _d 0
169 mlosch 1.8 RESID_HEAT(I,J) = 0.0 _d 0
170 mlosch 1.1 ENDDO
171     ENDDO
172     #ifdef ALLOW_AUTODIFF_TAMC
173 heimbach 1.2 CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
174     CADJ & key = iicekey, byte = isbyte
175     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
176     CADJ & key = iicekey, byte = isbyte
177 mlosch 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
178 heimbach 1.2 DO J=1,sNy
179     DO I=1,sNx
180 dimitri 1.6 C COMPUTE ACTUAL ICE THICKNESS AND PUT MINIMUM/MAXIMUM
181     C ON ICE THICKNESS FOR BUDGET COMPUTATION
182 mlosch 1.8 C The default of A22 = 0.15 is a common threshold for defining
183     C the ice edge. This ice concentration usually does not occur
184     C due to thermodynamics but due to advection.
185 mlosch 1.7 areaLoc = MAX(A22,AREA(I,J,2,bi,bj))
186     HICE(I,J) = HEFF(I,J,2,bi,bj)/areaLoc
187 mlosch 1.8 C Do we know what this is for?
188 dimitri 1.6 HICE(I,J) = MAX(HICE(I,J),0.05 _d +00)
189 mlosch 1.8 C Capping the actual ice thickness effectively enforces a
190     C minimum of heat flux through the ice and helps getting rid of
191     C very thick ice.
192 dimitri 1.6 HICE(I,J) = MIN(HICE(I,J),9.0 _d +00)
193 mlosch 1.7 hSnwLoc(I,J) = HSNOW(I,J,bi,bj)/areaLoc
194 heimbach 1.2 ENDDO
195     ENDDO
196 mlosch 1.1
197     C NOW DETERMINE MIXED LAYER TEMPERATURE
198     DO J=1,sNy
199     DO I=1,sNx
200     TMIX(I,J,bi,bj)=theta(I,J,kSurface,bi,bj)+273.16 _d +00
201     #ifdef SEAICE_DEBUG
202     TMIX(I,J,bi,bj)=MAX(TMIX(I,J,bi,bj),271.2 _d +00)
203     #endif
204     ENDDO
205     ENDDO
206    
207     C THERMAL WIND OF ATMOSPHERE
208     DO J=1,sNy
209     DO I=1,sNx
210     CML#ifdef SEAICE_EXTERNAL_FORCING
211     CMLC this seems to be more natural as we do compute the wind speed in
212     CMLC pkg/exf/exf_wind.F, but it changes the results
213     CML UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))
214     CML#else
215     SPEED_SQ = UWIND(I,J,bi,bj)**2 + VWIND(I,J,bi,bj)**2
216     IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
217     UG(I,J)=SEAICE_EPS
218     ELSE
219     UG(I,J)=SQRT(SPEED_SQ)
220     ENDIF
221     CML#endif /* SEAICE_EXTERNAL_FORCING */
222     ENDDO
223     ENDDO
224    
225    
226     #ifdef ALLOW_AUTODIFF_TAMC
227 heimbach 1.2 cphCADJ STORE heff = comlev1, key = ikey_dynamics
228     cphCADJ STORE hsnow = comlev1, key = ikey_dynamics
229     cphCADJ STORE uwind = comlev1, key = ikey_dynamics
230     cphCADJ STORE vwind = comlev1, key = ikey_dynamics
231     c
232 mlosch 1.1 CADJ STORE tice = comlev1, key = ikey_dynamics
233 mlosch 1.7 # ifdef SEAICE_MULTICATEGORY
234 mlosch 1.1 CADJ STORE tices = comlev1, key = ikey_dynamics
235     # endif
236     #endif /* ALLOW_AUTODIFF_TAMC */
237    
238     C NOW DETERMINE GROWTH RATES
239     C FIRST DO OPEN WATER
240     CALL SEAICE_BUDGET_OCEAN(
241     I UG,
242     U TMIX,
243     O QNETO, QSWO,
244 mlosch 1.11 I bi, bj, myThid )
245 mlosch 1.1
246     C NOW DO ICE
247 mlosch 1.7 #ifdef SEAICE_MULTICATEGORY
248     C-- Start loop over muli-categories
249 mlosch 1.1 DO IT=1,MULTDIM
250     #ifdef ALLOW_AUTODIFF_TAMC
251     ilockey = (iicekey-1)*MULTDIM + IT
252     CADJ STORE tices(:,:,it,bi,bj) = comlev1_multdim,
253     CADJ & key = ilockey, byte = isbyte
254     #endif /* ALLOW_AUTODIFF_TAMC */
255 mlosch 1.7 RK=REAL(IT)
256 mlosch 1.1 DO J=1,sNy
257     DO I=1,sNx
258 mlosch 1.7 HICEP(I,J)=(HICE(I,J)/MULTDIM)*((2.0 _d 0*RK)-1.0 _d 0)
259 mlosch 1.1 TICE(I,J,bi,bj)=TICES(I,J,IT,bi,bj)
260     ENDDO
261     ENDDO
262     CALL SEAICE_BUDGET_ICE(
263 mlosch 1.5 I UG, HICEP, hSnwLoc,
264 mlosch 1.1 U TICE,
265 mlosch 1.7 O FICEP, QSWIP,
266 mlosch 1.11 I bi, bj, myThid )
267 mlosch 1.1 DO J=1,sNy
268     DO I=1,sNx
269 mlosch 1.7 C average surface heat fluxes/growth rates
270     FICE (I,J) = FICE(I,J) + FICEP(I,J)/MULTDIM
271     QSWI (I,J) = QSWI(I,J) + QSWIP(I,J)/MULTDIM
272     TICES(I,J,IT,bi,bj) = TICE(I,J,bi,bj)
273 mlosch 1.1 ENDDO
274     ENDDO
275     ENDDO
276 mlosch 1.7 C-- End loop over multi-categories
277     #else /* SEAICE_MULTICATEGORY */
278 mlosch 1.1 CALL SEAICE_BUDGET_ICE(
279     I UG, HICE, hSnwLoc,
280     U TICE,
281     O FICE, QSWI,
282 mlosch 1.11 I bi, bj, myThid )
283 mlosch 1.7 #endif /* SEAICE_MULTICATEGORY */
284 mlosch 1.1
285 mlosch 1.3 #ifdef ALLOW_AUTODIFF_TAMC
286     CADJ STORE theta(:,:,:,bi,bj)= comlev1_bibj,
287     CADJ & key = iicekey, byte = isbyte
288     CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
289     CADJ & key = iicekey, byte = isbyte
290     #endif /* ALLOW_AUTODIFF_TAMC */
291 mlosch 1.8 C
292     C-- compute and apply ice growth due to oceanic heat flux from below
293     C
294 mlosch 1.3 DO J=1,sNy
295     DO I=1,sNx
296     C-- Create or melt sea-ice so that first-level oceanic temperature
297     C is approximately at the freezing point when there is sea-ice.
298 mlosch 1.8 C Initially the units of YNEG/availHeat are m of sea-ice.
299 mlosch 1.3 C The factor dRf(1)/72.0764, used to convert temperature
300     C change in deg K to m of sea-ice, is approximately:
301     C dRf(1) * (sea water heat capacity = 3996 J/kg/K)
302     C * (density of sea-water = 1026 kg/m^3)
303     C / (latent heat of fusion of sea-ice = 334000 J/kg)
304     C / (density of sea-ice = 910 kg/m^3)
305 mlosch 1.8 C Negative YNEG/availHeat leads to ice growth.
306     C Positive YNEG/availHeat leads to ice melting.
307 mlosch 1.3 IF ( .NOT. inAdMode ) THEN
308     #ifdef SEAICE_VARIABLE_FREEZING_POINT
309     TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) + 0.0901 _d 0
310     #endif /* SEAICE_VARIABLE_FREEZING_POINT */
311 mlosch 1.8 availHeat = (theta(I,J,kSurface,bi,bj)-TBC)
312 mlosch 1.3 & *dRf(1)/72.0764 _d 0
313     ELSE
314 mlosch 1.8 availHeat = 0.
315 mlosch 1.3 ENDIF
316 mlosch 1.8 C local copy of old effective ice thickness
317     hEffOld = HEFF(I,J,1,bi,bj)
318     C Melt (availHeat>0) or create (availHeat<0) sea ice
319     HEFF(I,J,1,bi,bj) = MAX(ZERO,HEFF(I,J,1,bi,bj)-availHeat)
320     C
321 dimitri 1.16 YNEG(I,J,bi,bj) = hEffOld - HEFF(I,J,1,bi,bj)
322 mlosch 1.8 C
323 dimitri 1.16 saltWtrIce(I,J) = saltWtrIce(I,J) - YNEG(I,J,bi,bj)
324     RESID_HEAT(I,J) = availHeat - YNEG(I,J,bi,bj)
325 mlosch 1.3 C YNEG now contains m of ice melted (>0) or created (<0)
326 dimitri 1.16 C saltWtrIce contains m of ice melted (<0) or created (>0)
327 mlosch 1.3 C RESID_HEAT is residual heat above freezing in equivalent m of ice
328     ENDDO
329     ENDDO
330    
331 mlosch 1.1 cph(
332     #ifdef ALLOW_AUTODIFF_TAMC
333     cphCADJ STORE heff = comlev1, key = ikey_dynamics
334     cphCADJ STORE hsnow = comlev1, key = ikey_dynamics
335     #endif
336     cph)
337     c
338     #ifdef ALLOW_AUTODIFF_TAMC
339     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
340     CADJ & key = iicekey, byte = isbyte
341     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
342     CADJ & key = iicekey, byte = isbyte
343 heimbach 1.2 CADJ STORE fice(:,:) = comlev1_bibj,
344 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
345     #endif /* ALLOW_AUTODIFF_TAMC */
346     cph)
347 mlosch 1.8 C
348     C-- compute and apply ice growth due to atmospheric fluxes from above
349     C
350 mlosch 1.1 DO J=1,sNy
351     DO I=1,sNx
352 mlosch 1.3 C NOW CALCULATE CORRECTED effective growth in J/m^2 (>0=melt)
353     GHEFF(I,J)=-SEAICE_deltaTtherm*FICE(I,J)*AREA(I,J,2,bi,bj)
354 mlosch 1.1 ENDDO
355     ENDDO
356 heimbach 1.2
357 mlosch 1.1 #ifdef ALLOW_AUTODIFF_TAMC
358 mlosch 1.3 CADJ STORE fice(:,:) = comlev1_bibj,
359     CADJ & key = iicekey, byte = isbyte
360 mlosch 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
361    
362     DO J=1,sNy
363     DO I=1,sNx
364     IF(FICE(I,J).LT.ZERO.AND.AREA(I,J,2,bi,bj).GT.ZERO) THEN
365 mlosch 1.8 C use FICE to melt snow and CALCULATE CORRECTED GROWTH
366     C effective snow thickness in J/m^2
367     snowEnergy=HSNOW(I,J,bi,bj)*QS
368     IF(GHEFF(I,J).LE.snowEnergy) THEN
369     C not enough heat to melt all snow; use up all heat flux FICE
370 mlosch 1.1 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)-GHEFF(I,J)/QS
371 mlosch 1.8 C SNOW CONVERTED INTO WATER AND THEN INTO equivalent m of ICE melt
372     C The factor 1/ICE2SNOW converts m of snow to m of sea-ice
373     frWtrIce(I,J) = frWtrIce(I,J) - GHEFF(I,J)/(QS*ICE2SNOW)
374     FICE (I,J) = ZERO
375 mlosch 1.1 ELSE
376 mlosch 1.8 C enought heat to melt snow completely;
377     C compute remaining heat flux that will melt ice
378     FICE(I,J)=-(GHEFF(I,J)-snowEnergy)/
379 mlosch 1.1 & SEAICE_deltaTtherm/AREA(I,J,2,bi,bj)
380     C convert all snow to melt water (fresh water flux)
381 mlosch 1.8 frWtrIce(I,J) = frWtrIce(I,J)
382     & -HSNOW(I,J,bi,bj)/ICE2SNOW
383 mlosch 1.1 HSNOW(I,J,bi,bj)=0.0
384     END IF
385     END IF
386 heimbach 1.2 ENDDO
387     ENDDO
388 mlosch 1.1
389 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
390 mlosch 1.3 CADJ STORE fice(:,:) = comlev1_bibj,
391     CADJ & key = iicekey, byte = isbyte
392 heimbach 1.2 #endif /* ALLOW_AUTODIFF_TAMC */
393    
394     DO J=1,sNy
395     DO I=1,sNx
396 mlosch 1.8 C now get cell averaged growth rate in W/m^2, >0 causes ice growth
397 mlosch 1.1 FHEFF(I,J)= FICE(I,J) * AREA(I,J,2,bi,bj)
398     & + QNETO(I,J) * (ONE-AREA(I,J,2,bi,bj))
399     ENDDO
400     ENDDO
401     cph(
402     #ifdef ALLOW_AUTODIFF_TAMC
403     CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
404     CADJ & key = iicekey, byte = isbyte
405     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
406     CADJ & key = iicekey, byte = isbyte
407 mlosch 1.3 CADJ STORE fice(:,:) = comlev1_bibj,
408 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
409 mlosch 1.3 CADJ STORE fheff(:,:) = comlev1_bibj,
410 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
411 mlosch 1.3 CADJ STORE qneto(:,:) = comlev1_bibj,
412 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
413 mlosch 1.3 CADJ STORE qswi(:,:) = comlev1_bibj,
414 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
415 mlosch 1.3 CADJ STORE qswo(:,:) = comlev1_bibj,
416 mlosch 1.1 CADJ & key = iicekey, byte = isbyte
417     #endif /* ALLOW_AUTODIFF_TAMC */
418     cph)
419     #ifdef ALLOW_AUTODIFF_TAMC
420     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
421     CADJ & key = iicekey, byte = isbyte
422     #endif
423 mlosch 1.8 C
424     C First update (freeze or melt) ice area
425     C
426 mlosch 1.1 DO J=1,sNy
427     DO I=1,sNx
428 mlosch 1.8 C negative growth in meters of ice (>0 for melting)
429     growthNeg = -SEAICE_deltaTtherm*FHEFF(I,J)*recip_QI
430     C negative growth must not exceed effective ice thickness (=volume)
431     C (that is, cannot melt more than all the ice)
432     growthHEFF = -ONE*MIN(HEFF(I,J,1,bi,bj),growthNeg)
433     C growthHEFF < 0 means melting
434     HCORR(I,J) = MIN(ZERO,growthHEFF)
435     C gain of new effective ice thickness over open water (>0 by definition)
436     GAREA(I,J) = MAX(ZERO,SEAICE_deltaTtherm*QNETO(I,J)*recip_QI)
437     CML removed these loops and moved TAMC store directive up
438     CML ENDDO
439     CML ENDDO
440     CML#ifdef ALLOW_AUTODIFF_TAMC
441     CMLCADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
442     CMLCADJ & key = iicekey, byte = isbyte
443     CML#endif
444     CML DO J=1,sNy
445     CML DO I=1,sNx
446     C Here we finally compute the new AREA
447     AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)+
448     & (ONE-AREA(I,J,2,bi,bj))*GAREA(I,J)/HO
449     & +HALF*HCORR(I,J)*AREA(I,J,2,bi,bj)
450     & /(HEFF(I,J,1,bi,bj)+.00001 _d 0)
451 mlosch 1.1 ENDDO
452     ENDDO
453     #ifdef ALLOW_AUTODIFF_TAMC
454     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
455     CADJ & key = iicekey, byte = isbyte
456     #endif
457 mlosch 1.8 C
458     C now update (freeze or melt) HEFF
459     C
460 mlosch 1.1 DO J=1,sNy
461     DO I=1,sNx
462 mlosch 1.8 C negative growth (>0 for melting) of existing ice in meters
463     growthNeg = -SEAICE_deltaTtherm*
464     & FICE(I,J)*recip_QI*AREA(I,J,2,bi,bj)
465     C negative growth must not exceed effective ice thickness (=volume)
466     C (that is, cannot melt more than all the ice)
467     growthHEFF = -ONE*MIN(HEFF(I,J,1,bi,bj),growthNeg)
468     C growthHEFF < 0 means melting
469     HEFF(I,J,1,bi,bj)= HEFF(I,J,1,bi,bj) + growthHEFF
470     C add effective growth to fresh water of ice
471 dimitri 1.16 saltWtrIce(I,J) = saltWtrIce(I,J) + growthHEFF
472 mlosch 1.8
473     C now calculate QNETI under ice (if any) as the difference between
474     C the available "heat flux" growthNeg and the actual growthHEFF;
475     C keep in mind that growthNeg and growthHEFF have different signs
476     C by construction
477     QNETI(I,J) = (growthHEFF + growthNeg)*QI/SEAICE_deltaTtherm
478 mlosch 1.1
479 mlosch 1.8 C now update other things
480 mlosch 1.1
481 jmc 1.13 #ifdef ALLOW_ATM_TEMP
482 mlosch 1.1 IF(FICE(I,J).GT.ZERO) THEN
483 mlosch 1.8 C freezing, add precip as snow
484 mlosch 1.3 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm*
485 mlosch 1.1 & PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*SDF
486     ELSE
487 mlosch 1.8 C add precip as rain, water converted into equivalent m of
488 mlosch 1.10 C ice by 1/ICE2WATR.
489 mlosch 1.9 C Do not get confused by the sign:
490     C precip > 0 for downward flux of fresh water
491     C frWtrIce > 0 for more ice (corresponds to an upward "fresh water flux"),
492     C so that here the rain is added *as if* it is melted ice (which is not
493     C true, but just a trick; physically the rain just runs as water
494     C through the ice into the ocean)
495 mlosch 1.8 frWtrIce(I,J) = frWtrIce(I,J)
496 mlosch 1.1 & -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*
497 mlosch 1.10 & SEAICE_deltaTtherm/ICE2WATR
498 mlosch 1.1 ENDIF
499 jmc 1.13 #else /* ALLOW_ATM_TEMP */
500     STOP 'ABNORMAL END: S/R THSICE_GROWTH: ATM_TEMP undef'
501     #endif /* ALLOW_ATM_TEMP */
502 mlosch 1.1
503 mlosch 1.8 C Now melt snow if there is residual heat left in surface level
504     C Note that units of YNEG and frWtrIce are m of ice
505 heimbach 1.4 cph( very sensitive bit here by JZ
506 mlosch 1.8 IF( RESID_HEAT(I,J) .GT. ZERO .AND.
507     & HSNOW(I,J,bi,bj) .GT. ZERO ) THEN
508 mlosch 1.10 GHEFF(I,J) = MIN( HSNOW(I,J,bi,bj)/SDF/ICE2WATR,
509 mlosch 1.3 & RESID_HEAT(I,J) )
510 mlosch 1.8 YNEG(I,J,bi,bj) = YNEG(I,J,bi,bj) +GHEFF(I,J)
511 mlosch 1.10 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-GHEFF(I,J)*SDF*ICE2WATR
512 mlosch 1.8 frWtrIce(I,J) = frWtrIce(I,J) -GHEFF(I,J)
513 mlosch 1.1 ENDIF
514 heimbach 1.4 cph)
515 mlosch 1.1
516 dimitri 1.16 #ifdef ALLOW_ATM_TEMP
517    
518 mlosch 1.1 C NOW GET FRESH WATER FLUX
519 mlosch 1.3 EmPmR(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
520 mlosch 1.9 & ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
521     & * ( ONE - AREA(I,J,2,bi,bj) )
522 jmc 1.14 #ifdef ALLOW_RUNOFF
523 mlosch 1.9 & - RUNOFF(I,J,bi,bj)
524 jmc 1.14 #endif
525 mlosch 1.10 & + frWtrIce(I,J)*ICE2WATR/SEAICE_deltaTtherm
526 dimitri 1.18 & + saltWtrIce(I,J)*ICE2WATR/SEAICE_deltaTtherm
527 mlosch 1.1 & )
528 jmc 1.15 #ifdef ALLOW_DIAGNOSTICS
529     frWtrAtm(I,J) = maskC(I,J,kSurface,bi,bj)*(
530     & PRECIP(I,J,bi,bj)
531     & - EVAP(I,J,bi,bj)
532     & *( ONE - AREA(I,J,2,bi,bj) )
533     & + RUNOFF(I,J,bi,bj)
534     & )
535     #endif
536 dimitri 1.16
537     C NOW GET SALT FLUX
538 dimitri 1.18 saltFlux(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*
539 dimitri 1.16 & saltWtrIce(I,J)*ICE2WATR/SEAICE_deltaTtherm*
540     & rhoConstFresh*SEAICE_salinity
541    
542 jmc 1.13 #else /* ALLOW_ATM_TEMP */
543     STOP 'ABNORMAL END: S/R THSICE_GROWTH: ATM_TEMP undef'
544     #endif /* ALLOW_ATM_TEMP */
545 mlosch 1.1
546     C NOW GET TOTAL QNET AND QSW
547 mlosch 1.3 QNET(I,J,bi,bj) = QNETI(I,J) * AREA(I,J,2,bi,bj)
548     & +QNETO(I,J) * (ONE-AREA(I,J,2,bi,bj))
549     QSW(I,J,bi,bj) = QSWI(I,J) * AREA(I,J,2,bi,bj)
550     & +QSWO(I,J) * (ONE-AREA(I,J,2,bi,bj))
551 mlosch 1.1
552     C Now convert YNEG back to deg K.
553 mlosch 1.3 YNEG(I,J,bi,bj) = YNEG(I,J,bi,bj)*recip_dRf(1)*72.0764 _d 0
554 mlosch 1.1
555     C Add YNEG contribution to QNET
556 mlosch 1.3 QNET(I,J,bi,bj) = QNET(I,J,bi,bj)
557 mlosch 1.1 & +YNEG(I,J,bi,bj)/SEAICE_deltaTtherm
558     & *maskC(I,J,kSurface,bi,bj)
559     & *HeatCapacity_Cp*recip_horiVertRatio*rhoConst
560     & *drF(kSurface)*hFacC(i,j,kSurface,bi,bj)
561    
562     ENDDO
563     ENDDO
564 jmc 1.15 #ifdef ALLOW_DIAGNOSTICS
565     IF ( useDiagnostics ) THEN
566     CALL DIAGNOSTICS_FILL(frWtrAtm,'SIatmFW ',0,1 ,2,bi,bj,myThid)
567     ENDIF
568     #endif /* ALLOW_DIAGNOSTICS */
569 mlosch 1.1
570     #ifdef SEAICE_DEBUG
571     c CALL PLOT_FIELD_XYRS( UWIND,'Current UWIND ', myIter, myThid )
572     c CALL PLOT_FIELD_XYRS( VWIND,'Current VWIND ', myIter, myThid )
573     CALL PLOT_FIELD_XYRS( GWATX,'Current GWATX ', myIter, myThid )
574     CALL PLOT_FIELD_XYRS( GWATY,'Current GWATY ', myIter, myThid )
575     CML CALL PLOT_FIELD_XYRL( FO,'Current FO ', myIter, myThid )
576     CML CALL PLOT_FIELD_XYRL( FHEFF,'Current FHEFF ', myIter, myThid )
577     CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
578     CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
579     CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
580 mlosch 1.8 CML DO j=1-OLy,sNy+OLy
581     CML DO i=1-OLx,sNx+OLx
582     CML GHEFF(I,J)=SQRT(UICE(I,J,1,bi,bj)**2+VICE(I,J,1,bi,bj)**2)
583     CML GAREA(I,J)=HEFF(I,J,1,bi,bj)
584     CML print*,'I J QNET:',I, J, QNET(i,j,bi,bj), QSW(I,J,bi,bj)
585     CML ENDDO
586     CML ENDDO
587     CML CALL PLOT_FIELD_XYRL( GHEFF,'Current UICE ', myIter, myThid )
588     CML CALL PLOT_FIELD_XYRL( GAREA,'Current HEFF ', myIter, myThid )
589 mlosch 1.1 DO j=1-OLy,sNy+OLy
590     DO i=1-OLx,sNx+OLx
591     if(HEFF(i,j,1,bi,bj).gt.1.) then
592     print '(A,2i4,3f10.2)','#### i j heff theta yneg',i,j,
593     & HEFF(i,j,1,bi,bj),theta(I,J,1,bi,bj),yneg(I,J,bi,bj)
594     print '(A,3f10.2)','QSW, QNET before/after correction',
595     & QSW(I,J,bi,bj),QNETI(I,J)*AREA(I,J,2,bi,bj)+
596     & (ONE-AREA(I,J,2,bi,bj))*QNETO(I,J), QNET(I,J,bi,bj)
597     endif
598     ENDDO
599     ENDDO
600     #endif /* SEAICE_DEBUG */
601    
602     crg Added by Ralf Giering: do we need DO_WE_NEED_THIS ?
603     #define DO_WE_NEED_THIS
604     C NOW ZERO OUTSIDE POINTS
605     #ifdef ALLOW_AUTODIFF_TAMC
606     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
607     CADJ & key = iicekey, byte = isbyte
608     CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
609     CADJ & key = iicekey, byte = isbyte
610     #endif /* ALLOW_AUTODIFF_TAMC */
611     DO J=1,sNy
612     DO I=1,sNx
613     C NOW SET AREA(I,J,1,bi,bj)=0 WHERE NO ICE IS
614     AREA(I,J,1,bi,bj)=MIN(AREA(I,J,1,bi,bj)
615     & ,HEFF(I,J,1,bi,bj)/.0001 _d 0)
616     ENDDO
617     ENDDO
618     #ifdef ALLOW_AUTODIFF_TAMC
619     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
620     CADJ & key = iicekey, byte = isbyte
621     #endif /* ALLOW_AUTODIFF_TAMC */
622     DO J=1,sNy
623     DO I=1,sNx
624     C NOW TRUNCATE AREA
625     #ifdef DO_WE_NEED_THIS
626     AREA(I,J,1,bi,bj)=MIN(ONE,AREA(I,J,1,bi,bj))
627     ENDDO
628     ENDDO
629     #ifdef ALLOW_AUTODIFF_TAMC
630     CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
631     CADJ & key = iicekey, byte = isbyte
632     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
633     CADJ & key = iicekey, byte = isbyte
634     #endif /* ALLOW_AUTODIFF_TAMC */
635     DO J=1,sNy
636     DO I=1,sNx
637 mlosch 1.3 AREA(I,J,1,bi,bj) = MAX(ZERO,AREA(I,J,1,bi,bj))
638     HSNOW(I,J,bi,bj) = MAX(ZERO,HSNOW(I,J,bi,bj))
639 mlosch 1.1 #endif
640 mlosch 1.3 AREA(I,J,1,bi,bj) = AREA(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)
641     HEFF(I,J,1,bi,bj) = HEFF(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)
642 mlosch 1.1 #ifdef DO_WE_NEED_THIS
643     c HEFF(I,J,1,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,1,bi,bj))
644     #endif
645 mlosch 1.3 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
646 mlosch 1.1 ENDDO
647     ENDDO
648    
649     #ifdef ALLOW_SEAICE_FLOODING
650     IF ( SEAICEuseFlooding ) THEN
651     C convert snow to ice if submerged
652 dimitri 1.16 C for time being, sea ice salinity is assumed constant, that is,
653     C contribution of snow flooding to freshening of sea ice is neglected
654 mlosch 1.1 DO J=1,sNy
655     DO I=1,sNx
656     hDraft = (HSNOW(I,J,bi,bj)*330. _d 0
657     & +HEFF(I,J,1,bi,bj)*SEAICE_rhoIce)/1000. _d 0
658     hFlood = hDraft - MIN(hDraft,HEFF(I,J,1,bi,bj))
659     HEFF(I,J,1,bi,bj) = HEFF(I,J,1,bi,bj) + hFlood
660 mlosch 1.10 HSNOW(I,J,bi,bj) = MAX(0. _d 0,
661     & HSNOW(I,J,bi,bj)-hFlood*ICE2SNOW)
662 mlosch 1.1 ENDDO
663     ENDDO
664     ENDIF
665     #endif /* ALLOW_SEAICE_FLOODING */
666    
667     IF ( useRealFreshWaterFlux ) THEN
668     DO J=1,sNy
669     DO I=1,sNx
670     sIceLoad(i,j,bi,bj) = HEFF(I,J,1,bi,bj)*SEAICE_rhoIce
671     & + HSNOW(I,J,bi,bj)* 330. _d 0
672     ENDDO
673     ENDDO
674     ENDIF
675    
676     ENDDO
677     ENDDO
678    
679     RETURN
680     END
  ViewVC Help
Powered by ViewVC 1.1.22