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 |