C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/ocean_inversion_project/code_bombC14/ini_tr1.F,v 1.3 2006/05/25 20:57:50 dimitri Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_TR1
C     !INTERFACE:
      SUBROUTINE INI_TR1( myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_TR1                                        
C     | o Set initial tracer 1 distribution.                      
C     *==========================================================*
C     | Passive tracers 1-N can be initialised so as to dye fluid.
C     | This routine is hardcoded for Nir's bomb C14 experiment.
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"
#ifdef ALLOW_PASSIVE_TRACER
#include "TR1.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid -  Number of this instance of INI_TR1
      INTEGER myThid

#ifdef ALLOW_PASSIVE_TRACER

C     !LOCAL VARIABLES:
C     == Local variables ==
C     iC, jC - Center of domain
C     iD, jD - Disitance from domain center.
C     rad    - Radius of initial patch
C     rD     - Radial displacement of point I,J
C     iG, jG - Global coordinate index
C     bi,bj  - Loop counters
C     I,J,K
      INTEGER iC, jC, iD, jD
      INTEGER iG, jG
      INTEGER bi, bj
      INTEGER I, J, K, localWarnings
      _RL     rad, rD
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

      _BARRIER

C--   Initialise tracer to inline.
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO K=1,Nr
         DO J=1,sNy
          DO I=1,sNx
            tr1(I,J,K,bi,bj) = 0. _d 0
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C     Set initial tendency terms
      localWarnings=0
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO K=1,Nr
         DO J=1,sNy
          DO I=1,sNx
           gtr1   (I,J,K,bi,bj) = 0. _d 0
           gtr1NM1(I,J,K,bi,bj) = 0. _d 0
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C
      _EXCH_XYZ_R8(tr1 , myThid )
      _EXCH_XYZ_R8(gtr1 , myThid )
      _EXCH_XYZ_R8(gtr1NM1 , myThid )

C--   Input file names
      FiceFile = 'FICE'
      XkwFile  = 'XKW'
      PatmFile = 'PATM'
      WOASssFile  = 'WOA01_S'
      WOASstFile  = 'WOA01_T'
      SurfDicPreindFile = 'surf_dic_preind'
      SurfDicAnthroFile = 'surf_dic_anthro'
      SurfNatD14CFile   = 'surface_naturalD14C'

C--   Compute numbers needed to determine DC14 and CO2_atm
      CALL read_dc14
      CALL read_CO2

      _BEGIN_MASTER(myThid)
#define OCMIP_USE_INTERP
C-- OCMIP_USE_INTERP option is useful for high-resolution integrations.
C   For low-resolution, less that 1-deg, it's best to generate files
C   separately because of sea-ice fraction.
C   Caution: OCMIP_USE_INTERP as used is not thread-safe.
#ifdef OCMIP_USE_INTERP
            CALL EXF_INTERP( SurfDicPreindFile,readBinaryPrec,
     &           SurfDicPreind,1,
     &           xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )
            CALL EXF_INTERP( SurfDicAnthroFile,readBinaryPrec,
     &           SurfDicAnthro,1,
     &           xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )
            CALL EXF_INTERP( SurfNatD14CFile,readBinaryPrec,
     &           SurfNatD14C,1,
     &           xC,yC,lon0,lon_inc,lat0,lat_inc,nlon,nlat,1,mythid )

#else
            CALL MDSREADFIELD ( SurfDicPreindFile, readBinaryPrec,
     &           'RS', 1, SurfDicPreind, 1,, myThid )
            CALL MDSREADFIELD ( SurfDicAnthroFile, readBinaryPrec,
     &           'RS', 1, SurfDicAnthro, 1,, myThid )
            CALL MDSREADFIELD ( SurfNatD14CFile, readBinaryPrec,
     &           'RS', 1, SurfNatD14C, 1,, myThid )
#endif
            _END_MASTER(myThid)
            _EXCH_XY_R4(SurfDicPreind , myThid )
            _EXCH_XY_R4(SurfDicAnthro , myThid )
            _EXCH_XY_R4(SurfNatD14C , myThid )

#endif /*  ALLOW_PASSIVE_TRACER */

      RETURN
      END

C====================================================================

      SUBROUTINE read_dc14

C     --------------------------------------------------------------------
C     Reads temporal history of atmospheric DC14 in permil
C
C     DC14_atm_year = year of measurement (mid-year)
C     DC14_atm_n    = DC14 in northern hemisphere
C     DC14_atm_eq   = DC14 at equator
C     DC14_atm_s    = DC14 in southern hemisphere
C     --------------------------------------------------------------------

      IMPLICIT NONE 
#include "SIZE.h"
#include "EEPARAMS.h"
#include "TR1.h"

      INTEGER n, iu
      PARAMETER (iu=10)

      CHARACTER*80 filen
C
C     OPEN FILE
C     ---------
      filen='DC14_atm'
      OPEN(UNIT=iu,FILE=filen,STATUS='old')
C
C     READ FILE
C     ---------
      DO n = 1, 240
        READ(iu,*)
     +        DC14_atm_year(n), DC14_atm_n(n), 
     +        DC14_atm_eq  (n), DC14_atm_s(n)
C       WRITE(6,100)
C    +        DC14_atm_year(n), DC14_atm_n(n), 
C    +        DC14_atm_eq  (n), DC14_atm_s(n)
      END DO 

 100  FORMAT(f7.2, 4f8.2)

      RETURN
      END 

C====================================================================

      SUBROUTINE read_CO2

C     --------------------------------------------------------------------
C     Reads temporal history of atmospheric CO2 in ppm
C
C     CO2_atm_year = year of measurement
C     CO2_atm      = atmospheric CO2
C     --------------------------------------------------------------------

      IMPLICIT NONE 
#include "SIZE.h"
#include "EEPARAMS.h"
#include "TR1.h"

      INTEGER n, iu
      PARAMETER (iu=10)

      CHARACTER*80 filen
C
C     OPEN FILE
C     ---------
      filen='CO2_atm'
      OPEN(UNIT=iu,FILE=filen,STATUS='old')
C
C     READ FILE
C     ---------
      DO n = 1, 478
        READ(iu,*)   CO2_atm_year(n), CO2_atm(n)
C       WRITE(6,100) CO2_atm_year(n), CO2_atm(n)
      END DO 

 100  FORMAT(f7.2, f8.2)

      RETURN
      END