C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/ocean_inversion_project/code/ptracers_init.F,v 1.9 2004/10/06 20:26:54 dimitri Exp $
C $Name:  $

#include "PTRACERS_OPTIONS.h"

CBOP
C !ROUTINE: PTRACERS_INIT

C !INTERFACE: ==========================================================
      SUBROUTINE PTRACERS_INIT( myThid )

C !DESCRIPTION:
C     Initialize PTRACERS data structures
C     This file is customized to compute CO2 perturbations from
C     30 ocean regions for Gruber's ocean inversion project.

C !USES: ===============================================================
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "EESUPPORT.h"
#include "PARAMS.h"
#include "GRID.h"
#include "PTRACERS.h"


C !INPUT PARAMETERS: ===================================================
C  myThid               :: thread number
      INTEGER myThid

C !OUTPUT PARAMETERS: ==================================================
C  none

#ifdef ALLOW_PTRACERS

C !LOCAL VARIABLES: ====================================================
C  i,j,k,bi,bj,iTracer  :: loop indices
      INTEGER i,j,k,bi,bj,iTracer,iMonth,iUnit,iRec
      CHARACTER*(10) suff
      _RL SumPtracer

#ifdef OCEAN_INVERSION_NETCDF
      Real*8 global(Nx,Ny)
      _RL    local(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      Real lon(Nx,Ny),bounds_lon(Nx,Ny,2,2)
      Real lat(Nx,Ny),bounds_lat(Nx,Ny,2,2)
      Real depth(Nr),bounds_depth(Nr,2)
      Real MASK(Nx,Ny,Nr)
      Real AREA(Nx,Ny)
      Real BATHY(Nx,Ny)
      Real REGION_MASK(Nx,Ny,PTRACERS_num)
      Real REGION_AREA(PTRACERS_num)
#endif /* OCEAN_INVERSION_NETCDF */

CEOP
        
C Loop over tracers
      DO iTracer = 1, PTRACERS_num

C Loop over tiles
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)

C Initialize arrays in common blocks :
         DO k=1,Nr
          DO j=1-Oly,sNy+OLy
           DO i=1-Olx,sNx+Olx
            pTracer(i,j,k,bi,bj,iTracer) = 0. _d 0
            gPtr(i,j,k,bi,bj,iTracer)    = 0. _d 0
            gPtrNM1(i,j,k,bi,bj,iTracer) = 0. _d 0
           ENDDO
          ENDDO
         ENDDO
         DO j=1-Oly,sNy+OLy
          DO i=1-Olx,sNx+Olx
           surfaceTendencyPtr(i,j,bi,bj,iTracer) = 0. _d 0
          ENDDO
         ENDDO

C end bi,bj loops
        ENDDO
       ENDDO

C end of Tracer loop
      ENDDO

C Read from a pickup file if nIter0 is not zero
      IF (nIter0.NE.0) THEN
C--     Suffix for pickup files
       IF (pickupSuff.EQ.' ') THEN
         WRITE(suff,'(I10.10)') nIter0
       ELSE
         WRITE(suff,'(A10)') pickupSuff
       ENDIF
       CALL PTRACERS_READ_CHECKPOINT( nIter0,myThid )
      ENDIF

C     Initialize pTracerMasks
      CALL PTRACERS_READ_MASK ( mythid )

C     Initialize pTracerTakahashi
      CALL PTRACERS_READ_Takahashi ( mythid )

