Back to home page

MITgcm

 
 

    


File indexing completed on 2022-11-23 06:09:41 UTC

view on githubraw file Latest commit 20dee616 on 2022-11-22 15:45:38 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 
113c43f200 Jean*0006 CBOP 0
                0007 C     !ROUTINE: EXF_MAPFIELDS
7f861c1808 Patr*0008 
113c43f200 Jean*0009 C     !INTERFACE:
                0010       SUBROUTINE EXF_MAPFIELDS( myTime, myIter, myThid )
7f861c1808 Patr*0011 
113c43f200 Jean*0012 C     !DESCRIPTION:
                0013 C     ==================================================================
                0014 C     SUBROUTINE EXF_MAPFIELDS
                0015 C     ==================================================================
                0016 C
                0017 C     o Map external forcing fields (ustress, vstress, hflux, sflux,
                0018 C       swflux, apressure, climsss, climsst, etc.) onto ocean model
                0019 C       arrays (fu, fv, Qnet, EmPmR, Qsw, pLoad, SSS, SST, etc.).
                0020 C       This routine is included to separate the ocean state estimation
                0021 C       tool as much as possible from the ocean model.  Unit and sign
                0022 C       conventions can be customized using variables exf_outscal_*,
                0023 C       which are set in exf_readparms.F.  See the header files
                0024 C       EXF_FIELDS.h and FFIELDS.h for definitions of the various input
                0025 C       and output fields and for default unit and sign convetions.
                0026 C
                0027 C     started: Christian Eckert eckert@mit.edu  09-Aug-1999
                0028 C
                0029 C     changed: Christian Eckert eckert@mit.edu  11-Jan-2000
                0030 C              - Restructured the code in order to create a package
                0031 C                for the MITgcmUV.
                0032 C
                0033 C              Christian Eckert eckert@mit.edu  12-Feb-2000
                0034 C              - Changed Routine names (package prefix: exf_)
                0035 C
                0036 C              Patrick Heimbach, heimbach@mit.edu  06-May-2000
                0037 C              - added and changed CPP flag structure for
                0038 C                ALLOW_BULKFORMULAE, ALLOW_ATM_TEMP
                0039 C
                0040 C              Patrick Heimbach, heimbach@mit.edu  23-May-2000
                0041 C              - sign change of ustress/vstress incorporated into
                0042 C                scaling factors exf_outscal_ust, exf_outscal_vst
                0043 C
                0044 C     mods for pkg/seaice: menemenlis@jpl.nasa.gov 20-Dec-2002
                0045 C
                0046 C     ==================================================================
                0047 C     SUBROUTINE EXF_MAPFIELDS
                0048 C     ==================================================================
7f861c1808 Patr*0049 
113c43f200 Jean*0050 C     !USES:
                0051       IMPLICIT NONE
7f861c1808 Patr*0052 
113c43f200 Jean*0053 C     == global variables ==
7f861c1808 Patr*0054 #include "EEPARAMS.h"
                0055 #include "SIZE.h"
ec80b41947 Patr*0056 #include "PARAMS.h"
7f861c1808 Patr*0057 #include "FFIELDS.h"
aca832c4c0 Patr*0058 #include "GRID.h"
60bf7ca903 Jean*0059 #include "DYNVARS.h"
aca832c4c0 Patr*0060 
082e18c36c Jean*0061 #include "EXF_PARAM.h"
                0062 #include "EXF_CONSTANTS.h"
                0063 #include "EXF_FIELDS.h"
a8fd6c497b Patr*0064 #ifdef ALLOW_AUTODIFF_TAMC
                0065 # include "tamc.h"
                0066 #endif
7f861c1808 Patr*0067 
113c43f200 Jean*0068 C     !INPUT/OUTPUT PARAMETERS:
                0069 C     myTime  :: Current time in simulation
                0070 C     myIter  :: Current iteration number
                0071 C     myThid  :: my Thread Id number
                0072       _RL     myTime
                0073       INTEGER myIter
                0074       INTEGER myThid
7f861c1808 Patr*0075 
113c43f200 Jean*0076 C     !LOCAL VARIABLES:
e1fb02e8f0 Jean*0077       INTEGER bi,bj
                0078       INTEGER i,j,ks
