69 |
C afx :: horizontal advective flux, x direction |
C afx :: horizontal advective flux, x direction |
70 |
C afy :: horizontal advective flux, y direction |
C afy :: horizontal advective flux, y direction |
71 |
C iceFrc :: (new) sea-ice fraction |
C iceFrc :: (new) sea-ice fraction |
|
C iceFld :: (new) effective sea-ice thickness |
|
72 |
C iceVol :: temporary array used in advection S/R |
C iceVol :: temporary array used in advection S/R |
73 |
C oldVol :: (old) sea-ice volume |
C oldVol :: (old) sea-ice volume |
74 |
C msgBuf :: Informational/error meesage buffer |
C msgBuf :: Informational/error meesage buffer |
84 |
_RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
85 |
_RS maskOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RS maskOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
86 |
_RL iceFrc (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL iceFrc (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
|
_RL iceFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
|
87 |
_RL iceVol (1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
_RL iceVol (1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
88 |
_RL oldVol (1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
_RL oldVol (1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
89 |
_RL minIcHeff, minIcArea, r_minArea |
_RL r_minArea |
90 |
_RL meanCellArea, areaEpsil, vol_Epsil |
_RL meanCellArea, areaEpsil, vol_Epsil |
91 |
#ifdef ALLOW_DIAGNOSTICS |
#ifdef ALLOW_DIAGNOSTICS |
92 |
CHARACTER*8 diagName |
CHARACTER*8 diagName |
118 |
areaEpsil = 1. _d -10 * meanCellArea |
areaEpsil = 1. _d -10 * meanCellArea |
119 |
vol_Epsil = 1. _d -15 * meanCellArea |
vol_Epsil = 1. _d -15 * meanCellArea |
120 |
|
|
|
minIcArea = iceMaskMin |
|
|
minIcHeff = iceMaskMin*himin |
|
121 |
r_minArea = 0. _d 0 |
r_minArea = 0. _d 0 |
122 |
IF ( minIcArea.GT.0. _d 0 ) r_minArea = 1. _d 0 / minIcArea |
IF ( iceMaskMin.GT.0. _d 0 ) r_minArea = 1. _d 0 / iceMaskMin |
123 |
|
|
124 |
thSIce_multiDimAdv = .TRUE. |
thSIce_multiDimAdv = .TRUE. |
125 |
#ifdef ALLOW_GENERIC_ADVDIFF |
#ifdef ALLOW_GENERIC_ADVDIFF |
362 |
ENDIF |
ENDIF |
363 |
#endif |
#endif |
364 |
|
|
365 |
C-- Update Ice Fraction, Ice thickness and snow thickness: |
C-- Update Ice Fraction: ensure that fraction is > iceMaskMin & < 1 |
366 |
C and adjust sea-ice state if not enough ice. |
C and adjust Ice thickness and snow thickness accordingly |
367 |
DO j=1,sNy |
DO j=1,sNy |
368 |
DO i=1,sNx |
DO i=1,sNx |
369 |
C- store new effective ice-thickness |
IF ( iceFrc(i,j) .GT. 1. _d 0 ) THEN |
|
iceFld(i,j) = iceHeight(i,j,bi,bj)*iceFrc(i,j) |
|
|
IF ( iceFld(i,j) .GE. minIcHeff ) THEN |
|
|
C- where there is enough ice, ensure that Ice fraction is > minIcArea & < 1 |
|
|
IF ( iceFrc(i,j) .GT. 1. _d 0 ) THEN |
|
370 |
iceMask(i,j,bi,bj) = 1. _d 0 |
iceMask(i,j,bi,bj) = 1. _d 0 |
371 |
iceHeight(i,j,bi,bj) = iceFld(i,j) |
iceHeight(i,j,bi,bj) = iceHeight(i,j,bi,bj) *iceFrc(i,j) |
372 |
snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj)*iceFrc(i,j) |
snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj)*iceFrc(i,j) |
373 |
ELSEIF ( iceFrc(i,j) .LT. minIcArea ) THEN |
ELSEIF ( iceFrc(i,j) .LT. iceMaskMin ) THEN |
374 |
iceMask(i,j,bi,bj) = minIcArea |
iceMask(i,j,bi,bj) = iceMaskMin |
375 |
iceHeight(i,j,bi,bj) = iceFld(i,j)*r_minArea |
iceHeight(i,j,bi,bj) = iceHeight(i,j,bi,bj) |
376 |
|
& *iceFrc(i,j)*r_minArea |
377 |
snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj) |
snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj) |
378 |
& *iceFrc(i,j)*r_minArea |
& *iceFrc(i,j)*r_minArea |
|
ELSE |
|
|
iceMask(i,j,bi,bj) = iceFrc(i,j) |
|
|
ENDIF |
|
379 |
ELSE |
ELSE |
380 |
|
iceMask(i,j,bi,bj) = iceFrc(i,j) |
381 |
|
ENDIF |
382 |
|
ENDDO |
383 |
|
ENDDO |
384 |
|
C- adjust sea-ice state if not enough ice. |
385 |
|
DO j=1,sNy |
386 |
|
DO i=1,sNx |
387 |
|
IF ( iceHeight(i,j,bi,bj).LT.himin ) THEN |
388 |
C- Not enough ice, melt the tiny amount of snow & ice: |
C- Not enough ice, melt the tiny amount of snow & ice: |
389 |
C and return frsh-water, salt & energy to the ocean (flx > 0 = into ocean) |
C and return fresh-water, salt & energy to the ocean (flx > 0 = into ocean) |
390 |
C- - Note: using 1rst.Order Upwind, I can get the same results as when |
C- - Note: using 1rst.Order Upwind, I can get the same results as when |
391 |
C using seaice_advdiff (with SEAICEadvScheme=1) providing I comment |
C using seaice_advdiff (with SEAICEadvScheme=1) providing I comment |
392 |
C out the following lines (and then loose conservation). |
C out the following lines (and then loose conservation). |
393 |
C- - |
C- - |
394 |
oceFWfx(i,j,bi,bj) = ( rhos*snowHeight(i,j,bi,bj) |
oceFWfx(i,j,bi,bj) = ( rhos*snowHeight(i,j,bi,bj) |
395 |
& +rhoi*iceHeight(i,j,bi,bj) ) |
& +rhoi*iceHeight(i,j,bi,bj) |
396 |
& *iceFrc(i,j)/thSIce_deltaT |
& )*iceMask(i,j,bi,bj)/thSIce_deltaT |
397 |
oceSflx(i,j,bi,bj) =saltice*rhoi*iceFld(i,j)/thSIce_deltaT |
oceSflx(i,j,bi,bj) = rhoi*iceHeight(i,j,bi,bj)*saltice |
398 |
oceQnet(i,j,bi,bj) = -qsnow*rhos*snowHeight(i,j,bi,bj) |
& *iceMask(i,j,bi,bj)/thSIce_deltaT |
399 |
& *iceFrc(i,j)/thSIce_deltaT |
oceQnet(i,j,bi,bj) = -( rhos*snowHeight(i,j,bi,bj)*qsnow |
400 |
& -( Qice1(i,j,bi,bj) |
& +rhoi*iceHeight(i,j,bi,bj) |
401 |
& +Qice2(i,j,bi,bj) )*0.5 _d 0 |
& *( Qice1(i,j,bi,bj) |
402 |
& *rhoi*iceFld(i,j)/thSIce_deltaT |
& +Qice2(i,j,bi,bj) )*0.5 _d 0 |
403 |
|
& )*iceMask(i,j,bi,bj)/thSIce_deltaT |
404 |
C- - |
C- - |
405 |
c flx2oc (i,j) = flx2oc (i,j) + |
c flx2oc (i,j) = flx2oc (i,j) + |
406 |
c frw2oc (i,j) = frw2oc (i,j) + |
c frw2oc (i,j) = frw2oc (i,j) + |