C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/exf/exf_mapfields.F,v 1.33 2017/10/06 00:03:56 jmc Exp $
C $Name:  $

#include "EXF_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif

CBOP 0
C     !ROUTINE: EXF_MAPFIELDS

C     !INTERFACE:
      SUBROUTINE EXF_MAPFIELDS( myTime, myIter, myThid )

C     !DESCRIPTION:
C     ==================================================================
C     SUBROUTINE EXF_MAPFIELDS
C     ==================================================================
C
C     o Map external forcing fields (ustress, vstress, hflux, sflux,
C       swflux, apressure, climsss, climsst, etc.) onto ocean model
C       arrays (fu, fv, Qnet, EmPmR, Qsw, pLoad, SSS, SST, etc.).
C       This routine is included to separate the ocean state estimation
C       tool as much as possible from the ocean model.  Unit and sign
C       conventions can be customized using variables exf_outscal_*,
C       which are set in exf_readparms.F.  See the header files
C       EXF_FIELDS.h and FFIELDS.h for definitions of the various input
C       and output fields and for default unit and sign convetions.
C
C     started: Christian Eckert eckert@mit.edu  09-Aug-1999
C
C     changed: Christian Eckert eckert@mit.edu  11-Jan-2000
C              - Restructured the code in order to create a package
C                for the MITgcmUV.
C
C              Christian Eckert eckert@mit.edu  12-Feb-2000
C              - Changed Routine names (package prefix: exf_)
C
C              Patrick Heimbach, heimbach@mit.edu  06-May-2000
C              - added and changed CPP flag structure for
C                ALLOW_BULKFORMULAE, ALLOW_ATM_TEMP
C
C              Patrick Heimbach, heimbach@mit.edu  23-May-2000
C              - sign change of ustress/vstress incorporated into
C                scaling factors exf_outscal_ust, exf_outscal_vst
C
C     mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
C
C     ==================================================================
C     SUBROUTINE EXF_MAPFIELDS
C     ==================================================================

C     !USES:
      IMPLICIT NONE

C     == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#include "GRID.h"
#include "DYNVARS.h"

#include "EXF_PARAM.h"
#include "EXF_CONSTANTS.h"
#include "EXF_FIELDS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     myTime  :: Current time in simulation
C     myIter  :: Current iteration number
C     myThid  :: my Thread Id number
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

C     !LOCAL VARIABLES:
      INTEGER bi,bj
      INTEGER i,j,ks
      INTEGER imin, imax
      INTEGER jmin, jmax
      PARAMETER ( imin = 1-OLx , imax = sNx+OLx )
      PARAMETER ( jmin = 1-OLy , jmax = sNy+OLy )
CEOP

C--   set surface level index:
      ks = 1

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

#ifdef ALLOW_AUTODIFF_TAMC
          act1 = bi - myBxLo(myThid)
          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
          act2 = bj - myByLo(myThid)
          max2 = myByHi(myThid) - myByLo(myThid) + 1
          act3 = myThid - 1
          max3 = nTx*nTy
          act4 = ikey_dynamics - 1
          ikey = (act1 + 1) + act2*max1
     &                      + act3*max1*max2
     &                      + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

C     Heat flux.
          DO j = jmin,jmax
            DO i = imin,imax
             Qnet(i,j,bi,bj) = exf_outscal_hflux*hflux(i,j,bi,bj)
            ENDDO
          ENDDO
          IF ( hfluxfile .EQ. ' ' ) THEN
           DO j = jmin,jmax
            DO i = imin,imax
             Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) -
     &            exf_outscal_hflux * ( hflux_exfremo_intercept +
     &            hflux_exfremo_slope*(myTime-startTime) )
            ENDDO
           ENDDO
          ENDIF

C     Freshwater flux.
          DO j = jmin,jmax
            DO i = imin,imax
             EmPmR(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)
     &                                          *rhoConstFresh
            ENDDO
          ENDDO
          IF ( sfluxfile .EQ. ' ' ) THEN
           DO j = jmin,jmax
            DO i = imin,imax
             EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj) - rhoConstFresh*
     &            exf_outscal_sflux * ( sflux_exfremo_intercept +
     &            sflux_exfremo_slope*(myTime-startTime) )
            ENDDO
           ENDDO
          ENDIF

