Back to home page

MITgcm
 
 

    

File indexing completed on 2022-04-21 05:09:20 UTC

view on githubraw file Latest commit 113b21cd on 2022-04-20 17:14:29 UTC
809c36b928 Patr*0001 #include "SEAICE_OPTIONS.h"
772b2ed80e Gael*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
09510da3bb Dimi*0005 
b7e868d5f4 Jean*0006 CBOP
                0007 C     !ROUTINE: ADVECT
                0008 C     !INTERFACE:
                0009       SUBROUTINE ADVECT(
                0010      I                   UI, VI,
                0011      U                   fld,
                0012      O                   fldNm1,
                0013      I                   iceMsk, myThid )
                0014 
                0015 C     !DESCRIPTION: \bv
3ee9cb5aa5 Jean*0016 C     *==========================================================*
b7e868d5f4 Jean*0017 C     | S/R ADVECT
                0018 C     | o Calculate advection and diffusion
                0019 C     |   and update the input ice-field
3ee9cb5aa5 Jean*0020 C     *==========================================================*
b7e868d5f4 Jean*0021 C     \ev
                0022 
                0023 C     !USES:
809c36b928 Patr*0024       IMPLICIT NONE
09510da3bb Dimi*0025 
b7e868d5f4 Jean*0026 C     == Global variables ===
809c36b928 Patr*0027 #include "SIZE.h"
                0028 #include "EEPARAMS.h"
                0029 #include "PARAMS.h"
4366d31d92 Mart*0030 #include "GRID.h"
7303eab4f2 Patr*0031 #include "SEAICE_SIZE.h"
809c36b928 Patr*0032 #include "SEAICE_PARAMS.h"
7109a141b2 Patr*0033 #ifdef ALLOW_AUTODIFF_TAMC
                0034 # include "tamc.h"
                0035 #endif
                0036 
b7e868d5f4 Jean*0037 C     !INPUT/OUTPUT PARAMETERS:
                0038 C     == Routine arguments ==
                0039 C     UI, VI     :: ice velocity components
                0040 C     fld        :: input and updated ice-field
                0041 C     fldNm1     :: copy of the input ice-field
                0042 C     iceMsk     :: Ocean/Land mask
                0043 C     myThid     :: my Thread Id. number
