/[MITgcm]/MITgcm/pkg/thsice/thsice_advdiff.F
ViewVC logotype

Diff of /MITgcm/pkg/thsice/thsice_advdiff.F

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

revision 1.2 by jmc, Thu Apr 5 22:06:43 2007 UTC revision 1.5 by heimbach, Mon Apr 16 22:38:24 2007 UTC
# Line 2  C $Header$ Line 2  C $Header$
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
# Line 32  C   oceQnet   :: heat flux to the ocean Line 35  C   oceQnet   :: heat flux to the ocean
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 ===
# Line 64  C     vTrIce    :: sea-ice volume transp Line 73  C     vTrIce    :: sea-ice volume transp
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
# Line 80  C     msgBuf    :: Informational/error m Line 88  C     msgBuf    :: Informational/error m
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
# Line 98  C     msgBuf    :: Informational/error m Line 105  C     msgBuf    :: Informational/error m
105  CEOP  CEOP
106    
107  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
108    
109    #ifdef ALLOW_AUTODIFF_TAMC
110              act1 = bi - myBxLo(myThid)
111              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
112              act2 = bj - myByLo(myThid)
113              max2 = myByHi(myThid) - myByLo(myThid) + 1
114              act3 = myThid - 1
115              max3 = nTx*nTy
116              act4 = ikey_dynamics - 1
117              iicekey = (act1 + 1) + act2*max1
118         &                         + act3*max1*max2
119         &                         + act4*max1*max2*max3
120    #endif /* ALLOW_AUTODIFF_TAMC */
121    
122  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:
123  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
124  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)
# Line 115  C     but what matter here is just to ha Line 136  C     but what matter here is just to ha
136        areaEpsil = 1. _d -10 * meanCellArea        areaEpsil = 1. _d -10 * meanCellArea
137        vol_Epsil = 1. _d -15 * meanCellArea        vol_Epsil = 1. _d -15 * meanCellArea
138    
       minIcArea = iceMaskMin  
       minIcHeff = iceMaskMin*himin  
