Back to home page

MITgcm
 
 

    

File indexing completed on 2023-05-28 05:10:13 UTC

view on githubraw file Latest commit b4daa243 on 2023-05-28 03:53:22 UTC
6d54cf9ca1 Ed H*0001 #include "EXF_OPTIONS.h"
3a255f48df Gael*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
7f861c1808 Patr*0005 
9fe4461eea Patr*0006 CBOI
                0007 C
                0008 C !TITLE: EXTERNAL FORCING
14a20b48a3 Jean*0009 C !AUTHORS: mitgcm developers ( mitgcm-support@mitgcm.org )
9fe4461eea Patr*0010 C !AFFILIATION: Massachussetts Institute of Technology
                0011 C !DATE:
                0012 C !INTRODUCTION: External forcing package
14a20b48a3 Jean*0013 C \bv
                0014 C * The external forcing package, in conjunction with the
                0015 C   calendar package (cal), enables the handling of realistic forcing
                0016 C   fields of differing temporal forcing patterns.
                0017 C * It comprises climatological restoring and relaxation
                0018 C * Bulk formulae are implemented to convert atmospheric fields
                0019 C   to surface fluxes.
                0020 C * An interpolation routine provides on-the-fly interpolation of
                0021 C   forcing fields an arbitrary grid onto the model grid.
                0022 C * A list of EXF variables and units is in EXF_FIELDS.h
                0023 C
9fe4461eea Patr*0024 C     !CALLING SEQUENCE:
14a20b48a3 Jean*0025 C ...
423768d890 Jean*0026 C  EXF_GETFORCING (TOP LEVEL ROUTINE)
14a20b48a3 Jean*0027 C  |
423768d890 Jean*0028 C  |-- EXF_GETCLIM (get climatological fields used e.g. for relax.)
14a20b48a3 Jean*0029 C  |   |--- exf_set_climtemp (relax. to 3-D temperature field)
                0030 C  |   |--- exf_set_climsalt (relax. to 3-D salinity field)
                0031 C  |   |--- exf_set_climsst  (relax. to 2-D SST field)
                0032 C  |   |--- exf_set_climsss  (relax. to 2-D SSS field)
                0033 C  |   o
                0034 C  |
423768d890 Jean*0035 C  |-- EXF_GETFFIELDS <- this one does almost everything
14a20b48a3 Jean*0036 C  |   |   1. reads in fields, either flux or atmos. state,
                0037 C  |   |      depending on CPP options (for each variable two fields
                0038 C  |   |      consecutive in time are read in and interpolated onto
                0039 C  |   |      current time step).
                0040 C  |   |   2. If forcing is atmos. state and control is atmos. state,
                0041 C  |   |      then the control variable anomalies are read here
                0042 C  |   |          * ctrl_getatemp
                0043 C  |   |          * ctrl_getaqh
                0044 C  |   |          * ctrl_getuwind
                0045 C  |   |          * ctrl_getvwind
                0046 C  |   |      If forcing and control are fluxes, then
                0047 C  |   |      controls are added later.
                0048 C  |   o
                0049 C  |
423768d890 Jean*0050 C  |-- EXF_CHECK_RANGE
                0051 C  |   |   Check whether read fields are within assumed range
                0052 C  |   |   (may capture mismatches in units)
14a20b48a3 Jean*0053 C  |   o
                0054 C  |
423768d890 Jean*0055 C  |-- EXF_RADIATION
                0056 C  |   |   Compute net or downwelling radiative fluxes via
                0057 C  |   |   Stefan-Boltzmann law in case only one is known.
                0058 C  |-- EXF_WIND
                0059 C  |   |   Compute air-sea wind-stress from winds (or the other way)
                0060 C  |-- EXF_BULKFORMULAE
                0061 C  |   |   Compute air-sea buoyancy fluxes from atmospheric
                0062 C  |   |   state following Large and Pond, JPO, 1981/82
14a20b48a3 Jean*0063 C  |   o
                0064 C  |
                0065 C  |-- < add time-mean river runoff here, if available >
                0066 C  |
                0067 C  |-- < update tile edges here >
                0068 C  |
423768d890 Jean*0069 C  |-- EXF_GETSURFACEFLUXES
                0070 C  |   |   If forcing and control are fluxes, then
                0071 C  |   |   control vector anomalies are added here.
                0072 C  |   |--- ctrl_get_gen
14a20b48a3 Jean*0073 C  |   o
                0074 C  |
                0075 C  |-- < treatment of hflux w.r.t. swflux >
                0076 C  |
