/[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.4 by jmc, Mon Apr 9 15:51:22 2007 UTC revision 1.7 by jmc, Mon Aug 27 13:20:57 2007 UTC
# Line 42  C   oceQnet   :: heat flux to the ocean Line 42  C   oceQnet   :: heat flux to the ocean
42  #ifdef ALLOW_GENERIC_ADVDIFF  #ifdef ALLOW_GENERIC_ADVDIFF
43  # include "GAD.h"  # include "GAD.h"
44  #endif  #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 92  C     msgBuf    :: Informational/error m Line 96  C     msgBuf    :: Informational/error m
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 101  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 146  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 175  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 195  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 208  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 234  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 362  c           WRITE(6,'(A,1P4E14.6)') 'ICE Line 406  c           WRITE(6,'(A,1P4E14.6)') 'ICE
406          ENDIF          ENDIF
407  #endif  #endif
408    
409    #ifdef ALLOW_AUTODIFF_TAMC
410    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  C--   Update Ice Fraction: ensure that fraction is > iceMaskMin & < 1
419  C      and adjust Ice thickness and snow thickness accordingly  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            IF ( iceFrc(i,j) .GT. 1. _d 0 ) THEN            IF ( iceFrc(i,j) .GT. 1. _d 0 ) THEN
423    c           IF ( dBug(i,j,bi,bj) )
424              iceMask(i,j,bi,bj)    = 1. _d 0              iceMask(i,j,bi,bj)    = 1. _d 0
425              iceHeight(i,j,bi,bj)  = iceHeight(i,j,bi,bj) *iceFrc(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. iceMaskMin ) THEN            ELSEIF ( iceFrc(i,j) .LT. iceMaskMin ) THEN
428    c           IF ( dBug(i,j,bi,bj) )
429              iceMask(i,j,bi,bj)    = iceMaskMin              iceMask(i,j,bi,bj)    = iceMaskMin
430              iceHeight(i,j,bi,bj)  = iceHeight(i,j,bi,bj)              iceHeight(i,j,bi,bj)  = iceHeight(i,j,bi,bj)
431       &                             *iceFrc(i,j)*r_minArea       &                             *iceFrc(i,j)*r_minArea
# Line 384  C      and adjust Ice thickness and snow Line 439  C      and adjust Ice thickness and snow
439  C-     adjust sea-ice state if not enough ice.  C-     adjust sea-ice state if not enough ice.
440          DO j=1,sNy          DO j=1,sNy
441           DO i=1,sNx           DO i=1,sNx
442            IF ( iceHeight(i,j,bi,bj).LT.himin ) THEN            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 fresh-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
# Line 394  C- - Line 449  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       &                            )*iceMask(i,j,bi,bj)/thSIce_deltaT       &                            )*iceMask(i,j,bi,bj)/thSIce_deltaT
452              oceSflx(i,j,bi,bj) =    rhoi*iceHeight(i,j,bi,bj)*saltice              oceSflx(i,j,bi,bj) =    rhoi*iceHeight(i,j,bi,bj)*saltIce
453       &                             *iceMask(i,j,bi,bj)/thSIce_deltaT       &                             *iceMask(i,j,bi,bj)/thSIce_deltaT
454              oceQnet(i,j,bi,bj) = -( rhos*snowHeight(i,j,bi,bj)*qsnow              oceQnet(i,j,bi,bj) = -( rhos*snowHeight(i,j,bi,bj)*qsnow
455       &                             +rhoi*iceHeight(i,j,bi,bj)       &                             +rhoi*iceHeight(i,j,bi,bj)
456       &                                  *( Qice1(i,j,bi,bj)       &                                  *( Qice1(i,j,bi,bj)
457       &                                    +Qice2(i,j,bi,bj) )*0.5 _d 0       &                                    +Qice2(i,j,bi,bj) )*0.5 _d 0
458       &                            )*iceMask(i,j,bi,bj)/thSIce_deltaT       &                            )*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 443  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.4  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22