#include "SBO_OPTIONS.h"
#ifdef ALLOW_SEAICE
# include "SEAICE_OPTIONS.h"
#endif

CBOP
C !ROUTINE: SBO_CALC

C !INTERFACE: ==========================================================
      SUBROUTINE SBO_CALC( myTime, myIter, myThid )

C !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SBO_CALC                                      |
C     | o Do SBO diagnostic output.                              |
C     *==========================================================*
C     | NOTE: The following subtleties are ignored for time      |
C     | being but may need revisiting at some point in time.     |
C     | 1) The model is volume-preserving and Boussinesq so      |
C     |    quantities like oceanic mass need to be interpreted   |
C     |    with some care.  We remove spurious mass variations   |
C     |    using the Greatbatch correction. Real freshwater      |
C     |    fluxes retained in mass load.                         |
C     | 2) The sea surface height variable etaN might lag other  |
C     |    prognostic variables by half a time step.  This lag   |
C     |    is ignored in SBO computations.                       |
C     | 3) OAM due to currents assumes constant density          |
C     |    (=rhoConst), rms differences using variable density   |
C     |    is less than 1%, assuming rhoConst is a good measure  |
C     |    of the actual mean density                            |
C     | 4) Seaice motion added to OAMC.  Seaice mass is in OAMP  |
C     |    and COM.   Net freshwater flux is between atmosphere  |
C     |    and liquid ocean plus seaice.  I.e. changes in seaice |
C     |    mass due to melt/freeze with liquid ocean do not      |
C     |    change net freshwater flux.                           |
C     *==========================================================*

C=======================================================================
C
C     Based on ftp://euler.jpl.nasa.gov/sbo/software/calc_sbo2.f
C     Written  by Richard Gross (Richard.Gross@jpl.nasa.gov)
C     Reference
C           Gross, R. S., F. O. Bryan, Y. Chao, J. O. Dickey, S. L. Marcus,
C           R. M. Ponte, and R. Tokmakian, The IERS Special Bureau for the
C           Oceans, in IERS Technical Note on the IERS Global Geophysical
C           Fluids Center, edited by B. Chao, in press, Observatoire de Paris,
C           Paris, France, 2002.
C
C     June 10, 2001: Modified for online computations in MIT GCM UV
C              by Dimitris Menemenlis (Menemenlis@jpl.nasa.gov)
C     Jan 7, 2014: Modified for real freshwater flux and coordinates other
C              than spherical polar by Katy Quinn (kquinn@aer.com)
C
C     Purpose
C           calc_sbo calculates the core products of the IERS Special Bureau
C           for the Oceans including oceanic mass, center-of-mass, and angular
C           momentum.
C
C=======================================================================
C \ev

C !USES: ===============================================================
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#include "SBO.h"
#ifdef ALLOW_SEAICE
# include "SEAICE_SIZE.h"
# include "SEAICE.h"
#endif

C !INPUT PARAMETERS: ===================================================
C     == Routine arguments ==
C     myTime  :: Current time of simulation ( s )
C     myIter  :: Iteration number
C     myThid  :: Number of this instance of SBO_CALC
      _RL     myTime
      INTEGER myIter, myThid

#ifdef ALLOW_SBO