423768d890 Jean*0077 C  |-- EXF_DIAGNOSTICS_FILL
                0078 C  |   |   Do EXF-related diagnostics output here.
                0079 C  |-- EXF_MONITOR
                0080 C  |   |   Monitor EXF-forcing fields
14a20b48a3 Jean*0081 C  |   o
                0082 C  |
423768d890 Jean*0083 C  |-- EXF_MAPFIELDS
                0084 C  |   |   Forcing fields from exf package are mapped onto
                0085 C  |   |   mitgcm forcing arrays.
                0086 C  |   |   Mapping enables a runtime rescaling of fields
14a20b48a3 Jean*0087 C  |   o
                0088 C
                0089 C \ev
9fe4461eea Patr*0090 CEOI
                0091 
                0092 CBOP
14a20b48a3 Jean*0093 C     !ROUTINE: EXF_GETFORCING
9fe4461eea Patr*0094 C     !INTERFACE:
14a20b48a3 Jean*0095       SUBROUTINE EXF_GETFORCING( myTime, myIter, myThid )
9fe4461eea Patr*0096 
                0097 C     !DESCRIPTION: \bv
14a20b48a3 Jean*0098 C     *=================================================================
                0099 C     | SUBROUTINE EXF_GETFORCING
                0100 C     *=================================================================
                0101 C     o Get the forcing fields for the current time step. The switches
                0102 C       for the inclusion of the individual forcing components have to
                0103 C       be set in EXF_OPTIONS.h (or ECCO_CPPOPTIONS.h).
                0104 C       A note on surface fluxes:
                0105 C       The MITgcm-UV vertical coordinate z is positive upward.
                0106 C       This implies that a positive flux is out of the ocean
                0107 C       model. However, the wind stress forcing is not treated
                0108 C       this way. A positive zonal wind stress accelerates the
                0109 C       model ocean towards the east.
                0110 C       started: eckert@mit.edu, heimbach@mit.edu, ralf@ocean.mit.edu
                0111 C       mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
                0112 C     *=================================================================
                0113 C     | SUBROUTINE EXF_GETFORCING
                0114 C     *=================================================================
9fe4461eea Patr*0115 C     \ev
7f861c1808 Patr*0116 
9fe4461eea Patr*0117 C     !USES:
14a20b48a3 Jean*0118       IMPLICIT NONE
7f861c1808 Patr*0119 
14a20b48a3 Jean*0120 C     == global variables ==
bdec91d862 Patr*0121 #include "EEPARAMS.h"
                0122 #include "SIZE.h"
                0123 #include "PARAMS.h"
                0124 #include "GRID.h"
                0125 
082e18c36c Jean*0126 #include "EXF_PARAM.h"
                0127 #include "EXF_FIELDS.h"
                0128 #include "EXF_CONSTANTS.h"
a0a3896567 Patr*0129 #ifdef ALLOW_AUTODIFF_TAMC
                0130 # include "tamc.h"
                0131 #endif
b4daa24319 Shre*0132 #ifdef ALLOW_TAPENADE
                0133 # include "EXF_INTERP_SIZE.h"
                0134 # include "EXF_INTERP_PARAM.h"
                0135 #endif /* ALLOW_TAPENADE */
bdec91d862 Patr*0136 
9fe4461eea Patr*0137 C     !INPUT/OUTPUT PARAMETERS:
14a20b48a3 Jean*0138 C     == routine arguments ==
                0139       _RL     myTime
                0140       INTEGER myIter
                0141       INTEGER myThid
7f861c1808 Patr*0142 
9fe4461eea Patr*0143 C     !LOCAL VARIABLES:
14a20b48a3 Jean*0144 C     == local variables ==
                0145       INTEGER bi,bj
0320e25227 Mart*0146       INTEGER i,j,ks
14a20b48a3 Jean*0147 C     == end of interface ==
9fe4461eea Patr*0148 CEOP
7f861c1808 Patr*0149 
14a20b48a3 Jean*0150 C     Get values of climatological fields.
423768d890 Jean*0151       CALL EXF_GETCLIM( myTime, myIter, myThid )
7f861c1808 Patr*0152 
14a20b48a3 Jean*0153 C     Get the surface forcing fields.
423768d890 Jean*0154       CALL EXF_GETFFIELDS( myTime, myIter, myThid )
358649780a Gael*0155       IF ( .NOT.useAtmWind ) THEN
423768d890 Jean*0156        IF ( stressIsOnCgrid .AND. ustressfile.NE.' '
                0157      &                      .AND. vstressfile.NE.' ' )
