/[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.7 by jmc, Mon Aug 27 13:20:57 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
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
# Line 98  C     msgBuf    :: Informational/error m Line 108  C     msgBuf    :: Informational/error m
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
# Line 115  C     but what matter here is just to ha Line 139  C     but what matter here is just to ha
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
# Line 143  C--   Initialisation (+ build oceanic ma Line 167  C--   Initialisation (+ build oceanic ma
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,
# Line 172  C--   Fractional area Line 202  C--   Fractional area
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,
# Line 192  C--   Snow thickness Line 227  C--   Snow thickness
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
# Line 205  C--   sea-ice Thickness Line 244  C--   sea-ice Thickness
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.
# Line 231  C-      Check that updated iceVol = iceF Line 277  C-      Check that updated iceVol = iceF
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
# Line 359  c           WRITE(6,'(A,1P4E14.6)') 'ICE Line 406  c           WRITE(6,'(A,1P4E14.6)') 'ICE
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) +
# Line 437  C---  if not multiDimAdvection Line 499  C---  if not multiDimAdvection
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
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22