/[MITgcm]/MITgcm_contrib/dgoldberg/shelfice/shelfice_thermodynamics.F
ViewVC logotype

Diff of /MITgcm_contrib/dgoldberg/shelfice/shelfice_thermodynamics.F

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

revision 1.1 by dgoldberg, Wed Aug 27 19:26:17 2014 UTC revision 1.3 by dgoldberg, Tue Apr 21 11:07:10 2015 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "SHELFICE_OPTIONS.h"  #include "SHELFICE_OPTIONS.h"
5  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_AUTODIFF
6  # include "STREAMICE_OPTIONS.h"  # include "AUTODIFF_OPTIONS.h"
7    #endif
8    #ifdef ALLOW_CTRL
9    # include "CTRL_OPTIONS.h"
10  #endif  #endif
11    
12  CBOP  CBOP
# Line 21  C     | Line 24  C     |
24  C     | stresses at the ice/water interface are computed in separate  C     | stresses at the ice/water interface are computed in separate
25  C     | routines that are called from mom_fluxform/mom_vecinv  C     | routines that are called from mom_fluxform/mom_vecinv
26  C     *=============================================================*  C     *=============================================================*
27    C     \ev
28    
29  C     !USES:  C     !USES:
30        IMPLICIT NONE        IMPLICIT NONE
# Line 33  C     === Global variables === Line 37  C     === Global variables ===
37  #include "DYNVARS.h"  #include "DYNVARS.h"
38  #include "FFIELDS.h"  #include "FFIELDS.h"
39  #include "SHELFICE.h"  #include "SHELFICE.h"
40    #include "SHELFICE_COST.h"
41  #ifdef ALLOW_AUTODIFF  #ifdef ALLOW_AUTODIFF
42  # include "CTRL_SIZE.h"  # include "CTRL_SIZE.h"
43  # include "ctrl.h"  # include "ctrl.h"
# Line 44  C     === Global variables === Line 49  C     === Global variables ===
49  #  include "tamc_keys.h"  #  include "tamc_keys.h"
50  # endif /* SHI_ALLOW_GAMMAFRICT */  # endif /* SHI_ALLOW_GAMMAFRICT */
51  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 #ifdef ALLOW_STREAMICE  
 # include "STREAMICE.h"  
 #endif  
   
52    
53  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
54  C     === Routine arguments ===  C     === Routine arguments ===
# Line 57  C     myThid :: thread number for this i Line 58  C     myThid :: thread number for this i
58        _RL  myTime        _RL  myTime
59        INTEGER myIter        INTEGER myIter
60        INTEGER myThid        INTEGER myThid
 CEOP  