9c3e24f78c Jean*0158      &  CALL EXCH_UV_XY_RL( ustress, vstress, .TRUE., myThid )
358649780a Gael*0159       ENDIF
7f861c1808 Patr*0160 
98238efc8b Patr*0161 #ifdef ALLOW_AUTODIFF_TAMC
9c41af81f6 Timo*0162 C     Store fields after reading them so that we do not need to save
                0163 C     their 0/1 levels to comlev1. Not all fields are required here (in
                0164 C     most cases only u/vwind or u/vstress, aqh, atemp, precip,
                0165 C     snowprecip, runoff), but we have directives for all potential
                0166 C     candidates here.
                0167 CADJ STORE ustress      = comlev1, key=ikey_dynamics, kind=isbyte
                0168 CADJ STORE vstress      = comlev1, key=ikey_dynamics, kind=isbyte
                0169 CADJ STORE uwind        = comlev1, key=ikey_dynamics, kind=isbyte
                0170 CADJ STORE vwind        = comlev1, key=ikey_dynamics, kind=isbyte
                0171 CADJ STORE wspeed       = comlev1, key=ikey_dynamics, kind=isbyte
                0172 # ifdef ALLOW_ATM_TEMP
                0173 CADJ STORE aqh        = comlev1, key=ikey_dynamics, kind=isbyte
                0174 CADJ STORE atemp      = comlev1, key=ikey_dynamics, kind=isbyte
                0175 CADJ STORE precip     = comlev1, key=ikey_dynamics, kind=isbyte
                0176 CADJ STORE lwflux     = comlev1, key=ikey_dynamics, kind=isbyte
                0177 CADJ STORE swflux     = comlev1, key=ikey_dynamics, kind=isbyte
                0178 CADJ STORE snowprecip = comlev1, key=ikey_dynamics, kind=isbyte
                0179 #  ifdef ALLOW_READ_TURBFLUXES
                0180 CADJ STORE hs         = comlev1, key=ikey_dynamics, kind=isbyte
                0181 CADJ STORE hl         = comlev1, key=ikey_dynamics, kind=isbyte
                0182 #  endif /* ALLOW_READ_TURBFLUXES */
                0183 #  ifdef EXF_READ_EVAP
                0184 CADJ STORE evap       = comlev1, key=ikey_dynamics, kind=isbyte
                0185 #  endif /* EXF_READ_EVAP */
                0186 #  ifdef ALLOW_DOWNWARD_RADIATION
                0187 CADJ STORE swdown     = comlev1, key=ikey_dynamics, kind=isbyte
                0188 CADJ STORE lwdown     = comlev1, key=ikey_dynamics, kind=isbyte
                0189 #  endif
                0190 # else /* ALLOW_ATM_TEMP undef */
                0191 #  ifdef SHORTWAVE_HEATING
                0192 CADJ STORE swflux     = comlev1, key=ikey_dynamics, kind=isbyte
                0193 #  endif
                0194 # endif /* ALLOW_ATM_TEMP */
                0195 # ifdef ATMOSPHERIC_LOADING
                0196 CADJ STORE apressure  = comlev1, key=ikey_dynamics, kind=isbyte
                0197 # endif
                0198 # ifdef ALLOW_RUNOFF
                0199 CADJ STORE runoff     = comlev1, key=ikey_dynamics, kind=isbyte
                0200 # endif
                0201 # ifdef ALLOW_SALTFLX
                0202 CADJ STORE saltflx    = comlev1, key=ikey_dynamics, kind=isbyte
                0203 # endif
                0204 # ifdef EXF_SEAICE_FRACTION
                0205 CADJ STORE areamask   = comlev1, key=ikey_dynamics, kind=isbyte
                0206 # endif
                0207 #endif /* ALLOW_AUTODIFF_TAMC */
                0208 
5d96c49326 Jean*0209 #ifdef ALLOW_AUTODIFF
14a20b48a3 Jean*0210 # ifdef ALLOW_AUTODIFF_MONITOR
98238efc8b Patr*0211         CALL EXF_ADJOINT_SNAPSHOTS( 2, myTime, myIter, myThid )
                0212 # endif
5d96c49326 Jean*0213 #endif /* ALLOW_AUTODIFF */
98238efc8b Patr*0214 
423768d890 Jean*0215 #ifdef ALLOW_DOWNWARD_RADIATION
14a20b48a3 Jean*0216 C     Set radiative fluxes
423768d890 Jean*0217       CALL EXF_RADIATION( myTime, myIter, myThid )
                0218 #endif
