--- MITgcm/pkg/mom_vecinv/mom_vecinv.F 2005/12/08 15:44:34 1.55 +++ MITgcm/pkg/mom_vecinv/mom_vecinv.F 2008/04/01 01:27:33 1.62 @@ -1,9 +1,9 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mom_vecinv/mom_vecinv.F,v 1.55 2005/12/08 15:44:34 heimbach Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mom_vecinv/mom_vecinv.F,v 1.62 2008/04/01 01:27:33 jmc Exp $ C $Name: $ #include "MOM_VECINV_OPTIONS.h" - SUBROUTINE MOM_VECINV( + SUBROUTINE MOM_VECINV( I bi,bj,iMin,iMax,jMin,jMax,k,kUp,kDown, I KappaRU, KappaRV, U fVerU, fVerV, @@ -38,6 +38,10 @@ #ifdef ALLOW_TIMEAVE #include "TIMEAVE_STATV.h" #endif +#ifdef ALLOW_AUTODIFF_TAMC +# include "tamc.h" +# include "tamc_keys.h" +#endif C == Routine arguments == C fVerU :: Flux of momentum in the vertical direction, out of the upper @@ -71,7 +75,6 @@ _RL vrF(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL vCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy) -c _RL mT (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS r_hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy) @@ -94,13 +97,14 @@ INTEGER i,j,k C xxxFac - On-off tracer parameters used for switching terms off. _RL ArDudrFac -c _RL mtFacU _RL ArDvdrFac -c _RL mtFacV _RL sideMaskFac LOGICAL bottomDragTerms LOGICAL writeDiag LOGICAL harmonic,biharmonic,useVariableViscosity +#ifdef ALLOW_AUTODIFF_TAMC + INTEGER imomkey +#endif #ifdef ALLOW_MNC INTEGER offsets(9) @@ -116,6 +120,23 @@ fVerV(1,1,kUp) = fVerV(1,1,kUp) #endif +#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 */ + writeDiag = DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) #ifdef ALLOW_MNC @@ -135,18 +156,17 @@ offsets(i) = 0 ENDDO offsets(3) = k -C write(*,*) 'offsets = ',(offsets(i),i=1,9) +c write(*,*) 'offsets = ',(offsets(i),i=1,9) ENDIF #endif /* ALLOW_MNC */ -C Initialise intermediate terms - DO J=1-OLy,sNy+OLy - DO I=1-OLx,sNx+OLx +C-- Initialise intermediate terms + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx vF(i,j) = 0. vrF(i,j) = 0. uCf(i,j) = 0. vCf(i,j) = 0. -c mT(i,j) = 0. del2u(i,j) = 0. del2v(i,j) = 0. dStar(i,j) = 0. @@ -156,14 +176,16 @@ vort3(i,j) = 0. omega3(i,j)= 0. KE(i,j) = 0. +C- need to initialise hDiv for MOM_VI_DEL2UV(call FILL_CS_CORNER_TR_RL) + hDiv(i,j) = 0. viscAh_Z(i,j) = 0. viscAh_D(i,j) = 0. viscA4_Z(i,j) = 0. viscA4_D(i,j) = 0. -#ifdef ALLOW_AUTODIFF_TAMC strain(i,j) = 0. _d 0 tension(i,j) = 0. _d 0 +#ifdef ALLOW_AUTODIFF_TAMC hFacZ(i,j) = 0. _d 0 #endif ENDDO @@ -172,10 +194,8 @@ C-- Term by term tracer parmeters C o U momentum equation ArDudrFac = vfFacMom*1. -c mTFacU = mtFacMom*1. C o V momentum equation ArDvdrFac = vfFacMom*1. -c mTFacV = mtFacMom*1. C note: using standard stencil (no mask) results in under-estimating C vorticity at a no-slip boundary by a factor of 2 = sideDragFactor @@ -213,9 +233,14 @@ CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid) IF (momViscosity) THEN -C-- For viscous term, compute horizontal divergence, tension & strain +C-- For viscous term, compute horizontal divergence, tension & strain C and mask relative vorticity (free-slip case): +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE vort3(:,:) = +CADJ & comlev1_bibj_k, key = imomkey, byte = isbyte +#endif + CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid) CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid) @@ -291,7 +316,7 @@ I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D, I harmonic,biharmonic,useVariableViscosity, O guDiss,gvDiss, - & myThid) + & myThid) ENDIF C-- if (momViscosity) end of block. ENDIF @@ -327,7 +352,7 @@ ENDDO 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 (momViscosity.AND.no_slip_sides) THEN C- No-slip BCs impose a drag at walls... CALL MOM_U_SIDEDRAG( @@ -352,6 +377,17 @@ ENDDO ENDDO ENDIF +#ifdef ALLOW_SHELFICE + IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) 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--- Other dissipation terms in Meridional momentum equation @@ -381,7 +417,7 @@ ENDDO 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 (momViscosity.AND.no_slip_sides) THEN C- No-slip BCs impose a drag at walls... CALL MOM_V_SIDEDRAG( @@ -406,6 +442,17 @@ ENDDO ENDDO ENDIF +#ifdef ALLOW_SHELFICE + IF (useShelfIce.AND.momViscosity.AND.bottomDragTerms) 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- Vorticity diagnostics: IF ( writeDiag ) THEN @@ -459,7 +506,6 @@ gV(i,j,k,bi,bj) = vCf(i,j) ENDDO ENDDO - IF ( writeDiag ) THEN IF (snapshot_mdsio) THEN CALL WRITE_LOCAL_RL('fV','I10',1,uCf,bi,bj,k,myIter,myThid) @@ -480,7 +526,6 @@ CALL DIAGNOSTICS_FILL(vCf,'Vm_Cori ',k,1,2,bi,bj,myThid) ENDIF #endif /* ALLOW_DIAGNOSTICS */ - ELSE DO j=jMin,jMax DO i=iMin,iMax @@ -492,13 +537,14 @@ IF (momAdvection) THEN C-- Horizontal advection of relative (or absolute) vorticity - IF (highOrderVorticity.AND.useAbsVorticity) THEN + IF ( (highOrderVorticity.OR.upwindVorticity) + & .AND.useAbsVorticity ) THEN CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,omega3,r_hFacZ, & uCf,myThid) - ELSEIF (highOrderVorticity) THEN + ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,vFld,vort3, r_hFacZ, & uCf,myThid) - ELSEIF (useAbsVorticity) THEN + ELSEIF ( useAbsVorticity ) THEN CALL MOM_VI_U_CORIOLIS(bi,bj,K,vFld,omega3,hFacZ,r_hFacZ, & uCf,myThid) ELSE @@ -510,13 +556,14 @@ gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j) ENDDO ENDDO - IF (highOrderVorticity.AND.useAbsVorticity) THEN + IF ( (highOrderVorticity.OR.upwindVorticity) + & .AND.useAbsVorticity ) THEN CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,omega3,r_hFacZ, & vCf,myThid) - ELSEIF (highOrderVorticity) THEN + ELSEIF ( (highOrderVorticity.OR.upwindVorticity) ) THEN CALL MOM_VI_V_CORIOLIS_C4(bi,bj,K,uFld,vort3, r_hFacZ, & vCf,myThid) - ELSEIF (useAbsVorticity) THEN + ELSEIF ( useAbsVorticity ) THEN CALL MOM_VI_V_CORIOLIS(bi,bj,K,uFld,omega3,hFacZ,r_hFacZ, & vCf,myThid) ELSE @@ -612,22 +659,40 @@ C-- end if momAdvection ENDIF -C-- Metric terms for curvilinear grid systems -c IF (usingSphericalPolarMTerms) THEN -C o Spherical polar grid metric terms -c CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,mT,myThid) -c DO j=jMin,jMax -c DO i=iMin,iMax -c gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j) -c ENDDO -c ENDDO -c CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,mT,myThid) -c DO j=jMin,jMax -c DO i=iMin,iMax -c gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j) -c ENDDO -c ENDDO -c ENDIF +C-- 3.D Coriolis term (horizontal momentum, Eastward component: -f'*w) + IF ( use3dCoriolis ) THEN + CALL MOM_U_CORIOLIS_NH(bi,bj,k,wVel,uCf,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j) + 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,vCf,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j) + ENDDO + ENDDO + ENDIF + ENDIF + +C-- Non-Hydrostatic (spherical) metric terms + IF ( useNHMTerms ) THEN + CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,uCf,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j) + ENDDO + ENDDO + CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,vCf,myThid) + DO j=jMin,jMax + DO i=iMin,iMax + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j) + ENDDO + ENDDO + ENDIF C-- Set du/dt & dv/dt on boundaries to zero DO j=jMin,jMax