Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:38:18 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
138482fdf6 Ed H*0001 #include "CD_CODE_OPTIONS.h"
d2f0d8534f Ed H*0002 
                0003 CBOP
c0d4d5aecc Patr*0004 C !ROUTINE: CD_CODE_SCHEME
d2f0d8534f Ed H*0005 
                0006 C !INTERFACE: ==========================================================
0aa8fbef89 Jean*0007       SUBROUTINE CD_CODE_SCHEME(
d2f0d8534f Ed H*0008      I        bi,bj,k, dPhiHydX,dPhiHydY, guFld,gvFld,
                0009      O        guCor,gvCor,
                0010      I        myTime, myIter, myThid)
                0011 
                0012 C !DESCRIPTION:
                0013 C The C-D scheme. The less said the better :-)
                0014 
                0015 C !USES: ===============================================================
                0016 C     == Global variables ==
                0017       IMPLICIT NONE
                0018 #include "SIZE.h"
                0019 #include "DYNVARS.h"
a5bd9cbb61 Ed H*0020 #ifdef ALLOW_CD_CODE
                0021 #include "CD_CODE_VARS.h"
                0022 #endif
d2f0d8534f Ed H*0023 #include "EEPARAMS.h"
                0024 #include "PARAMS.h"
                0025 #include "GRID.h"
                0026 #include "SURFACE.h"
                0027 
                0028 C !INPUT PARAMETERS: ===================================================
                0029 C  bi,bj                :: tile indices
                0030 C  k                    :: vertical level
                0031 C     dPhiHydX,Y        :: Gradient (X & Y dir.) of Hydrostatic Potential
                0032 C  guFld,gvFld          :: Acceleration (U & V compon.) from the C grid
                0033 C  guCor,gvCor          :: Coriolis terms (2 compon.) computed on C grid
                0034 C  myTime               :: current time
                0035 C  myIter               :: current time-step number
                0036 C  myThid               :: thread number
                0037 
                0038       INTEGER bi,bj,k
a30e0ce2b4 Jean*0039       _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0040       _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0041       _RL    guFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0042       _RL    gvFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0043       _RL    guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0044       _RL    gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d2f0d8534f Ed H*0045       _RL     myTime
                0046       INTEGER myIter
                0047       INTEGER myThid
                0048 
                0049 
                0050 C !LOCAL VARIABLES: ====================================================
                0051 #ifdef ALLOW_CD_CODE
                0052 C  i,j                  :: loop indices
                0053 C  pF                   :: pressure gradient
                0054 C  vF                   :: work space
                0055 C  aF                   :: work space
0aa8fbef89 Jean*0056       INTEGER i,j
d2f0d8534f Ed H*0057       _RL pF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0058       _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0059       _RL aF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0060       _RL ab15,ab05
                0061       _RL phxFac, phyFac
0aa8fbef89 Jean*0062 C     Define ranges
                0063       INTEGER iMin,iMax, jMin,jMax
a30e0ce2b4 Jean*0064       PARAMETER( iMin = 1-OLx+1 , iMax = sNx+OLx-1 )
                0065       PARAMETER( jMin = 1-OLy+1 , jMax = sNy+OLy-1 )
d2f0d8534f Ed H*0066 CEOP
                0067 
                0068 C     Adams-Bashforth weighting factors
2434c713b2 Jean*0069       IF ( myIter.EQ.0 ) THEN
                0070         ab15   =  1. _d 0
                0071         ab05   = -0. _d 0
                0072       ELSE
                0073         ab15   =  1.5 _d 0 + epsAB_CD
                0074         ab05   = -0.5 _d 0 - epsAB_CD
                0075       ENDIF
d2f0d8534f Ed H*0076 
                0077 C-- stagger time stepping: grad Phi_Hyp is not in gU,gV and needs to be added:
                0078       IF (staggerTimeStep) THEN
                0079         phxFac = pfFacMom
                0080         phyFac = pfFacMom
                0081       ELSE
                0082         phxFac = 0.
                0083         phyFac = 0.
                0084       ENDIF
                0085 
                0086 C-    Initialize output (dummy) arrays:
a30e0ce2b4 Jean*0087 c     DO j=1-OLy,sNy+OLy
                0088 c      DO i=1-OLx,sNx+OLx
d2f0d8534f Ed H*0089 c       guCor(i,j) = 0. _d 0
                0090 c       gvCor(i,j) = 0. _d 0
                0091 c      ENDDO
                0092 c     ENDDO
                0093 
                0094 C     Pressure extrapolated forward in time
a30e0ce2b4 Jean*0095       DO j=1-OLy,sNy+OLy
                0096        DO i=1-OLx,sNx+OLx
