C $Header: /u/gcmpack/MITgcm/pkg/cfc/cfc12_forcing.F,v 1.12 2013/06/10 02:52:57 jmc Exp $
C $Name:  $

C modified for external_forcing_DIC.F  August 1999
C
C modified swd Oct 01 and Feb 02, for use as package for c40_patch1
C modified to use with c44 and ptracers: swd May 2002
C modified to have carbonate and biological influences: swd June 2002
C modified for cfc: swd Sep 2003

#include "GCHEM_OPTIONS.h"
#define OCMIP_GRAD
#undef STEPH_GRAD

CBOP
C     !ROUTINE: CFC12_FORCING
C     !INTERFACE:
      SUBROUTINE CFC12_FORCING(
     I                          pTr_CFC12,
     U                          gCFC12,
     I                          bi, bj, iMin, iMax, jMin, jMax,
     I                          myTime, myIter, myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | SUBROUTINE CFC12_FORCING
C     | o Calculate the changes to CFC12 through air-sea  fluxes
C     *==========================================================*

C     !USES:
      IMPLICIT NONE

C     == GLobal variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "CFC.h"
#include "CFC_ATMOS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     pTr_CFC12  :: ocean CFC12 concentration
C     gCFC12     :: CFC12 tendency
C     bi, bj     :: current tile indices
C     iMin,iMax  :: computation domain, 1rst index bounds
C     jMin,jMax  :: computation domain, 2nd  index bounds
C     myTime     :: current time in simulation
C     myIter     :: current iteration number
C     myThid     :: my Thread Id number
      _RL  pTr_CFC12(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL  gCFC12(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      INTEGER bi, bj
      INTEGER iMin, iMax, jMin, jMax
      _RL  myTime
      INTEGER myIter
      INTEGER myThid
CEOP

#ifdef ALLOW_PTRACERS
#ifdef ALLOW_CFC
C     !FUNCTIONS:
      LOGICAL  DIFFERENT_MULTIPLE
      EXTERNAL DIFFERENT_MULTIPLE

C     !LOCAL VARIABLES:
C     AtmosCFC12 :: atmospheric CFC12 field
C     fluxCFC12  :: air-sea CFC12 fluxes
C     msgBuf     :: message buffer
      _RL  fluxCFC12(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL  AtmosCFC12(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER i, j
      INTEGER intimeP, intime0, intime1, iRec0, iRec1
      _RL cfcTime, aWght, bWght
      _RL ACFC12north, ACFC12south
      _RL recip_dLat, weight
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef STEPH_GRAD
       _RL a1, a2
#endif

C--   Find atmospheric CFC :
C     assume that cfcTime=0 corresponds to the beginning of the 1rst record
C     time-period. This is consistent with 1rst record value = time-averaged
C     atmos-CFC over time period: cfcTime= 0 to cfcTime= 1 x atmCFC_recSepTime
C---------------------------
       cfcTime = myTime + atmCFC_timeOffset
       CALL GET_PERIODIC_INTERVAL(
     O                   intimeP, intime0, intime1, bWght, aWght,
     I                   zeroRL, atmCFC_recSepTime,
     I                   deltaTclock, cfcTime, myThid )
       iRec0 = MAX( 1, MIN( ACFCnRec, intime0 ) )
       iRec1 = MAX( 1, MIN( ACFCnRec, intime1 ) )
       ACFC12north = ACFC12( iRec0, 1 )*bWght
     &             + ACFC12( iRec1, 1 )*aWght
       ACFC12south = ACFC12( iRec0, 2 )*bWght
     &             + ACFC12( iRec1, 2 )*aWght

C-    Print to check:
       IF ( DIFFERENT_MULTIPLE( CFC_monFreq, myTime, deltaTClock )
     &      .AND. bi*bj.EQ.1 ) THEN
         WRITE(msgBuf,'(A,6X,I10,I6,F9.4,F7.1)')
     &    'CFC12_FORCING: iter,rec0,w0,yr0 =', myIter,
     &        intime0, bWght, ACFCyear(iRec0)
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                       SQUEEZE_RIGHT, myThid )
         WRITE(msgBuf,'(A,1PE16.7,I6,0PF9.4,F7.1)')
     &    'CFC12_FORCING: cfcT,rec1,w1,yr1 =', cfcTime,
     &        intime1, aWght, ACFCyear(iRec1)
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                       SQUEEZE_RIGHT, myThid )
         WRITE(msgBuf,'(2(A,F14.6))')
     &    'CFC12_FORCING: aCFC12_N =', ACFC12north,
     &                ' , aCFC12_S =', ACFC12south
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                       SQUEEZE_RIGHT, myThid )
       ENDIF

C--   Provide gradient between N and S values
#ifdef STEPH_GRAD
C STEPH S INITIAL VERSION
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
          if ((j.gt.int(sNy/2)+3.and.j.le.sNy).or.j.lt.1) then
             AtmosCFC12(i,j)=ACFC12north
          endif
          if (j.ge.int(sNy/2)-3.and.j.le.int(sNy/2)+3) then
             a1=(float(j-int(sNy/2)+3)+.5)/7
             a2=1.d0-a1
             AtmosCFC12(i,j)=a1*ACFC12south +
     &                       a2*ACFC12north
          endif
          if ((j.lt.int(sNy/2)-3.and.j.gt.0).or.j.gt.sNy) then
             AtmosCFC12(i,j)=ACFC12south
          endif
        ENDDO
       ENDDO
#endif
#ifdef OCMIP_GRAD
C-    OCMIP VERSION
C     between N & S lat boundaries, do linear interpolation ; and
C     beyond N or S lat boundaries, just take the hemispheric value
       recip_dLat = 1. _d 0 / ( atmCFC_yNorthBnd - atmCFC_ySouthBnd )
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
          weight = ( yC(i,j,bi,bj) - atmCFC_ySouthBnd )*recip_dLat
          weight = MAX( zeroRL, MIN( oneRL, weight ) )
          AtmosCFC12(i,j)= weight * ACFC12north
     &         + ( oneRL - weight )*ACFC12south

        ENDDO
c         print*,'QQ cfc12', j, ATMOSCFC12(1,j,bi,bj)
       ENDDO
#endif
C--   cfc12 air-sea fluxes
       CALL CFC12_SURFFORCING(
     I                    pTr_CFC12, AtmosCFC12,
     O                    fluxCFC12,
     I                    bi, bj, iMin, iMax, jMin, jMax,
     I                    myTime, myIter, myThid )

C--   update surface tendencies
       DO j=jMin,jMax
        DO i=iMin,iMax
          gCFC12(i,j,1) = gCFC12(i,j,1)
c    &     + fluxCFC12(i,j)*recip_drF(1)*maskC(i,j,1,bi,bj)
     &     + fluxCFC12(i,j)*recip_drF(1)*recip_hFacC(i,j,1,bi,bj)
        ENDDO
       ENDDO

#endif /* ALLOW_CFC */
#endif /* ALLOW_PTRACERS */

       RETURN
       END