3752238fd8 Patr*0219 
14a20b48a3 Jean*0220 C     Set wind fields
423768d890 Jean*0221       CALL EXF_WIND( myTime, myIter, myThid )
                0222 
                0223 #ifdef ALLOW_ATM_TEMP
                0224 # ifdef ALLOW_BULKFORMULAE
                0225 #  ifdef ALLOW_AUTODIFF_TAMC
9c41af81f6 Timo*0226 C     Here we probably only need to store uwind, vwind, wstress but we
                0227 C     keep the other fields for the paranoid AD-modeller
358649780a Gael*0228 CADJ STORE uwind        = comlev1, key=ikey_dynamics, kind=isbyte
                0229 CADJ STORE vwind        = comlev1, key=ikey_dynamics, kind=isbyte
                0230 CADJ STORE wspeed       = comlev1, key=ikey_dynamics, kind=isbyte
9c41af81f6 Timo*0231 CADJ STORE ustress      = comlev1, key=ikey_dynamics, kind=isbyte
                0232 CADJ STORE vstress      = comlev1, key=ikey_dynamics, kind=isbyte
                0233 CADJ STORE wstress      = comlev1, key=ikey_dynamics, kind=isbyte
423768d890 Jean*0234 #  endif
14a20b48a3 Jean*0235 C     Compute turbulent fluxes (and surface stress) from bulk formulae
423768d890 Jean*0236       CALL EXF_BULKFORMULAE( myTime, myIter, myThid )
9c41af81f6 Timo*0237 #  ifdef ALLOW_AUTODIFF_TAMC
                0238 CADJ STORE evap         = comlev1, key=ikey_dynamics, kind=isbyte
                0239 #  endif
423768d890 Jean*0240 # endif /* ALLOW_BULKFORMULAE */
                0241 #endif /* ALLOW_ATM_TEMP */
bdec91d862 Patr*0242 
14a20b48a3 Jean*0243       DO bj = myByLo(myThid), myByHi(myThid)
                0244        DO bi = myBxLo(myThid), myBxHi(myThid)
423768d890 Jean*0245 
                0246 #ifdef ALLOW_ATM_TEMP
                0247 C     compute hflux & sflux from multiple components
0320e25227 Mart*0248         DO j = 1,sNy
                0249          DO i = 1,sNx
                0250 C     Net surface heat flux.
                0251           hflux(i,j,bi,bj) =
                0252      &         - hs(i,j,bi,bj)
                0253      &         - hl(i,j,bi,bj)
                0254      &         + lwflux(i,j,bi,bj)
3752238fd8 Patr*0255 #ifndef SHORTWAVE_HEATING
0320e25227 Mart*0256      &         + swflux(i,j,bi,bj)
3752238fd8 Patr*0257 #endif
14a20b48a3 Jean*0258 C             fresh-water flux from Precipitation and Evaporation.
0320e25227 Mart*0259           sflux(i,j,bi,bj) = evap(i,j,bi,bj) - precip(i,j,bi,bj)
                0260          ENDDO
                0261         ENDDO
3752238fd8 Patr*0262 #endif /* ALLOW_ATM_TEMP */
423768d890 Jean*0263 
                0264 C     Apply runoff, masks and exchanges
0320e25227 Mart*0265         ks = 1
                0266         IF ( usingPCoords ) ks = Nr
                0267         DO j = 1,sNy
                0268          DO i = 1,sNx
bdec91d862 Patr*0269 #ifdef ALLOW_RUNOFF
0320e25227 Mart*0270           sflux(i,j,bi,bj) = sflux(i,j,bi,bj) - runoff(i,j,bi,bj)
bdec91d862 Patr*0271 #endif
0320e25227 Mart*0272           hflux(i,j,bi,bj) = hflux(i,j,bi,bj)*maskC(i,j,ks,bi,bj)
                0273           sflux(i,j,bi,bj) = sflux(i,j,bi,bj)*maskC(i,j,ks,bi,bj)
                0274          ENDDO
14a20b48a3 Jean*0275         ENDDO
0320e25227 Mart*0276 
                0277        ENDDO
14a20b48a3 Jean*0278       ENDDO
bdec91d862 Patr*0279 
14a20b48a3 Jean*0280 C     Update the tile edges: needed for some EXF fields involved in horizontal
                0281 C     averaging, e.g., wind-stress; fields used by main model or other pkgs
                0282 C     are exchanged in EXF_MAPFIELDS.
                0283 c     _EXCH_XY_RL(hflux,   myThid)
                0284 c     _EXCH_XY_RL(sflux,   myThid)
