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

Annotation 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 - (hide 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 mmazloff 1.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