9c3e24f78c Jean*0079       INTEGER imin, imax
                0080       INTEGER jmin, jmax
                0081       PARAMETER ( imin = 1-OLx , imax = sNx+OLx )
                0082       PARAMETER ( jmin = 1-OLy , jmax = sNy+OLy )
7c50f07931 Mart*0083 #ifdef ALLOW_AUTODIFF_TAMC
                0084       INTEGER ikey
                0085 #endif
113c43f200 Jean*0086 CEOP
7f861c1808 Patr*0087 
e1fb02e8f0 Jean*0088 C--   set surface level index:
                0089       ks = 1
0320e25227 Mart*0090       IF ( usingPCoords ) ks = Nr
e1fb02e8f0 Jean*0091 
9c3e24f78c Jean*0092       DO bj = myByLo(myThid),myByHi(myThid)
113c43f200 Jean*0093        DO bi = myBxLo(myThid),myBxHi(myThid)
a8fd6c497b Patr*0094 
                0095 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0096           ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
a8fd6c497b Patr*0097 #endif /* ALLOW_AUTODIFF_TAMC */
                0098 
113c43f200 Jean*0099 C     Heat flux.
                0100           DO j = jmin,jmax
                0101             DO i = imin,imax
                0102              Qnet(i,j,bi,bj) = exf_outscal_hflux*hflux(i,j,bi,bj)
                0103             ENDDO
                0104           ENDDO
                0105           IF ( hfluxfile .EQ. ' ' ) THEN
                0106            DO j = jmin,jmax
                0107             DO i = imin,imax
                0108              Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) -
ec80b41947 Patr*0109      &            exf_outscal_hflux * ( hflux_exfremo_intercept +
113c43f200 Jean*0110      &            hflux_exfremo_slope*(myTime-startTime) )
                0111             ENDDO
                0112            ENDDO
                0113           ENDIF
a8fd6c497b Patr*0114 
a66aad0124 Gael*0115 C     Freshwater flux.
113c43f200 Jean*0116           DO j = jmin,jmax
                0117             DO i = imin,imax
6206cdb986 Jean*0118              EmPmR(i,j,bi,bj)= exf_outscal_sflux*sflux(i,j,bi,bj)
                0119      &                                          *rhoConstFresh
113c43f200 Jean*0120             ENDDO
                0121           ENDDO
                0122           IF ( sfluxfile .EQ. ' ' ) THEN
                0123            DO j = jmin,jmax
                0124             DO i = imin,imax
                0125              EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj) - rhoConstFresh*
ec80b41947 Patr*0126      &            exf_outscal_sflux * ( sflux_exfremo_intercept +
113c43f200 Jean*0127      &            sflux_exfremo_slope*(myTime-startTime) )
                0128             ENDDO
                0129            ENDDO
                0130           ENDIF
7f861c1808 Patr*0131 
60bf7ca903 Jean*0132 #ifdef ALLOW_ATM_TEMP
                0133 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0134           IF ( temp_EvPrRn .NE. UNSET_RL ) THEN
                0135 C--   Account for energy content of Precip + RunOff & Evap. Assumes:
                0136 C     1) Rain has same temp as Air
                0137 C     2) Snow has no heat capacity (consistent with seaice & thsice pkgs)
                0138 C     3) No distinction between sea-water Cp and fresh-water Cp
e603dbf008 Dimi*0139 C     4) By default, RunOff comes at the temp of surface water (with same Cp);
                0140 C        ifdef ALLOW_RUNOFTEMP, RunOff temp can be specified in runoftempfile.
