C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_getforcing.F,v 1.46 2014/06/05 15:37:46 jmc Exp $
C $Name: $
#include "EXF_OPTIONS.h"
CBOI
C
C !TITLE: EXTERNAL FORCING
C !AUTHORS: mitgcm developers ( mitgcm-support@mitgcm.org )
C !AFFILIATION: Massachussetts Institute of Technology
C !DATE:
C !INTRODUCTION: External forcing package
C \bv
C * The external forcing package, in conjunction with the
C calendar package (cal), enables the handling of realistic forcing
C fields of differing temporal forcing patterns.
C * It comprises climatological restoring and relaxation
C * Bulk formulae are implemented to convert atmospheric fields
C to surface fluxes.
C * An interpolation routine provides on-the-fly interpolation of
C forcing fields an arbitrary grid onto the model grid.
C * A list of EXF variables and units is in EXF_FIELDS.h
C
C !CALLING SEQUENCE:
C ...
C exf_getforcing (TOP LEVEL ROUTINE)
C |
C |-- exf_getclim (get climatological fields used e.g. for relax.)
C | |--- exf_set_climtemp (relax. to 3-D temperature field)
C | |--- exf_set_climsalt (relax. to 3-D salinity field)
C | |--- exf_set_climsst (relax. to 2-D SST field)
C | |--- exf_set_climsss (relax. to 2-D SSS field)
C | o
C |
C |-- exf_getffields <- this one does almost everything
C | | 1. reads in fields, either flux or atmos. state,
C | | depending on CPP options (for each variable two fields
C | | consecutive in time are read in and interpolated onto
C | | current time step).
C | | 2. If forcing is atmos. state and control is atmos. state,
C | | then the control variable anomalies are read here
C | | * ctrl_getatemp
C | | * ctrl_getaqh
C | | * ctrl_getuwind
C | | * ctrl_getvwind
C | | If forcing and control are fluxes, then
C | | controls are added later.
C | o
C |
C |-- exf_check_range
C | | 1. Check whether read fields are within assumed range
C | | (may capture mismatches in units)
C | o
C |
C |-- exf_bulkformulae
C | | 1. Compute net or downwelling radiative fluxes via
C | | Stefan-Boltzmann law in case only one is known.
C | | 2. Compute air-sea momentum and buoyancy fluxes from
C | | atmospheric state following Large and Pond, JPO, 1981/82
C | o
C |
C |-- < add time-mean river runoff here, if available >
C |
C |-- < update tile edges here >
C |
C |-- exf_getsurfacefluxes
C | | 1. If forcing and control are fluxes, then
C | | controls are added here.
C | o
C |
C |-- < treatment of hflux w.r.t. swflux >
C |
C |-- exf_diagnostics_fill
C | | 1. Do EXF-related diagnostics output here.
C | o
C |
C |-- exf_mapfields
C | | 1. Map the EXF variables onto the core MITgcm
C | | forcing fields.
C | o
C |
C |-- exf_bulkformulae
C | If ALLOW_BULKFORMULAE, compute fluxes via bulkformulae
C |
C |-- exf_getsurfacefluxes
C | If forcing and control is flux, then the
C | control vector anomalies are read here
C | * ctrl_getheatflux
C | * ctrl_getsaltflux
C | * ctrl_getzonstress
C | * CALL ctrl_getmerstress
C |
C |-- exf_mapfields
C | Forcing fields from exf package are mapped onto
C | mitgcm forcing arrays.
C | Mapping enables a runtime rescaling of fields
C
C \ev
CEOI
CBOP
C !ROUTINE: EXF_GETFORCING
C !INTERFACE:
SUBROUTINE EXF_GETFORCING( myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *=================================================================
C | SUBROUTINE EXF_GETFORCING
C *=================================================================
C o Get the forcing fields for the current time step. The switches
C for the inclusion of the individual forcing components have to
C be set in EXF_OPTIONS.h (or ECCO_CPPOPTIONS.h).
C A note on surface fluxes:
C The MITgcm-UV vertical coordinate z is positive upward.
C This implies that a positive flux is out of the ocean
C model. However, the wind stress forcing is not treated
C this way. A positive zonal wind stress accelerates the
C model ocean towards the east.
C started: eckert@mit.edu, heimbach@mit.edu, ralf@ocean.mit.edu
C mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
C *=================================================================
C | SUBROUTINE EXF_GETFORCING
C *=================================================================
C \ev
C !USES:
IMPLICIT NONE
C == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#include "EXF_PARAM.h"
#include "EXF_FIELDS.h"
#include "EXF_CONSTANTS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C == routine arguments ==
_RL myTime
INTEGER myIter
INTEGER myThid
C !LOCAL VARIABLES:
C == local variables ==
INTEGER bi,bj
INTEGER i,j,k
C == end of interface ==
CEOP
C Get values of climatological fields.
CALL EXF_GETCLIM( myTime, myIter, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
# ifdef ALLOW_ATM_TEMP
CADJ STORE precip0 = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE precip1 = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE snowprecip0 = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE snowprecip1 = comlev1, key=ikey_dynamics, kind=isbyte
# endif
#endif
C Get the surface forcing fields.
CALL EXF_GETFFIELDS( myTime, myIter, myThid )
IF ( .NOT.useAtmWind ) THEN
IF ( stressIsOnCgrid .AND. ustressfile.NE.' '
& .AND. vstressfile.NE.' ' )
& CALL EXCH_UV_XY_RL( ustress, vstress, .TRUE., myThid )
ENDIF
#ifdef ALLOW_AUTODIFF_TAMC
# ifdef ALLOW_AUTODIFF_MONITOR
CALL EXF_ADJOINT_SNAPSHOTS( 2, myTime, myIter, myThid )
# endif
#endif
#ifdef ALLOW_BULKFORMULAE
C Set radiative fluxes
CALL EXF_RADIATION( myTime, myIter, myThid )
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ustress = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE vstress = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE uwind = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE vwind = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE wspeed = comlev1, key=ikey_dynamics, kind=isbyte
# endif
C Set wind fields
CALL EXF_WIND( myTime, myIter, myThid )
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ustress = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE vstress = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE uwind = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE vwind = comlev1, key=ikey_dynamics, kind=isbyte
CADJ STORE wspeed = comlev1, key=ikey_dynamics, kind=isbyte
# endif
C Compute turbulent fluxes (and surface stress) from bulk formulae
CALL EXF_BULKFORMULAE( myTime, myIter, myThid )
#endif
C Apply runoff, masks and exchanges
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
k = 1
DO j = 1,sNy
DO i = 1,sNx
#ifdef ALLOW_ATM_TEMP
C Net surface heat flux.
hflux(i,j,bi,bj) =
& - hs(i,j,bi,bj)
& - hl(i,j,bi,bj)
& + lwflux(i,j,bi,bj)
#ifndef SHORTWAVE_HEATING
& + swflux(i,j,bi,bj)
#endif
C fresh-water flux from Precipitation and Evaporation.
sflux(i,j,bi,bj) = evap(i,j,bi,bj) - precip(i,j,bi,bj)
#endif /* ALLOW_ATM_TEMP */
#ifdef ALLOW_RUNOFF
sflux(i,j,bi,bj) = sflux(i,j,bi,bj) - runoff(i,j,bi,bj)
#endif
hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskC(i,j,1,bi,bj)
sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskC(i,j,1,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
C Update the tile edges: needed for some EXF fields involved in horizontal
C averaging, e.g., wind-stress; fields used by main model or other pkgs
C are exchanged in EXF_MAPFIELDS.
c _EXCH_XY_RL(hflux, myThid)
c _EXCH_XY_RL(sflux, myThid)
IF ( stressIsOnCgrid ) THEN
CALL EXCH_UV_XY_RL( ustress, vstress, .TRUE., myThid )
ELSE
CALL EXCH_UV_AGRID_3D_RL(ustress, vstress, .TRUE., 1, myThid)
ENDIF
#ifdef SHORTWAVE_HEATING
c _EXCH_XY_RL(swflux, myThid)
#endif
#ifdef ATMOSPHERIC_LOADING
c _EXCH_XY_RL(apressure, myThid)
#endif
#ifdef EXF_SEAICE_FRACTION
c _EXCH_XY_RL(areamask, myThid)
#endif
C Get values of the surface flux anomalies.
CALL EXF_GETSURFACEFLUXES( myTime, myIter, myThid )
IF ( useExfCheckRange .AND.
& ( myIter.EQ.nIter0 .OR. exf_debugLev.GE.debLevC ) ) THEN
CALL EXF_CHECK_RANGE( myTime, myIter, myThid )
ENDIF
#ifdef ALLOW_AUTODIFF_TAMC
# ifdef ALLOW_AUTODIFF_MONITOR
CALL EXF_ADJOINT_SNAPSHOTS( 1, myTime, myIter, myThid )
# endif
#endif
#ifdef SHORTWAVE_HEATING
C Treatment of qnet
C The location of te summation of Qnet in exf_mapfields is unfortunate.
C For backward compatibility issues we want it to happen after
C applying control variables, but before exf_diagnostics_fill.
C Therefore, we DO it exactly here:
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO j = 1-oLy,sNy+oLy
DO i = 1-oLx,sNx+oLx
hflux(i,j,bi,bj) = hflux(i,j,bi,bj) + swflux(i,j,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
#endif
C Diagnostics output
CALL EXF_DIAGNOSTICS_FILL( myTime, myIter, myThid )
C Monitor output
CALL EXF_MONITOR( myTime, myIter, myThid )
C Map the forcing fields onto the corresponding model fields.
CALL EXF_MAPFIELDS( myTime, myIter, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
# ifdef ALLOW_AUTODIFF_MONITOR
IF ( .NOT. useSEAICE )
& CALL EXF_ADJOINT_SNAPSHOTS( 3, myTime, myIter, myThid )
# endif
#endif
RETURN
END