pkg/seaice/diffus.F

809c36b928 Patr*0001 #include "SEAICE_OPTIONS.h"
37e58bbaf7 Jean*0002
b7e868d5f4 Jean*0003 CBOP
                0004 C     !ROUTINE: DIFFUS
                0005 C     !INTERFACE:
                0006       SUBROUTINE DIFFUS(
                0007      U                   fld,
                0008      I                   DIFFA, iceMsk, myThid )
                0009
                0010 C     !DESCRIPTION: \bv
37e58bbaf7 Jean*0011 C     *==========================================================*
b7e868d5f4 Jean*0012 C     | S/R DIFFUS
                0013 C     | o Calculate laplacian of input field
                0014 C     |   and return the result in the same array
37e58bbaf7 Jean*0015 C     *==========================================================*
b7e868d5f4 Jean*0016 C     \ev
                0017
                0018 C     !USES:
809c36b928 Patr*0019       IMPLICIT NONE
37e58bbaf7 Jean*0020
b7e868d5f4 Jean*0021 C     == Global variables ===
809c36b928 Patr*0022 #include "SIZE.h"
                0023 #include "EEPARAMS.h"
4366d31d92 Mart*0024 #include "GRID.h"
7303eab4f2 Patr*0025 #include "SEAICE_SIZE.h"
9cf7ff39d8 Mart*0026 #include "SEAICE_PARAMS.h"
37e58bbaf7 Jean*0027
b7e868d5f4 Jean*0028 C     !INPUT/OUTPUT PARAMETERS:
                0029 C     == Routine Arguments ==
                0030 C     fld        :: In: input ice-field ; Out: laplacian of input field
                0031 C     DIFFA      :: grid length scale
                0032 C     iceMsk     :: Ocean/Land mask
                0033 C     myThid     :: my Thread Id. number
