C $Header: /u/gcmpack/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