C     Normalize pTracerTakahashi so that 1e18 mol/yr is released
C     from each model region defined in pTracerMasks.
C     It is assumed that each year is 365.24167 days (31556880 s)
C     long and that each month is 2629740 s.

      DO iTracer = 1, PTRACERS_num
       SumPtracer = 0.
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           DO iMonth=1,12
            SumPtracer = SumPtracer + 2629740. * rA(I,J,bi,bj) *
     &             pTracerMasks    (I,J,iTracer,bi,bj) *
     &             pTracerTakahashi(I,J,iMonth ,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       _GLOBAL_SUM_R8( SumPtracer, myThid )
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           DO iMonth=1,12
            IF ( pTracerMasks(I,J,iTracer,bi,bj) .eq. 1. )
     &             pTracerTakahashi(I,J,iMonth ,bi,bj) = 1.e18 *
     &             pTracerTakahashi(I,J,iMonth ,bi,bj) / SumPtracer
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#ifdef OCEAN_INVERSION_PROJECT_TIME_DEPENDENT
C     Initialize atmospheric CO2 array
      call mdsfindunit( iUnit, mythid)
      open(iUnit,file='splco2_cis92a.dat',status='old',form="FORMATTED")
      do iRec=1,pTracerAtm_Nrec
         read(iUnit,*) pTracerAtmYear(iRec), pTracerAtmCO2(iRec)
cdb      print*, '###', iRec, pTracerAtmYear(iRec), pTracerAtmCO2(iRec)
      enddo
      close(iUnit)
#endif /* OCEAN_INVERSION_PROJECT_TIME_DEPENDENT */

#ifdef OCEAN_INVERSION_NETCDF

C--   lon
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          local (i,j,bi,bj) = xC (i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      call gather_2d ( global, local, myThid )
      DO j=1,Ny
       DO i=1,Nx
        if ( global(i,j) .gt. 180. ) then
         lon (i,j) = global(i,j) - 360.
        else
         lon (i,j) = global(i,j)
        endif
       ENDDO
      ENDDO

C--   bounds_lon
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          local (i,j,bi,bj) = xG (i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      call gather_2d ( global, local, myThid )
      DO j=1,Ny
       DO i=1,Nx
        if ( global(i,j) .gt. 180. ) then
         bounds_lon (i,j,1,1) = global(i,j) - 360.
         bounds_lon (i,j,1,2) = global(i,j) - 360.
        else
         bounds_lon (i,j,1,1) = global(i,j)
         bounds_lon (i,j,1,2) = global(i,j)
        endif
       ENDDO
      ENDDO
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          local (i,j,bi,bj) = xG (i+1,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      call gather_2d ( global, local, myThid )
      DO j=1,Ny
       DO i=1,Nx
        if ( global(i,j) .gt. 180. ) then
         bounds_lon (i,j,2,1) = global(i,j) - 360.
         bounds_lon (i,j,2,2) = global(i,j) - 360.
        else
         bounds_lon (i,j,2,1) = global(i,j)
         bounds_lon (i,j,2,2) = global(i,j)
        endif
       ENDDO
      ENDDO

C--   lat
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          local (i,j,bi,bj) = yC (i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      call gather_2d ( global, local, myThid )
      DO j=1,Ny
       DO i=1,Nx
        lat (i,j) = global(i,j)
       ENDDO
      ENDDO

C--   bounds_lat
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          local (i,j,bi,bj) = yG (i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      call gather_2d ( global, local, myThid )
      DO j=1,Ny
       DO i=1,Nx
        bounds_lat (i,j,1,1) = global(i,j)
        bounds_lat (i,j,2,1) = global(i,j)
       ENDDO
      ENDDO
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          local (i,j,bi,bj) = yG (i,j+1,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      call gather_2d ( global, local, myThid )
      DO j=1,Ny
       DO i=1,Nx
        bounds_lat (i,j,1,2) = global(i,j)
        bounds_lat (i,j,2,2) = global(i,j)
       ENDDO
      ENDDO

C--   AREA
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          local (i,j,bi,bj) = rA (i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      call gather_2d ( global, local, myThid )
      DO j=1,Ny
       DO i=1,Nx
        AREA (i,j) = global(i,j)
       ENDDO
      ENDDO

C--   BATHY
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          local (i,j,bi,bj) = R_low (i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      call gather_2d ( global, local, myThid )
      DO j=1,Ny
       DO i=1,Nx
        BATHY (i,j) = abs ( global(i,j) )
       ENDDO
      ENDDO

C--   MASK
      DO k=1,Nr
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           local (i,j,bi,bj) = hFacC (i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       call gather_2d ( global, local, myThid )
       DO j=1,Ny
        DO i=1,Nx
         MASK (i,j,k) = global(i,j)
        ENDDO
       ENDDO
      ENDDO

C--   REGION_MASK
      DO iTracer = 1, PTRACERS_num
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           local (i,j,bi,bj) = pTracerMasks(i,j,iTracer,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       call gather_2d ( global, local, myThid )
       DO j=1,Ny
        DO i=1,Nx
         REGION_MASK (i,j,iTracer) = global(i,j)
        ENDDO
       ENDDO
      ENDDO

C--   depth, bounds_depth
      DO k=1,Nr
       depth        (k  ) = abs ( rC (k)   )
       bounds_depth (k,1) = abs ( rF (k+1) )
       bounds_depth (k,2) = abs ( rF (k)   )
      ENDDO

      DO iTracer = 1, PTRACERS_num
       REGION_AREA (iTracer) = 0.
       DO j=1,Ny
        DO i=1,Nx
         REGION_AREA (iTracer) = REGION_AREA (iTracer) +
     &          REGION_MASK (i,j,iTracer) * AREA (i,j)
        ENDDO
       ENDDO
      ENDDO

C--   Master thread of process 0 writes netcdf output file
      _BEGIN_MASTER( myThid )
#ifdef ALLOW_USE_MPI
       IF( mpiMyId .EQ. 0 ) THEN
#endif
        call WRITE_NC_MaskAreaBathy(
     &         'ECCO','MIT GCM Release 1',
     &         Nx,Ny,Nr,PTRACERS_num,
     &         lon,bounds_lon,lat,bounds_lat,depth,bounds_depth,
     &         MASK,AREA,BATHY,REGION_MASK,REGION_AREA)
#ifdef ALLOW_USE_MPI
       ENDIF
#endif
      _END_MASTER( myThid )

#endif /* OCEAN_INVERSION_NETCDF */

#endif /* ALLOW_PTRACERS */

      RETURN
      END