139        r_minArea = 0. _d 0        r_minArea = 0. _d 0
140        IF ( minIcArea.GT.0. _d 0 ) r_minArea = 1. _d 0 / minIcArea        IF ( iceMaskMin.GT.0. _d 0 ) r_minArea = 1. _d 0 / iceMaskMin
141    
142        thSIce_multiDimAdv = .TRUE.        thSIce_multiDimAdv = .TRUE.
143    #ifdef ALLOW_GENERIC_ADVDIFF
144        IF ( thSIceAdvScheme.EQ.ENUM_CENTERED_2ND        IF ( thSIceAdvScheme.EQ.ENUM_CENTERED_2ND
145       & .OR.thSIceAdvScheme.EQ.ENUM_UPWIND_3RD       & .OR.thSIceAdvScheme.EQ.ENUM_UPWIND_3RD
146       & .OR.thSIceAdvScheme.EQ.ENUM_CENTERED_4TH ) THEN       & .OR.thSIceAdvScheme.EQ.ENUM_CENTERED_4TH ) THEN
147         thSIce_multiDimAdv = .FALSE.         thSIce_multiDimAdv = .FALSE.
148        ENDIF        ENDIF
149    #endif /* ALLOW_GENERIC_ADVDIFF */
150    
151  C--   Initialisation (+ build oceanic mask)  C--   Initialisation (+ build oceanic mask)
152        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
# Line 143  C--   Initialisation (+ build oceanic ma Line 164  C--   Initialisation (+ build oceanic ma
164         ENDDO         ENDDO
165        ENDDO        ENDDO
166    
167    #ifdef ALLOW_AUTODIFF_TAMC
168    CADJ STORE iceHeight(i,j,bi,bj)  
169    CADJ &     = comlev1_bibj, key=iicekey, byte=isbyte
170    CADJ STORE snowHeight(i,j,bi,bj)
171    CADJ &     = comlev1_bibj, key=iicekey, byte=isbyte
172    #endif
173        IF ( thSIce_diffK .GT. 0. ) THEN        IF ( thSIce_diffK .GT. 0. ) THEN
174          CALL THSICE_DIFFUSION(          CALL THSICE_DIFFUSION(
175       I              maskOce,       I              maskOce,
# Line 172  C--   Fractional area Line 199  C--   Fractional area
199            iceFrc(i,j) = iceMask(i,j,bi,bj)            iceFrc(i,j) = iceMask(i,j,bi,bj)
200           ENDDO           ENDDO
201          ENDDO          ENDDO
202    #ifdef ALLOW_AUTODIFF_TAMC
203    CADJ STORE icevol(i,j) = comlev1_bibj, key=iicekey, byte=isbyte
204    CADJ STORE utrice(i,j) = comlev1_bibj, key=iicekey, byte=isbyte
205    CADJ STORE vtrice(i,j) = comlev1_bibj, key=iicekey, byte=isbyte
206    #endif
207          CALL THSICE_ADVECTION(          CALL THSICE_ADVECTION(
208       I       GAD_SI_FRAC,  thSIceAdvScheme, .TRUE.,       I       GAD_SI_FRAC,  thSIceAdvScheme, .TRUE.,
209       I       uTrIce, vTrIce, maskOce, thSIce_deltaT, areaEpsil,       I       uTrIce, vTrIce, maskOce, thSIce_deltaT, areaEpsil,
# Line 192  C--   Snow thickness Line 224  C--   Snow thickness
224       O       afx, afy,       O       afx, afy,
225       I       bi, bj, myTime, myIter, myThid )       I       bi, bj, myTime, myIter, myThid )
226    
227    #ifdef ALLOW_AUTODIFF_TAMC
228    CADJ STORE iceHeight(i,j,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
229    CADJ STORE iceMask(i,j,bi,bj)   = comlev1_bibj, key=iicekey, byte=isbyte
230    #endif
231  C--   sea-ice Thickness  C--   sea-ice Thickness
232          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
233           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
# Line 205  C--   sea-ice Thickness Line 241  C--   sea-ice Thickness
241       U       iceVol, iceHeight(1-Olx,1-Oly,bi,bj),       U       iceVol, iceHeight(1-Olx,1-Oly,bi,bj),
242       O       uTrIce, vTrIce,       O       uTrIce, vTrIce,
243       I       bi, bj, myTime, myIter, myThid )       I       bi, bj, myTime, myIter, myThid )
244    
245    #ifdef ALLOW_AUTODIFF_TAMC
246    CADJ STORE qice2(i,j,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
247    CADJ STORE utrice(i,j)      = comlev1_bibj, key=iicekey, byte=isbyte
248    CADJ STORE vtrice(i,j)      = comlev1_bibj, key=iicekey, byte=isbyte
249    #endif
250    
251  #ifdef ALLOW_DBUG_THSICE  #ifdef ALLOW_DBUG_THSICE
252          IF ( dBugFlag ) THEN          IF ( dBugFlag ) THEN
253           sumVar1 = 0.           sumVar1 = 0.
# Line 359  c           WRITE(6,'(A,1P4E14.6)') 'ICE Line 402  c           WRITE(6,'(A,1P4E14.6)') 'ICE
402          ENDIF          ENDIF
403  #endif  #endif
404    
405  C--   Update Ice Fraction, Ice thickness and snow thickness:  #ifdef ALLOW_AUTODIFF_TAMC
406  C      and adjust sea-ice state if not enough ice.  CADJ STORE iceHeight(:,:,bi,bj) =
407    CADJ &     comlev1_bibj, key=iicekey, byte=isbyte
408    CADJ STORE snowHeight(:,:,bi,bj) =
409    CADJ &     comlev1_bibj, key=iicekey, byte=isbyte
410    CADJ STORE iceFrc(:,:) =
411    CADJ &     comlev1_bibj, key=iicekey, byte=isbyte
412    #endif
413    
414    C--   Update Ice Fraction: ensure that fraction is > iceMaskMin & < 1
415    C      and adjust Ice thickness and snow thickness accordingly
416          DO j=1,sNy          DO j=1,sNy
417           DO i=1,sNx           DO i=1,sNx
418  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  
419              iceMask(i,j,bi,bj)    = 1. _d 0              iceMask(i,j,bi,bj)    = 1. _d 0
420              iceHeight(i,j,bi,bj)  = iceFld(i,j)              iceHeight(i,j,bi,bj)  = iceHeight(i,j,bi,bj) *iceFrc(i,j)
421              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)
422             ELSEIF ( iceFrc(i,j) .LT. minIcArea ) THEN            ELSEIF ( iceFrc(i,j) .LT. iceMaskMin ) THEN
423              iceMask(i,j,bi,bj)    = minIcArea              iceMask(i,j,bi,bj)    = iceMaskMin
424              iceHeight(i,j,bi,bj)  = iceFld(i,j)*r_minArea              iceHeight(i,j,bi,bj)  = iceHeight(i,j,bi,bj)
425         &                             *iceFrc(i,j)*r_minArea
426              snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj)              snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj)
427       &                             *iceFrc(i,j)*r_minArea       &                             *iceFrc(i,j)*r_minArea
            ELSE  
             iceMask(i,j,bi,bj)    = iceFrc(i,j)  
            ENDIF  
