Back to home page

MITgcm
 
 

    

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

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
2e2b97e510 Jean*0001 #include "CPP_OPTIONS.h"
                0002 
                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-|--+----|
                0011 CBOP
                0012 C     !ROUTINE: APPLY_FORCING_U
                0013 C     !INTERFACE:
                0014       SUBROUTINE APPLY_FORCING_U(
                0015      U                     gU_arr,
                0016      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0017      I                     myTime, myIter, myThid )
                0018 C     !DESCRIPTION: \bv
                0019 C     *==========================================================*
                0020 C     | S/R APPLY_FORCING_U
                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
                0027 
                0028 C     !USES:
                0029       IMPLICIT NONE
                0030 C     == Global data ==
                0031 #include "SIZE.h"
                0032 #include "EEPARAMS.h"
                0033 #include "PARAMS.h"
                0034 #include "GRID.h"
                0035 #include "SURFACE.h"
                0036 #include "DYNVARS.h"
                0037 #include "FFIELDS.h"
                0038 
                0039 C     !INPUT/OUTPUT PARAMETERS:
                0040 C     gU_arr    :: the tendency array
                0041 C     iMin,iMax :: Working range of x-index for applying forcing.
                0042 C     jMin,jMax :: Working range of y-index for applying forcing.
                0043 C     k         :: Current vertical level index
                0044 C     bi,bj     :: Current tile indices
                0045 C     myTime    :: Current time in simulation
                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
                0053       INTEGER myThid
                0054 
                0055 C     !LOCAL VARIABLES:
                0056 C     i,j       :: Loop counters
                0057       INTEGER i, j
                0058 CEOP
                0059       _RL recip_P0g, termP, rFullDepth
                0060       _RL kV, kF, sigma_b
                0061 
                0062 C--   Forcing term
                0063       kF = 1. _d 0/86400. _d 0
                0064       sigma_b = 0.7 _d 0
                0065       rFullDepth = rF(1)-rF(Nr+1)
                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
                0070         IF ( maskW(i,j,k,bi,bj).EQ.oneRS ) THEN
                0071          IF ( selectSigmaCoord.EQ.0 ) THEN
                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
                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:
                0085           termP = aHybSigmC(k)*rFullDepth
                0086 #ifdef NONLIN_FRSURF
                0087      &          /(etaHw(i,j,bi,bj)+rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
                0088 #else
                0089      &          /(rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj))
                0090 #endif
                0091      &          + bHybSigmC(k)
                0092          ENDIF
                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)
                0096         ENDIF
                0097        ENDDO
                0098       ENDDO
                0099 
                0100       RETURN
                0101       END
                0102 
                0103 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0104 CBOP
                0105 C     !ROUTINE: APPLY_FORCING_V
                0106 C     !INTERFACE:
                0107       SUBROUTINE APPLY_FORCING_V(
                0108      U                     gV_arr,
                0109      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0110      I                     myTime, myIter, myThid )
                0111 C     !DESCRIPTION: \bv
                0112 C     *==========================================================*
                0113 C     | S/R APPLY_FORCING_V
                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
                0120 
                0121 C     !USES:
                0122       IMPLICIT NONE
                0123 C     == Global data ==
                0124 #include "SIZE.h"
                0125 #include "EEPARAMS.h"
                0126 #include "PARAMS.h"
                0127 #include "GRID.h"
                0128 #include "SURFACE.h"
                0129 #include "DYNVARS.h"
                0130 #include "FFIELDS.h"
                0131 
                0132 C     !INPUT/OUTPUT PARAMETERS:
                0133 C     gV_arr    :: the tendency array
                0134 C     iMin,iMax :: Working range of x-index for applying forcing.
                0135 C     jMin,jMax :: Working range of y-index for applying forcing.
                0136 C     k         :: Current vertical level index
                0137 C     bi,bj     :: Current tile indices
                0138 C     myTime    :: Current time in simulation
                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
                0146       INTEGER myThid
                0147 
                0148 C     !LOCAL VARIABLES:
                0149 C     i,j       :: Loop counters
                0150       INTEGER i, j
                0151 CEOP
                0152       _RL recip_P0g, termP, rFullDepth
                0153       _RL kV, kF, sigma_b
                0154 
                0155 C--   Forcing term
                0156       kF = 1. _d 0/86400. _d 0
                0157       sigma_b = 0.7 _d 0
                0158       rFullDepth = rF(1)-rF(Nr+1)
                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
                0163         IF ( maskS(i,j,k,bi,bj).EQ.oneRS ) THEN
                0164          IF ( selectSigmaCoord.EQ.0 ) THEN
                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
                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:
                0178           termP = aHybSigmC(k)*rFullDepth
                0179 #ifdef NONLIN_FRSURF
                0180      &          /(etaHs(i,j,bi,bj)+rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
                0181 #else
                0182      &          /(rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj))
                0183 #endif
                0184      &          + bHybSigmC(k)
                0185          ENDIF
                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)
                0189         ENDIF
                0190        ENDDO
                0191       ENDDO
                0192 
                0193       RETURN
                0194       END
                0195 
                0196 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0197 CBOP
                0198 C     !ROUTINE: APPLY_FORCING_T
                0199 C     !INTERFACE:
                0200       SUBROUTINE APPLY_FORCING_T(
                0201      U                     gT_arr,
                0202      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0203      I                     myTime, myIter, myThid )
                0204 C     !DESCRIPTION: \bv
                0205 C     *==========================================================*
                0206 C     | S/R APPLY_FORCING_T
                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
                0213 
                0214 C     !USES:
                0215       IMPLICIT NONE
                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 
                0224 C     !INPUT/OUTPUT PARAMETERS:
                0225 C     gT_arr    :: the tendency array
                0226 C     iMin,iMax :: Working range of x-index for applying forcing.
                0227 C     jMin,jMax :: Working range of y-index for applying forcing.
                0228 C     k         :: Current vertical level index
                0229 C     bi,bj     :: Current tile indices
                0230 C     myTime    :: Current time in simulation
                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
                0238       INTEGER myThid
                0239 
                0240 C     !LOCAL VARIABLES:
                0241 C     i,j       :: Loop counters
                0242 C     kSurface  :: index of surface level
                0243       INTEGER i, j
                0244 CEOP
                0245       _RL thetaLim, kT, ka, ks, sigma_b, term1, term2, thetaEq
                0246       _RL termP, rFullDepth
                0247 
                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)
                0251       sigma_b = 0.7 _d 0
                0252       rFullDepth = rF(1)-rF(Nr+1)
                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)
                0258      &            *(COS(yC(i,j,bi,bj)*deg2rad)**2)
                0259          thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa)
                0260          thetaEq = 315. _d 0 - term1 - term2
                0261          thetaEq = MAX(thetaLim,thetaEq)
                0262          IF ( selectSigmaCoord.EQ.0 ) THEN
                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)
                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:
                0275           termP = aHybSigmC(k)*rFullDepth
                0276 #ifdef NONLIN_FRSURF
                0277      &          /(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
                0278 #else
                0279      &          /(Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
                0280 #endif
                0281      &          + bHybSigmC(k)
                0282          ENDIF
                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)
                0289        ENDDO
                0290       ENDDO
                0291 
                0292       RETURN
                0293       END
                0294 
                0295 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0296 CBOP
                0297 C     !ROUTINE: APPLY_FORCING_S
                0298 C     !INTERFACE:
                0299       SUBROUTINE APPLY_FORCING_S(
                0300      U                     gS_arr,
                0301      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0302      I                     myTime, myIter, myThid )
                0303 C     !DESCRIPTION: \bv
                0304 C     *==========================================================*
                0305 C     | S/R APPLY_FORCING_S
                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
                0312 
                0313 C     !USES:
                0314       IMPLICIT NONE
                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"
                0322 #include "SURFACE.h"
                0323 
                0324 C     !INPUT/OUTPUT PARAMETERS:
                0325 C     gS_arr    :: the tendency array
                0326 C     iMin,iMax :: Working range of x-index for applying forcing.
                0327 C     jMin,jMax :: Working range of y-index for applying forcing.
                0328 C     k         :: Current vertical level index
                0329 C     bi,bj     :: Current tile indices
                0330 C     myTime    :: Current time in simulation
                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
                0338       INTEGER myThid
                0339 
                0340 C     !LOCAL VARIABLES:
                0341 C     i,j       :: Loop counters
                0342 c     INTEGER i, j
                0343 CEOP
                0344 
                0345 C--   Forcing term
                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!