#ifdef ALLOW_ATM_TEMP
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
          IF ( temp_EvPrRn .NE. UNSET_RL ) THEN
C--   Account for energy content of Precip + RunOff & Evap. Assumes:
C     1) Rain has same temp as Air
C     2) Snow has no heat capacity (consistent with seaice & thsice pkgs)
C     3) No distinction between sea-water Cp and fresh-water Cp
C     4) By default, RunOff comes at the temp of surface water (with same Cp);
C        ifdef ALLOW_RUNOFTEMP, RunOff temp can be specified in runoftempfile.
C     5) Evap is released to the Atmos @ surf-temp (=SST); should be using
C        the water-vapor heat capacity here and consistently in Bulk-Formulae;
C        Could also be put directly into Latent Heat flux.
           IF ( snowPrecipFile .NE. ' ' ) THEN
C--   Melt snow (if provided) into the ocean and account for rain-temp
            DO j = 1, sNy
             DO i = 1, sNx
              Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
     &              + flami*snowPrecip(i,j,bi,bj)*rhoConstFresh
     &              - HeatCapacity_Cp
     &               *( atemp(i,j,bi,bj) - cen2kel - temp_EvPrRn )
     &               *( precip(i,j,bi,bj)- snowPrecip(i,j,bi,bj) )
     &               *rhoConstFresh
             ENDDO
            ENDDO
           ELSE
C--   Make snow (according to Air Temp) and melt it in the ocean
C     note: here we just use the same criteria as over seaice but would be
C           better to consider a higher altitude air temp, e.g., 850.mb
            DO j = 1, sNy
             DO i = 1, sNx
              IF ( atemp(i,j,bi,bj).LT.cen2kel ) THEN
               Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
     &              + flami*precip(i,j,bi,bj)*rhoConstFresh
              ELSE
C--   Account for rain-temp
               Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
     &              - HeatCapacity_Cp
     &               *( atemp(i,j,bi,bj) - cen2kel - temp_EvPrRn )
     &               *precip(i,j,bi,bj)*rhoConstFresh
              ENDIF
             ENDDO
            ENDDO
           ENDIF
#ifdef ALLOW_RUNOFF
C--   Account for energy content of RunOff:
           DO j = 1, sNy
            DO i = 1, sNx
              Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
     &              - HeatCapacity_Cp
     &               *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
     &               *runoff(i,j,bi,bj)*rhoConstFresh
            ENDDO
           ENDDO
#endif
C--   Account for energy content of Evap:
           DO j = 1, sNy
            DO i = 1, sNx
              Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
     &              + HeatCapacity_Cp
     &               *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
     &               *evap(i,j,bi,bj)*rhoConstFresh
              Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)*maskC(i,j,ks,bi,bj)
            ENDDO
           ENDDO
          ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#endif /* ALLOW_ATM_TEMP */
#if defined(ALLOW_RUNOFF) && defined(ALLOW_RUNOFTEMP)
          IF ( runoftempfile .NE. ' ' ) THEN
C--   Add energy content of RunOff
           DO j = 1, sNy
            DO i = 1, sNx
               Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
     &              + HeatCapacity_Cp
     &              *( theta(i,j,ks,bi,bj) - runoftemp(i,j,bi,bj) )
     &              *runoff(i,j,bi,bj)*rhoConstFresh
            ENDDO
           ENDDO
          ENDIF
#endif

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
#endif
          DO j = jmin,jmax
            DO i = imin,imax
C             Zonal wind stress.
              IF (ustress(i,j,bi,bj).GT.windstressmax) THEN
                ustress(i,j,bi,bj)=windstressmax
              ENDIF
            ENDDO
          ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
#endif
          DO j = jmin,jmax
            DO i = imin,imax
              IF (ustress(i,j,bi,bj).LT.-windstressmax) THEN
                ustress(i,j,bi,bj)=-windstressmax
              ENDIF
            ENDDO
          ENDDO
          IF ( stressIsOnCgrid ) THEN
           DO j = jmin,jmax
            DO i = imin+1,imax
              fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)
            ENDDO
           ENDDO
          ELSE
           DO j = jmin,jmax
            DO i = imin+1,imax
C     Shift wind stresses calculated at Grid-center to W/S points
              fu(i,j,bi,bj) = exf_outscal_ustress*
     &              (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))
     &              *exf_half*maskW(i,j,ks,bi,bj)
            ENDDO
           ENDDO
          ENDIF

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
#endif
          DO j = jmin,jmax
            DO i = imin,imax
