C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_getforcing.F,v 1.49 2017/01/27 17:22:55 jmc Exp $ C $Name: $ #include "EXF_OPTIONS.h" #ifdef ALLOW_AUTODIFF # include "AUTODIFF_OPTIONS.h" #endif 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 | | Check whether read fields are within assumed range C | | (may capture mismatches in units) C | o C | C |-- EXF_RADIATION C | | Compute net or downwelling radiative fluxes via C | | Stefan-Boltzmann law in case only one is known. C |-- EXF_WIND C | | Compute air-sea wind-stress from winds (or the other way) C |-- EXF_BULKFORMULAE C | | Compute air-sea buoyancy fluxes from atmospheric C | | 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 | | If forcing and control are fluxes, then C | | control vector anomalies are added here. C | |--- ctrl_get_gen C | o C | C |-- < treatment of hflux w.r.t. swflux > C | C |-- EXF_DIAGNOSTICS_FILL C | | Do EXF-related diagnostics output here. C |-- EXF_MONITOR C | | Monitor EXF-forcing fields C | o 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 | o 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_DOWNWARD_RADIATION C Set radiative fluxes CALL EXF_RADIATION( myTime, myIter, myThid ) #endif #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_ATM_TEMP # ifdef ALLOW_BULKFORMULAE # 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 /* ALLOW_BULKFORMULAE */ #endif /* ALLOW_ATM_TEMP */ DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) #ifdef ALLOW_ATM_TEMP C compute hflux & sflux from multiple components DO j = 1,sNy DO i = 1,sNx 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) ENDDO ENDDO #endif /* ALLOW_ATM_TEMP */ C Apply runoff, masks and exchanges k = 1 DO j = 1,sNy DO i = 1,sNx #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,k,bi,bj) sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskC(i,j,k,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 # ifdef ALLOW_AUTODIFF_MONITOR CALL EXF_ADJOINT_SNAPSHOTS( 1, myTime, myIter, myThid ) # endif #endif /* ALLOW_AUTODIFF */ #ifdef ALLOW_ATM_TEMP # ifdef SHORTWAVE_HEATING C Treatment of qnet C The location of the 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 /* SHORTWAVE_HEATING */ #endif /* ALLOW_ATM_TEMP */ 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 # ifdef ALLOW_AUTODIFF_MONITOR IF ( .NOT. useSEAICE ) & CALL EXF_ADJOINT_SNAPSHOTS( 3, myTime, myIter, myThid ) # endif #endif /* ALLOW_AUTODIFF */ RETURN END