61    
62  #ifdef ALLOW_SHELFICE  #ifdef ALLOW_SHELFICE
63  C     !LOCAL VARIABLES :  C     !LOCAL VARIABLES :
# Line 72  C                         (negative dens Line 72  C                         (negative dens
72  C     convertFW2SaltLoc:: local copy of convertFW2Salt  C     convertFW2SaltLoc:: local copy of convertFW2Salt
73  C     cFac             :: 1 for conservative form, 0, otherwise  C     cFac             :: 1 for conservative form, 0, otherwise
74  C     rFac             :: realFreshWaterFlux factor  C     rFac             :: realFreshWaterFlux factor
75  C     dFac             :: 0 for diffusive heat flux (Holland and Jenkins, 1999, eq21)  C     dFac             :: 0 for diffusive heat flux (Holland and Jenkins, 1999,
76    C                           eq21)
77  C                         1 for advective and diffusive heat flux (eq22, 26, 31)  C                         1 for advective and diffusive heat flux (eq22, 26, 31)
78  C     fwflxFac         :: only effective for dFac=1, 1 if we expect a melting fresh  C     fwflxFac         :: only effective for dFac=1, 1 if we expect a melting
79  C                         water flux, 0 otherwise  C                         fresh water flux, 0 otherwise
80  C     auxiliary variables and abbreviations:  C     auxiliary variables and abbreviations:
81  C     a0, a1, a2, b, c0  C     a0, a1, a2, b, c0
82  C     eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8  C     eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8
# Line 95  C     drKp1, recip_drLoc Line 96  C     drKp1, recip_drLoc
96        _RL cFac, rFac, dFac, fwflxFac        _RL cFac, rFac, dFac, fwflxFac
97        _RL aqe, bqe, cqe, discrim, recip_aqe        _RL aqe, bqe, cqe, discrim, recip_aqe
98        _RL drKp1, recip_drLoc        _RL drKp1, recip_drLoc
99          _RL recip_latentHeat
100        _RL tmpFac        _RL tmpFac
101    
102  #ifdef SHI_ALLOW_GAMMAFRICT  #ifdef SHI_ALLOW_GAMMAFRICT
# Line 102  C     drKp1, recip_drLoc Line 104  C     drKp1, recip_drLoc
104        _RL gammaTmoleT, gammaTmoleS, gammaTurb, gammaTurbConst        _RL gammaTmoleT, gammaTmoleS, gammaTurb, gammaTurbConst
105        _RL ustar, ustarSq, etastar        _RL ustar, ustarSq, etastar
106        PARAMETER ( shiTwoThirds = 0.66666666666666666666666666667D0 )        PARAMETER ( shiTwoThirds = 0.66666666666666666666666666667D0 )
107    #ifdef ALLOW_DIAGNOSTICS
108          _RL uStarDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
109    #endif /* ALLOW_DIAGNOSTICS */
110  #endif  #endif
111    
112  #ifndef ALLOW_OPENAD  #ifndef ALLOW_OPENAD
# Line 112  C     drKp1, recip_drLoc Line 117  C     drKp1, recip_drLoc
117  #ifdef ALLOW_SHIFWFLX_CONTROL  #ifdef ALLOW_SHIFWFLX_CONTROL
118        _RL xx_shifwflx_loc(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL xx_shifwflx_loc(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
119  #endif  #endif
120    CEOP
121  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
122    
123  #ifdef SHI_ALLOW_GAMMAFRICT  #ifdef SHI_ALLOW_GAMMAFRICT
124  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
125  C     re-initialize here again, curtesy to TAF  C     re-initialize here again, curtesy to TAF
126        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
127         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 127  C     re-initialize here again, curtesy Line 133  C     re-initialize here again, curtesy
133          ENDDO          ENDDO
134         ENDDO         ENDDO
135        ENDDO        ENDDO
136  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
137        IF ( SHELFICEuseGammaFrict ) THEN        IF ( SHELFICEuseGammaFrict ) THEN
138  C     Implement friction velocity-dependent transfer coefficient  C     Implement friction velocity-dependent transfer coefficient
139  C     of Holland and Jenkins, JPO, 1999  C     of Holland and Jenkins, JPO, 1999
# Line 143  C     instead of etastar = sqrt(1+zetaN* Line 149  C     instead of etastar = sqrt(1+zetaN*
149         etastar = 1. _d 0         etastar = 1. _d 0
150         gammaTurbConst  = 1. _d 0 / (2. _d 0 * shiZetaN*etastar)         gammaTurbConst  = 1. _d 0 / (2. _d 0 * shiZetaN*etastar)
151       &      - recip_shiKarman       &      - recip_shiKarman
152  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
153         DO bj = myByLo(myThid), myByHi(myThid)         DO bj = myByLo(myThid), myByHi(myThid)
154          DO bi = myBxLo(myThid), myBxHi(myThid)          DO bi = myBxLo(myThid), myBxHi(myThid)
155           DO J = 1-OLy,sNy+OLy           DO J = 1-OLy,sNy+OLy
# Line 154  C     instead of etastar = sqrt(1+zetaN* Line 160  C     instead of etastar = sqrt(1+zetaN*
160           ENDDO           ENDDO
161          ENDDO          ENDDO
162         ENDDO         ENDDO
163  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
164        ENDIF        ENDIF
165  #endif /* SHI_ALLOW_GAMMAFRICT */  #endif /* SHI_ALLOW_GAMMAFRICT */
166    
167          recip_latentHeat = 0. _d 0
168          IF ( SHELFICElatentHeat .NE. 0. _d 0 )
169         &     recip_latentHeat = 1. _d 0/SHELFICElatentHeat
170  C     are we doing the conservative form of Jenkins et al. (2001)?  C     are we doing the conservative form of Jenkins et al. (2001)?
171        recip_Cp = 1. _d 0 / HeatCapacity_Cp        recip_Cp = 1. _d 0 / HeatCapacity_Cp
172        cFac = 0. _d 0        cFac = 0. _d 0
# Line 200  C     where this value is part of the pr Line 209  C     where this value is part of the pr
209            shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0            shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0
210            shelficeForcingT      (I,J,bi,bj) = 0. _d 0            shelficeForcingT      (I,J,bi,bj) = 0. _d 0
211            shelficeForcingS      (I,J,bi,bj) = 0. _d 0            shelficeForcingS      (I,J,bi,bj) = 0. _d 0
212    #if (defined SHI_ALLOW_GAMMAFRICT && defined ALLOW_DIAGNOSTICS)
213              uStarDiag             (I,J,bi,bj) = 0. _d 0
214    #endif /* SHI_ALLOW_GAMMAFRICT and ALLOW_DIAGNOSTICS */
215           ENDDO           ENDDO
216          ENDDO          ENDDO
217         ENDDO         ENDDO
# Line 214  C     where this value is part of the pr Line 226  C     where this value is part of the pr
226          ENDDO          ENDDO
227         ENDDO         ENDDO
228        ENDDO        ENDDO
229        CALL CTRL_GET_GEN (  #ifdef ALLOW_CTRL
230          if (useCTRL) CALL CTRL_GET_GEN (
231       &     xx_shifwflx_file, xx_shifwflxstartdate, xx_shifwflxperiod,       &     xx_shifwflx_file, xx_shifwflxstartdate, xx_shifwflxperiod,
232       &     maskSHI, xx_shifwflx_loc, xx_shifwflx0, xx_shifwflx1,       &     maskSHI, xx_shifwflx_loc, xx_shifwflx0, xx_shifwflx1,
233       &     xx_shifwflx_dummy,       &     xx_shifwflx_dummy,
234       &     xx_shifwflx_remo_intercept, xx_shifwflx_remo_slope,       &     xx_shifwflx_remo_intercept, xx_shifwflx_remo_slope,
235         &     wshifwflx,
236       &     myTime, myIter, myThid )       &     myTime, myIter, myThid )
237    #endif
238  #endif /* ALLOW_SHIFWFLX_CONTROL */  #endif /* ALLOW_SHIFWFLX_CONTROL */
239        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
240         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 265  C--   overlap into lower cell Line 280  C--   overlap into lower cell
280              drKp1 = drF(K)*( 1. _d 0 - _hFacC(I,J,K,bi,bj) )              drKp1 = drF(K)*( 1. _d 0 - _hFacC(I,J,K,bi,bj) )
281  C--   lower cell may not be as thick as required  C--   lower cell may not be as thick as required
282              drKp1 = MIN( drKp1, drF(Kp1) * _hFacC(I,J,Kp1,bi,bj) )              drKp1 = MIN( drKp1, drF(Kp1) * _hFacC(I,J,Kp1,bi,bj) )
283                drKp1 = MAX( drKp1, 0. _d 0 )
284              recip_drLoc = 1. _d 0 /              recip_drLoc = 1. _d 0 /
285       &           ( drF(K)*_hFacC(I,J,K,bi,bj) + drKp1 )       &           ( drF(K)*_hFacC(I,J,K,bi,bj) + drKp1 )
286              tLoc(I,J) = ( tLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)              tLoc(I,J) = ( tLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj)
# Line 309  C--   to the surface Line 325  C--   to the surface
325              ustarSq = shiCdrag * MAX( 1.D-6,              ustarSq = shiCdrag * MAX( 1.D-6,
326       &           0.25 _d 0 *(uLoc(I,J)*uLoc(I,J)+vLoc(I,J)*vLoc(I,J)) )       &           0.25 _d 0 *(uLoc(I,J)*uLoc(I,J)+vLoc(I,J)*vLoc(I,J)) )
327              ustar   = SQRT(ustarSq)              ustar   = SQRT(ustarSq)
328    #ifdef ALLOW_DIAGNOSTICS
329                uStarDiag(I,J,bi,bj) = ustar
330    #endif /* ALLOW_DIAGNOSTICS */
331  C     instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc))  C     instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc))
332  C           etastar = 1. _d 0  C           etastar = 1. _d 0
333  C           gammaTurbConst  = 1. _d 0 / (2. _d 0 * shiZetaN*etastar)  C           gammaTurbConst  = 1. _d 0 / (2. _d 0 * shiZetaN*etastar)
# Line 365  C     thus a downward (negative upward) Line 384  C     thus a downward (negative upward)
384  C     and vice versa  C     and vice versa
385              shelfIceFreshWaterFlux(I,J,bi,bj) =              shelfIceFreshWaterFlux(I,J,bi,bj) =
386       &           - shelfIceHeatFlux(I,J,bi,bj)       &           - shelfIceHeatFlux(I,J,bi,bj)
387       &           *recip_SHELFICElatentHeat       &           *recip_latentHeat
388  C--   compute surface tendencies  C--   compute surface tendencies
389              shelficeForcingT(i,j,bi,bj) =              shelficeForcingT(i,j,bi,bj) =
390       &           - shelfIceHeatFlux(I,J,bi,bj)       &           - shelfIceHeatFlux(I,J,bi,bj)
# Line 517  C     endif (not) useISOMIPTD Line 536  C     endif (not) useISOMIPTD
536         ENDDO         ENDDO
537        ENDDO        ENDDO
538    
539        IF (SHELFICEallowThinIceMass) then        IF (SHELFICEMassStepping) THEN
 #ifdef ALLOW_STREAMICE  
        DO bj = myByLo(myThid), myByHi(myThid)  
         DO bi = myBxLo(myThid), myBxHi(myThid)  
          DO j=1-Oly,sNy+Oly-1  
           DO i=1-Olx+1,sNx+Olx-1  
   
            if (streamice_hmask(i,j,bi,bj).eq.1 .or.  
      &         streamice_hmask(i,j,bi,bj).eq.2) then  
   
             shelficeMass(i,j,bi,bj) =  
      &       H_streamice(I,J,bi,bj) * streamice_density  
   
            endif  
   
           ENDDO  
          ENDDO  
         ENDDO  
        ENDDO  
 #else  
