/[MITgcm]/MITgcm_contrib/SOSE/code_ad/seaice_growth_if.F
ViewVC logotype

Contents of /MITgcm_contrib/SOSE/code_ad/seaice_growth_if.F

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


Revision 1.1 - (show annotations) (download)
Fri Apr 23 19:55:13 2010 UTC (15 years, 3 months ago) by mmazloff
Branch: MAIN
CVS Tags: HEAD
original files

1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_growth_if.F,v 1.8 2009/08/02 19:45:42 jmc Exp $
2 C $Name: $
3
4 #include "SEAICE_OPTIONS.h"
5
6 C StartOfInterface
7 SUBROUTINE SEAICE_GROWTH_IF( myTime, myIter, myThid )
8 C /==========================================================\
9 C | SUBROUTINE seaice_growth_if |
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 #ifdef ALLOW_EXF
25 # include "EXF_OPTIONS.h"
26 # include "EXF_FIELDS.h"
27 # include "EXF_PARAM.h"
28 #endif
29 #ifdef ALLOW_SALT_PLUME
30 # include "SALT_PLUME.h"
31 #endif
32 #ifdef ALLOW_AUTODIFF_TAMC
33 # include "tamc.h"
34 #endif
35 C === Routine arguments ===
36 C myTime - Simulation time
37 C myIter - Simulation timestep number
38 C myThid - Thread no. that called this routine.
39 _RL myTime
40 INTEGER myIter, myThid
41 C EndOfInterface(global-font-lock-mode 1)
42
43 C === Local variables ===
44 C i,j,bi,bj - Loop counters
45
46 INTEGER i, j, bi, bj
47 C number of surface interface layer
48 INTEGER kSurface
49
50 C constants
51 _RL TBC, salinity_ice, SDF, ICE2SNOW,TMELT
52
53 #ifdef ALLOW_SEAICE_FLOODING
54 _RL hDraft, hFlood
55 #endif /* ALLOW_SEAICE_FLOODING */
56
57 C QNETI - net surface heat flux under ice in W/m^2
58 C QSWO - short wave heat flux over ocean in W/m^2
59 C QSWI - short wave heat flux under ice in W/m^2
60
61 _RL QNETI (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62 _RL QSWO (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63 _RL QSWI (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64
65 _RL QSWO_IN_FIRST_LAYER
66 & (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67 _RL QSWO_BELOW_FIRST_LAYER
68 & (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69
70 _RL QSW_absorb_in_ML (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71 _RL QSW_absorb_below_ML (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72
73 C actual ice thickness with upper and lower limit
74 _RL HICE_ACTUAL (1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
75
76 C actual snow thickness
77 _RL HSNOW_ACTUAL(1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
78
79 C wind speed
80 _RL UG (1-OLx:sNx+OLx, 1-OLy:sNy+OLy)
81 _RL SPEED_SQ
82
83 C IAN
84 _RL RHOI, RHOFW,CPW,LI,QI,QS,GAMMAT,GAMMA,RHOSW,RHOSN
85 _RL FL_C1,FL_C2,FL_C3,FL_C4,deltaHS,deltaHI
86
87 _RL NetExistingIceGrowthRate (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88 _RL IceGrowthRateUnderExistingIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89 _RL IceGrowthRateFromSurface (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90 _RL IceGrowthRateOpenWater (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91 _RL IceGrowthRateMixedLayer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92 _RL S_a_from_IGROW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93
94 _RL PredTempChange
95 & (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96 _RL PredTempChangeFromQSW
97 & (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98 _RL PredTempChangeFromOA_MQNET
99 & (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100 _RL PredTempChangeFromFIA
101 & (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102 _RL PredTempChangeFromNewIceVol
103 & (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104 _RL PredTempChangeFromF_IA_NET
105 & (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106 _RL PredTempChangeFromF_IO_NET
107 & (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108
109 _RL ExpectedIceVolumeChange (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110 _RL ExpectedSnowVolumeChange (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111 _RL ActualNewTotalVolumeChange(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112 _RL ActualNewTotalSnowMelt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113
114 _RL EnergyInNewTotalIceVolume (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115 _RL NetEnergyFluxOutOfSystem (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116
117 _RL ResidualHeatOutOfSystem (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118
119 _RL SnowAccRateOverIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120 _RL SmowAccOverIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121 _RL PrecipRateOverIceSurfaceToSea (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122
123 _RL PotSnowMeltRateFromSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124 _RL PotSnowMeltFromSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125 _RL SnowMeltFromSurface (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126 _RL SnowMeltRateFromSurface (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127
128 _RL FreshwaterContribFromSnowMelt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129 _RL FreshwaterContribFromIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130
131 _RL SurfHeatFluxConvergToSnowMelt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132 _RL EnergyToMeltSnowAndIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133 _RL EnergyToMeltSnowAndIce2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134
135 C dA/dt = S_a
136 _RL S_a (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137 C dh/dt = S_h
138 _RL S_h (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139 _RL S_hsnow (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140 _RL HSNOW_ORIG (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141
142 C F_ia - heat flux from ice to atmosphere (W/m^2)
143 C >0 causes ice growth, <0 causes snow and sea ice melt
144 _RL F_ia (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145 _RL F_ia_net (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146 _RL F_ia_net_before_snow (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147 _RL F_io_net (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148
149 C F_ao - heat flux from atmosphere to ocean (W/m^2)
150 _RL F_ao (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151
152 C F_mi - heat flux from mixed layer to ice (W/m^2)
153 _RL F_mi (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154
155 c the theta to use for the calculation of mixed layer-> ice heat fluxes
156 _RL surf_theta
157
158 _RL FLUX_TO_DELTA_TEMP,ENERGY_TO_DELTA_TEMP
159
160 if ( buoyancyRelation .eq. 'OCEANICP' ) then
161 kSurface = Nr
162 else
163 kSurface = 1
164 endif
165
166 FLUX_TO_DELTA_TEMP = SEAICE_deltaTtherm*
167 & recip_Cp*recip_rhoConst * recip_drF(1)
168
169 ENERGY_TO_DELTA_TEMP = recip_Cp*recip_rhoConst*recip_drF(1)
170
171 C ICE SALINITY (g/kg)
172 salinity_ice = 4.0
173
174 C FREEZING TEMP. OF SEA WATER (deg C)
175 TBC = SEAICE_freeze
176
177 C FREEZING POINT OF FRESHWATER
178 TMELT = 273.15
179
180 C IAN
181
182 c Sea ice density (kg m^-3)
183 RHOI = 917.0
184
185 c Seawater density (kg m^-3)
186 RHOSW = 1026.0
187
188 c Freshwater density (KG M^-3)
189 RHOFW = 1000.0
190
191 C Snow density
192 RHOSN = SEAICE_rhoSnow
193
194 C Heat capacity of seawater (J m^-3 K^-1)
195 CPW = 4010.0
196
197 c latent heat of fusion for ice (J kg^-1)
198 LI = 3.340e5
199 c conversion between Joules and m^3 of ice (m^3)
200 QI = 1/rhoi/Li
201 QS = 1/RHOSN/Li
202
203 c FOR FLOODING
204 FL_C2 = RHOI/RHOSW
205 CMM4IF FL_C2 = RHOSN/RHOSW
206 FL_C3 = (RHOSW-RHOI)/RHOSN
207 FL_C4 = RHOSN/RHOI
208
209 c Timescale for melting of ice from a warm ML (3 days in seconds)
210 c Damping term for mixed layer heat to melt existing ice
211 GAMMA = dRf(1)/SEAICE_gamma_t
212
213 DO bj=myByLo(myThid),myByHi(myThid)
214 DO bi=myBxLo(myThid),myBxHi(myThid)
215 c
216 #ifdef ALLOW_AUTODIFF_TAMC
217 act1 = bi - myBxLo(myThid)
218 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
219 act2 = bj - myByLo(myThid)
220 max2 = myByHi(myThid) - myByLo(myThid) + 1
221 act3 = myThid - 1
222 max3 = nTx*nTy
223 act4 = ikey_dynamics - 1
224 iicekey = (act1 + 1) + act2*max1
225 & + act3*max1*max2
226 & + act4*max1*max2*max3
227 #endif /* ALLOW_AUTODIFF_TAMC */
228 C
229 C initialise a few fields
230 C
231 #ifdef ALLOW_AUTODIFF_TAMC
232 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
233 CADJ & key = iicekey, byte = isbyte
234 CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj,
235 CADJ & key = iicekey, byte = isbyte
236 CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj,
237 CADJ & key = iicekey, byte = isbyte
238 #endif /* ALLOW_AUTODIFF_TAMC */
239 DO J=1,sNy
240 DO I=1,sNx
241 F_ia_net (I,J) = 0.0
242 F_ia_net_before_snow(I,J) = 0.0
243 F_io_net (I,J) = 0.0
244
245 F_ia (I,J) = 0.0
246 F_ao (I,J) = 0.0
247 F_mi (I,J) = 0.0
248
249 QNETI(I,J) = 0.0
250 QSWO (I,J) = 0.0
251 QSWI (I,J) = 0.0
252
253 QSWO_BELOW_FIRST_LAYER (I,J) = 0.0
254 QSWO_IN_FIRST_LAYER (I,J) = 0.0
255
256 S_a (I,J) = 0.0
257 S_h (I,J) = 0.0
258
259 IceGrowthRateUnderExistingIce (I,J) = 0.0
260 IceGrowthRateFromSurface (I,J) = 0.0
261 NetExistingIceGrowthRate (I,J) = 0.0
262 S_a_from_IGROW (I,J) = 0.0
263
264 PredTempChange (I,J) = 0.0
265 PredTempChangeFromQSW (I,J) = 0.0
266 PredTempChangeFromOA_MQNET (I,J) = 0.0
267 PredTempChangeFromFIA (I,J) = 0.0
268 PredTempChangeFromF_IA_NET (I,J) = 0.0
269 PredTempChangeFromF_IO_NET (I,J) = 0.0
270 PredTempChangeFromNewIceVol (I,J) = 0.0
271
272 IceGrowthRateOpenWater (I,J) = 0.0
273 IceGrowthRateMixedLayer (I,J) = 0.0
274
275 ExpectedIceVolumeChange (I,J) = 0.0
276 ExpectedSnowVolumeChange (I,J) = 0.0
277 ActualNewTotalVolumeChange (I,J) = 0.0
278 ActualNewTotalSnowMelt (I,J) = 0.0
279
280 EnergyInNewTotalIceVolume (I,J) = 0.0
281 NetEnergyFluxOutOfSystem (I,J) = 0.0
282 ResidualHeatOutOfSystem (I,J) = 0.0
283 QSW_absorb_in_ML (I,J) = 0.0
284 QSW_absorb_below_ML (I,J) = 0.0
285
286 SnowAccRateOverIce (I,J) = 0.0
287 SmowAccOverIce (I,J) = 0.0
288 PrecipRateOverIceSurfaceToSea (I,J) = 0.0
289
290 PotSnowMeltRateFromSurf (I,J) = 0.0
291 PotSnowMeltFromSurf (I,J) = 0.0
292 SnowMeltFromSurface (I,J) = 0.0
293 SnowMeltRateFromSurface (I,J) = 0.0
294 SurfHeatFluxConvergToSnowMelt (I,J) = 0.0
295
296 FreshwaterContribFromSnowMelt (I,J) = 0.0
297 FreshwaterContribFromIce (I,J) = 0.0
298
299 c the post sea ice advection and diffusion ice state are in time level 1.
300 c move these to the time level 2 before thermo. after this routine
301 c the updated ice state will be in time level 1 again. (except for snow
302 c which does not have 3 time levels for some reason)
303 HEFFNm1(I,J,bi,bj) = HEFF(I,J,bi,bj)
304 AREANm1(I,J,bi,bj) = AREA(I,J,bi,bj)
305
306 ENDDO
307 ENDDO
308
309 #ifdef ALLOW_AUTODIFF_TAMC
310 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
311 CADJ & key = iicekey, byte = isbyte
312 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
313 CADJ & key = iicekey, byte = isbyte
314 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
315 CADJ & key = iicekey, byte = isbyte
316 CADJ STORE tice(:,:,bi,bj) = comlev1_bibj,
317 CADJ & key = iicekey, byte = isbyte
318 CADJ STORE precip(:,:,bi,bj) = comlev1_bibj,
319 CADJ & key = iicekey, byte = isbyte
320 #endif /* ALLOW_AUTODIFF_TAMC */
321
322 DO J=1,sNy
323 DO I=1,sNx
324 C WE HAVE TO BE CAREFUL HERE SINCE ADVECTION/DIFFUSION COULD HAVE
325 C MAKE EITHER (BUT NOT BOTH) HEFF OR AREA ZERO OR NEGATIVE
326 C HSNOW COULD ALSO BECOME NEGATIVE
327 HEFFNm1(I,J,bi,bj) = MAX(0. _d 0,HEFFNm1(I,J,bi,bj))
328 HSNOW(I,J,bi,bj) = MAX(0. _d 0,HSNOW(I,J,bi,bj) )
329 AREANm1(I,J,bi,bj) = MAX(0. _d 0,AREANm1(I,J,bi,bj))
330 cif this is hack to prevent negative precip. somehow negative precips
331 cif escapes my exf_checkrange hack
332 cph-checkthis
333 IF (PRECIP(I,J,bi,bj) .LT. 0.0 _d 0) THEN
334 PRECIP(I,J,bi,bj) = 0.0 _d 0
335 ENDIF
336 ENDDO
337 ENDDO
338
339 #ifdef ALLOW_AUTODIFF_TAMC
340 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
341 CADJ & key = iicekey, byte = isbyte
342 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
343 CADJ & key = iicekey, byte = isbyte
344 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
345 CADJ & key = iicekey, byte = isbyte
346 CADJ STORE precip(:,:,bi,bj) = comlev1_bibj,
347 CADJ & key = iicekey, byte = isbyte
348 #endif
349 DO J=1,sNy
350 DO I=1,sNx
351
352 IF (HEFFNm1(I,J,bi,bj) .EQ. 0.0) THEN
353 AREANm1(I,J,bi,bj) = 0.0 _d 0
354 HSNOW(I,J,bi,bj) = 0.0 _d 0
355 ENDIF
356
357 IF (AREANm1(I,J,bi,bj) .EQ. 0.0) THEN
358 HEFFNm1(I,J,bi,bj) = 0.0 _d 0
359 HSNOW(I,J,bi,bj) = 0.0 _d 0
360 ENDIF
361
362 C PROCEED ONLY IF WE ARE CERTAIN TO HAVE ICE (AREA > 0)
363
364 IF (AREANm1(I,J,bi,bj) .GT. 0.) THEN
365 HICE_ACTUAL(I,J) =
366 & HEFFNm1(I,J,bi,bj)/AREANm1(I,J,bi,bj)
367
368 HSNOW_ACTUAL(I,J) = HSNOW(I,J,bi,bj)/
369 & AREANm1(I,J,bi,bj)
370
371 c ACCUMULATE SNOW
372 c Is the ice/surface below freezing or at the freezing
373 c point (melting). If it is freezing the precip is
374 c felt as snow and will accumulate over the ice. Else,
375 c precip makes its way, like all things in time, to the sea.
376 IF (TICE(I,J,bi,bj) .LT. TMELT) THEN
377 c Snow falls onto freezing surface remaining as snow
378 SnowAccRateOverIce(I,J) =
379 & PRECIP(I,J,bi,bj)*RHOFW/RHOSN
380
381 c None of the precipitation falls into the sea
382 PrecipRateOverIceSurfaceToSea(I,J) = 0.0
383
384 ELSE
385 c The snow melts on impact is is considered
386 c nothing more than rain. Since meltponds are
387 c not explicitly represented,this rain runs
388 c immediately into the sea
389
390 SnowAccRateOverIce(I,J) = 0.0
391 C The rate of rainfall over melting ice.
392 PrecipRateOverIceSurfaceToSea(I,J)=
393 & PRECIP(I,J,bi,bj)
394 ENDIF
395
396 c In m of mean snow thickness.
397 SmowAccOverIce(I,J) =
398 & SnowAccRateOverIce(I,J)
399 & *SEAICE_deltaTtherm*AreaNm1(I,J,bi,bj)
400
401 ELSE
402 HEFFNm1(I,J,bi,bj) = 0.0
403 HICE_ACTUAL(I,J) = 0.0
404 HSNOW_ACTUAL(I,J) = 0.0
405 HSNOW(I,J,bi,bj) = 0.0
406 ENDIF
407 HSNOW_ORIG(I,J) = HSNOW(I,J,bi,bj)
408 ENDDO
409 ENDDO
410
411 C FIND ATM. WIND SPEED
412 DO J=1,sNy
413 DO I=1,sNx
414 #ifdef SEAICE_EXTERNAL_FORCING
415 c USE EXF PACKAGE
416 UG(I,J) = MAX(SEAICE_EPS,wspeed(I,J,bi,bj))
417 #else
418 C CALCULATE IT HERE
419 SPEED_SQ = UWIND(I,J,bi,bj)**2 + VWIND(I,J,bi,bj)**2
420 IF ( SPEED_SQ .LE. SEAICE_EPS_SQ ) THEN
421 UG(I,J)=SEAICE_EPS
422 ELSE
423 UG(I,J)=SQRT(SPEED_SQ)
424 ENDIF
425 #endif /* SEAICE_EXTERNAL_FORCING */
426 ENDDO
427 ENDDO
428
429 #ifdef ALLOW_AUTODIFF_TAMC
430 cphCADJ STORE heff = comlev1, key = ikey_dynamics
431 cphCADJ STORE hsnow = comlev1, key = ikey_dynamics
432 cphCADJ STORE uwind = comlev1, key = ikey_dynamics
433 cphCADJ STORE vwind = comlev1, key = ikey_dynamics
434 CADJ STORE tice = comlev1, key = ikey_dynamics
435 #endif /* ALLOW_AUTODIFF_TAMC */
436
437 C SET LAYER TEMPERATURE IN KELVIN
438 DO J=1,sNy
439 DO I=1,sNx
440 TMIX(I,J,bi,bj)=
441 & theta(I,J,kSurface,bi,bj) + TMELT
442 ENDDO
443 ENDDO
444
445 C NOW DO ICE
446
447 CALL SEAICE_BUDGET_ICE_IF(
448 I UG, HICE_ACTUAL, HSNOW_ACTUAL,
449 U TICE,
450 O F_io_net,F_ia_net,F_ia, QSWI,
451 I bi, bj)
452
453 C Sometimes it's nice to have a setup without ice-atmosphere heat
454 C fluxes. This flag turns those fluxes to zero but leaves the
455 C Ice ocean fluxes intact. Thus, the first oceanic cell can transfer
456 C heat to the ice leading to melting in F_ml and it can release
457 C heat to the atmosphere through leads and open area thus growing it in
458 C F_ao
459
460 #ifdef FORBID_ICE_SURFACE_ATMOSPHERE_HEAT_FLUXES
461 DO J=1,sNy
462 DO I=1,sNx
463 F_ia_net (I,J) = 0.0
464 F_ia (I,J) = 0.0
465 F_io_net(I,J) = 0.0
466 ENDDO
467 ENDDO
468 #endif
469
470 C-- NET HEAT FLUX TO ICE FROM MIXED LAYER (POSITIVE MEANS NET OUT)
471 DO J=1,sNy
472 DO I=1,sNx
473
474 #ifdef SEAICE_DEBUG
475 IF ( (I .EQ. SEAICE_debugPointX) .and.
476 & (J .EQ. SEAICE_debugPointY) ) THEN
477
478 print *,'sig: I,J,F_ia,F_ia_net',I,J,F_ia(I,J),
479 & F_ia_net(I,J)
480
481 ENDIF
482 #endif
483
484 F_ia_net_before_snow(I,J) = F_ia_net(I,J)
485
486 IF (AreaNm1(I,J,bi,bj)*HEFFNm1(I,J,bi,bj).LE.0.) THEN
487 IceGrowthRateUnderExistingIce(I,J) = 0.0
488 IceGrowthRateFromSurface(I,J) = 0.0
489 NetExistingIceGrowthRate(I,J) = 0.0
490 ELSE
491 c The growth rate under existing ice is given by the upward
492 c ocean-ice conductive flux, F_io_net, and QI, which converts
493 c Joules to meters of ice. This quantity has units of meters
494 c of ice per second.
495 IceGrowthRateUnderExistingIce(I,J)=F_io_net(I,J)*QI
496
497 c Snow/Ice surface heat convergence is first used to melt
498 c snow. If all of this heat convergence went into melting
499 c snow, this is the rate at which it would do it
500 c F_ia_net must be negative, -> PSMRFW is positive for melting
501 PotSnowMeltRateFromSurf(I,J)= - F_ia_net(I,J)*QS
502
503 c This is the depth of snow that would be melted at this rate
504 c and the seaice delta t. In meters of snow.
505 PotSnowMeltFromSurf(I,J) =
506 & PotSnowMeltRateFromSurf(I,J)* SEAICE_deltaTtherm
507
508 c If we can melt MORE than is actually there, then we will
509 c reduce the melt rate so that only that which is there
510 c is melted in one time step. In this case not all of the
511 c heat flux convergence at the surface is used to melt snow,
512 c The leftover energy is going to melt ice.
513 c SurfHeatFluxConvergToSnowMelt is the part of the total heat
514 c flux convergence going to melt snow.
515
516 IF (PotSnowMeltFromSurf(I,J) .GE.
517 & HSNOW_ACTUAL(I,J)) THEN
518 c Snow melt and melt rate in actual snow thickness.
519 SnowMeltFromSurface(I,J) = HSNOW_ACTUAL(I,J)
520
521 SnowMeltRateFromSurface(I,J) =
522 & SnowMeltFromSurface(I,J)/ SEAICE_deltaTtherm
523
524 c Since F_ia_net is focused only over ice, its reduction
525 c requires knowing how much snow is actually melted
526 SurfHeatFluxConvergToSnowMelt(I,J) =
527 & -HSNOW_ACTUAL(I,J)/QS/SEAICE_deltaTtherm
528 ELSE
529 c In this case there will be snow remaining after melting.
530 c All of the surface heat convergence will be redirected to
531 c this effort.
532 SnowMeltFromSurface(I,J)=PotSnowMeltFromSurf(I,J)
533
534 SnowMeltRateFromSurface(I,J) =
535 & PotSnowMeltRateFromSurf(I,J)
536
537 SurfHeatFluxConvergToSnowMelt(I,J) =F_ia_net(I,J)
538 ENDIF
539
540 c Reduce the heat flux convergence available to melt surface
541 c ice by the amount used to melt snow
542 F_ia_net(I,J) =
543 & F_ia_net(I,J)-SurfHeatFluxConvergToSnowMelt(I,J)
544
545 IceGrowthRateFromSurface(I,J) = F_ia_net(I,J)*QI
546
547 NetExistingIceGrowthRate(I,J) =
548 & IceGrowthRateUnderExistingIce(I,J) +
549 & IceGrowthRateFromSurface(I,J)
550 ENDIF
551 ENDDO
552 ENDDO
553
554 c HERE WE WILL MELT SNOW AND ADJUST NET EXISTING ICE GROWTH RATE
555 C TO REFLECT REDUCTION IN SEA ICE MELT.
556
557 C NOW DETERMINE GROWTH RATES
558 C FIRST DO OPEN WATER
559 CALL SEAICE_BUDGET_OCEAN_IF(
560 I UG,
561 U TMIX,
562 O F_ao, QSWO,
563 I bi, bj, myThid )
564
565 #ifdef SEAICE_DEBUG
566 print *,'myiter', myIter
567 print '(A,2i4,2(1x,1P2E15.3))',
568 & 'ifice sigr, dbgx,dby, (netHF, SWHeatFlux)',
569 & SEAICE_debugPointX, SEAICE_debugPointY,
570 & F_ao(SEAICE_debugPointX, SEAICE_debugPointY),
571 & QSWO(SEAICE_debugPointX, SEAICE_debugPointY)
572 #endif
573
574
575 C-- NET HEAT FLUX TO ICE FROM MIXED LAYER (POSITIVE MEANS NET OUT)
576 c-- not all of the sw radiation is absorbed in the first layer, only that
577 c which is absorbed melts ice. SWFRACB is calculated in seaice_init_vari.F
578 DO J=1,sNy
579 DO I=1,sNx
580
581 c The contribution of shortwave heating is
582 c not included without SHORTWAVE_HEATING
583 #ifdef SHORTWAVE_HEATING
584 QSWO_BELOW_FIRST_LAYER(i,j)= QSWO(I,J)*SWFRACB
585 QSWO_IN_FIRST_LAYER(I,J) = QSWO(I,J)*(1.0 - SWFRACB)
586 #else
587 QSWO_BELOW_FIRST_LAYER(i,j)= 0. _d 0
588 QSWO_IN_FIRST_LAYER(I,J) = 0. _d 0
589 #endif
590 IceGrowthRateOpenWater(I,J)= QI*
591 & (F_ao(I,J) - QSWO(I,J) + QSWO_IN_FIRST_LAYER(I,J))
592
593 ENDDO
594 ENDDO
595
596
597 #ifdef ALLOW_AUTODIFF_TAMC
598 CADJ STORE theta(:,:,:,bi,bj)= comlev1_bibj,
599 CADJ & key = iicekey, byte = isbyte
600 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
601 CADJ & key = iicekey, byte = isbyte
602 #endif /* ALLOW_AUTODIFF_TAMC */
603
604
605 C-- NET HEAT FLUX TO ICE FROM MIXED LAYER (POSITIVE MEANS FLUX INTO ICE
606 C AND MELTING)
607 DO J=1,sNy
608 DO I=1,sNx
609
610 C FIND THE FREEZING POINT OF SEAWATER IN THIS CELL
611 #ifdef SEAICE_VARIABLE_FREEZING_POINT
612 TBC = -0.0575 _d 0*salt(I,J,kSurface,bi,bj) +
613 & 0.0901 _d 0
614 #endif /* SEAICE_VARIABLE_FREEZING_POINT */
615
616 c example: theta(i,j,ksurf) = 0, tbc = -2,
617 c fmi = -gamm*rhocpw * (0-(-2)) = - 2 * gamm * rhocpw,
618 c a NEGATIVE number. Heat flux INTO ice.
619
620 c It is fantastic that the model frequently generates thetas less
621 c then the freezing point. Just fantastic. When this happens,
622 c throw your hands up into the air, shut off the mixed layer
623 c heat flux, and hope for the best.
624 surf_theta = max(theta(I,J,kSurface,bi,bj), TBC)
625
626 F_mi(I,J) = -GAMMA*RHOSW*CPW *(surf_theta - TBC)
627
628 IceGrowthRateMixedLayer(I,J) = F_mi(I,J)*QI
629 ENDDO
630 ENDDO
631
632 #ifdef ALLOW_AUTODIFF_TAMC
633 CADJ STORE S_h(:,:) = comlev1_bibj,
634 CADJ & key = iicekey, byte = isbyte
635 #endif /* ALLOW_AUTODIFF_TAMC */
636
637 C CALCULATE THICKNESS DERIVATIVE (S_h)
638 DO J=1,sNy
639 DO I=1,sNx
640 S_h(I,J) =
641 & NetExistingIceGrowthRate(I,J)*AREANm1(I,J,bi,bj)+
642 & (1. -AREANm1(I,J,bi,bj))*
643 & IceGrowthRateOpenWater(I,J) +
644 & IceGrowthRateMixedLayer(I,J)
645
646 c Both the accumulation and melt rates are in terms
647 c of actual snow thickness. As with ice, multiplying
648 c with area converts to mean snow thickness.
649 S_hsnow(I,J) = AREANm1(I,J,bi,bj)* (
650 & SnowAccRateOverIce(I,J) -
651 & SnowMeltRateFromSurface(I,J) )
652 ENDDO
653 ENDDO
654
655 #ifdef ALLOW_AUTODIFF_TAMC
656 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
657 CADJ & key = iicekey, byte = isbyte
658 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
659 CADJ & key = iicekey, byte = isbyte
660 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
661 CADJ & key = iicekey, byte = isbyte
662 CADJ STORE S_h(:,:) = comlev1_bibj,
663 CADJ & key = iicekey, byte = isbyte
664 CADJ STORE S_hsnow(:,:) = comlev1_bibj,
665 CADJ & key = iicekey, byte = isbyte
666 #endif /* ALLOW_AUTODIFF_TAMC */
667
668 DO J=1,sNy
669 DO I=1,sNx
670 S_a(I,J) = 0.0
671 C IF THE OPEN WATER GROWTH RATE IS POSITIVE
672 C THEN EXTEND ICE AREAL COVER, S_a > 0
673
674 C TWO CASES, IF THERE IS ALREADY ICE PRESENT THEN EXTEND THE AREA USING THE
675 C OPEN WATER GROWTH RATE. IF THERE IS NO ICE PRESENT DO NOT EXTEND THE ICE
676 C UNTIL THE NET ICE THICKNESS RATE IS POSITIVE. I.E. IF THE MIXED LAYER
677 C HEAT FLUX INTO THE NEW ICE IS ENOUGH TO IMMEDIATELY MELT IT, DO NOT GROW
678 C IT.
679 IF (IceGrowthRateOpenWater(I,J) .GT. 0) THEN
680 IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
681 S_a_from_IGROW(I,J) = (ONE-AREANm1(I,J,bi,bj))*
682 & IceGrowthRateOpenWater(I,J)/HO_south
683 ELSE
684 S_a_from_IGROW(I,J) = (ONE-AREANm1(I,J,bi,bj))*
685 & IceGrowthRateOpenWater(I,J)/HO
686 ENDIF
687
688 IF (AREANm1(I,J,bi,bj) .GT. 0.) THEN
689 S_a(I,J) = S_a(I,J) + S_a_from_IGROW(I,J)
690 ELSE
691 IF (S_h(I,J) .GT. 0) THEN
692 S_a(I,J) = S_a(I,J) + S_a_from_IGROW(I,J)
693 ENDIF
694 ENDIF
695 ENDIF
696
697 C REDUCE THE ICE COVER IF ICE IS PRESENT
698 IF ( (S_h(I,J) .LT. 0.) .AND.
699 & (AREANm1(I,J,bi,bj).GT. 0.) .AND.
700 & (HEFFNm1(I,J,bi,bj).NE. 0.) ) THEN
701
702 S_a(I,J) = S_a(I,J)
703 & + AREANm1(I,J,bi,bj)/(2.0*HEFFNm1(I,J,bi,bj))*
704 & IceGrowthRateOpenWater(I,J)*
705 & (1-AREANm1(I,J,bi,bj))
706 ELSE
707 S_a(I,J) = S_a(I,J) + 0.0
708 ENDIF
709
710 C REDUCE THE ICE COVER IF ICE IS PRESENT
711 IF ( (IceGrowthRateMixedLayer(I,J) .LE. 0.) .AND.
712 & (AREANm1(I,J,bi,bj).GT. 0.) .AND.
713 & (HEFFNm1(I,J,bi,bj).NE. 0.) ) THEN
714
715 S_a(I,J) = S_a(I,J)
716 & + AREANm1(I,J,bi,bj)/(2.0*HEFFNm1(I,J,bi,bj))*
717 & IceGrowthRateMixedLayer(I,J)
718
719 ELSE
720 S_a(I,J) = S_a(I,J) + 0.0
721 ENDIF
722
723 C REDUCE THE ICE COVER IF ICE IS PRESENT
724 IF ( (NetExistingIceGrowthRate(I,J) .LE. 0.) .AND.
725 & (AREANm1(I,J,bi,bj).GT. 0.) .AND.
726 & (HEFFNm1(I,J,bi,bj).NE. 0.) ) THEN
727
728 S_a(I,J) = S_a(I,J)
729 & + AREANm1(I,J,bi,bj)/(2.0*HEFFNm1(I,J,bi,bj))*
730 & NetExistingIceGrowthRate(I,J)*AREANm1(I,J,bi,bj)
731
732 ELSE
733 S_a(I,J) = S_a(I,J) + 0.0
734 ENDIF
735
736 ENDDO
737 ENDDO
738
739
740 #ifdef ALLOW_AUTODIFF_TAMC
741 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
742 CADJ & key = iicekey, byte = isbyte
743 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
744 CADJ & key = iicekey, byte = isbyte
745 CADJ STORE S_a(:,:) = comlev1_bibj,
746 CADJ & key = iicekey, byte = isbyte
747 CADJ STORE S_h(:,:) = comlev1_bibj,
748 CADJ & key = iicekey, byte = isbyte
749 CADJ STORE f_ao(:,:) = comlev1_bibj,
750 CADJ & key = iicekey, byte = isbyte
751 CADJ STORE qswi(:,:) = comlev1_bibj,
752 CADJ & key = iicekey, byte = isbyte
753 CADJ STORE qswo(:,:) = comlev1_bibj,
754 CADJ & key = iicekey, byte = isbyte
755 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
756 CADJ & key = iicekey, byte = isbyte
757 #endif
758
759 C ACTUALLY CHANGE THE AREA AND THICKNESS
760 DO J=1,sNy
761 DO I=1,sNx
762 AREA(I,J,bi,bj) = AREANm1(I,J,bi,bj) +
763 & SEAICE_deltaTtherm * S_a(I,J)
764 HEFF(I,J,bi,bj) = HEFFNm1(I,J,bi,bj) +
765 & SEAICE_deltaTTherm * S_h(I,J)
766 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) +
767 & SEAICE_deltaTTherm * S_hsnow(I,J)
768 ENDDO
769 ENDDO
770
771 #ifdef ALLOW_AUTODIFF_TAMC
772 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
773 CADJ & key = iicekey, byte = isbyte
774 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
775 CADJ & key = iicekey, byte = isbyte
776 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
777 CADJ & key = iicekey, byte = isbyte
778 #endif
779
780 DO J=1,sNy
781 DO I=1,sNx
782 C SET LIMIT ON AREA etc.
783 AREA(I,J,bi,bj) = MIN(1. _d 0,AREA(I,J,bi,bj))
784 AREA(I,J,bi,bj) = MAX(0. _d 0,AREA(I,J,bi,bj))
785 HEFF(I,J,bi,bj) = MAX(0. _d 0, HEFF(I,J,bi,bj))
786 HSNOW(I,J,bi,bj) = MAX(0. _d 0, HSNOW(I,J,bi,bj))
787 ENDDO
788 ENDDO
789
790 #ifdef ALLOW_AUTODIFF_TAMC
791 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
792 CADJ & key = iicekey, byte = isbyte
793 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
794 CADJ & key = iicekey, byte = isbyte
795 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
796 CADJ & key = iicekey, byte = isbyte
797 #endif
798
799 DO J=1,sNy
800 DO I=1,sNx
801 IF (AREA(I,J,bi,bj) .GT. 0.0) THEN
802 HICE_ACTUAL(I,J) =
803 & HEFF(I,J,bi,bj)/AREA(I,J,bi,bj)
804 HSNOW_ACTUAL(I,J) =
805 & HSNOW(I,J,bi,bj)/AREA(I,J,bi,bj)
806 ELSE
807 HICE_ACTUAL(I,J) = 0.0
808 HSNOW_ACTUAL(I,J) = 0.0
809 ENDIF
810 ENDDO
811 ENDDO
812
813 #ifdef ALLOW_AUTODIFF_TAMC
814 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
815 CADJ & key = iicekey, byte = isbyte
816 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
817 CADJ & key = iicekey, byte = isbyte
818 #endif /* ALLOW_AUTODIFF_TAMC */
819
820 c constrain area is no thickness and vice versa.
821 DO J=1,sNy
822 DO I=1,sNx
823 IF (HEFF(I,J,bi,bj) .LE. 0.0 .OR.
824 & AREA(I,J,bi,bj) .LE. 0.0) THEN
825
826 AREA(I,J,bi,bj) = 0.0
827 HEFF(I,J,bi,bj) = 0.0
828 HICE_ACTUAL(I,J) = 0.0
829 HSNOW(I,J,bi,bj) = 0.0
830 HSNOW_ACTUAL(I,J) = 0.0
831 ENDIF
832 ENDDO
833 ENDDO
834
835 #ifdef ALLOW_AUTODIFF_TAMC
836 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
837 CADJ & key = iicekey, byte = isbyte
838 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
839 CADJ & key = iicekey, byte = isbyte
840 #endif /* ALLOW_AUTODIFF_TAMC */
841
842 DO J=1,sNy
843 DO I=1,sNx
844
845 c The amount of new mean thickness we expect to grow
846 ExpectedIceVolumeChange(I,J) = S_h(I,J) *
847 & SEAICE_deltaTtherm
848
849 ExpectedSnowVolumeChange(I,J) = S_hsnow(I,J)*
850 & SEAICE_deltaTtherm
851
852 c THE EFFECTIVE SHORTWAVE HEATING RATE
853 #ifdef SHORTWAVE_HEATING
854 QSW(I,J,bi,bj) =
855 & QSWI(I,J) * ( AREANm1(I,J,bi,bj)) +
856 & QSWO(I,J) * (1. - AREANm1(I,J,bi,bj))
857 #else
858 QSW(I,J,bi,bj) = 0. _d 0
859 #endif
860
861 ActualNewTotalVolumeChange(I,J) =
862 & HEFF(I,J,bi,bj) - HEFFNm1(I,J,bi,bj)
863
864 c The net average snow thickness melt that is actually realized. e.g.
865 c hsnow_orig = 0.25 m (e.g. 1 m of ice over a cell 1/4 covered in snow)
866 c hsnow_new = 0.20 m
867 c snow accum = 0.05 m
868 c melt = 0.25 + 0.05 - 0.2 = 0.1 m
869
870 c since this is in mean snow thickness it might have been 0.4 of actual
871 c snow thickness over the 1/4 of the cell which is ice covered.
872 ActualNewTotalSnowMelt(I,J) =
873 & HSNOW_ORIG(I,J) +
874 & SmowAccOverIce(I,J) -
875 & HSNOW(I,J,bi,bj)
876
877 c The latent heat of fusion of the new ice
878 EnergyInNewTotalIceVolume(I,J) =
879 & ActualNewTotalVolumeChange(I,J)/QI
880
881 c This is the net energy flux out of the ice+ocean system
882 c Remember -----
883 c F_ia_net : 0 if under freezing conditions (F_c < 0)
884 c The sum of the non-conductive surfice ice fluxes otherwise
885 c
886 c F_io_net : The conductive fluxes under freezing conditions (F_c < 0)
887 c 0 under melting conditions (no energy flux from ice to
888 c ocean)
889 c
890 c So if we are freezing, F_io_net is the conductive flux and there
891 c is energy balance at ice surface, F_ia_net =0. If we are melting
892 c There is a convergence of energy into the ice from above
893 NetEnergyFluxOutOfSystem(I,J) = SEAICE_deltaTtherm *
894 & (AREANm1(I,J,bi,bj) *
895 & (F_ia_net(I,J) + F_io_net(I,J) + QSWI(I,J))
896 & + (1.0 - AREANm1(I,J,bi,bj)) *
897 & F_ao(I,J))
898
899 c THE QUANTITY OF HEAT WHICH IS THE RESIDUAL TO THE QUANTITY OF
900 c ML temperature. If the net energy flux is exactly balanced by the
901 c latent energy of fusion in the new ice created then we won't
902 c change the ML temperature at all.
903
904 ResidualHeatOutOfSystem(I,J) =
905 & NetEnergyFluxOutOfSystem(I,J) -
906 & EnergyInNewTotalIceVolume(I,J)
907
908 C NOW FORMULATE QNET, which time LEVEL, ORIG 2.
909 C THIS QNET WILL DETERMINE THE TEMPERATURE CHANGE OF THE MIXED LAYER
910 C QNET IS A DEPTH AVERAGED HEAT FLUX FOR THE OCEAN COLUMN
911 C BECAUSE OF THE
912 QNET(I,J,bi,bj) =
913 & ResidualHeatOutOfSystem(I,J) / SEAICE_deltaTtherm
914
915
916 c Like snow melt, if there is melting, this quantity is positive.
917 c The change of freshwater content is per unit area over the entire
918 c cell, not just over the ice covered bits.
919 FreshwaterContribFromIce(I,J) =
920 & -ActualNewTotalVolumeChange(I,J)*RHOI/RHOFW
921
922 c The freshwater contribution from snow comes only in the form of melt
923 c unlike ice, which takes freshwater upon growth and yields freshwater
924 c upon melt. This is why the the actual new average snow melt was determined.
925 c In m/m^2 over the entire cell.
926 FreshwaterContribFromSnowMelt(I,J) =
927 & ActualNewTotalSnowMelt(I,J)*RHOSN/RHOFW
928
929 c This seems to be in m/s, original time level 2 for area
930 c Only the precip and evap need to be area weighted. The runoff
931 c and freshwater contribs from ice and snow melt are already mean
932 c weighted
933 EmPmR(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
934 & ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
935 & * ( ONE - AREANm1(I,J,bi,bj) )
936 & - PrecipRateOverIceSurfaceToSea(I,J)*
937 & AREANm1(I,J,bi,bj)
938 #ifdef ALLOW_RUNOFF
939 & - RUNOFF(I,J,bi,bj)
940 #endif
941 & - (FreshwaterContribFromIce(I,J) +
942 & FreshwaterContribFromSnowMelt(I,J))/
943 & SEAICE_deltaTtherm )*rhoConstFresh
944
945 C DO SOME DEBUGGING CALCULATIONS. MAKE SURE SUMS ALL ADD UP.
946 #ifdef SEAICE_DEBUG
947
948 C THE SHORTWAVE ENERGY FLUX ABSORBED IN THE SURFACE LAYER
949 #ifdef SHORTWAVE_HEATING
950 QSW_absorb_in_ML(I,J) = QSW(I,J,bi,bj)*
951 & (1.0 - SWFRACB)
952 #else
953
954 QSW_absorb_in_ML(I,J) = 0. _d 0
955 #endif
956
957 C THE SHORTWAVE ENERGY FLUX PENETRATING BELOW THE SURFACE LAYER
958 QSW_absorb_below_ML(I,J) =
959 & QSW(I,J,bi,bj) - QSW_absorb_in_ML(I,J);
960
961 PredTempChangeFromQSW(I,J) =
962 & - QSW_absorb_in_ML(I,J) * FLUX_TO_DELTA_TEMP
963
964 PredTempChangeFromOA_MQNET(I,J) =
965 & -(QNET(I,J,bi,bj)-QSWO(I,J))*(1. -AREANm1(I,J,bi,bj))
966 & * FLUX_TO_DELTA_TEMP
967
968 PredTempChangeFromF_IO_NET(I,J) =
969 & -F_io_net(I,J)*AREANm1(I,J,bi,bj)*FLUX_TO_DELTA_TEMP
970
971 PredTempChangeFromF_IA_NET(I,J) =
972 & -F_ia_net(I,J)*AREANm1(I,J,bi,bj)*FLUX_TO_DELTA_TEMP
973
974 PredTempChangeFromNewIceVol(I,J) =
975 & EnergyInNewTotalIceVolume(I,J)*ENERGY_TO_DELTA_TEMP
976
977 PredTempChange(I,J) =
978 & PredTempChangeFromQSW(I,J) +
979 & PredTempChangeFromOA_MQNET(I,J) +
980 & PredTempChangeFromF_IO_NET(I,J) +
981 & PredTempChangeFromF_IA_NET(I,J) +
982 & PredTempChangeFromNewIceVol(I,J)
983 #endif
984
985 ENDDO
986 ENDDO
987
988
989 #ifdef ALLOW_AUTODIFF_TAMC
990 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
991 CADJ & key = iicekey, byte = isbyte
992 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
993 CADJ & key = iicekey, byte = isbyte
994 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
995 CADJ & key = iicekey, byte = isbyte
996 #endif /* ALLOW_AUTODIFF_TAMC */
997
998 DO J=1,sNy
999 DO I=1,sNx
1000 AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1001 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1002 HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
1003 ENDDO
1004 ENDDO
1005
1006 #ifdef ALLOW_AUTODIFF_TAMC
1007 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
1008 CADJ & key = iicekey, byte = isbyte
1009 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
1010 CADJ & key = iicekey, byte = isbyte
1011 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
1012 CADJ & key = iicekey, byte = isbyte
1013 #endif /* ALLOW_AUTODIFF_TAMC */
1014
1015
1016 #ifdef ALLOW_SEAICE_FLOODING
1017 IF(SEAICEuseFlooding) THEN
1018
1019 DO J = 1,sNy
1020 DO I = 1,sNx
1021 EnergyToMeltSnowAndIce(I,J) =
1022 & HEFF(I,J,bi,bj)/QI +
1023 & HSNOW(I,J,bi,bj)/QS
1024
1025 deltaHS = FL_C2*( HSNOW_ACTUAL(I,J) -
1026 & HICE_ACTUAL(I,J)*FL_C3 )
1027
1028 IF (deltaHS .GT. 0.0) THEN
1029 deltaHI = FL_C4*deltaHS
1030
1031 HICE_ACTUAL(I,J) = HICE_ACTUAL(I,J)
1032 & + deltaHI
1033
1034 HSNOW_ACTUAL(I,J)= HSNOW_ACTUAL(I,J)
1035 & - deltaHS
1036
1037 HEFF(I,J,bi,bj)= HICE_ACTUAL(I,J) *
1038 & AREA(I,J,bi,bj)
1039
1040 HSNOW(I,J,bi,bj) = HSNOW_ACTUAL(I,J)*
1041 & AREA(I,J,bi,bj)
1042
1043 EnergyToMeltSnowAndIce2(I,J) =
1044 & HEFF(I,J,bi,bj)/QI +
1045 & HSNOW(I,J,bi,bj)/QS
1046
1047 #ifdef SEAICE_DEBUG
1048 IF ( (I .EQ. SEAICE_debugPointX) .and.
1049 & (J .EQ. SEAICE_debugPointY) ) THEN
1050
1051 print *,'Energy to melt snow+ice: pre,post,delta',
1052 & EnergyToMeltSnowAndIce(I,J),
1053 & EnergyToMeltSnowAndIce2(I,J),
1054 & EnergyToMeltSnowAndIce(I,J) -
1055 & EnergyToMeltSnowAndIce2(I,J)
1056 ENDIF
1057 c SEAICE DEBUG
1058 #endif
1059 c there is any flooding to be had
1060 ENDIF
1061 ENDDO
1062 ENDDO
1063
1064 c SEAICEuseFlooding
1065 ENDIF
1066 c ALLOW_SEAICE_FLOODING
1067 #endif
1068
1069 #ifdef ALLOW_AUTODIFF_TAMC
1070 CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
1071 CADJ & key = iicekey, byte = isbyte
1072 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
1073 CADJ & key = iicekey, byte = isbyte
1074 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
1075 CADJ & key = iicekey, byte = isbyte
1076 #endif /* ALLOW_AUTODIFF_TAMC */
1077
1078
1079 #ifdef ATMOSPHERIC_LOADING
1080 IF ( useRealFreshWaterFlux ) THEN
1081 DO J=1,sNy
1082 DO I=1,sNx
1083 sIceLoad(i,j,bi,bj) = HEFF(I,J,bi,bj)*
1084 & SEAICE_rhoIce + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
1085 ENDDO
1086 ENDDO
1087 ENDIF
1088 #endif
1089
1090 #ifdef SEAICE_DEBUG
1091 DO j=1-OLy,sNy+OLy
1092 DO i=1-OLx,sNx+OLx
1093
1094 IF ( (i .EQ. SEAICE_debugPointX) .and.
1095 & (j .EQ. SEAICE_debugPointY) ) THEN
1096
1097 print *,'ifsig: myTime,myIter:',myTime,myIter
1098
1099 print '(A,2i4,2(1x,1P3E15.4))',
1100 & 'ifice i j -------------- ',i,j
1101
1102 print '(A,2i4,2(1x,1P3E15.4))',
1103 & 'ifice i j IGR(ML OW ICE) ',i,j,
1104 & IceGrowthRateMixedLayer(i,j),
1105 & IceGrowthRateOpenWater(i,j),
1106 & NetExistingIceGrowthRate(i,j),
1107 & SEAICE_deltaTtherm
1108
1109 print '(A,2i4,2(1x,1P3E15.4))',
1110 & 'ifice i j F(mi ao) ',
1111 & i,j,F_mi(i,j), F_ao(i,j)
1112
1113 print '(A,2i4,2(1x,1P3E15.4))',
1114 & 'ifice i j Fi(a,ant2/1 ont)',
1115 & i,j,F_ia(i,j),
1116 & F_ia_net_before_snow(i,j),
1117 & F_ia_net(i,j),
1118 & F_io_net(i,j)
1119
1120 print '(A,2i4,2(1x,1P3E15.4))',
1121 & 'ifice i j AREA2/1 HEFF2/1 ',i,j,
1122 & AREANm1(I,J,bi,bj),
1123 & AREA(i,j,bi,bj),
1124 & HEFFNm1(I,J,bi,bj),
1125 & HEFF(i,j,bi,bj)
1126
1127 print '(A,2i4,2(1x,1P3E15.4))',
1128 & 'ifice i j HSNOW2/1 TMX TBC',i,j,
1129 & HSNOW_ORIG(I,J),
1130 & HSNOW(I,J,bi,bj),
1131 & TMIX(i,j,bi,bj)- TMELT,
1132 & TBC
1133
1134 print '(A,2i4,2(1x,1P3E15.4))',
1135 & 'ifice i j TI ATP LWD ',i,j,
1136 & TICE(i,j,bi,bj) - TMELT,
1137 & ATEMP(i,j,bi,bj) -TMELT,
1138 & LWDOWN(i,j,bi,bj)
1139
1140
1141 print '(A,2i4,2(1x,1P3E15.4))',
1142 & 'ifice i j S_a S_h S_hsnow ',i,j,
1143 & S_a(i,j),
1144 & S_h(i,j),
1145 & S_hsnow(i,j)
1146
1147 print '(A,2i4,2(1x,1P3E15.4))',
1148 & 'ifice i j IVC(E A ENIN) ',i,j,
1149 & ExpectedIceVolumeChange(i,j),
1150 & ActualNewTotalVolumeChange(i,j),
1151 & EnergyInNewTotalIceVolume(i,j)
1152
1153 print '(A,2i4,2(1x,1P3E15.4))',
1154 & 'ifice i j EF(NOS RE) QNET ',i,j,
1155 & NetEnergyFluxOutOfSystem(i,j),
1156 & ResidualHeatOutOfSystem(i,j),
1157 & QNET(I,J,bi,bj)
1158
1159 print '(A,2i4,3(1x,1P3E15.4))',
1160 & 'ifice i j QSW QSWO QSWI ',i,j,
1161 & QSW(i,j,bi,bj),
1162 & QSWO(i,j),
1163 & QSWI(i,j)
1164
1165 print '(A,2i4,3(1x,1P3E15.4))',
1166 & 'ifice i j SW(BML IML SW) ',i,j,
1167 & QSW_absorb_below_ML(i,j),
1168 & QSW_absorb_in_ML(i,j),
1169 & SWFRACB
1170
1171 print '(A,2i4,3(1x,1P3E15.4))',
1172 & 'ifice i j ptc(to, qsw, oa)',i,j,
1173 & PredTempChange(i,j),
1174 & PredTempChangeFromQSW (i,j),
1175 & PredTempChangeFromOA_MQNET(i,j)
1176
1177
1178 print '(A,2i4,3(1x,1P3E15.4))',
1179 & 'ifice i j ptc(fion,ian,ia)',i,j,
1180 & PredTempChangeFromF_IO_NET(i,j),
1181 & PredTempChangeFromF_IA_NET(i,j),
1182 & PredTempChangeFromFIA(i,j)
1183
1184 print '(A,2i4,3(1x,1P3E15.4))',
1185 & 'ifice i j ptc(niv) ',i,j,
1186 & PredTempChangeFromNewIceVol(i,j)
1187
1188
1189 print '(A,2i4,3(1x,1P3E15.4))',
1190 & 'ifice i j EmPmR EVP PRE RU',i,j,
1191 & EmPmR(I,J,bi,bj),
1192 & EVAP(I,J,bi,bj),
1193 & PRECIP(I,J,bi,bj),
1194 & RUNOFF(I,J,bi,bj)
1195
1196 print '(A,2i4,3(1x,1P3E15.4))',
1197 & 'ifice i j PRROIS,SAOI(R .)',i,j,
1198 & PrecipRateOverIceSurfaceToSea(I,J),
1199 & SnowAccRateOverIce(I,J),
1200 & SmowAccOverIce(I,J)
1201
1202 print '(A,2i4,4(1x,1P3E15.4))',
1203 & 'ifice i j SM(PM PMR . .R) ',i,j,
1204 & PotSnowMeltFromSurf(I,J),
1205 & PotSnowMeltRateFromSurf(I,J),
1206 & SnowMeltFromSurface(I,J),
1207 & SnowMeltRateFromSurface(I,J)
1208
1209 print '(A,2i4,4(1x,1P3E15.4))',
1210 & 'ifice i j TotSnwMlt ExSnVC',i,j,
1211 & ActualNewTotalSnowMelt(I,J),
1212 & ExpectedSnowVolumeChange(I,J)
1213
1214
1215 print '(A,2i4,4(1x,1P3E15.4))',
1216 & 'ifice i j fw(CFICE, CFSM) ',i,j,
1217 & FreshwaterContribFromIce(I,J),
1218 & FreshwaterContribFromSnowMelt(I,J)
1219
1220 print '(A,2i4,2(1x,1P3E15.4))',
1221 & 'ifice i j -------------- ',i,j
1222
1223 ENDIF
1224 ENDDO
1225 ENDDO
1226 #endif /* SEAICE_DEBUG */
1227
1228
1229 C BI,BJ'S
1230 ENDDO
1231 ENDDO
1232
1233 RETURN
1234 END

  ViewVC Help
Powered by ViewVC 1.1.22