SUBROUTINE INTERNALROT
!
!  INTERNAL computes the following terms of the
!  moment equations:
!
!   -(g*D**2/RHO0)*Integral from sigma to 0 of
!                  ((dRHO/dx) - (sigma/D)*(dsigma/dx)*(dRHO/dsigma) )
!
!   -(g*D**2/RHO0)*Integral from sigma to 0 of
!                  ((dRHO/dy) - (sigma/D)*(dsigma/dy)*(dRHO/dsigma)
!
!  where g is gravity, D ia dynamic depth
!  and RHO0 reference density,
!  with an explicit technique.
!
!  The results are stored in DELTAU and DELTAV.
!
!  The method is based on first estimating the gradients in a coordinate-system
!  rotated 45 degrees, then rotation to X-Y coord, and then interpolation.
!
  USE STATE
  IMPLICIT NONE
  REAL, DIMENSION(IM,JM,KB) :: DRHODX,DRHODY,DR1,DR2
  REAL, DIMENSION(KB) :: ZZ2(KB)
  REAL REFDENS,XFACTOR,YFACTOR, DIFF
  INTEGER I,J,K
  REAL DRHO
  REAL DR1NEW,DR2NEW,SQRT2HALV

  LOGICAL :: INITROT=.TRUE.
  
!
  SQRT2HALV = 0.50*SQRT(2.0)
!
  RHO = RHO - RMEAN
!
  REFDENS = 1.0 / RHO0
  ZZ2 = 0.
!$omp parallel do
  DO K = 2,KB-1
     ZZ2(K) = 0.5*(ZZ(K-1)+ZZ(K))
 END DO


 IF(INITROT)THEN
     DCORNER = 0.
     DCORNER(1,1) = H(1,1)*FSM(1,1)
     DO I = 2,IM
         IF (FSM(I-1,1).GT.0. .AND. FSM(I,J).GT.0.) THEN
             DCORNER(I,1) = 0.5*(H(I-1,1)+H(I,1))
         ENDIF
     END DO
     DO J = 2,JM
         IF (FSM(1,J-1).GT.0. .AND. FSM(1,J).GT.0.) THEN
             DCORNER(1,J) = 0.5*(H(1,J-1)+H(1,J))
         ENDIF
     END DO
     DO I = 2,IM
         DO J = 2,JM
             IF (FSM(I,J-1).GT.0. .AND. FSM(I-1,J).GT.0. .AND. &
                  &     FSM(I-1,J-1).GT.0. .AND. FSM(I,J).GT.0.) THEN
                 DCORNER(I,J) = 0.25*(H(I,J-1)+H(I-1,J)+H(I-1,J-1)+H(I,J))
             ENDIF
         END DO
     END DO
     INITROT=.FALSE.
 END IF
 
 
!
!    COMPONENT OF BAROCLINIC PRESSURE GRADIENT in 45 Degree direction
!
  DRHODX = 0.
  DRHODY = 0.
  DR1 = 0.
  DR2 = 0.
  DO J=2,JM-1
     DO I=2,IM-1
        DR1(I,J,1)=GRAV*(-ZZ(1))*DCORNER(I,J)*(RHO(I,J,1)-RHO(I-1,J-1,1))
     END DO
  END DO
  DO K=2,KB-1
     DO J=2,JM-1
        DO I=2,IM-1
           DR1(I,J,K)=DR1(I,J,K-1) + GRAV*.5E0*DZZ(K-1)*DCORNER(I,J) &
                &      *(RHO(I,J,K)-RHO(I-1,J-1,K)+RHO(I,J,K-1)-RHO(I-1,J-1,K-1)) &
                &      +GRAV*.5E0*ZZ2(K)*(D(I,J)-D(I-1,J-1)) &
                &      *(RHO(I,J,K)+RHO(I-1,J-1,K)-RHO(I,J,K-1)-RHO(I-1,J-1,K-1))
        END DO
     END DO
  END DO
!
!  COMPONENT OF BAROCLINIC PRESSURE GRADIENT in 135 Degree direction
!
  DO J=2,JM-1
     DO I=2,IM-1
        DR2(I,J,1)=GRAV*(-ZZ(1))*DCORNER(I,J)*(RHO(I-1,J,1)-RHO(I,J-1,1))
     END DO
  END DO
  DO K=2,KB-1
     DO J=2,JM-1
        DO I=2,IM-1
           DR2(I,J,K)=DR2(I,J,K-1) +GRAV*.5E0*DZZ(K-1)*DCORNER(I,J) &
                &        *(RHO(I-1,J,K)-RHO(I,J-1,K)+RHO(I-1,J,K-1)-RHO(I,J-1,K-1)) &
                &        +GRAV*.5E0*ZZ2(K)*(D(I-1,J)-D(I,J-1)) &
                &        *(RHO(I-1,J,K)+RHO(I,J-1,K)-RHO(I-1,J,K-1)-RHO(I,J-1,K-1)) 
        END DO
     END DO
  END DO
!
!  Rotate DR1 and DR2 45 degrees to get aligned with X and Y.
!
  DO K = 1,KB-1
    DO J = 2,JM-1
      DO I = 2,IM-1
        DR1NEW = (DR1(I,J,K)-DR2(I,J,K))*SQRT2HALV
        DR2NEW = (DR1(I,J,K)+DR2(I,J,K))*SQRT2HALV
        DR1(I,J,K) = DR1NEW
        DR2(I,J,K) = DR2NEW
      END DO
    END DO
  END DO
!
!  Interpolate to DRHODX and DRHODY
!
  DO K = 1,KB-1
    DO J = 2,JM-1
      DO I = 2,IM-1
         DRHODX(I,J,K) = 0.5*(DR1(I,J,K)+DR1(I,J+1,K))
         DRHODY(I,J,K) = 0.5*(DR2(I,J,K)+DR2(I+1,J,K))
      END DO
    END DO
  END DO
        
!
!  Include effect of internal pressure.
!
  XFACTOR = REFDENS/(SQRT(2.0)*DX)
  DELTAU = 0.
  DELTAV = 0.
  YFACTOR = REFDENS/(SQRT(2.0)*DY)
  DO K = 1,KB-1
     DO J = 2,JM-1
        DO I = 2,IM-1
           DELTAU(I,J,K) = - XFACTOR*DRHODX(I,J,K)*DU(I,J)*DUM(I,J)
           DELTAV(I,J,K) = - YFACTOR*DRHODY(I,J,K)*DV(I,J)*DVM(I,J)
        END DO
     END DO
  END DO
!
  RHO = RHO + RMEAN
!
  RETURN
END SUBROUTINE INTERNALROT