428            ELSE            ELSE
429                iceMask(i,j,bi,bj)    = iceFrc(i,j)
430              ENDIF
431             ENDDO
432            ENDDO
433    C-     adjust sea-ice state if not enough ice.
434            DO j=1,sNy
435             DO i=1,sNx
436              IF ( iceHeight(i,j,bi,bj).LT.himin ) THEN
437  C-    Not enough ice, melt the tiny amount of snow & ice:  C-    Not enough ice, melt the tiny amount of snow & ice:
438  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)
439  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
440  C     using seaice_advdiff (with SEAICEadvScheme=1) providing I comment  C     using seaice_advdiff (with SEAICEadvScheme=1) providing I comment
441  C     out the following lines (and then loose conservation).  C     out the following lines (and then loose conservation).
442  C- -  C- -
443              oceFWfx(i,j,bi,bj) = ( rhos*snowHeight(i,j,bi,bj)              oceFWfx(i,j,bi,bj) =  ( rhos*snowHeight(i,j,bi,bj)
444       &                            +rhoi*iceHeight(i,j,bi,bj) )       &                             +rhoi*iceHeight(i,j,bi,bj)
445       &                                      *iceFrc(i,j)/thSIce_deltaT       &                            )*iceMask(i,j,bi,bj)/thSIce_deltaT
446              oceSflx(i,j,bi,bj) =saltice*rhoi*iceFld(i,j)/thSIce_deltaT              oceSflx(i,j,bi,bj) =    rhoi*iceHeight(i,j,bi,bj)*saltice
447              oceQnet(i,j,bi,bj) = -qsnow*rhos*snowHeight(i,j,bi,bj)       &                             *iceMask(i,j,bi,bj)/thSIce_deltaT
448       &                                      *iceFrc(i,j)/thSIce_deltaT              oceQnet(i,j,bi,bj) = -( rhos*snowHeight(i,j,bi,bj)*qsnow
449       &                         -( Qice1(i,j,bi,bj)       &                             +rhoi*iceHeight(i,j,bi,bj)
450       &                           +Qice2(i,j,bi,bj) )*0.5 _d 0       &                                  *( Qice1(i,j,bi,bj)
451       &                                 *rhoi*iceFld(i,j)/thSIce_deltaT       &                                    +Qice2(i,j,bi,bj) )*0.5 _d 0
452         &                            )*iceMask(i,j,bi,bj)/thSIce_deltaT
453  C- -  C- -
454  c           flx2oc (i,j) = flx2oc (i,j) +  c           flx2oc (i,j) = flx2oc (i,j) +
455  c           frw2oc (i,j) = frw2oc (i,j) +  c           frw2oc (i,j) = frw2oc (i,j) +
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22