da4112fbea Jean*0097 #ifdef CD_CODE_NO_AB_CORIOLIS
                0098 C     Keep this just to reproduce old results (to get same truncation)
0aa8fbef89 Jean*0099         pf(i,j) =
d2f0d8534f Ed H*0100      &   ab15*(  etaN(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
                0101      &  +ab05*(etaNm1(i,j,bi,bj)*Bo_surf(i,j,bi,bj) )
da4112fbea Jean*0102 #else /* CD_CODE_NO_AB_CORIOLIS */
                0103         pf(i,j) = Bo_surf(i,j,bi,bj)
                0104      &          *( ab15*etaN(i,j,bi,bj) + ab05*etaNm1(i,j,bi,bj) )
                0105 #endif /* CD_CODE_NO_AB_CORIOLIS */
d2f0d8534f Ed H*0106        ENDDO
                0107       ENDDO
                0108 
                0109 C--   Zonal velocity coriolis term
                0110 C     Note. As coded here, coriolis will not work with "thin walls"
                0111 C--   Coriolis with CD scheme allowed
                0112 C     grady(p) + gV
a30e0ce2b4 Jean*0113       DO j=1-OLy+1,sNy+OLy
                0114        DO i=1-OLx,sNx+OLx
dd55623897 Jean*0115         af(i,j) =
                0116      &        (  gvFld(i,j)
                0117      &          -( _recip_dyC(i,j,bi,bj)*(pf(i,j)-pf(i,j-1))
                0118      &            +phyFac*dPhiHydY(i,j) )
                0119      &        )*_maskS(i,j,k,bi,bj)
a30e0ce2b4 Jean*0120 #ifdef ALLOW_OBCS
                0121      &         *maskInC(i,j-1,bi,bj)*maskInC(i,j,bi,bj)
                0122 #endif
d2f0d8534f Ed H*0123        ENDDO
                0124       ENDDO
                0125 C     Average to Vd point and add coriolis
                0126       DO j=jMin,jMax
                0127        DO i=iMin,iMax
0aa8fbef89 Jean*0128         vf(i,j) =
                0129      &           ( (af(i,j)+af(i-1,j+1))
                0130      &            +(af(i-1,j)+af(i,j+1)) )*0.25 _d 0
                0131      &           *_maskW(i,j,k,bi,bj)
                0132      &          -( _fCori( i, j,bi,bj)
                0133      &            +_fCori(i-1,j,bi,bj) )*0.5 _d 0
da4112fbea Jean*0134 #ifdef CD_CODE_NO_AB_CORIOLIS
0aa8fbef89 Jean*0135      &           *uVel(i,j,k,bi,bj)
da4112fbea Jean*0136 #else /* CD_CODE_NO_AB_CORIOLIS */
                0137      &           *( ab15*uVel(i,j,k,bi,bj) + ab05*uNM1(i,j,k,bi,bj) )
                0138 #endif /* CD_CODE_NO_AB_CORIOLIS */
d2f0d8534f Ed H*0139        ENDDO
                0140       ENDDO
                0141 C     Step forward Vd
                0142       DO j=jMin,jMax
                0143        DO i=iMin,iMax
                0144         vVelD(i,j,k,bi,bj) = vVelD(i,j,k,bi,bj) + deltaTmom*vf(i,j)
                0145        ENDDO
                0146       ENDDO
                0147 C     Relax D grid V to C grid V
                0148       DO j=jMin,jMax
                0149        DO i=iMin,iMax
                0150          vVelD(i,j,k,bi,bj) = ( rCD*vVelD(i,j,k,bi,bj)
0aa8fbef89 Jean*0151      &                         +(1. _d 0 - rCD)
                0152      &         *( ab15*(
                0153      &                  (vVel(i,j,k,bi,bj)+vVel(i-1,j+1,k,bi,bj))
                0154      &                 +(vVel(i-1,j,k,bi,bj)+vVel(i,j+1,k,bi,bj))
                0155      &                 )*0.25 _d 0
                0156      &           +ab05*(
                0157      &                  (vNM1(i,j,k,bi,bj)+vNM1(i-1,j+1,k,bi,bj))
                0158      &                 +(vNM1(i-1,j,k,bi,bj)+vNM1(i,j+1,k,bi,bj))
                0159      &                 )*0.25 _d 0
                0160      &          )             )*_maskW(i,j,k,bi,bj)
d2f0d8534f Ed H*0161        ENDDO
                0162       ENDDO
                0163 C     Calculate coriolis force on U
                0164       DO j=jMin,jMax
                0165        DO i=iMin,iMax
0aa8fbef89 Jean*0166         guCor(i,j) =
                0167      &              ( _fCori( i, j,bi,bj)
                0168      &               +_fCori(i-1,j,bi,bj) )*0.5 _d 0
                0169      &              *vVelD(i,j,k,bi,bj)*cfFacMom
d2f0d8534f Ed H*0170        ENDDO
                0171       ENDDO
                0172 
                0173 C--   Meridional velocity coriolis term
                0174 C     gradx(p)+gU
a30e0ce2b4 Jean*0175       DO j=1-OLy,sNy+OLy
                0176        DO i=1-OLx+1,sNx+OLx
dd55623897 Jean*0177         af(i,j) =
                0178      &        (  guFld(i,j)
                0179      &          -( _recip_dxC(i,j,bi,bj)*(pf(i,j)-pf(i-1,j))
                0180      &            +phxFac*dPhiHydX(i,j) )
                0181      &        )*_maskW(i,j,k,bi,bj)
a30e0ce2b4 Jean*0182 #ifdef ALLOW_OBCS
                0183      &         *maskInC(i-1,j,bi,bj)*maskInC(i,j,bi,bj)
                0184 #endif
d2f0d8534f Ed H*0185        ENDDO
                0186       ENDDO
                0187 C     Average to Ud point and add coriolis
                0188       DO j=jMin,jMax
                0189        DO i=iMin,iMax
                0190         vf(i,j) =
0aa8fbef89 Jean*0191      &           ( (af(i,j)+af(i+1,j-1))
                0192      &            +(af(i+1,j)+af(i,j-1)) )*0.25 _d 0
                0193      &           *_maskS(i,j,k,bi,bj)
                0194      &          +( _fCori(i, j, bi,bj)
                0195      &            +_fCori(i,j-1,bi,bj) )*0.5 _d 0
da4112fbea Jean*0196 #ifdef CD_CODE_NO_AB_CORIOLIS
0aa8fbef89 Jean*0197      &           *vVel(i,j,k,bi,bj)
da4112fbea Jean*0198 #else /* CD_CODE_NO_AB_CORIOLIS */
                0199      &           *( ab15*vVel(i,j,k,bi,bj) + ab05*vNM1(i,j,k,bi,bj) )
                0200 #endif /* CD_CODE_NO_AB_CORIOLIS */
d2f0d8534f Ed H*0201        ENDDO
                0202       ENDDO
                0203 C     Step forward Ud
                0204       DO j=jMin,jMax
                0205        DO i=iMin,iMax
                0206         uVelD(i,j,k,bi,bj) = uVelD(i,j,k,bi,bj) + deltaTmom*vf(i,j)
                0207        ENDDO
                0208       ENDDO
                0209 C     Relax D grid U to C grid U
                0210       DO j=jMin,jMax
                0211        DO i=iMin,iMax
0aa8fbef89 Jean*0212          uVelD(i,j,k,bi,bj) = ( rCD*uVelD(i,j,k,bi,bj)
                0213      &                         +(1. _d 0 - rCD)
                0214      &         *( ab15*(
                0215      &                  (uVel(i,j,k,bi,bj)+uVel(i+1,j-1,k,bi,bj))
                0216      &                 +(uVel(i,j-1,k,bi,bj)+uVel(i+1,j,k,bi,bj))
                0217      &                 )*0.25 _d 0
                0218      &           +ab05*(
                0219      &                  (uNM1(i,j,k,bi,bj)+uNM1(i+1,j-1,k,bi,bj))
                0220      &                 +(uNM1(i,j-1,k,bi,bj)+uNM1(i+1,j,k,bi,bj))
                0221      &                 )*0.25 _d 0
                0222      &          )             )*_maskS(i,j,k,bi,bj)
d2f0d8534f Ed H*0223        ENDDO
                0224       ENDDO
                0225 C     Calculate coriolis force on V
                0226       DO j=jMin,jMax
                0227        DO i=iMin,iMax
                0228         gvCor(i,j) =
0aa8fbef89 Jean*0229      &             -( _fCori(i, j, bi,bj)
                0230      &               +_fCori(i,j-1,bi,bj) )*0.5 _d 0
                0231      &              *uVelD(i,j,k,bi,bj)*cfFacMom
d2f0d8534f Ed H*0232        ENDDO
                0233       ENDDO
                0234 
                0235 C--   Save "previous time level" variables
                0236       DO j=1-OLy,sNy+OLy
                0237        DO i=1-OLx,sNx+OLx
0aa8fbef89 Jean*0238          uNM1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
                0239          vNM1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
d2f0d8534f Ed H*0240        ENDDO
                0241       ENDDO
                0242 
                0243 #endif /* ALLOW_CD_CODE */
                0244 
                0245       RETURN
                0246       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!