60bf7ca903 Jean*0141 C     5) Evap is released to the Atmos @ surf-temp (=SST); should be using
                0142 C        the water-vapor heat capacity here and consistently in Bulk-Formulae;
                0143 C        Could also be put directly into Latent Heat flux.
                0144            IF ( snowPrecipFile .NE. ' ' ) THEN
                0145 C--   Melt snow (if provided) into the ocean and account for rain-temp
                0146             DO j = 1, sNy
                0147              DO i = 1, sNx
                0148               Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
                0149      &              + flami*snowPrecip(i,j,bi,bj)*rhoConstFresh
                0150      &              - HeatCapacity_Cp
                0151      &               *( atemp(i,j,bi,bj) - cen2kel - temp_EvPrRn )
                0152      &               *( precip(i,j,bi,bj)- snowPrecip(i,j,bi,bj) )
                0153      &               *rhoConstFresh
                0154              ENDDO
                0155             ENDDO
                0156            ELSE
                0157 C--   Make snow (according to Air Temp) and melt it in the ocean
                0158 C     note: here we just use the same criteria as over seaice but would be
                0159 C           better to consider a higher altitude air temp, e.g., 850.mb
                0160             DO j = 1, sNy
                0161              DO i = 1, sNx
                0162               IF ( atemp(i,j,bi,bj).LT.cen2kel ) THEN
                0163                Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
                0164      &              + flami*precip(i,j,bi,bj)*rhoConstFresh
                0165               ELSE
                0166 C--   Account for rain-temp
                0167                Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
                0168      &              - HeatCapacity_Cp
                0169      &               *( atemp(i,j,bi,bj) - cen2kel - temp_EvPrRn )
                0170      &               *precip(i,j,bi,bj)*rhoConstFresh
                0171               ENDIF
                0172              ENDDO
                0173             ENDDO
                0174            ENDIF
27ebbc08f6 Dimi*0175 #ifdef ALLOW_RUNOFF
                0176 C--   Account for energy content of RunOff:
                0177            DO j = 1, sNy
                0178             DO i = 1, sNx
                0179               Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
                0180      &              - HeatCapacity_Cp
                0181      &               *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
                0182      &               *runoff(i,j,bi,bj)*rhoConstFresh
                0183             ENDDO
                0184            ENDDO
                0185 #endif
60bf7ca903 Jean*0186 C--   Account for energy content of Evap:
                0187            DO j = 1, sNy
                0188             DO i = 1, sNx
                0189               Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
                0190      &              + HeatCapacity_Cp
                0191      &               *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
                0192      &               *evap(i,j,bi,bj)*rhoConstFresh
                0193               Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)*maskC(i,j,ks,bi,bj)
                0194             ENDDO
                0195            ENDDO
                0196           ENDIF
27ebbc08f6 Dimi*0197 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0198 #endif /* ALLOW_ATM_TEMP */
                0199 #if defined(ALLOW_RUNOFF) && defined(ALLOW_RUNOFTEMP)
                0200           IF ( runoftempfile .NE. ' ' ) THEN
                0201 C--   Add energy content of RunOff
e603dbf008 Dimi*0202            DO j = 1, sNy
                0203             DO i = 1, sNx
                0204                Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj)
27ebbc08f6 Dimi*0205      &              + HeatCapacity_Cp
                0206      &              *( theta(i,j,ks,bi,bj) - runoftemp(i,j,bi,bj) )
e603dbf008 Dimi*0207      &              *runoff(i,j,bi,bj)*rhoConstFresh
                0208             ENDDO
                0209            ENDDO
                0210           ENDIF
27ebbc08f6 Dimi*0211 #endif
60bf7ca903 Jean*0212 
a8fd6c497b Patr*0213 #ifdef ALLOW_AUTODIFF_TAMC
                0214 CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
                0215 #endif
113c43f200 Jean*0216           DO j = jmin,jmax
                0217             DO i = imin,imax
                0218 C             Zonal wind stress.
                0219               IF (ustress(i,j,bi,bj).GT.windstressmax) THEN
df66ce88b0 Mart*0220                 ustress(i,j,bi,bj)=windstressmax
113c43f200 Jean*0221               ENDIF
                0222             ENDDO
                0223           ENDDO
a8fd6c497b Patr*0224 #ifdef ALLOW_AUTODIFF_TAMC
                0225 CADJ STORE ustress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
                0226 #endif
113c43f200 Jean*0227           DO j = jmin,jmax
                0228             DO i = imin,imax
                0229               IF (ustress(i,j,bi,bj).LT.-windstressmax) THEN
df66ce88b0 Mart*0230                 ustress(i,j,bi,bj)=-windstressmax
113c43f200 Jean*0231               ENDIF
                0232             ENDDO
                0233           ENDDO
9c3e24f78c Jean*0234           IF ( stressIsOnCgrid ) THEN
113c43f200 Jean*0235            DO j = jmin,jmax
                0236             DO i = imin+1,imax
