--- MITgcm/pkg/mom_fluxform/mom_fluxform.F 2005/09/23 15:19:38 1.26 +++ MITgcm/pkg/mom_fluxform/mom_fluxform.F 2006/11/23 00:45:18 1.38 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.26 2005/09/23 15:19:38 jmc Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.38 2006/11/23 00:45:18 jmc Exp $ C $Name: $ CBOI @@ -31,7 +31,7 @@ C !ROUTINE: MOM_FLUXFORM C !INTERFACE: ========================================================== - SUBROUTINE MOM_FLUXFORM( + SUBROUTINE MOM_FLUXFORM( I bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown, I KappaRU, KappaRV, U fVerU, fVerV, @@ -52,6 +52,11 @@ #include "PARAMS.h" #include "GRID.h" #include "SURFACE.h" +#ifdef ALLOW_AUTODIFF_TAMC +# include "tamc.h" +# include "tamc_keys.h" +# include "MOM_FLUXFORM.h" +#endif C !INPUT PARAMETERS: =================================================== C bi,bj :: tile indices @@ -93,6 +98,9 @@ C fMer :: meridional fluxes C fVrUp,fVrDw :: vertical viscous fluxes at interface k-1 & k INTEGER i,j +#ifdef ALLOW_AUTODIFF_TAMC + INTEGER imomkey +#endif _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL cF(1-OLx:sNx+OLx,1-OLy:sNy+OLy) @@ -101,14 +109,14 @@ _RL fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL fVrUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL fVrDw(1-OLx:sNx+OLx,1-OLy:sNy+OLy) -C afFacMom - Tracer parameters for turning terms -C vfFacMom on and off. -C pfFacMom afFacMom - Advective terms +C afFacMom :: Tracer parameters for turning terms on and off. +C vfFacMom +C pfFacMom afFacMom - Advective terms C cfFacMom vfFacMom - Eddy viscosity terms -C mTFacMom pfFacMom - Pressure terms +C mtFacMom pfFacMom - Pressure terms C cfFacMom - Coriolis terms C foFacMom - Forcing -C mTFacMom - Metric term +C mtFacMom - Metric term C uDudxFac, AhDudxFac, etc ... individual term parameters for switching terms off _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy) @@ -137,6 +145,7 @@ _RL ArDudrFac _RL fuFac _RL mtFacU + _RL mtNHFacU _RL uDvdxFac _RL AhDvdxFac _RL vDvdyFac @@ -145,9 +154,28 @@ _RL ArDvdrFac _RL fvFac _RL mtFacV + _RL mtNHFacV + _RL sideMaskFac LOGICAL bottomDragTerms,harmonic,biharmonic,useVariableViscosity CEOP +#ifdef ALLOW_AUTODIFF_TAMC + act0 = k - 1 + max0 = Nr + act1 = bi - myBxLo(myThid) + max1 = myBxHi(myThid) - myBxLo(myThid) + 1 + act2 = bj - myByLo(myThid) + max2 = myByHi(myThid) - myByLo(myThid) + 1 + act3 = myThid - 1 + max3 = nTx*nTy + act4 = ikey_dynamics - 1 + imomkey = (act0 + 1) + & + act1*max0 + & + act2*max0*max1 + & + act3*max0*max1*max2 + & + act4*max0*max1*max2*max3 +#endif /* ALLOW_AUTODIFF_TAMC */ + C Initialise intermediate terms DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx @@ -161,6 +189,9 @@ fVrDw(i,j)= 0. rTransU(i,j)= 0. rTransV(i,j)= 0. +c KE(i,j) = 0. +c hDiv(i,j) = 0. + vort3(i,j) = 0. strain(i,j) = 0. tension(i,j)= 0. guDiss(i,j) = 0. @@ -176,7 +207,8 @@ AhDudyFac = vfFacMom*1. rVelDudrFac = afFacMom*1. ArDudrFac = vfFacMom*1. - mTFacU = mtFacMom*1. + mtFacU = mtFacMom*1. + mtNHFacU = 1. fuFac = cfFacMom*1. C o V momentum equation uDvdxFac = afFacMom*1. @@ -185,7 +217,8 @@ AhDvdyFac = vfFacMom*1. rVelDvdrFac = afFacMom*1. ArDvdrFac = vfFacMom*1. - mTFacV = mtFacMom*1. + mtFacV = mtFacMom*1. + mtNHFacV = 1. fvFac = cfFacMom*1. IF (implicitViscosity) THEN @@ -193,6 +226,14 @@ ArDvdrFac = 0. ENDIF +C note: using standard stencil (no mask) results in under-estimating +C vorticity at a no-slip boundary by a factor of 2 = sideDragFactor + IF ( no_slip_sides ) THEN + sideMaskFac = sideDragFactor + ELSE + sideMaskFac = 0. _d 0 + ENDIF + IF ( no_slip_bottom & .OR. bottomDragQuadratic.NE.0. & .OR. bottomDragLinear.NE.0.) THEN @@ -232,15 +273,46 @@ ENDDO CALL MOM_CALC_KE(bi,bj,k,2,uFld,vFld,KE,myThid) - CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid) - CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid) - CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid) - CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid) + IF ( momViscosity) THEN + CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid) + CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid) + CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid) + CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid) + DO j=1-Oly,sNy+Oly + DO i=1-Olx,sNx+Olx + IF ( hFacZ(i,j).EQ.0. ) THEN + vort3(i,j) = sideMaskFac*vort3(i,j) + strain(i,j) = sideMaskFac*strain(i,j) + ENDIF + ENDDO + ENDDO +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + CALL DIAGNOSTICS_FILL(hDiv, 'momHDiv ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(vort3, 'momVort3',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(tension,'Tension ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(strain, 'Strain ',k,1,2,bi,bj,myThid) + ENDIF +#endif + ENDIF C--- First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp) IF (momAdvection.AND.k.EQ.1) THEN C- Calculate vertical transports above U & V points (West & South face): + +#ifdef ALLOW_AUTODIFF_TAMC +# ifdef NONLIN_FRSURF +# ifndef DISABLE_RSTAR_CODE +CADJ STORE dwtransc(:,:,bi,bj) = +CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte +CADJ STORE dwtransu(:,:,bi,bj) = +CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte +CADJ STORE dwtransv(:,:,bi,bj) = +CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte +# endif +# endif /* NONLIN_FRSURF */ +#endif /* ALLOW_AUTODIFF_TAMC */ CALL MOM_CALC_RTRANS( k, bi, bj, O rTransU, rTransV, I myTime, myIter, myThid) @@ -322,6 +394,7 @@ #ifdef NONLIN_FRSURF C-- account for 3.D divergence of the flow in rStar coordinate: +# ifndef DISABLE_RSTAR_CODE IF ( select_rStar.GT.0 ) THEN DO j=jMin,jMax DO i=iMin,iMax @@ -339,6 +412,7 @@ ENDDO ENDDO ENDIF +# endif /* DISABLE_RSTAR_CODE */ #endif /* NONLIN_FRSURF */ ELSE @@ -356,16 +430,16 @@ C--- Calculate eddy fluxes (dissipation) between cells for zonal flow. C Bi-harmonic term del^2 U -> v4F - IF (biharmonic) + IF (biharmonic) & CALL MOM_U_DEL2U(bi,bj,k,uFld,hFacZ,v4f,myThid) C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon CALL MOM_U_XVISCFLUX(bi,bj,k,uFld,v4F,fZon, - I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid) + I viscAh_D,viscA4_D,myThid) C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,fMer, - I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid) + I viscAh_Z,viscA4_Z,myThid) C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw IF (.NOT.implicitViscosity) THEN @@ -400,10 +474,16 @@ ENDIF #endif -C-- No-slip and drag BCs appear as body forces in cell abutting topography +C-- No-slip and drag BCs appear as body forces in cell abutting topography IF (no_slip_sides) THEN C- No-slip BCs impose a drag at walls... - CALL MOM_U_SIDEDRAG(bi,bj,k,uFld,v4F,hFacZ,vF,myThid) + CALL MOM_U_SIDEDRAG( + I bi,bj,k, + I uFld, v4f, hFacZ, + I viscAh_Z,viscA4_Z, + I harmonic,biharmonic,useVariableViscosity, + O vF, + I myThid) DO j=jMin,jMax DO i=iMin,iMax gUdiss(i,j) = gUdiss(i,j) + vF(i,j) @@ -420,6 +500,17 @@ ENDDO ENDIF +#ifdef ALLOW_SHELFICE + IF (useShelfIce) THEN + CALL SHELFICE_U_DRAG(bi,bj,k,uFld,KE,KappaRU,vF,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + gUdiss(i,j) = gUdiss(i,j) + vF(i,j) + ENDDO + ENDDO + ENDIF +#endif /* ALLOW_SHELFICE */ + C- endif momViscosity ENDIF @@ -431,28 +522,30 @@ C-- Metric terms for curvilinear grid systems IF (useNHMTerms) THEN -C o Non-hydrosatic metric terms +C o Non-Hydrostatic (spherical) metric terms CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid) DO j=jMin,jMax DO i=iMin,iMax - gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j) + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtNHFacU*mT(i,j) ENDDO ENDDO ENDIF - IF (usingSphericalPolarMTerms) THEN + IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN +C o Spherical polar grid metric terms CALL MOM_U_METRIC_SPHERE(bi,bj,k,uFld,vFld,mT,myThid) DO j=jMin,jMax DO i=iMin,iMax - gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j) + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j) ENDDO ENDDO ENDIF - IF (usingCylindricalGrid) THEN - CALL MOM_U_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid) - DO j=jMin,jMax - DO i=iMin,iMax - gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j) - ENDDO + IF ( usingCylindricalGrid .AND. metricTerms ) THEN +C o Cylindrical grid metric terms + CALL MOM_U_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j) + ENDDO ENDDO ENDIF @@ -504,6 +597,7 @@ #ifdef NONLIN_FRSURF C-- account for 3.D divergence of the flow in rStar coordinate: +# ifndef DISABLE_RSTAR_CODE IF ( select_rStar.GT.0 ) THEN DO j=jMin,jMax DO i=iMin,iMax @@ -521,6 +615,7 @@ ENDDO ENDDO ENDIF +# endif /* DISABLE_RSTAR_CODE */ #endif /* NONLIN_FRSURF */ ELSE @@ -537,16 +632,16 @@ IF (momViscosity) THEN C--- Calculate eddy fluxes (dissipation) between cells for meridional flow. C Bi-harmonic term del^2 V -> v4F - IF (biharmonic) + IF (biharmonic) & CALL MOM_V_DEL2V(bi,bj,k,vFld,hFacZ,v4f,myThid) C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon CALL MOM_V_XVISCFLUX(bi,bj,k,vFld,v4f,hFacZ,fZon, - I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid) + I viscAh_Z,viscA4_Z,myThid) C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer CALL MOM_V_YVISCFLUX(bi,bj,k,vFld,v4f,fMer, - I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid) + I viscAh_D,viscA4_D,myThid) C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw IF (.NOT.implicitViscosity) THEN @@ -581,10 +676,16 @@ ENDIF #endif -C-- No-slip and drag BCs appear as body forces in cell abutting topography - IF (no_slip_sides) THEN +C-- No-slip and drag BCs appear as body forces in cell abutting topography + IF (no_slip_sides) THEN C- No-slip BCs impose a drag at walls... - CALL MOM_V_SIDEDRAG(bi,bj,k,vFld,v4F,hFacZ,vF,myThid) + CALL MOM_V_SIDEDRAG( + I bi,bj,k, + I vFld, v4f, hFacZ, + I viscAh_Z,viscA4_Z, + I harmonic,biharmonic,useVariableViscosity, + O vF, + I myThid) DO j=jMin,jMax DO i=iMin,iMax gvDiss(i,j) = gvDiss(i,j) + vF(i,j) @@ -601,6 +702,17 @@ ENDDO ENDIF +#ifdef ALLOW_SHELFICE + IF (useShelfIce) THEN + CALL SHELFICE_V_DRAG(bi,bj,k,vFld,KE,KappaRU,vF,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + gvDiss(i,j) = gvDiss(i,j) + vF(i,j) + ENDDO + ENDDO + ENDIF +#endif /* ALLOW_SHELFICE */ + C- endif momViscosity ENDIF @@ -612,29 +724,31 @@ C-- Metric terms for curvilinear grid systems IF (useNHMTerms) THEN -C o Spherical polar grid metric terms +C o Non-Hydrostatic (spherical) metric terms CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid) DO j=jMin,jMax DO i=iMin,iMax - gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j) + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtNHFacV*mT(i,j) ENDDO ENDDO ENDIF - IF (usingSphericalPolarMTerms) THEN + IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN +C o Spherical polar grid metric terms CALL MOM_V_METRIC_SPHERE(bi,bj,k,uFld,mT,myThid) DO j=jMin,jMax DO i=iMin,iMax - gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j) + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j) ENDDO ENDDO ENDIF - IF (usingCylindricalGrid) THEN - CALL MOM_V_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid) - DO j=jMin,jMax - DO i=iMin,iMax - gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j) - ENDDO - ENDDO + IF ( usingCylindricalGrid .AND. metricTerms ) THEN +C o Cylindrical grid metric terms + CALL MOM_V_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j) + ENDDO + ENDDO ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| @@ -667,13 +781,23 @@ #endif ENDIF - IF (nonHydrostatic.OR.quasiHydrostatic) THEN - CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid) - DO j=jMin,jMax - DO i=iMin,iMax - gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j) +C-- 3.D Coriolis term (horizontal momentum, Eastward component: -f'*w) + IF ( use3dCoriolis ) THEN + CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+fuFac*cf(i,j) + ENDDO ENDDO - ENDDO + IF ( usingCurvilinearGrid ) THEN +C- presently, non zero angleSinC array only supported with Curvilinear-Grid + CALL MOM_V_CORIOLIS_NH(bi,bj,k,wVel,cf,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+fvFac*cf(i,j) + ENDDO + ENDDO + ENDIF ENDIF C-- Set du/dt & dv/dt on boundaries to zero @@ -688,8 +812,7 @@ #ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN - IF (bottomDragTerms) - & CALL DIAGNOSTICS_FILL(KE, 'momKE ',k,1,2,bi,bj,myThid) + CALL DIAGNOSTICS_FILL(KE, 'momKE ',k,1,2,bi,bj,myThid) CALL DIAGNOSTICS_FILL(gU(1-Olx,1-Oly,k,bi,bj), & 'Um_Advec',k,1,2,bi,bj,myThid) CALL DIAGNOSTICS_FILL(gV(1-Olx,1-Oly,k,bi,bj),