8744e57318 Mart*0044       _RL UI     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0045       _RL VI     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0046       _RL fld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0047       _RL fldNm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0048       _RL iceMsk (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
809c36b928 Patr*0049       INTEGER myThid
b7e868d5f4 Jean*0050 CEOP
809c36b928 Patr*0051 
b7e868d5f4 Jean*0052 C     !LOCAL VARIABLES:
                0053 C     == Local variables ==
                0054 C     i,j,k,bi,bj :: Loop counters
baa476eeba Dimi*0055       INTEGER i, j, bi, bj
8744e57318 Mart*0056       INTEGER k
3cb5d2b11d Mart*0057       _RL DELTT
8744e57318 Mart*0058       _RL DIFFA  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0059       _RL tmpFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
04a836ae4f Mart*0060       _RL afx    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0061       _RL afy    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
809c36b928 Patr*0062 
8744e57318 Mart*0063       DELTT=SEAICE_deltaTtherm
                0064 C     save fld from previous time step
809c36b928 Patr*0065       DO bj=myByLo(myThid),myByHi(myThid)
                0066        DO bi=myBxLo(myThid),myBxHi(myThid)
                0067         DO j=1-OLy,sNy+OLy
                0068          DO i=1-OLx,sNx+OLx
b7e868d5f4 Jean*0069           fldNm1(i,j,bi,bj) = fld(i,j,bi,bj)
809c36b928 Patr*0070          ENDDO
                0071         ENDDO
                0072        ENDDO
                0073       ENDDO
                0074 
113b21cd45 Mart*0075 #ifdef ALLOW_AUTODIFF_TAMC
                0076 CADJ INIT loctape_advect = COMMON, 2
                0077 #endif
8744e57318 Mart*0078       DO k=1,2
113b21cd45 Mart*0079 #ifdef ALLOW_AUTODIFF_TAMC
                0080 CADJ STORE fld = loctape_advect, key = k
                0081 #endif
8744e57318 Mart*0082 C     Backward Euler correction step
                0083         DO bj=myByLo(myThid),myByHi(myThid)
                0084          DO bi=myBxLo(myThid),myBxHi(myThid)
                0085           DO j=1-OLy,sNy+OLy
                0086            DO i=1-OLx,sNx+OLx
b7e868d5f4 Jean*0087             tmpFld(i,j,bi,bj)=HALF*(fld(i,j,bi,bj)
8744e57318 Mart*0088      &           +fldNm1(i,j,bi,bj))
baa476eeba Dimi*0089            ENDDO
                0090           ENDDO
                0091          ENDDO
8744e57318 Mart*0092         ENDDO
                0093 
113b21cd45 Mart*0094 #ifdef ALLOW_AUTODIFF
                0095        DO j=1-OLy,sNy+OLy
                0096         DO i=1-OLx,sNx+OLx
b7e868d5f4 Jean*0097          afx(i,j) = 0. _d 0
                0098          afy(i,j) = 0. _d 0
04a836ae4f Mart*0099         ENDDO
                0100        ENDDO
113b21cd45 Mart*0101 #endif /* ALLOW_AUTODIFF */
baa476eeba Dimi*0102 
                0103 C NOW GO THROUGH STANDARD CONSERVATIVE ADVECTION
8744e57318 Mart*0104        IF ( .NOT. SEAICEuseFluxForm ) THEN
                0105         DO bj=myByLo(myThid),myByHi(myThid)
                0106          DO bi=myBxLo(myThid),myBxHi(myThid)
b7e868d5f4 Jean*0107           DO j=0,sNy+1
                0108            DO i=0,sNx+1
4366d31d92 Mart*0109 CML   This formulation gives the same result as the original code on a
                0110 CML   lat-lon-grid, but may not be accurate on irregular grids
b7e868d5f4 Jean*0111             fld(i,j,bi,bj)=fldNm1(i,j,bi,bj)
8744e57318 Mart*0112      &           -DELTT*(
b7e868d5f4 Jean*0113      &           ( tmpFld(i  ,j  ,bi,bj)+tmpFld(i+1,j  ,bi,bj))
                0114      &           *   UI(i+1,j,  bi,bj) -
                0115      &           ( tmpFld(i  ,j  ,bi,bj)+tmpFld(i-1,j  ,bi,bj))
                0116      &           *   UI(i  ,j,  bi,bj) )*maskInC(i,j,bi,bj)
                0117      &           *(HALF * _recip_dxF(i,j,bi,bj))
8744e57318 Mart*0118      &           -DELTT*(
b7e868d5f4 Jean*0119      &           ( tmpFld(i  ,j  ,bi,bj)+tmpFld(i  ,j+1,bi,bj))
                0120      &           *   VI(i  ,j+1,  bi,bj)
                0121      &           * _dxG(i  ,j+1,bi,bj) -
                0122      &           ( tmpFld(i  ,j  ,bi,bj)+tmpFld(i  ,j-1,bi,bj))
                0123      &           *   VI(i  ,j  ,  bi,bj)
                0124      &           * _dxG(i,j,bi,bj))*maskInC(i,j,bi,bj)
                0125      &           *(HALF * _recip_dyF(i,j,bi,bj) * _recip_dxF(i,j,bi,bj))
9cf7ff39d8 Mart*0126            ENDDO
                0127           ENDDO
                0128          ENDDO
8744e57318 Mart*0129         ENDDO
                0130        ELSE
04a836ae4f Mart*0131 C--   Use flux form for MITgcm compliance, unfortunately changes results
8744e57318 Mart*0132         DO bj=myByLo(myThid),myByHi(myThid)
                0133          DO bi=myBxLo(myThid),myBxHi(myThid)
04a836ae4f Mart*0134 C--   first compute fluxes across cell faces
b7e868d5f4 Jean*0135           DO j=1,sNy+1
                0136            DO i=1,sNx+1
                0137             afx(i,j) = _dyG(i,j,bi,bj) * UI(i,j,bi,bj)
                0138      &           * 0.5 _d 0 * (tmpFld(i,j,bi,bj)+tmpFld(i-1,j,bi,bj))
                0139             afy(i,j) = _dxG(i,j,bi,bj) * VI(i,j,bi,bj)
                0140      &           * 0.5 _d 0 * (tmpFld(i,j,bi,bj)+tmpFld(i,j-1,bi,bj))
04a836ae4f Mart*0141            ENDDO
                0142           ENDDO
b7e868d5f4 Jean*0143           DO j=1,sNy
                0144            DO i=1,sNx
                0145             fld(i,j,bi,bj)=fldNm1(i,j,bi,bj)
70baca3a4f Jean*0146      &           -DELTT * (
b7e868d5f4 Jean*0147      &             afx(i+1,j) - afx(i,j)
                0148      &           + afy(i,j+1) - afy(i,j)
                0149      &           )*recip_rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
baa476eeba Dimi*0150            ENDDO
                0151           ENDDO
                0152          ENDDO
8744e57318 Mart*0153         ENDDO
                0154        ENDIF
809c36b928 Patr*0155 
d2c42c7d63 Jean*0156        CALL EXCH_XY_RL( fld, myThid )
baa476eeba Dimi*0157 
809c36b928 Patr*0158       ENDDO
                0159 
b7e868d5f4 Jean*0160       IF ( DIFF1.GT.0. _d 0 ) THEN
70baca3a4f Jean*0161 C NOW DO DIFFUSION
b7e868d5f4 Jean*0162 
8744e57318 Mart*0163 C     make a working copy of field from last time step
b7e868d5f4 Jean*0164        DO bj=myByLo(myThid),myByHi(myThid)
                0165         DO bi=myBxLo(myThid),myBxHi(myThid)
                0166          DO j=1-OLy,sNy+OLy
                0167           DO i=1-OLx,sNx+OLx
                0168            tmpFld(i,j,bi,bj) = fldNm1(i,j,bi,bj)
                0169           ENDDO
809c36b928 Patr*0170          ENDDO
                0171         ENDDO
                0172        ENDDO
b7e868d5f4 Jean*0173 
809c36b928 Patr*0174 C NOW CALCULATE DIFFUSION COEF ROUGHLY
70baca3a4f Jean*0175 C  1rst pass: compute changes due to harmonic diffusion and add it to ice-field
b7e868d5f4 Jean*0176        DO bj=myByLo(myThid),myByHi(myThid)
                0177         DO bi=myBxLo(myThid),myBxHi(myThid)
8744e57318 Mart*0178           DO j=1-OLy,sNy+OLy
                0179            DO i=1-OLx,sNx+OLx
b7e868d5f4 Jean*0180             DIFFA(i,j,bi,bj) = MIN( _dxF(i,j,bi,bj), _dyF(i,j,bi,bj) )
8744e57318 Mart*0181            ENDDO
                0182           ENDDO
b7e868d5f4 Jean*0183         ENDDO
                0184        ENDDO
                0185 C-     Compute laplacian of ice-field; return result in same array
                0186        CALL DIFFUS( tmpFld, DIFFA, iceMsk, myThid )
                0187 
                0188        DO bj=myByLo(myThid),myByHi(myThid)
                0189         DO bi=myBxLo(myThid),myBxHi(myThid)
                0190          DO j=1-OLy,sNy+OLy
                0191           DO i=1-OLx,sNx+OLx
                0192            fld(i,j,bi,bj) = ( fld(i,j,bi,bj)
                0193      &                       +tmpFld(i,j,bi,bj)*DIFF1*DELTT
                0194      &                      )*iceMsk(i,j,bi,bj)
                0195           ENDDO
809c36b928 Patr*0196          ENDDO
                0197         ENDDO
b7e868d5f4 Jean*0198        ENDDO
                0199 
                0200 c     IF ( useBiHarmonic ) THEN
                0201 C  2nd  pass: compute changes due to biharmonic diffusion and add it to ice-field
                0202        _EXCH_XY_RL( tmpFld, myThid )
70baca3a4f Jean*0203 C     use some strange quadratic form for the second time around
b7e868d5f4 Jean*0204        DO bj=myByLo(myThid),myByHi(myThid)
                0205         DO bi=myBxLo(myThid),myBxHi(myThid)
8744e57318 Mart*0206           DO j=1-OLy,sNy+OLy
                0207            DO i=1-OLx,sNx+OLx
b7e868d5f4 Jean*0208 #ifdef ALLOW_AUTODIFF_TAMC
                0209 C      to avoid recomputations when there was a k loop; not needed anymore
                0210 c           DIFFA(i,j,bi,bj) = MIN( _dxF(i,j,bi,bj), _dyF(i,j,bi,bj) )
                0211 #endif
                0212             DIFFA(i,j,bi,bj) = - DIFFA(i,j,bi,bj)*DIFFA(i,j,bi,bj)
8744e57318 Mart*0213            ENDDO
                0214           ENDDO
                0215         ENDDO
b7e868d5f4 Jean*0216        ENDDO
                0217 C-     Compute bilaplacian (laplacian of laplacian); return result in same array
                0218        CALL DIFFUS( tmpFld, DIFFA, iceMsk, myThid )
8744e57318 Mart*0219 
                0220        DO bj=myByLo(myThid),myByHi(myThid)
                0221         DO bi=myBxLo(myThid),myBxHi(myThid)
                0222          DO j=1-OLy,sNy+OLy
                0223           DO i=1-OLx,sNx+OLx
b7e868d5f4 Jean*0224            fld(i,j,bi,bj) = ( fld(i,j,bi,bj)
                0225      &                       +tmpFld(i,j,bi,bj)*DIFF1*DELTT
                0226      &                      )*iceMsk(i,j,bi,bj)
8744e57318 Mart*0227           ENDDO
809c36b928 Patr*0228          ENDDO
                0229         ENDDO
                0230        ENDDO
b7e868d5f4 Jean*0231 C--   end biharmonic block
                0232 c     ENDIF
70baca3a4f Jean*0233 
b7e868d5f4 Jean*0234 C--   end DIFF1 > 0 block
                0235       ENDIF
809c36b928 Patr*0236 
                0237       RETURN
                0238       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!