--- MITgcm_contrib/dgoldberg/shelfice/shelfice_thermodynamics.F 2014/08/27 19:26:17 1.1 +++ MITgcm_contrib/dgoldberg/shelfice/shelfice_thermodynamics.F 2015/04/21 11:07:10 1.3 @@ -1,9 +1,12 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/dgoldberg/shelfice/shelfice_thermodynamics.F,v 1.1 2014/08/27 19:26:17 dgoldberg Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/dgoldberg/shelfice/shelfice_thermodynamics.F,v 1.3 2015/04/21 11:07:10 dgoldberg Exp $ C $Name: $ #include "SHELFICE_OPTIONS.h" -#ifdef ALLOW_STREAMICE -# include "STREAMICE_OPTIONS.h" +#ifdef ALLOW_AUTODIFF +# include "AUTODIFF_OPTIONS.h" +#endif +#ifdef ALLOW_CTRL +# include "CTRL_OPTIONS.h" #endif CBOP @@ -21,6 +24,7 @@ C | stresses at the ice/water interface are computed in separate C | routines that are called from mom_fluxform/mom_vecinv C *=============================================================* +C \ev C !USES: IMPLICIT NONE @@ -33,6 +37,7 @@ #include "DYNVARS.h" #include "FFIELDS.h" #include "SHELFICE.h" +#include "SHELFICE_COST.h" #ifdef ALLOW_AUTODIFF # include "CTRL_SIZE.h" # include "ctrl.h" @@ -44,10 +49,6 @@ # include "tamc_keys.h" # endif /* SHI_ALLOW_GAMMAFRICT */ #endif /* ALLOW_AUTODIFF_TAMC */ -#ifdef ALLOW_STREAMICE -# include "STREAMICE.h" -#endif - C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === @@ -57,7 +58,6 @@ _RL myTime INTEGER myIter INTEGER myThid -CEOP #ifdef ALLOW_SHELFICE C !LOCAL VARIABLES : @@ -72,10 +72,11 @@ C convertFW2SaltLoc:: local copy of convertFW2Salt C cFac :: 1 for conservative form, 0, otherwise C rFac :: realFreshWaterFlux factor -C dFac :: 0 for diffusive heat flux (Holland and Jenkins, 1999, eq21) +C dFac :: 0 for diffusive heat flux (Holland and Jenkins, 1999, +C eq21) C 1 for advective and diffusive heat flux (eq22, 26, 31) -C fwflxFac :: only effective for dFac=1, 1 if we expect a melting fresh -C water flux, 0 otherwise +C fwflxFac :: only effective for dFac=1, 1 if we expect a melting +C fresh water flux, 0 otherwise C auxiliary variables and abbreviations: C a0, a1, a2, b, c0 C eps1, eps2, eps3, eps3a, eps4, eps5, eps6, eps7, eps8 @@ -95,6 +96,7 @@ _RL cFac, rFac, dFac, fwflxFac _RL aqe, bqe, cqe, discrim, recip_aqe _RL drKp1, recip_drLoc + _RL recip_latentHeat _RL tmpFac #ifdef SHI_ALLOW_GAMMAFRICT @@ -102,6 +104,9 @@ _RL gammaTmoleT, gammaTmoleS, gammaTurb, gammaTurbConst _RL ustar, ustarSq, etastar PARAMETER ( shiTwoThirds = 0.66666666666666666666666666667D0 ) +#ifdef ALLOW_DIAGNOSTICS + _RL uStarDiag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +#endif /* ALLOW_DIAGNOSTICS */ #endif #ifndef ALLOW_OPENAD @@ -112,10 +117,11 @@ #ifdef ALLOW_SHIFWFLX_CONTROL _RL xx_shifwflx_loc(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy) #endif +CEOP C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef SHI_ALLOW_GAMMAFRICT -#ifdef ALLOW_AUTODIFF_TAMC +#ifdef ALLOW_AUTODIFF C re-initialize here again, curtesy to TAF DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) @@ -127,7 +133,7 @@ ENDDO ENDDO ENDDO -#endif /* ALLOW_AUTODIFF_TAMC */ +#endif /* ALLOW_AUTODIFF */ IF ( SHELFICEuseGammaFrict ) THEN C Implement friction velocity-dependent transfer coefficient C of Holland and Jenkins, JPO, 1999 @@ -143,7 +149,7 @@ etastar = 1. _d 0 gammaTurbConst = 1. _d 0 / (2. _d 0 * shiZetaN*etastar) & - recip_shiKarman -#ifdef ALLOW_AUTODIFF_TAMC +#ifdef ALLOW_AUTODIFF DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO J = 1-OLy,sNy+OLy @@ -154,10 +160,13 @@ ENDDO ENDDO ENDDO -#endif /* ALLOW_AUTODIFF_TAMC */ +#endif /* ALLOW_AUTODIFF */ ENDIF #endif /* SHI_ALLOW_GAMMAFRICT */ + recip_latentHeat = 0. _d 0 + IF ( SHELFICElatentHeat .NE. 0. _d 0 ) + & recip_latentHeat = 1. _d 0/SHELFICElatentHeat C are we doing the conservative form of Jenkins et al. (2001)? recip_Cp = 1. _d 0 / HeatCapacity_Cp cFac = 0. _d 0 @@ -200,6 +209,9 @@ shelfIceFreshWaterFlux(I,J,bi,bj) = 0. _d 0 shelficeForcingT (I,J,bi,bj) = 0. _d 0 shelficeForcingS (I,J,bi,bj) = 0. _d 0 +#if (defined SHI_ALLOW_GAMMAFRICT && defined ALLOW_DIAGNOSTICS) + uStarDiag (I,J,bi,bj) = 0. _d 0 +#endif /* SHI_ALLOW_GAMMAFRICT and ALLOW_DIAGNOSTICS */ ENDDO ENDDO ENDDO @@ -214,12 +226,15 @@ ENDDO ENDDO ENDDO - CALL CTRL_GET_GEN ( +#ifdef ALLOW_CTRL + if (useCTRL) CALL CTRL_GET_GEN ( & xx_shifwflx_file, xx_shifwflxstartdate, xx_shifwflxperiod, & maskSHI, xx_shifwflx_loc, xx_shifwflx0, xx_shifwflx1, & xx_shifwflx_dummy, & xx_shifwflx_remo_intercept, xx_shifwflx_remo_slope, + & wshifwflx, & myTime, myIter, myThid ) +#endif #endif /* ALLOW_SHIFWFLX_CONTROL */ DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) @@ -265,6 +280,7 @@ drKp1 = drF(K)*( 1. _d 0 - _hFacC(I,J,K,bi,bj) ) C-- lower cell may not be as thick as required drKp1 = MIN( drKp1, drF(Kp1) * _hFacC(I,J,Kp1,bi,bj) ) + drKp1 = MAX( drKp1, 0. _d 0 ) recip_drLoc = 1. _d 0 / & ( drF(K)*_hFacC(I,J,K,bi,bj) + drKp1 ) tLoc(I,J) = ( tLoc(I,J) * drF(K)*_hFacC(I,J,K,bi,bj) @@ -309,6 +325,9 @@ ustarSq = shiCdrag * MAX( 1.D-6, & 0.25 _d 0 *(uLoc(I,J)*uLoc(I,J)+vLoc(I,J)*vLoc(I,J)) ) ustar = SQRT(ustarSq) +#ifdef ALLOW_DIAGNOSTICS + uStarDiag(I,J,bi,bj) = ustar +#endif /* ALLOW_DIAGNOSTICS */ C instead of etastar = sqrt(1+zetaN*ustar./(f*Lo*Rc)) C etastar = 1. _d 0 C gammaTurbConst = 1. _d 0 / (2. _d 0 * shiZetaN*etastar) @@ -365,7 +384,7 @@ C and vice versa shelfIceFreshWaterFlux(I,J,bi,bj) = & - shelfIceHeatFlux(I,J,bi,bj) - & *recip_SHELFICElatentHeat + & *recip_latentHeat C-- compute surface tendencies shelficeForcingT(i,j,bi,bj) = & - shelfIceHeatFlux(I,J,bi,bj) @@ -517,28 +536,8 @@ ENDDO ENDDO - IF (SHELFICEallowThinIceMass) 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 + IF (SHELFICEMassStepping) THEN CALL SHELFICE_STEP_ICEMASS( myTime, myIter, myThid ) -#endif ENDIF C-- Calculate new loading anomaly (in case the ice-shelf mass was updated) @@ -556,10 +555,9 @@ ENDDO c ENDIF #endif /* ndef ALLOW_AUTODIFF */ - print *, "GOT HERE 0" + #ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN - print *, "GOT HERE 1" CALL DIAGNOSTICS_FILL_RS(shelfIceFreshWaterFlux,'SHIfwFlx', & 0,1,0,1,1,myThid) CALL DIAGNOSTICS_FILL_RS(shelfIceHeatFlux, 'SHIhtFlx', @@ -577,7 +575,11 @@ & 0,1,0,1,1,myThid) CALL DIAGNOSTICS_FILL(shiTransCoeffS,'SHIgammS', & 0,1,0,1,1,myThid) - print *, "GOT HERE 2" +C Friction velocity +#ifdef SHI_ALLOW_GAMMAFRICT + IF ( SHELFICEuseGammaFrict ) + & CALL DIAGNOSTICS_FILL(uStarDiag,'SHIuStar',0,1,0,1,1,myThid) +#endif /* SHI_ALLOW_GAMMAFRICT */ ENDIF #endif /* ALLOW_DIAGNOSTICS */