Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:45:40 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
aea29c8517 Alis*0001 #include "CPP_OPTIONS.h"
                0002 
1feb54070b Jean*0003 C--  File apply_forcing.F:
                0004 C--   Contents
                0005 C--   o APPLY_FORCING_U
                0006 C--   o APPLY_FORCING_V
                0007 C--   o APPLY_FORCING_T
                0008 C--   o APPLY_FORCING_S
                0009 
                0010 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
b4656da4c6 Jean*0011 CBOP
1feb54070b Jean*0012 C     !ROUTINE: APPLY_FORCING_U
b4656da4c6 Jean*0013 C     !INTERFACE:
1feb54070b Jean*0014       SUBROUTINE APPLY_FORCING_U(
                0015      U                     gU_arr,
                0016      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0017      I                     myTime, myIter, myThid )
b4656da4c6 Jean*0018 C     !DESCRIPTION: \bv
                0019 C     *==========================================================*
1feb54070b Jean*0020 C     | S/R APPLY_FORCING_U
b4656da4c6 Jean*0021 C     | o Contains problem specific forcing for zonal velocity.
                0022 C     *==========================================================*
                0023 C     | Adds terms to gU for forcing by external sources
                0024 C     | e.g. wind stress, bottom friction etc ...
                0025 C     *==========================================================*
                0026 C     \ev
aea29c8517 Alis*0027 
b4656da4c6 Jean*0028 C     !USES:
                0029       IMPLICIT NONE
aea29c8517 Alis*0030 C     == Global data ==
                0031 #include "SIZE.h"
                0032 #include "EEPARAMS.h"
                0033 #include "PARAMS.h"
                0034 #include "GRID.h"
fe959cba9b Jean*0035 #include "SURFACE.h"
aea29c8517 Alis*0036 #include "DYNVARS.h"
                0037 #include "FFIELDS.h"
                0038 
b4656da4c6 Jean*0039 C     !INPUT/OUTPUT PARAMETERS:
1feb54070b Jean*0040 C     gU_arr    :: the tendency array
b4656da4c6 Jean*0041 C     iMin,iMax :: Working range of x-index for applying forcing.
                0042 C     jMin,jMax :: Working range of y-index for applying forcing.
1feb54070b Jean*0043 C     k         :: Current vertical level index
b4656da4c6 Jean*0044 C     bi,bj     :: Current tile indices
                0045 C     myTime    :: Current time in simulation
1feb54070b Jean*0046 C     myIter    :: Current iteration number
                0047 C     myThid    :: my Thread Id number
                0048       _RL     gU_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0049       INTEGER iMin, iMax, jMin, jMax
                0050       INTEGER k, bi, bj
                0051       _RL     myTime
                0052       INTEGER myIter
aea29c8517 Alis*0053       INTEGER myThid
                0054 
b4656da4c6 Jean*0055 C     !LOCAL VARIABLES:
                0056 C     i,j       :: Loop counters
                0057       INTEGER i, j
                0058 CEOP
fe959cba9b Jean*0059       _RL recip_P0g, termP, rFullDepth
                0060       _RL kV, kF, sigma_b
aea29c8517 Alis*0061 
1feb54070b Jean*0062 C--   Forcing term
                0063       kF = 1. _d 0/86400. _d 0
5375376566 Jean*0064       sigma_b = 0.7 _d 0
fe959cba9b Jean*0065       rFullDepth = rF(1)-rF(Nr+1)
b4656da4c6 Jean*0066 c     DO j=1,sNy
                0067 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
                0068       DO j=0,sNy+1
                0069        DO i=1,sNx+1
1feb54070b Jean*0070         IF ( maskW(i,j,k,bi,bj).EQ.oneRS ) THEN
fe959cba9b Jean*0071          IF ( selectSigmaCoord.EQ.0 ) THEN
1feb54070b Jean*0072           recip_P0g = MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i-1,j,bi,bj))
                0073           termP = 0.5 _d 0*( MIN( rF(k)*recip_P0g, oneRL )
                0074      &                      +rF(k+1)*recip_P0g )
                0075 c         termP = 0.5 _d 0*( rF(k) + rF(k+1) )*recip_P0g