aca832c4c0 Patr*0237               fu(i,j,bi,bj) = exf_outscal_ustress*ustress(i,j,bi,bj)
113c43f200 Jean*0238             ENDDO
                0239            ENDDO
9c3e24f78c Jean*0240           ELSE
113c43f200 Jean*0241            DO j = jmin,jmax
                0242             DO i = imin+1,imax
                0243 C     Shift wind stresses calculated at Grid-center to W/S points
9c3e24f78c Jean*0244               fu(i,j,bi,bj) = exf_outscal_ustress*
                0245      &              (ustress(i,j,bi,bj)+ustress(i-1,j,bi,bj))
e1fb02e8f0 Jean*0246      &              *exf_half*maskW(i,j,ks,bi,bj)
113c43f200 Jean*0247             ENDDO
                0248            ENDDO
9c3e24f78c Jean*0249           ENDIF
7f861c1808 Patr*0250 
a8fd6c497b Patr*0251 #ifdef ALLOW_AUTODIFF_TAMC
                0252 CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
                0253 #endif
113c43f200 Jean*0254           DO j = jmin,jmax
                0255             DO i = imin,imax
                0256 C             Meridional wind stress.
                0257               IF (vstress(i,j,bi,bj).GT.windstressmax) THEN
df66ce88b0 Mart*0258                 vstress(i,j,bi,bj)=windstressmax
113c43f200 Jean*0259               ENDIF
                0260             ENDDO
                0261           ENDDO