9c3e24f78c Jean*0285       IF ( stressIsOnCgrid ) THEN
0320e25227 Mart*0286        CALL EXCH_UV_XY_RL( ustress, vstress, .TRUE., myThid )
9c3e24f78c Jean*0287       ELSE
0320e25227 Mart*0288        CALL EXCH_UV_AGRID_3D_RL(ustress, vstress, .TRUE., 1, myThid)
9c3e24f78c Jean*0289       ENDIF
bdec91d862 Patr*0290 #ifdef SHORTWAVE_HEATING
14a20b48a3 Jean*0291 c     _EXCH_XY_RL(swflux, myThid)
bdec91d862 Patr*0292 #endif
                0293 #ifdef ATMOSPHERIC_LOADING
14a20b48a3 Jean*0294 c     _EXCH_XY_RL(apressure, myThid)
bdec91d862 Patr*0295 #endif
24da7525ba Jean*0296 #ifdef EXF_SEAICE_FRACTION
14a20b48a3 Jean*0297 c     _EXCH_XY_RL(areamask, myThid)
8f277f2728 Gael*0298 #endif
bdec91d862 Patr*0299 
14a20b48a3 Jean*0300 C     Get values of the surface flux anomalies.
423768d890 Jean*0301       CALL EXF_GETSURFACEFLUXES( myTime, myIter, myThid )
7f861c1808 Patr*0302 
14a20b48a3 Jean*0303       IF ( useExfCheckRange .AND.
9f46642c85 Jean*0304      &     ( myIter.EQ.nIter0 .OR. exf_debugLev.GE.debLevC ) ) THEN
0320e25227 Mart*0305        CALL EXF_CHECK_RANGE( myTime, myIter, myThid )
14a20b48a3 Jean*0306       ENDIF
d7ee8fe52e Patr*0307 
423768d890 Jean*0308 #ifdef ALLOW_AUTODIFF
14a20b48a3 Jean*0309 # ifdef ALLOW_AUTODIFF_MONITOR
0320e25227 Mart*0310       CALL EXF_ADJOINT_SNAPSHOTS( 1, myTime, myIter, myThid )
98238efc8b Patr*0311 # endif
423768d890 Jean*0312 #endif /* ALLOW_AUTODIFF */
98238efc8b Patr*0313 
da754645e6 Jean*0314 #ifdef ALLOW_ATM_TEMP
                0315 # ifdef SHORTWAVE_HEATING
14a20b48a3 Jean*0316 C     Treatment of qnet
da754645e6 Jean*0317 C     The location of the summation of Qnet in exf_mapfields is unfortunate:
14a20b48a3 Jean*0318 C     For backward compatibility issues we want it to happen after
                0319 C     applying control variables, but before exf_diagnostics_fill.
                0320 C     Therefore, we DO it exactly here:
                0321       DO bj = myByLo(myThid), myByHi(myThid)
                0322        DO bi = myBxLo(myThid), myBxHi(myThid)
0320e25227 Mart*0323         DO j = 1-OLy,sNy+OLy
                0324          DO i = 1-OLx,sNx+OLx
92db8eb413 Patr*0325           hflux(i,j,bi,bj) = hflux(i,j,bi,bj) + swflux(i,j,bi,bj)
14a20b48a3 Jean*0326          ENDDO
                0327         ENDDO
                0328        ENDDO
                0329       ENDDO
da754645e6 Jean*0330 # endif /* SHORTWAVE_HEATING */
                0331 #endif /* ALLOW_ATM_TEMP */
b7bd8a1e5a Patr*0332 
14a20b48a3 Jean*0333 C     Diagnostics output
423768d890 Jean*0334       CALL EXF_DIAGNOSTICS_FILL( myTime, myIter, myThid )
b1e3781773 Patr*0335 
14a20b48a3 Jean*0336 C     Monitor output
423768d890 Jean*0337       CALL EXF_MONITOR( myTime, myIter, myThid )
d9263fe447 Jean*0338 
14a20b48a3 Jean*0339 C     Map the forcing fields onto the corresponding model fields.
423768d890 Jean*0340       CALL EXF_MAPFIELDS( myTime, myIter, myThid )
7f861c1808 Patr*0341 
423768d890 Jean*0342 #ifdef ALLOW_AUTODIFF
14a20b48a3 Jean*0343 # ifdef ALLOW_AUTODIFF_MONITOR
                0344       IF ( .NOT. useSEAICE )
7e9378bd0f Patr*0345      &     CALL EXF_ADJOINT_SNAPSHOTS( 3, myTime, myIter, myThid )
                0346 # endif
423768d890 Jean*0347 #endif /* ALLOW_AUTODIFF */
7e9378bd0f Patr*0348 
e1fb02e8f0 Jean*0349       RETURN
                0350       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!