C !LOCAL VARIABLES: ====================================================
C     bi, bj    :: array indices
C     i         :: index over longitude grid points
C     j         :: index over latitude  grid points
C     k         :: index over layers
C     lat       :: latitude  of grid point (radians)
C     lon       :: longitude of grid point (radians)
C     darea     :: element of surface area (m**2)
C     dvolume   :: element of volume (m**3)
C     ae        :: mean radius of Earth (m) (PREM value)
C     sbo_omega :: mean angular velocity of Earth (rad/s)
C     UE,VN     :: geographic (east,north) ocean velocities at cell centers (m/s)
C     UEice,VNice  :: geographic (east,north) seaice velocities at cell centers (m/s)
C     Mload     :: total mass load (kg/m**2)
C     GCload    :: mass load for Greatbatch correction (kg/m**2)
C     FWLoad    :: real freshwater flux mass load (kg/m**2)
      integer bi, bj, i, j, k
      _RL lat, lon, darea, dvolume
      _RL ae, sbo_omega
      PARAMETER ( ae        = 6.3710 _d 6    )
      PARAMETER ( sbo_omega = 7.292115 _d -5 )
      _RL UE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL VN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL UEice(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL VNice(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL Mload(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL GCload, FWload
C Tiled global sums
      _RL tile_FWload(nSx,nSy)
      _RL tile_sboarea(nSx,nSy)
      _RL tile_GCload(nSx,nSy)
      _RL tile_mass(nSx,nSy)
      _RL tile_xcom(nSx,nSy)
      _RL tile_ycom(nSx,nSy)
      _RL tile_zcom(nSx,nSy)
      _RL tile_xoamc(nSx,nSy)
      _RL tile_yoamc(nSx,nSy)
      _RL tile_zoamc(nSx,nSy)
      _RL tile_xoamp(nSx,nSy)
      _RL tile_yoamp(nSx,nSy)
      _RL tile_zoamp(nSx,nSy)
      _RL tile_xoamc_si(nSx,nSy)
      _RL tile_yoamc_si(nSx,nSy)
      _RL tile_zoamc_si(nSx,nSy)
      _RL tile_mass_si(nSx,nSy)
      _RL tile_mass_fw(nSx,nSy)
      _RL tile_xcom_fw(nSx,nSy)
      _RL tile_ycom_fw(nSx,nSy)
      _RL tile_zcom_fw(nSx,nSy)
      _RL tile_xoamp_fw(nSx,nSy)
      _RL tile_yoamp_fw(nSx,nSy)
      _RL tile_zoamp_fw(nSx,nSy)
      _RL tile_mass_gc (nSx,nSy)
C Pre-computed cos(lat), sin(lat), cos(lon), sin(lon)
      _RL COSlat(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL SINlat(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL COSlon(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL SINlon(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
CEOP

C Initialize variables to be computed---------------------------------
C-    note: only done once (by master thread) for var in common block
      _BEGIN_MASTER(myThid)

      xoamc = 0.0
      yoamc = 0.0
      zoamc = 0.0
      xoamp = 0.0
      yoamp = 0.0
      zoamp = 0.0
      mass  = 0.0
      xcom  = 0.0
      ycom  = 0.0
      zcom  = 0.0
      sboarea = 0.0

      xoamc_si = 0.0
      yoamc_si = 0.0
      zoamc_si = 0.0
      mass_si  = 0.0

      xoamp_fw = 0.0
      yoamp_fw = 0.0
      zoamp_fw = 0.0
      mass_fw  = 0.0
      xcom_fw  = 0.0
      ycom_fw  = 0.0
      zcom_fw  = 0.0

      mass_gc  = 0.0

      _END_MASTER(myThid)

C Get geographic (East,North) velocities------------------------------

      CALL ROTATE_UV2EN_RL(
     U          uVel, vVel,
     U          UE, VN,
     I          .TRUE., .TRUE., .FALSE., Nr, mythid )

#ifdef ALLOW_SEAICE
      IF ( useSEAICE ) THEN
        CALL ROTATE_UV2EN_RL(
     U          UICE, VICE,
     U          UEice, VNice,
     I          .TRUE., .TRUE., .FALSE., 1, mythid )
      ELSE
#else /* ALLOW_SEAICE */
      IF ( .TRUE. ) THEN
#endif /* ALLOW_SEAICE */
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             UEice(i,j,bi,bj) = 0.
             VNice(i,j,bi,bj) = 0.
           ENDDO
          ENDDO
         ENDDO
        ENDDO
      ENDIF

C Calculate mass load-------------------------------------------------
C     Calculate mass load (Mload), Greatbatch correction for
C     spurious mass but spatial mean freshwater flux retained.
C     Mload *needs* to be total mass (for center of mass), so add
C     back missing time invariant term: -R_low*rhoConst

c Calculate freshwater load
C Calculate Greatbatch correction load over global ocean volume
c Note: no halo regions in i,j loops, do not want to double book sums
      FWload = 0.0
      GCload = 0.0
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        tile_FWload(bi,bj) = 0.0
        tile_GCload(bi,bj) = 0.0
        tile_sboarea(bi,bj) = 0.0
        DO j = 1, sNy
         DO i = 1, sNx
            darea = rA(i,j,bi,bj)*maskC(i,j,1,bi,bj)
            tile_sboarea(bi,bj) = tile_sboarea(bi,bj) + darea
            tile_FWload(bi,bj) =  tile_FWload(bi,bj) +
     &           rhoConst*etaN(i,j,bi,bj)*darea +
     &           sIceLoad(i,j,bi,bj)*darea
            DO k = 1, Nr
               dvolume = rA(i,j,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
               tile_GCload(bi,bj) = tile_GCload(bi,bj) +
     &              rhoInSitu(i,j,k,bi,bj) * dvolume
            ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      CALL GLOBAL_SUM_TILE_RL( tile_FWload  , FWload  , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_sboarea , sboarea , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_GCload  , GCload  , myThid )
      FWload = FWload/sboarea
      GCload = -1.0 * GCload/sboarea

c Total mass load with freshwater flux
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j = 1-OLy, sNy+OLy
         DO i = 1-OLx, sNx+OLx
            Mload(i,j,bi,bj) =  rhoConst*etaN(i,j,bi,bj) +
     &           sIceLoad(i,j,bi,bj) +
     &           GCload - R_low(i,j,bi,bj)*rhoConst
            DO k = 1, Nr
               Mload(i,j,bi,bj) = Mload(i,j,bi,bj) +
     &           rhoInSitu(i,j,k,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
            ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

c Pre-compute cos(lat), sin(lat), cos(lon), sin(lon)
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j = 1-OLy, sNy+OLy
         DO i = 1-OLx, sNx+OLx
            lat = yC(i,j,bi,bj) * deg2rad
            lon = xC(i,j,bi,bj) * deg2rad
            COSlat(i,j,bi,bj) = COS(lat)
            SINlat(i,j,bi,bj) = SIN(lat)
            COSlon(i,j,bi,bj) = COS(lon)
            SINlon(i,j,bi,bj) = SIN(lon)
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C Main loops----------------------------------------------------------
C     loop over all grid points, accumulating mass, com, oam
C     Note: no halo regions in i,j loops, do not want to double book sums

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)

C Initialize tile sums---------------------------------
        tile_xoamc(bi,bj) = 0.0
        tile_yoamc(bi,bj) = 0.0
        tile_zoamc(bi,bj) = 0.0
        tile_xoamp(bi,bj) = 0.0
        tile_yoamp(bi,bj) = 0.0
        tile_zoamp(bi,bj) = 0.0
        tile_mass(bi,bj) = 0.0
        tile_xcom(bi,bj) = 0.0
        tile_ycom(bi,bj) = 0.0
        tile_zcom(bi,bj) = 0.0

        tile_xoamc_si(bi,bj) = 0.0
        tile_yoamc_si(bi,bj) = 0.0
        tile_zoamc_si(bi,bj) = 0.0
        tile_mass_si(bi,bj) = 0.0

        tile_xoamp_fw(bi,bj) = 0.0
        tile_yoamp_fw(bi,bj) = 0.0
        tile_zoamp_fw(bi,bj) = 0.0
        tile_mass_fw(bi,bj) = 0.0
        tile_xcom_fw(bi,bj) = 0.0
        tile_ycom_fw(bi,bj) = 0.0
        tile_zcom_fw(bi,bj) = 0.0

        tile_mass_gc(bi,bj) = 0.0

        DO j = 1, sNy
         DO i = 1, sNx

           IF ( maskC(i,j,1,bi,bj) .NE. 0. ) THEN

C     horizontal area
              darea = rA(i,j,bi,bj)*maskC(i,j,1,bi,bj)

C     accumulate mass of oceans, Greatbatch correction, seaice
              tile_mass(bi,bj) = tile_mass(bi,bj) +
     &          Mload(i,j,bi,bj)*darea
              tile_mass_gc(bi,bj) = tile_mass_gc(bi,bj) +
     &          GCload*darea
              tile_mass_si(bi,bj) = tile_mass_si(bi,bj) +
     &          sIceLoad(i,j,bi,bj)*darea

C     accumulate center-of-mass of oceans (need to divide by total mass at end)
              tile_xcom(bi,bj) = tile_xcom(bi,bj) +
     &          Mload(i,j,bi,bj)*COSlat(i,j,bi,bj)*COSlon(i,j,bi,bj)
     &          * ae * darea
              tile_ycom(bi,bj) = tile_ycom(bi,bj) +
     &          Mload(i,j,bi,bj)*COSlat(i,j,bi,bj)*SINlon(i,j,bi,bj)
     &          * ae * darea
              tile_zcom(bi,bj) = tile_zcom(bi,bj) +
     &          Mload(i,j,bi,bj)*SINlat(i,j,bi,bj)
     &          * ae * darea

C     accumulate oceanic angular momentum due to currents (need depth integral too)
C     Note: depth integral goes from k=1,Nr.  hFacC takes care of R_low and etaN (as per JMC)
              DO k = 1, Nr
                 dvolume = rA(i,j,bi,bj)*drF(k)
     &             * maskC(i,j,k,bi,bj)*hFacC(i,j,k,bi,bj)
                 tile_xoamc(bi,bj) = tile_xoamc(bi,bj) +
     &             ( VN(i,j,k,bi,bj)*SINlon(i,j,bi,bj) -
     &             UE(i,j,k,bi,bj)*
     &             SINlat(i,j,bi,bj)*COSlon(i,j,bi,bj) )
     &             * rhoConst * ae * dvolume
                 tile_yoamc(bi,bj) = tile_yoamc(bi,bj) +
     &             (-VN(i,j,k,bi,bj)*COSlon(i,j,bi,bj) -
     &             UE(i,j,k,bi,bj)*
     &             SINlat(i,j,bi,bj)*SINlon(i,j,bi,bj) )
     &             * rhoConst * ae * dvolume
                 tile_zoamc(bi,bj) = tile_zoamc(bi,bj) +
     &             UE(i,j,k,bi,bj)*COSlat(i,j,bi,bj)
     &             * rhoConst * ae * dvolume
              ENDDO

C     accumulate sea angular momentum due to motion (one layer, so no depth integral needed)
              tile_xoamc_si(bi,bj) = tile_xoamc_si(bi,bj) +
     &          ( VNice(i,j,bi,bj)*SINlon(i,j,bi,bj) -
     &          UEice(i,j,bi,bj)*
     &          SINlat(i,j,bi,bj)*COSlon(i,j,bi,bj) )
     &          * sIceLoad(i,j,bi,bj) * ae * darea
              tile_yoamc_si(bi,bj) = tile_yoamc_si(bi,bj) +
     &          (-VNice(i,j,bi,bj)*COSlon(i,j,bi,bj) -
     &          UEice(i,j,bi,bj)*
     &          SINlat(i,j,bi,bj)*SINlon(i,j,bi,bj) )
     &          * sIceLoad(i,j,bi,bj) * ae * darea
              tile_zoamc_si(bi,bj) = tile_zoamc_si(bi,bj) +
     &          UEice(i,j,bi,bj)*COSlat(i,j,bi,bj)
     &          * sIceLoad(i,j,bi,bj) * ae * darea

C     accumulate oceanic angular momentum due to pressure
              tile_xoamp(bi,bj) = tile_xoamp(bi,bj) -
     &          SINlat(i,j,bi,bj)*COSlat(i,j,bi,bj)*COSlon(i,j,bi,bj)
     &          * sbo_omega * Mload(i,j,bi,bj) * ae*ae * darea
              tile_yoamp(bi,bj) = tile_yoamp(bi,bj) -
     &          SINlat(i,j,bi,bj)*COSlat(i,j,bi,bj)*SINlon(i,j,bi,bj)
     &          * sbo_omega * Mload(i,j,bi,bj) * ae*ae * darea
              tile_zoamp(bi,bj) = tile_zoamp(bi,bj) +
     &          COSlat(i,j,bi,bj) * COSlat(i,j,bi,bj)
     &          * sbo_omega * Mload(i,j,bi,bj) * ae*ae * darea

C     accumulate mass of real freshwater flux
              tile_mass_fw(bi,bj) = tile_mass_fw(bi,bj) +
     &          FWload * darea

C     accumulate center-of-mass of real freshwater flux (need to divide by total FW mass at end)
              tile_xcom_fw(bi,bj) = tile_xcom_fw(bi,bj) +
     &          FWload * COSlat(i,j,bi,bj) * COSlon(i,j,bi,bj)
     &          * ae * darea
              tile_ycom_fw(bi,bj) = tile_ycom_fw(bi,bj) +
     &          FWload * COSlat(i,j,bi,bj) * SINlon(i,j,bi,bj)
     &          * ae * darea
              tile_zcom_fw(bi,bj) = tile_zcom_fw(bi,bj) +
     &          FWload * SINlat(i,j,bi,bj)
     &          * ae * darea

C     accumulate oceanic angular momentum due to real freshwater flux
              tile_xoamp_fw(bi,bj) = tile_xoamp_fw(bi,bj) -
     &          SINlat(i,j,bi,bj)*COSlat(i,j,bi,bj)*COSlon(i,j,bi,bj)
     &          * sbo_omega * FWload * ae*ae * darea
              tile_yoamp_fw(bi,bj) = tile_yoamp_fw(bi,bj) -
     &          SINlat(i,j,bi,bj)*COSlat(i,j,bi,bj)*SINlon(i,j,bi,bj)
     &          * sbo_omega * FWload * ae*ae * darea
              tile_zoamp_fw(bi,bj) = tile_zoamp_fw(bi,bj) +
     &          COSlat(i,j,bi,bj) * COSlat(i,j,bi,bj)
     &          * sbo_omega * FWload * ae*ae * darea

C     end if over ocean
           ENDIF
C     end loop over i,j
         ENDDO
        ENDDO

C     end loop over bi,bj
       ENDDO
      ENDDO

C     sum all global values across model tiles
      CALL GLOBAL_SUM_TILE_RL( tile_mass  , mass  , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_xcom  , xcom  , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_ycom  , ycom  , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_zcom  , zcom  , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_xoamc , xoamc , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_yoamc , yoamc , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_zoamc , zoamc , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_xoamp , xoamp , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_yoamp , yoamp , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_zoamp , zoamp , myThid )

      CALL GLOBAL_SUM_TILE_RL( tile_xoamc_si , xoamc_si , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_yoamc_si , yoamc_si , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_zoamc_si , zoamc_si , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_mass_si , mass_si , myThid )

      CALL GLOBAL_SUM_TILE_RL( tile_mass_fw  , mass_fw  , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_xcom_fw  , xcom_fw  , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_ycom_fw  , ycom_fw  , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_zcom_fw  , zcom_fw  , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_xoamp_fw , xoamp_fw , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_yoamp_fw , yoamp_fw , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_zoamp_fw , zoamp_fw , myThid )
      CALL GLOBAL_SUM_TILE_RL( tile_mass_gc  , mass_gc  , myThid )

C     finish calculating center-of-mass of oceans
C-    note: only master thread updates/modifies var in common block
      _BEGIN_MASTER(myThid)

      IF ( mass.NE.zeroRL ) THEN
        xcom = xcom / mass
        ycom = ycom / mass
        zcom = zcom / mass
      ENDIF

      IF ( mass_fw.NE.zeroRL ) THEN
        xcom_fw = xcom_fw / mass_fw
        ycom_fw = ycom_fw / mass_fw
        zcom_fw = zcom_fw / mass_fw
      ENDIF

C     Add seaice OAMC to total OAMC
      xoamc = xoamc + xoamc_si
      yoamc = yoamc + yoamc_si
      zoamc = zoamc + zoamc_si

      _END_MASTER(myThid)

#endif /* ALLOW_SBO */

      RETURN
      END
