2 |
C $Name$ |
C $Name$ |
3 |
|
|
4 |
#include "THSICE_OPTIONS.h" |
#include "THSICE_OPTIONS.h" |
5 |
|
#ifdef ALLOW_GENERIC_ADVDIFF |
6 |
|
# include "GAD_OPTIONS.h" |
7 |
|
#endif |
8 |
|
|
9 |
CBOP |
CBOP |
10 |
C !ROUTINE: THSICE_ADVDIFF |
C !ROUTINE: THSICE_ADVDIFF |
35 |
#include "EEPARAMS.h" |
#include "EEPARAMS.h" |
36 |
#include "PARAMS.h" |
#include "PARAMS.h" |
37 |
#include "GRID.h" |
#include "GRID.h" |
|
#include "GAD.h" |
|
38 |
#include "THSICE_SIZE.h" |
#include "THSICE_SIZE.h" |
39 |
#include "THSICE_PARAMS.h" |
#include "THSICE_PARAMS.h" |
40 |
#include "THSICE_VARS.h" |
#include "THSICE_VARS.h" |
41 |
#include "THSICE_2DYN.h" |
#include "THSICE_2DYN.h" |
42 |
|
#ifdef ALLOW_GENERIC_ADVDIFF |
43 |
|
# include "GAD.h" |
44 |
|
#endif |
45 |
|
#ifdef ALLOW_AUTODIFF_TAMC |
46 |
|
# include "tamc.h" |
47 |
|
# include "tamc_keys.h" |
48 |
|
#endif |
49 |
|
|
50 |
C !INPUT PARAMETERS: =================================================== |
C !INPUT PARAMETERS: =================================================== |
51 |
C === Routine arguments === |
C === Routine arguments === |
73 |
C afx :: horizontal advective flux, x direction |
C afx :: horizontal advective flux, x direction |
74 |
C afy :: horizontal advective flux, y direction |
C afy :: horizontal advective flux, y direction |
75 |
C iceFrc :: (new) sea-ice fraction |
C iceFrc :: (new) sea-ice fraction |
|
C iceFld :: (new) effective sea-ice thickness |
|
76 |
C iceVol :: temporary array used in advection S/R |
C iceVol :: temporary array used in advection S/R |
77 |
C oldVol :: (old) sea-ice volume |
C oldVol :: (old) sea-ice volume |
78 |
C msgBuf :: Informational/error meesage buffer |
C msgBuf :: Informational/error meesage buffer |
88 |
_RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
89 |
_RS maskOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RS maskOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
90 |
_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) |
|
91 |
_RL iceVol (1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
_RL iceVol (1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
92 |
_RL oldVol (1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
_RL oldVol (1-Olx:sNx+Olx,1-Oly:sNy+Oly) |
93 |
_RL minIcHeff, minIcArea, r_minArea |
_RL r_minArea |
94 |
_RL meanCellArea, areaEpsil, vol_Epsil |
_RL meanCellArea, areaEpsil, vol_Epsil |
95 |
#ifdef ALLOW_DIAGNOSTICS |
#ifdef ALLOW_DIAGNOSTICS |
96 |
CHARACTER*8 diagName |
CHARACTER*8 diagName |
97 |
CHARACTER*4 THSICE_DIAG_SUFX, diagSufx |
CHARACTER*4 THSICE_DIAG_SUFX, diagSufx |
98 |
EXTERNAL THSICE_DIAG_SUFX |
EXTERNAL THSICE_DIAG_SUFX |
99 |
|
LOGICAL DIAGNOSTICS_IS_ON |
100 |
|
EXTERNAL DIAGNOSTICS_IS_ON |
101 |
|
_RL tmpFac |
102 |
#endif |
#endif |
103 |
#ifdef ALLOW_DBUG_THSICE |
#ifdef ALLOW_DBUG_THSICE |
104 |
_RL tmpVar, sumVar1, sumVar2 |
_RL tmpVar, sumVar1, sumVar2 |
108 |
CEOP |
CEOP |
109 |
|
|
110 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
111 |
|
|
112 |
|
#ifdef ALLOW_AUTODIFF_TAMC |
113 |
|
act1 = bi - myBxLo(myThid) |
114 |
|
max1 = myBxHi(myThid) - myBxLo(myThid) + 1 |
115 |
|
act2 = bj - myByLo(myThid) |
116 |
|
max2 = myByHi(myThid) - myByLo(myThid) + 1 |
117 |
|
act3 = myThid - 1 |
118 |
|
max3 = nTx*nTy |
119 |
|
act4 = ikey_dynamics - 1 |
120 |
|
iicekey = (act1 + 1) + act2*max1 |
121 |
|
& + act3*max1*max2 |
122 |
|
& + act4*max1*max2*max3 |
123 |
|
#endif /* ALLOW_AUTODIFF_TAMC */ |
124 |
|
|
125 |
C areaEpsil, vol_Epsil are 2 small numbers for ice area & ice volume: |
C areaEpsil, vol_Epsil are 2 small numbers for ice area & ice volume: |
126 |
C if ice area (=ice fraction * grid-cell area) or ice volume (= effective |
C if ice area (=ice fraction * grid-cell area) or ice volume (= effective |
127 |
C thickness * grid-cell area) are too small (i.e.: < areaEpsil,vol_Epsil) |
C thickness * grid-cell area) are too small (i.e.: < areaEpsil,vol_Epsil) |
128 |
C will assume that ice is gone, and will loose mass or energy. |
C will assume that ice is gone, and will loose mass or energy. |
129 |
C However, if areaEpsil,vol_Epsil are much smaller than minimun ice area |
C However, if areaEpsil,vol_Epsil are much smaller than minimun ice area |
130 |
C (iceMaskMin*rAc) and minimum ice volume (iceMaskMin*himin*rAc), |
C (iceMaskMin*rAc) and minimum ice volume (iceMaskMin*hIceMin*rAc), |
131 |
C good chance that this will never happen within 1 time step. |
C good chance that this will never happen within 1 time step. |
132 |
|
|
133 |
dBugFlag = debugLevel.GE.debLevB |
dBugFlag = debugLevel.GE.debLevB |
139 |
areaEpsil = 1. _d -10 * meanCellArea |
areaEpsil = 1. _d -10 * meanCellArea |
140 |
vol_Epsil = 1. _d -15 * meanCellArea |
vol_Epsil = 1. _d -15 * meanCellArea |
141 |
|
|
|
minIcArea = iceMaskMin |
|
|
minIcHeff = iceMaskMin*himin |
|
142 |
r_minArea = 0. _d 0 |
r_minArea = 0. _d 0 |
143 |
IF ( minIcArea.GT.0. _d 0 ) r_minArea = 1. _d 0 / minIcArea |
IF ( iceMaskMin.GT.0. _d 0 ) r_minArea = 1. _d 0 / iceMaskMin |
144 |
|
|
145 |
thSIce_multiDimAdv = .TRUE. |
thSIce_multiDimAdv = .TRUE. |
146 |
|
#ifdef ALLOW_GENERIC_ADVDIFF |
147 |
IF ( thSIceAdvScheme.EQ.ENUM_CENTERED_2ND |
IF ( thSIceAdvScheme.EQ.ENUM_CENTERED_2ND |
148 |
& .OR.thSIceAdvScheme.EQ.ENUM_UPWIND_3RD |
& .OR.thSIceAdvScheme.EQ.ENUM_UPWIND_3RD |
149 |
& .OR.thSIceAdvScheme.EQ.ENUM_CENTERED_4TH ) THEN |
& .OR.thSIceAdvScheme.EQ.ENUM_CENTERED_4TH ) THEN |
150 |
thSIce_multiDimAdv = .FALSE. |
thSIce_multiDimAdv = .FALSE. |
151 |
ENDIF |
ENDIF |
152 |
|
#endif /* ALLOW_GENERIC_ADVDIFF */ |
153 |
|
|
154 |
C-- Initialisation (+ build oceanic mask) |
C-- Initialisation (+ build oceanic mask) |
155 |
DO j=1-OLy,sNy+OLy |
DO j=1-OLy,sNy+OLy |
167 |
ENDDO |
ENDDO |
168 |
ENDDO |
ENDDO |
169 |
|
|
170 |
|
#ifdef ALLOW_AUTODIFF_TAMC |
171 |
|
CADJ STORE iceHeight(i,j,bi,bj) |
172 |
|
CADJ & = comlev1_bibj, key=iicekey, byte=isbyte |
173 |
|
CADJ STORE snowHeight(i,j,bi,bj) |
174 |
|
CADJ & = comlev1_bibj, key=iicekey, byte=isbyte |
175 |
|
#endif |
176 |
IF ( thSIce_diffK .GT. 0. ) THEN |
IF ( thSIce_diffK .GT. 0. ) THEN |
177 |
CALL THSICE_DIFFUSION( |
CALL THSICE_DIFFUSION( |
178 |
I maskOce, |
I maskOce, |
202 |
iceFrc(i,j) = iceMask(i,j,bi,bj) |
iceFrc(i,j) = iceMask(i,j,bi,bj) |
203 |
ENDDO |
ENDDO |
204 |
ENDDO |
ENDDO |
205 |
|
#ifdef ALLOW_AUTODIFF_TAMC |
206 |
|
CADJ STORE icevol(i,j) = comlev1_bibj, key=iicekey, byte=isbyte |
207 |
|
CADJ STORE utrice(i,j) = comlev1_bibj, key=iicekey, byte=isbyte |
208 |
|
CADJ STORE vtrice(i,j) = comlev1_bibj, key=iicekey, byte=isbyte |
209 |
|
#endif |
210 |
CALL THSICE_ADVECTION( |
CALL THSICE_ADVECTION( |
211 |
I GAD_SI_FRAC, thSIceAdvScheme, .TRUE., |
I GAD_SI_FRAC, thSIceAdvScheme, .TRUE., |
212 |
I uTrIce, vTrIce, maskOce, thSIce_deltaT, areaEpsil, |
I uTrIce, vTrIce, maskOce, thSIce_deltaT, areaEpsil, |
227 |
O afx, afy, |
O afx, afy, |
228 |
I bi, bj, myTime, myIter, myThid ) |
I bi, bj, myTime, myIter, myThid ) |
229 |
|
|
230 |
|
#ifdef ALLOW_AUTODIFF_TAMC |
231 |
|
CADJ STORE iceHeight(i,j,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte |
232 |
|
CADJ STORE iceMask(i,j,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte |
233 |
|
#endif |
234 |
C-- sea-ice Thickness |
C-- sea-ice Thickness |
235 |
DO j=1-Oly,sNy+Oly |
DO j=1-Oly,sNy+Oly |
236 |
DO i=1-Olx,sNx+Olx |
DO i=1-Olx,sNx+Olx |
244 |
U iceVol, iceHeight(1-Olx,1-Oly,bi,bj), |
U iceVol, iceHeight(1-Olx,1-Oly,bi,bj), |
245 |
O uTrIce, vTrIce, |
O uTrIce, vTrIce, |
246 |
I bi, bj, myTime, myIter, myThid ) |
I bi, bj, myTime, myIter, myThid ) |
247 |
|
|
248 |
|
#ifdef ALLOW_AUTODIFF_TAMC |
249 |
|
CADJ STORE qice2(i,j,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte |
250 |
|
CADJ STORE utrice(i,j) = comlev1_bibj, key=iicekey, byte=isbyte |
251 |
|
CADJ STORE vtrice(i,j) = comlev1_bibj, key=iicekey, byte=isbyte |
252 |
|
#endif |
253 |
|
|
254 |
#ifdef ALLOW_DBUG_THSICE |
#ifdef ALLOW_DBUG_THSICE |
255 |
IF ( dBugFlag ) THEN |
IF ( dBugFlag ) THEN |
256 |
sumVar1 = 0. |
sumVar1 = 0. |
277 |
& 'ICE_ADV: uIce=', uIce(i,j), uIce(i+1,j), |
& 'ICE_ADV: uIce=', uIce(i,j), uIce(i+1,j), |
278 |
& ' , vIce=', vIce(i,j), vIce(i,j+1) |
& ' , vIce=', vIce(i,j), vIce(i,j+1) |
279 |
WRITE(6,'(2(A,1P2E14.6))') |
WRITE(6,'(2(A,1P2E14.6))') |
280 |
c & 'ICE_ADV: heff_b,a=', HEFF(i,j,2,bi,bj),HEFF(i,j,1,bi,bj) |
& 'ICE_ADV: area_b,a=', iceMask(i,j,bi,bj), iceFrc(i,j), |
281 |
c WRITE(6,'(A,1P4E14.6)') 'ICE_ADV: mFx=', gFld(i,j) |
& ' , Heff_b,a=', oldVol(i,j)*recip_rA(i,j,bi,bj), |
282 |
|
& iceHeight(i,j,bi,bj)*iceFrc(i,j) |
283 |
ENDIF |
ENDIF |
284 |
ENDDO |
ENDDO |
285 |
ENDDO |
ENDDO |
406 |
ENDIF |
ENDIF |
407 |
#endif |
#endif |
408 |
|
|
409 |
C-- Update Ice Fraction, Ice thickness and snow thickness: |
#ifdef ALLOW_AUTODIFF_TAMC |
410 |
C and adjust sea-ice state if not enough ice. |
CADJ STORE iceHeight(:,:,bi,bj) = |
411 |
|
CADJ & comlev1_bibj, key=iicekey, byte=isbyte |
412 |
|
CADJ STORE snowHeight(:,:,bi,bj) = |
413 |
|
CADJ & comlev1_bibj, key=iicekey, byte=isbyte |
414 |
|
CADJ STORE iceFrc(:,:) = |
415 |
|
CADJ & comlev1_bibj, key=iicekey, byte=isbyte |
416 |
|
#endif |
417 |
|
|
418 |
|
C-- Update Ice Fraction: ensure that fraction is > iceMaskMin & < 1 |
419 |
|
C and adjust Ice thickness and snow thickness accordingly |
420 |
DO j=1,sNy |
DO j=1,sNy |
421 |
DO i=1,sNx |
DO i=1,sNx |
422 |
C- store new effective ice-thickness |
IF ( iceFrc(i,j) .GT. 1. _d 0 ) THEN |
423 |
iceFld(i,j) = iceHeight(i,j,bi,bj)*iceFrc(i,j) |
c IF ( dBug(i,j,bi,bj) ) |
|
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 |
|
424 |
iceMask(i,j,bi,bj) = 1. _d 0 |
iceMask(i,j,bi,bj) = 1. _d 0 |
425 |
iceHeight(i,j,bi,bj) = iceFld(i,j) |
iceHeight(i,j,bi,bj) = iceHeight(i,j,bi,bj) *iceFrc(i,j) |
426 |
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) |
427 |
ELSEIF ( iceFrc(i,j) .LT. minIcArea ) THEN |
ELSEIF ( iceFrc(i,j) .LT. iceMaskMin ) THEN |
428 |
iceMask(i,j,bi,bj) = minIcArea |
c IF ( dBug(i,j,bi,bj) ) |
429 |
iceHeight(i,j,bi,bj) = iceFld(i,j)*r_minArea |
iceMask(i,j,bi,bj) = iceMaskMin |
430 |
|
iceHeight(i,j,bi,bj) = iceHeight(i,j,bi,bj) |
431 |
|
& *iceFrc(i,j)*r_minArea |
432 |
snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj) |
snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj) |
433 |
& *iceFrc(i,j)*r_minArea |
& *iceFrc(i,j)*r_minArea |
|
ELSE |
|
|
iceMask(i,j,bi,bj) = iceFrc(i,j) |
|
|
ENDIF |
|
434 |
ELSE |
ELSE |
435 |
|
iceMask(i,j,bi,bj) = iceFrc(i,j) |
436 |
|
ENDIF |
437 |
|
ENDDO |
438 |
|
ENDDO |
439 |
|
C- adjust sea-ice state if not enough ice. |
440 |
|
DO j=1,sNy |
441 |
|
DO i=1,sNx |
442 |
|
IF ( iceHeight(i,j,bi,bj).LT.hIceMin ) THEN |
443 |
C- Not enough ice, melt the tiny amount of snow & ice: |
C- Not enough ice, melt the tiny amount of snow & ice: |
444 |
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) |
445 |
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 |
446 |
C using seaice_advdiff (with SEAICEadvScheme=1) providing I comment |
C using seaice_advdiff (with SEAICEadvScheme=1) providing I comment |
447 |
C out the following lines (and then loose conservation). |
C out the following lines (and then loose conservation). |
448 |
C- - |
C- - |
449 |
oceFWfx(i,j,bi,bj) = ( rhos*snowHeight(i,j,bi,bj) |
oceFWfx(i,j,bi,bj) = ( rhos*snowHeight(i,j,bi,bj) |
450 |
& +rhoi*iceHeight(i,j,bi,bj) ) |
& +rhoi*iceHeight(i,j,bi,bj) |
451 |
& *iceFrc(i,j)/thSIce_deltaT |
& )*iceMask(i,j,bi,bj)/thSIce_deltaT |
452 |
oceSflx(i,j,bi,bj) =saltice*rhoi*iceFld(i,j)/thSIce_deltaT |
oceSflx(i,j,bi,bj) = rhoi*iceHeight(i,j,bi,bj)*saltIce |
453 |
oceQnet(i,j,bi,bj) = -qsnow*rhos*snowHeight(i,j,bi,bj) |
& *iceMask(i,j,bi,bj)/thSIce_deltaT |
454 |
& *iceFrc(i,j)/thSIce_deltaT |
oceQnet(i,j,bi,bj) = -( rhos*snowHeight(i,j,bi,bj)*qsnow |
455 |
& -( Qice1(i,j,bi,bj) |
& +rhoi*iceHeight(i,j,bi,bj) |
456 |
& +Qice2(i,j,bi,bj) )*0.5 _d 0 |
& *( Qice1(i,j,bi,bj) |
457 |
& *rhoi*iceFld(i,j)/thSIce_deltaT |
& +Qice2(i,j,bi,bj) )*0.5 _d 0 |
458 |
|
& )*iceMask(i,j,bi,bj)/thSIce_deltaT |
459 |
|
c IF ( dBug(i,j,bi,bj) ) |
460 |
C- - |
C- - |
461 |
c flx2oc (i,j) = flx2oc (i,j) + |
c flx2oc (i,j) = flx2oc (i,j) + |
462 |
c frw2oc (i,j) = frw2oc (i,j) + |
c frw2oc (i,j) = frw2oc (i,j) + |
499 |
C--- end if multiDimAdvection |
C--- end if multiDimAdvection |
500 |
ENDIF |
ENDIF |
501 |
|
|
502 |
|
#ifdef ALLOW_DIAGNOSTICS |
503 |
|
IF ( useDiagnostics ) THEN |
504 |
|
CALL DIAGNOSTICS_FILL(iceMask,'SI_AdvFr',0,1,1,bi,bj,myThid) |
505 |
|
C- Ice-fraction weighted quantities: |
506 |
|
tmpFac = 1. _d 0 |
507 |
|
CALL DIAGNOSTICS_FRACT_FILL( |
508 |
|
I iceHeight, iceMask,tmpFac,1,'SI_AdvHi', |
509 |
|
I 0,1,1,bi,bj,myThid) |
510 |
|
CALL DIAGNOSTICS_FRACT_FILL( |
511 |
|
I snowHeight,iceMask,tmpFac,1,'SI_AdvHs', |
512 |
|
I 0,1,1,bi,bj,myThid) |
513 |
|
C- Ice-Volume weighted quantities: |
514 |
|
IF ( DIAGNOSTICS_IS_ON('SI_AdvQ1',myThid) .OR. |
515 |
|
& DIAGNOSTICS_IS_ON('SI_AdvQ2',myThid) ) THEN |
516 |
|
DO j=1,sNy |
517 |
|
DO i=1,sNx |
518 |
|
iceVol(i,j) = iceMask(i,j,bi,bj)*iceHeight(i,j,bi,bj) |
519 |
|
ENDDO |
520 |
|
ENDDO |
521 |
|
CALL DIAGNOSTICS_FRACT_FILL( |
522 |
|
I Qice1(1-OLx,1-OLy,bi,bj), |
523 |
|
I iceVol,tmpFac,1,'SI_AdvQ1', |
524 |
|
I 0,1,2,bi,bj,myThid) |
525 |
|
CALL DIAGNOSTICS_FRACT_FILL( |
526 |
|
I Qice2(1-OLx,1-OLy,bi,bj), |
527 |
|
I iceVol,tmpFac,1,'SI_AdvQ2', |
528 |
|
I 0,1,2,bi,bj,myThid) |
529 |
|
ENDIF |
530 |
|
ENDIF |
531 |
|
#endif /* ALLOW_DIAGNOSTICS */ |
532 |
|
|
533 |
#endif /* ALLOW_THSICE */ |
#endif /* ALLOW_THSICE */ |
534 |
|
|
535 |
RETURN |
RETURN |