fe959cba9b Jean*0076          ELSE
                0077 C-- Pressure at U.point :
                0078 c         midP = rLowW(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
                0079 c    &         + bHybSigmC(k)
                0080 c    &          *(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
                0081 C-- Sigma at U.point :
                0082 c         termP = ( midP - rLowW(i,j,bi,bj))
                0083 c    &          /(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
                0084 C-  which simplifies to:
1feb54070b Jean*0085           termP = aHybSigmC(k)*rFullDepth
                0086 #ifdef NONLIN_FRSURF
fe959cba9b Jean*0087      &          /(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
1feb54070b Jean*0088 #else
                0089      &          /(rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
                0090 #endif
                0091      &          + bHybSigmC(k)
fe959cba9b Jean*0092          ENDIF
1feb54070b Jean*0093          kV = kF*MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )
                0094          gU_arr(i,j) = gU_arr(i,j)
                0095      &               - kV*uVel(i,j,k,bi,bj)
aea29c8517 Alis*0096         ENDIF
                0097        ENDDO
                0098       ENDDO
                0099 
                0100       RETURN
                0101       END
b4656da4c6 Jean*0102 
                0103 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0104 CBOP
1feb54070b Jean*0105 C     !ROUTINE: APPLY_FORCING_V
b4656da4c6 Jean*0106 C     !INTERFACE:
1feb54070b Jean*0107       SUBROUTINE APPLY_FORCING_V(
                0108      U                     gV_arr,
                0109      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0110      I                     myTime, myIter, myThid )
b4656da4c6 Jean*0111 C     !DESCRIPTION: \bv
                0112 C     *==========================================================*
1feb54070b Jean*0113 C     | S/R APPLY_FORCING_V
b4656da4c6 Jean*0114 C     | o Contains problem specific forcing for merid velocity.
                0115 C     *==========================================================*
                0116 C     | Adds terms to gV for forcing by external sources
                0117 C     | e.g. wind stress, bottom friction etc ...
                0118 C     *==========================================================*
                0119 C     \ev
aea29c8517 Alis*0120 
b4656da4c6 Jean*0121 C     !USES:
                0122       IMPLICIT NONE
aea29c8517 Alis*0123 C     == Global data ==
                0124 #include "SIZE.h"
                0125 #include "EEPARAMS.h"
                0126 #include "PARAMS.h"
                0127 #include "GRID.h"
fe959cba9b Jean*0128 #include "SURFACE.h"
aea29c8517 Alis*0129 #include "DYNVARS.h"
                0130 #include "FFIELDS.h"
                0131 
b4656da4c6 Jean*0132 C     !INPUT/OUTPUT PARAMETERS:
1feb54070b Jean*0133 C     gV_arr    :: the tendency array
b4656da4c6 Jean*0134 C     iMin,iMax :: Working range of x-index for applying forcing.
                0135 C     jMin,jMax :: Working range of y-index for applying forcing.
1feb54070b Jean*0136 C     k         :: Current vertical level index
b4656da4c6 Jean*0137 C     bi,bj     :: Current tile indices
                0138 C     myTime    :: Current time in simulation
1feb54070b Jean*0139 C     myIter    :: Current iteration number
                0140 C     myThid    :: my Thread Id number
                0141       _RL     gV_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0142       INTEGER iMin, iMax, jMin, jMax
                0143       INTEGER k, bi, bj
                0144       _RL     myTime
                0145       INTEGER myIter
aea29c8517 Alis*0146       INTEGER myThid
b4656da4c6 Jean*0147 
                0148 C     !LOCAL VARIABLES:
                0149 C     i,j       :: Loop counters
                0150       INTEGER i, j
                0151 CEOP
fe959cba9b Jean*0152       _RL recip_P0g, termP, rFullDepth
                0153       _RL kV, kF, sigma_b
aea29c8517 Alis*0154 
1feb54070b Jean*0155 C--   Forcing term
                0156       kF = 1. _d 0/86400. _d 0
5375376566 Jean*0157       sigma_b = 0.7 _d 0
fe959cba9b Jean*0158       rFullDepth = rF(1)-rF(Nr+1)
b4656da4c6 Jean*0159       DO j=1,sNy+1
                0160 c      DO i=1,sNx
                0161 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
                0162        DO i=0,sNx+1
1feb54070b Jean*0163         IF ( maskS(i,j,k,bi,bj).EQ.oneRS ) THEN
fe959cba9b Jean*0164          IF ( selectSigmaCoord.EQ.0 ) THEN
1feb54070b Jean*0165           recip_P0g = MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i,j-1,bi,bj))
                0166           termP = 0.5 _d 0*( MIN( rF(k)*recip_P0g, oneRL )
                0167      &                      +rF(k+1)*recip_P0g )
                0168 c         termP = 0.5 _d 0*( rF(k) + rF(k+1) )*recip_P0g
fe959cba9b Jean*0169          ELSE
                0170 C-- Pressure at V.point :
                0171 c         midP = rLowS(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
                0172 c    &         + bHybSigmC(k)
                0173 c    &          *(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
                0174 C-- Sigma at V.point :
                0175 c         termP = ( midP - rLowS(i,j,bi,bj))
                0176 c    &          /(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
                0177 C-  which simplifies to:
1feb54070b Jean*0178           termP = aHybSigmC(k)*rFullDepth
                0179 #ifdef NONLIN_FRSURF
fe959cba9b Jean*0180      &          /(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
1feb54070b Jean*0181 #else
                0182      &          /(rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
                0183 #endif
                0184      &          + bHybSigmC(k)
fe959cba9b Jean*0185          ENDIF
1feb54070b Jean*0186          kV = kF*MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )
                0187          gV_arr(i,j) = gV_arr(i,j)
                0188      &               - kV*vVel(i,j,k,bi,bj)
aea29c8517 Alis*0189         ENDIF
                0190        ENDDO
                0191       ENDDO
                0192 
                0193       RETURN
                0194       END
b4656da4c6 Jean*0195 
                0196 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0197 CBOP
1feb54070b Jean*0198 C     !ROUTINE: APPLY_FORCING_T
b4656da4c6 Jean*0199 C     !INTERFACE:
1feb54070b Jean*0200       SUBROUTINE APPLY_FORCING_T(
                0201      U                     gT_arr,
                0202      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0203      I                     myTime, myIter, myThid )
b4656da4c6 Jean*0204 C     !DESCRIPTION: \bv
                0205 C     *==========================================================*
1feb54070b Jean*0206 C     | S/R APPLY_FORCING_T
b4656da4c6 Jean*0207 C     | o Contains problem specific forcing for temperature.
                0208 C     *==========================================================*
                0209 C     | Adds terms to gT for forcing by external sources
                0210 C     | e.g. heat flux, climatalogical relaxation, etc ...
                0211 C     *==========================================================*
                0212 C     \ev
aea29c8517 Alis*0213 
b4656da4c6 Jean*0214 C     !USES:
                0215       IMPLICIT NONE
aea29c8517 Alis*0216 C     == Global data ==
                0217 #include "SIZE.h"
                0218 #include "EEPARAMS.h"
                0219 #include "PARAMS.h"
                0220 #include "GRID.h"
                0221 #include "DYNVARS.h"
                0222 #include "FFIELDS.h"
                0223 
b4656da4c6 Jean*0224 C     !INPUT/OUTPUT PARAMETERS:
1feb54070b Jean*0225 C     gT_arr    :: the tendency array
b4656da4c6 Jean*0226 C     iMin,iMax :: Working range of x-index for applying forcing.
                0227 C     jMin,jMax :: Working range of y-index for applying forcing.
1feb54070b Jean*0228 C     k         :: Current vertical level index
b4656da4c6 Jean*0229 C     bi,bj     :: Current tile indices
                0230 C     myTime    :: Current time in simulation
1feb54070b Jean*0231 C     myIter    :: Current iteration number
                0232 C     myThid    :: my Thread Id number
                0233       _RL     gT_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0234       INTEGER iMin, iMax, jMin, jMax
                0235       INTEGER k, bi, bj
                0236       _RL     myTime
                0237       INTEGER myIter
aea29c8517 Alis*0238       INTEGER myThid
                0239 
b4656da4c6 Jean*0240 C     !LOCAL VARIABLES:
                0241 C     i,j       :: Loop counters
1feb54070b Jean*0242 C     kSurface  :: index of surface level
b4656da4c6 Jean*0243       INTEGER i, j
                0244 CEOP
1feb54070b Jean*0245       _RL thetaLim, kT, ka, ks, sigma_b, term1, term2, thetaEq
fe959cba9b Jean*0246       _RL termP, rFullDepth
aea29c8517 Alis*0247 
1feb54070b Jean*0248 C--   Forcing term
                0249       ka = 1. _d 0/(40. _d 0*86400. _d 0)
                0250       ks = 1. _d 0/(4. _d 0 *86400. _d 0)
5375376566 Jean*0251       sigma_b = 0.7 _d 0
fe959cba9b Jean*0252       rFullDepth = rF(1)-rF(Nr+1)
1feb54070b Jean*0253       DO j=0,sNy+1
                0254        DO i=0,sNx+1
                0255          term1 = 60. _d 0*(SIN(yC(i,j,bi,bj)*deg2rad)**2)
                0256          termP = 0.5 _d 0*( rF(k) + rF(k+1) )
                0257          term2 = 10. _d 0*LOG(termP/atm_po)
b4656da4c6 Jean*0258      &            *(COS(yC(i,j,bi,bj)*deg2rad)**2)
5375376566 Jean*0259          thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)
1feb54070b Jean*0260          thetaEq = 315. _d 0 - term1 - term2
                0261          thetaEq = MAX(thetaLim,thetaEq)
fe959cba9b Jean*0262          IF ( selectSigmaCoord.EQ.0 ) THEN
1feb54070b Jean*0263           termP = 0.5 _d 0*( MIN(rF(k),Ro_surf(i,j,bi,bj))
                0264      &                     + rF(k+1) )
                0265      &                    *recip_Rcol(i,j,bi,bj)
fe959cba9b Jean*0266          ELSE
                0267 C-- Pressure at T.point :
                0268 c         midP = R_low(i,j,bi,bj) + aHybSigmC(k)*rFullDepth
                0269 c    &         + bHybSigmC(k)
                0270 c    &          *(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
                0271 C-- Sigma at T.point :
                0272 c         termP = ( midP - R_low(i,j,bi,bj))
                0273 c    &          /(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
                0274 C-  which simplifies to:
1feb54070b Jean*0275           termP = aHybSigmC(k)*rFullDepth
                0276 #ifdef NONLIN_FRSURF
fe959cba9b Jean*0277      &          /(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
1feb54070b Jean*0278 #else
                0279      &          /(Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
                0280 #endif
                0281      &          + bHybSigmC(k)
fe959cba9b Jean*0282          ENDIF
1feb54070b Jean*0283          kT = ka+(ks-ka)
                0284      &      *MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) )
                0285      &      *COS((yC(i,j,bi,bj)*deg2rad))**4
                0286          gT_arr(i,j) = gT_arr(i,j)
                0287      &               - kT*( theta(i,j,k,bi,bj)-thetaEq )
                0288      &                *maskC(i,j,k,bi,bj)
aea29c8517 Alis*0289        ENDDO
                0290       ENDDO
                0291 
                0292       RETURN
                0293       END
b4656da4c6 Jean*0294 
                0295 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0296 CBOP
1feb54070b Jean*0297 C     !ROUTINE: APPLY_FORCING_S
b4656da4c6 Jean*0298 C     !INTERFACE:
1feb54070b Jean*0299       SUBROUTINE APPLY_FORCING_S(
                0300      U                     gS_arr,
                0301      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0302      I                     myTime, myIter, myThid )
b4656da4c6 Jean*0303 C     !DESCRIPTION: \bv
                0304 C     *==========================================================*
1feb54070b Jean*0305 C     | S/R APPLY_FORCING_S
b4656da4c6 Jean*0306 C     | o Contains problem specific forcing for merid velocity.
                0307 C     *==========================================================*
                0308 C     | Adds terms to gS for forcing by external sources
                0309 C     | e.g. fresh-water flux, climatalogical relaxation, etc ...
                0310 C     *==========================================================*
                0311 C     \ev
aea29c8517 Alis*0312 
b4656da4c6 Jean*0313 C     !USES:
                0314       IMPLICIT NONE
aea29c8517 Alis*0315 C     == Global data ==
                0316 #include "SIZE.h"
                0317 #include "EEPARAMS.h"
                0318 #include "PARAMS.h"
                0319 #include "GRID.h"
                0320 #include "DYNVARS.h"
                0321 #include "FFIELDS.h"
1feb54070b Jean*0322 #include "SURFACE.h"
aea29c8517 Alis*0323 
b4656da4c6 Jean*0324 C     !INPUT/OUTPUT PARAMETERS:
1feb54070b Jean*0325 C     gS_arr    :: the tendency array
b4656da4c6 Jean*0326 C     iMin,iMax :: Working range of x-index for applying forcing.
                0327 C     jMin,jMax :: Working range of y-index for applying forcing.
1feb54070b Jean*0328 C     k         :: Current vertical level index
b4656da4c6 Jean*0329 C     bi,bj     :: Current tile indices
                0330 C     myTime    :: Current time in simulation
1feb54070b Jean*0331 C     myIter    :: Current iteration number
                0332 C     myThid    :: my Thread Id number
                0333       _RL     gS_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0334       INTEGER iMin, iMax, jMin, jMax
                0335       INTEGER k, bi, bj
                0336       _RL     myTime
                0337       INTEGER myIter
aea29c8517 Alis*0338       INTEGER myThid
                0339 
b4656da4c6 Jean*0340 C     !LOCAL VARIABLES:
                0341 C     i,j       :: Loop counters
                0342 c     INTEGER i, j
                0343 CEOP
aea29c8517 Alis*0344 
1feb54070b Jean*0345 C--   Forcing term
aea29c8517 Alis*0346 
                0347       RETURN
                0348       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!