C             Meridional wind stress.
              IF (vstress(i,j,bi,bj).GT.windstressmax) THEN
                vstress(i,j,bi,bj)=windstressmax
              ENDIF
            ENDDO
          ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
#endif
          DO j = jmin,jmax
            DO i = imin,imax
              IF (vstress(i,j,bi,bj).LT.-windstressmax) THEN
                vstress(i,j,bi,bj)=-windstressmax
              ENDIF
            ENDDO
          ENDDO
          IF ( stressIsOnCgrid ) THEN
           DO j = jmin+1,jmax
            DO i = imin,imax
              fv(i,j,bi,bj) = exf_outscal_vstress*vstress(i,j,bi,bj)
            ENDDO
           ENDDO
          ELSE
           DO j = jmin+1,jmax
            DO i = imin,imax
C     Shift wind stresses calculated at C-points to W/S points
              fv(i,j,bi,bj) = exf_outscal_vstress*
     &              (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))
     &              *exf_half*maskS(i,j,ks,bi,bj)
            ENDDO
           ENDDO
          ENDIF

#if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
C             Short wave radiative flux.
          DO j = jmin,jmax
            DO i = imin,imax
             Qsw(i,j,bi,bj)  = exf_outscal_swflux*swflux(i,j,bi,bj)
            ENDDO
          ENDDO
#endif

#ifdef ALLOW_CLIMSST_RELAXATION
          DO j = jmin,jmax
            DO i = imin,imax
             SST(i,j,bi,bj)  = exf_outscal_sst*climsst(i,j,bi,bj)
            ENDDO
          ENDDO
#endif

#ifdef ALLOW_CLIMSSS_RELAXATION
          DO j = jmin,jmax
            DO i = imin,imax
             SSS(i,j,bi,bj)  = exf_outscal_sss*climsss(i,j,bi,bj)
            ENDDO
          ENDDO
#endif

#ifdef ATMOSPHERIC_LOADING
          DO j = jmin,jmax
            DO i = imin,imax
             pLoad(i,j,bi,bj)=exf_outscal_apressure*apressure(i,j,bi,bj)
            ENDDO
          ENDDO
#endif

#ifdef EXF_ALLOW_TIDES
          DO j = jmin,jmax
            DO i = imin,imax
             phiTide2d(i,j,bi,bj)=exf_outscal_tidePot*tidePot(i,j,bi,bj)
            ENDDO
          ENDDO
#endif /* EXF_ALLOW_TIDES */

#ifdef ALLOW_SALTFLX
          DO j = jmin,jmax
            DO i = imin,imax
              saltFlux(I,J,bi,bj) = saltflx(I,J,bi,bj)
            ENDDO
          ENDDO
#endif

#ifdef EXF_SEAICE_FRACTION
          DO j = jmin,jmax
            DO i = imin,imax
              exf_iceFraction(i,j,bi,bj) =
     &           exf_outscal_areamask*areamask(i,j,bi,bj)
              exf_iceFraction(I,J,bi,bj) =
     &           MIN( MAX(exf_iceFraction(I,J,bi,bj),zeroRS), oneRS )
            ENDDO
          ENDDO
#endif

       ENDDO
      ENDDO

C--   Update the tile edges.
      _EXCH_XY_RS(  Qnet, myThid )
      _EXCH_XY_RS( EmPmR, myThid )
       CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
c#if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
#ifdef SHORTWAVE_HEATING
C     Qsw used in SHORTWAVE_HEATING code & for diagnostics (<- EXCH not needed)
      _EXCH_XY_RS(   Qsw, myThid )
#endif
#ifdef ALLOW_CLIMSST_RELAXATION
      _EXCH_XY_RS(   SST, myThid )
#endif
#ifdef ALLOW_CLIMSSS_RELAXATION
      _EXCH_XY_RS(   SSS, myThid )
#endif
#ifdef ATMOSPHERIC_LOADING
      _EXCH_XY_RS( pLoad, myThid )
#endif
#ifdef EXF_ALLOW_TIDES
      _EXCH_XY_RS( phiTide2d, myThid )
#endif
#ifdef EXF_SEAICE_FRACTION
      _EXCH_XY_RS( exf_iceFraction, myThid )
#endif

      RETURN
      END