C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/ocean_inversion_project/code/ptracers_forcing_surf.F,v 1.4 2003/12/18 06:11:14 dimitri Exp $
C $Name:  $

#include "PTRACERS_OPTIONS.h"

CBOP
C !ROUTINE: PTRACERS_FORCING

C !INTERFACE: ==========================================================
      SUBROUTINE PTRACERS_FORCING_SURF(
     I                            bi, bj, iMin, iMax, jMin, jMax,
     I                            myTime,myIter,myThid )

C !DESCRIPTION:
C     Precomputes surface forcing term for pkg/ptracers.
C     Precomputation is needed because of non-local KPP transport term,
C     routine KPP_TRANSPORT_PTR.
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 "PTRACERS.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#include "DYNVARS.h"
#include "GRID.h"

C !INPUT PARAMETERS: ===================================================
C  bi,bj                :: tile indices
C  myTime               :: model time
C  myIter               :: time-step number
C  myThid               :: thread number
      INTEGER bi, bj, iMin, iMax, jMin, jMax
      _RL myTime
      INTEGER myIter
      INTEGER myThid

#ifdef ALLOW_PTRACERS

C !LOCAL VARIABLES: ====================================================
C  i,j                  :: loop indices
      INTEGER i, j, iTracer
      integer fldstartdate(4), mydate(4), difftime(4)
      integer count0,  count1
      _RL repeatPeriod, fldperiod, frac, fldsecs, fldsecs0, fldsecs1
#ifdef OCEAN_INVERSION_PROJECT_TIME_DEPENDENT
      _RL     AtmCO2frac
      integer AtmCO2rec0, AtmCO2rec1 
#endif /* OCEAN_INVERSION_PROJECT_TIME_DEPENDENT */

CEOP

C     Compute time interpolation indices and factors.
C     It is assumed that each year is 365.24167 days (31556880 s)
C     long and that each month is 2629740 s.
      repeatPeriod = 31556880
      fldperiod    = 2629740

C     Takahashi climatology start date is January 15, 5:15am
      fldstartdate(1) = 17650115
      fldstartdate(2) = 51500
      fldstartdate(3) = 1
      fldstartdate(4) = 5

C     Determine offset in seconds from beginning of input data
C     to current date.
      call cal_GetDate( myiter, mytime, mydate, mythid )
      call cal_TimePassed( fldstartdate, mydate, difftime, mythid )
      call cal_ToSeconds( difftime, fldsecs, mythid )

C     Determine the flux records just before and after mycurrentdate.
      do while ( fldsecs     .lt. 0 )
         fldsecs  = fldsecs     + repeatPeriod
      enddo
      fldsecs0    = mod(fldsecs,repeatPeriod)
      count0      = int((fldsecs0+0.5)/fldperiod) + 1
      fldsecs1    = mod(fldsecs+fldperiod,repeatPeriod)
      count1      = int((fldsecs1+0.5)/fldperiod) + 1
      fldsecs     = fldsecs0-int((fldsecs0+0.5)/fldperiod)*fldperiod

C     Weight belonging to irec for linear interpolation purposes.
C     Note: The weight as chosen here is 1. - frac of the "old"
C           MITgcm's estimation program.
      frac = 1. - fldsecs/fldperiod

#ifdef OCEAN_INVERSION_PROJECT_TIME_DEPENDENT
C     Atmospheric CO2 flux start date is January 1, 1765
      repeatPeriod = repeatPeriod / 2
      fldstartdate(1) = 17650101
      fldstartdate(2) = 00000
      fldstartdate(3) = 1
      fldstartdate(4) = 5
      call cal_GetDate( myiter, mytime, mydate, mythid )
      call cal_TimePassed( fldstartdate, mydate, difftime, mythid )
      call cal_ToSeconds( difftime, fldsecs, mythid )
      AtmCO2rec0 = int(fldsecs/repeatPeriod)+1
      AtmCO2rec1 = int(fldsecs/repeatPeriod)+2
      AtmCO2frac = ( fldsecs - repeatPeriod *
     &     int(fldsecs/repeatPeriod) ) / repeatPeriod
      AtmCO2frac = (1-AtmCO2frac) * pTracerAtmCO2(AtmCO2rec0) +
     &     AtmCO2frac * pTracerAtmCO2(AtmCO2rec1) - pTracerAtmCO2(1)
cdb   print*,'###',myiter,mytime,AtmCO2rec0,AtmCO2rec1,AtmCO2frac
#endif /* OCEAN_INVERSION_PROJECT_TIME_DEPENDENT */

      DO iTracer=1,PTRACERS_numInUse
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         surfaceTendencyPtr(i,j,bi,bj,iTracer) =
     &          pTracerMasks(i,j,iTracer,bi,bj) * (
     &          pTracerTakahashi(I,J,count0,bi,bj) * frac +
     &          pTracerTakahashi(I,J,count1,bi,bj) * (1-frac) ) *
     &          recip_drF(1) * recip_hFacC(i,j,1,bi,bj)
#ifdef OCEAN_INVERSION_PROJECT_TIME_DEPENDENT
     &        * AtmCO2frac
#endif /* OCEAN_INVERSION_PROJECT_TIME_DEPENDENT */
        ENDDO
       ENDDO
      ENDDO

#endif /* ALLOW_PTRACERS */

      RETURN
      END