a8fd6c497b Patr*0262 #ifdef ALLOW_AUTODIFF_TAMC
                0263 CADJ STORE vstress(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
                0264 #endif
113c43f200 Jean*0265           DO j = jmin,jmax
                0266             DO i = imin,imax
                0267               IF (vstress(i,j,bi,bj).LT.-windstressmax) THEN
df66ce88b0 Mart*0268                 vstress(i,j,bi,bj)=-windstressmax
113c43f200 Jean*0269               ENDIF
                0270             ENDDO
                0271           ENDDO
9c3e24f78c Jean*0272           IF ( stressIsOnCgrid ) THEN
113c43f200 Jean*0273            DO j = jmin+1,jmax
                0274             DO i = imin,imax
9c3e24f78c Jean*0275               fv(i,j,bi,bj) = exf_outscal_vstress*vstress(i,j,bi,bj)
113c43f200 Jean*0276             ENDDO
                0277            ENDDO
9c3e24f78c Jean*0278           ELSE
113c43f200 Jean*0279            DO j = jmin+1,jmax
                0280             DO i = imin,imax
                0281 C     Shift wind stresses calculated at C-points to W/S points
aca832c4c0 Patr*0282               fv(i,j,bi,bj) = exf_outscal_vstress*
9c3e24f78c Jean*0283      &              (vstress(i,j,bi,bj)+vstress(i,j-1,bi,bj))
e1fb02e8f0 Jean*0284      &              *exf_half*maskS(i,j,ks,bi,bj)
113c43f200 Jean*0285             ENDDO
                0286            ENDDO
9c3e24f78c Jean*0287           ENDIF
7f861c1808 Patr*0288 
113c43f200 Jean*0289 #if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
                0290 C             Short wave radiative flux.
                0291           DO j = jmin,jmax
                0292             DO i = imin,imax
                0293              Qsw(i,j,bi,bj)  = exf_outscal_swflux*swflux(i,j,bi,bj)
                0294             ENDDO
                0295           ENDDO
7f861c1808 Patr*0296 #endif
                0297 
                0298 #ifdef ALLOW_CLIMSST_RELAXATION
113c43f200 Jean*0299           DO j = jmin,jmax
                0300             DO i = imin,imax
                0301              SST(i,j,bi,bj)  = exf_outscal_sst*climsst(i,j,bi,bj)
                0302             ENDDO
                0303           ENDDO
7f861c1808 Patr*0304 #endif
                0305 
                0306 #ifdef ALLOW_CLIMSSS_RELAXATION
113c43f200 Jean*0307           DO j = jmin,jmax
                0308             DO i = imin,imax
                0309              SSS(i,j,bi,bj)  = exf_outscal_sss*climsss(i,j,bi,bj)
                0310             ENDDO
                0311           ENDDO
7f861c1808 Patr*0312 #endif
                0313 
a8fd6c497b Patr*0314 #ifdef ATMOSPHERIC_LOADING
c49e52d441 Jean*0315 C-    subtract "surf_pRef" to fill in atmos. pressure anomaly "pLoad"
113c43f200 Jean*0316           DO j = jmin,jmax
                0317             DO i = imin,imax
c49e52d441 Jean*0318              pLoad(i,j,bi,bj) =
                0319      &          exf_outscal_apressure*apressure(i,j,bi,bj) - surf_pRef
113c43f200 Jean*0320             ENDDO
                0321           ENDDO
a8fd6c497b Patr*0322 #endif
                0323 
497d85062c Jean*0324 #ifdef EXF_ALLOW_TIDES
                0325           DO j = jmin,jmax
                0326             DO i = imin,imax
                0327              phiTide2d(i,j,bi,bj)=exf_outscal_tidePot*tidePot(i,j,bi,bj)
                0328             ENDDO
                0329           ENDDO
                0330 #endif /* EXF_ALLOW_TIDES */
                0331 
a66aad0124 Gael*0332 #ifdef ALLOW_SALTFLX
                0333           DO j = jmin,jmax
                0334             DO i = imin,imax
0320e25227 Mart*0335               saltFlux(i,j,bi,bj) = saltflx(i,j,bi,bj)
a66aad0124 Gael*0336             ENDDO
                0337           ENDDO
                0338 #endif
                0339 
24da7525ba Jean*0340 #ifdef EXF_SEAICE_FRACTION
113c43f200 Jean*0341           DO j = jmin,jmax
                0342             DO i = imin,imax
24da7525ba Jean*0343               exf_iceFraction(i,j,bi,bj) =
d877a5eaeb Patr*0344      &           exf_outscal_areamask*areamask(i,j,bi,bj)
0320e25227 Mart*0345               exf_iceFraction(i,j,bi,bj) =
                0346      &           MIN( MAX(exf_iceFraction(i,j,bi,bj),zeroRS), oneRS )
113c43f200 Jean*0347             ENDDO
                0348           ENDDO
d877a5eaeb Patr*0349 #endif
                0350 
113c43f200 Jean*0351        ENDDO
9c3e24f78c Jean*0352       ENDDO
7f861c1808 Patr*0353 
113c43f200 Jean*0354 C--   Update the tile edges.
                0355       _EXCH_XY_RS(  Qnet, myThid )
                0356       _EXCH_XY_RS( EmPmR, myThid )
8193892475 Curt*0357        CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
113c43f200 Jean*0358 c#if defined(ALLOW_ATM_TEMP) || defined(SHORTWAVE_HEATING)
6060ec2938 Dimi*0359 #ifdef SHORTWAVE_HEATING
113c43f200 Jean*0360 C     Qsw used in SHORTWAVE_HEATING code & for diagnostics (<- EXCH not needed)
                0361       _EXCH_XY_RS(   Qsw, myThid )
7f861c1808 Patr*0362 #endif
                0363 #ifdef ALLOW_CLIMSST_RELAXATION
113c43f200 Jean*0364       _EXCH_XY_RS(   SST, myThid )
7f861c1808 Patr*0365 #endif
                0366 #ifdef ALLOW_CLIMSSS_RELAXATION
113c43f200 Jean*0367       _EXCH_XY_RS(   SSS, myThid )
7f861c1808 Patr*0368 #endif
a8fd6c497b Patr*0369 #ifdef ATMOSPHERIC_LOADING
113c43f200 Jean*0370       _EXCH_XY_RS( pLoad, myThid )
a8fd6c497b Patr*0371 #endif
497d85062c Jean*0372 #ifdef EXF_ALLOW_TIDES
                0373       _EXCH_XY_RS( phiTide2d, myThid )
                0374 #endif
24da7525ba Jean*0375 #ifdef EXF_SEAICE_FRACTION
                0376       _EXCH_XY_RS( exf_iceFraction, myThid )
d877a5eaeb Patr*0377 #endif
7f861c1808 Patr*0378 
9c3e24f78c Jean*0379       RETURN
                0380       END