1a11ffcc11 Mart*0034       _RL fld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0035       _RL DIFFA  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
b7e868d5f4 Jean*0036       _RL iceMsk (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
809c36b928 Patr*0037       INTEGER myThid
b7e868d5f4 Jean*0038 CEOP
809c36b928 Patr*0039
b7e868d5f4 Jean*0040 C     !LOCAL VARIABLES:
                0041 C     == Local variables ==
                0042 C     i,j,bi,bj :: Loop counters
cee16b76ae Dimi*0043       INTEGER i, j, bi, bj
b7e868d5f4 Jean*0044       _RL DELTXX, DELTYY
                0045       _RL tmpFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
04a836ae4f Mart*0046       _RL dfx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0047       _RL dfy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
809c36b928 Patr*0048
b7e868d5f4 Jean*0049 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0050
809c36b928 Patr*0051       DO bj=myByLo(myThid),myByHi(myThid)
                0052        DO bi=myBxLo(myThid),myBxHi(myThid)
b7e868d5f4 Jean*0053         IF ( SEAICEuseFluxForm ) THEN
                0054 C--   Use flux form for MITgcm compliance, unfortunately changes results
809c36b928 Patr*0055
b7e868d5f4 Jean*0056           DO j=1-Oly,sNy+Oly
                0057            DO i=1-Olx,sNx+Olx
                0058             dfx(i,j) = 0. _d 0
                0059             dfy(i,j) = 0. _d 0
                0060            ENDDO
                0061           ENDDO
                0062 C--   first compute fluxes across cell faces
                0063           DO j=1,sNy+1
                0064            DO i=1,sNx+1
                0065             dfx(i,j) = _dyG(i,j,bi,bj) * _recip_dxC(i,j,bi,bj)
                0066      &         * (fld(i,j,bi,bj)-fld(i-1,j,bi,bj))
                0067      &         * cosFacU(j,bi,bj)
                0068      &         * iceMsk(i,j,bi,bj)*iceMsk(i-1,j,bi,bj)
                0069      &         * ( DIFFA(i,j,bi,bj)+DIFFA(i-1,j,bi,bj) )*HALF
036cdc98c8 Jean*0070 #ifdef ALLOW_OBCS
b7e868d5f4 Jean*0071      &         * maskInW(i,j,bi,bj)
036cdc98c8 Jean*0072 #endif
b7e868d5f4 Jean*0073             dfy(i,j) = _dxG(i,j,bi,bj) * _recip_dyC(i,j,bi,bj)
                0074      &         * (fld(i,j,bi,bj)-fld(i,j-1,bi,bj))
                0075 #ifdef ISOTROPIC_COS_SCALING
                0076      &         * cosFacV(j,bi,bj)
036cdc98c8 Jean*0077 #endif
b7e868d5f4 Jean*0078      &         * iceMsk(i,j,bi,bj)*iceMsk(i,j-1,bi,bj)
                0079      &         * ( DIFFA(i,j,bi,bj)+DIFFA(i,j-1,bi,bj) )*HALF
036cdc98c8 Jean*0080 #ifdef ALLOW_OBCS
b7e868d5f4 Jean*0081      &         * maskInS(i,j,bi,bj)
036cdc98c8 Jean*0082 #endif
b7e868d5f4 Jean*0083            ENDDO
                0084           ENDDO
                0085           DO j=1-OLy,sNy+OLy
                0086            DO i=1-OLx,sNx+OLx
                0087             fld(i,j,bi,bj) = 0. _d 0
                0088            ENDDO
                0089           ENDDO
                0090 C--   compute Laplacian as flux divergence
                0091           DO j=1,sNy
                0092            DO i=1,sNx
                0093             fld(i,j,bi,bj) = (
                0094      &                         ( dfx(i+1,j) - dfx(i,j) )
                0095      &                       + ( dfy(i,j+1) - dfy(i,j) )
                0096      &                       ) * recip_rA(i,j,bi,bj)
                0097            ENDDO
                0098           ENDDO
                0099
                0100         ELSE
                0101 C NOW DO DIFFUSION WITH NUIT CONVERSION
                0102
                0103           DO j=1-OLy,sNy+OLy
                0104            DO i=1-OLx,sNx+OLx
                0105             tmpFld(i,j) = 0.0 _d 0
                0106            ENDDO
                0107           ENDDO
                0108
                0109           DO j=1,sNy
                0110            DO i=1,sNx
                0111             DELTXX = DIFFA(i,j,bi,bj)
                0112      &             * _recip_dxF(i,j,bi,bj)*_recip_dxF(i,j,bi,bj)
                0113             DELTYY = DIFFA(i,j,bi,bj)
                0114      &             * _recip_dyF(i,j,bi,bj)*_recip_dyF(i,j,bi,bj)
                0115      &             * _recip_dxF(i,j,bi,bj)
                0116             tmpFld(i,j) =
                0117      &        DELTXX*(
                0118      &                 (fld(i+1,j,bi,bj)-fld(i,  j,bi,bj))
                0119      &                 *iceMsk(i+1,j,bi,bj)
036cdc98c8 Jean*0120 #ifdef ALLOW_OBCS
b7e868d5f4 Jean*0121      &                 *maskInW(i+1,j,bi,bj)
036cdc98c8 Jean*0122 #endif
b7e868d5f4 Jean*0123      &                -(fld(i,  j,bi,bj)-fld(i-1,j,bi,bj))
                0124      &                 *iceMsk(i-1,j,bi,bj)
036cdc98c8 Jean*0125 #ifdef ALLOW_OBCS
b7e868d5f4 Jean*0126      &                 *maskInW(i,j,bi,bj)
036cdc98c8 Jean*0127 #endif
b7e868d5f4 Jean*0128      &               )
                0129      &       +DELTYY*(
                0130      &                 (fld(i,j+1,bi,bj)-fld(i,j,  bi,bj))
                0131      &                 * _dxG(i,j+1,bi,bj)*iceMsk(i,j+1,bi,bj)
                0132 #ifdef ALLOW_OBCS
                0133      &                 *maskInS(i,j+1,bi,bj)
37e58bbaf7 Jean*0134 #endif
b7e868d5f4 Jean*0135      &                -(fld(i,j,  bi,bj)-fld(i,j-1,bi,bj))
                0136      &                 * _dxG(i,j,  bi,bj)*iceMsk(i,j-1,bi,bj)
036cdc98c8 Jean*0137 #ifdef ALLOW_OBCS
b7e868d5f4 Jean*0138      &                 *maskInS(i,j,bi,bj)
036cdc98c8 Jean*0139 #endif
b7e868d5f4 Jean*0140      &               )
                0141            ENDDO
                0142           ENDDO
                0143           DO j=1-OLy,sNy+OLy
                0144            DO i=1-OLx,sNx+OLx
                0145             fld(i,j,bi,bj) = tmpFld(i,j)
                0146            ENDDO
                0147           ENDDO
809c36b928 Patr*0148
b7e868d5f4 Jean*0149 C--  end flux-form / non flux-form
                0150         ENDIF
809c36b928 Patr*0151        ENDDO
                0152       ENDDO
                0153
                0154       RETURN
                0155       END