540         CALL SHELFICE_STEP_ICEMASS( myTime, myIter, myThid )         CALL SHELFICE_STEP_ICEMASS( myTime, myIter, myThid )
 #endif  
541        ENDIF        ENDIF
542    
543  C--  Calculate new loading anomaly (in case the ice-shelf mass was updated)  C--  Calculate new loading anomaly (in case the ice-shelf mass was updated)
# Line 556  c     IF ( SHELFICEloadAnomalyFile .EQ. Line 555  c     IF ( SHELFICEloadAnomalyFile .EQ.
555         ENDDO         ENDDO
556  c     ENDIF  c     ENDIF
557  #endif /* ndef ALLOW_AUTODIFF */  #endif /* ndef ALLOW_AUTODIFF */
558        print *, "GOT HERE 0"  
559  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
560        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
         print *, "GOT HERE 1"  
561         CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx',         CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx',
562       &      0,1,0,1,1,myThid)       &      0,1,0,1,1,myThid)
563         CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux,      'SHIhtFlx',         CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux,      'SHIhtFlx',
# Line 577  C     Transfer coefficients Line 575  C     Transfer coefficients
575       &      0,1,0,1,1,myThid)       &      0,1,0,1,1,myThid)
576         CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS',         CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS',
577       &      0,1,0,1,1,myThid)       &      0,1,0,1,1,myThid)
578         print *, "GOT HERE 2"  C     Friction velocity
579    #ifdef SHI_ALLOW_GAMMAFRICT
580           IF ( SHELFICEuseGammaFrict )
581         &  CALL DIAGNOSTICS_FILL(uStarDiag,'SHIuStar',0,1,0,1,1,myThid)
582    #endif /* SHI_ALLOW_GAMMAFRICT */
583        ENDIF        ENDIF
584  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
585    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22