pkg/layers/layers_diapycnal.F

f0f170d54b Mart*0001 #include "LAYERS_OPTIONS.h"
                0002
                0003 CBOP 0
                0004 C     !ROUTINE: LAYERS_DIAPYCNAL
                0005 C     !INTERFACE:
                0006       SUBROUTINE LAYERS_DIAPYCNAL(
                0007      I                  tracer, iLa,
                0008      O                  TtendSurf, TtendDiffh, TtendDiffr,
                0009      O                  TtendAdvh, TtendAdvr, TtendTot,
                0010      O                  StendSurf, StendDiffh, StendDiffr,
                0011      O                  StendAdvh, StendAdvr, StendTot,
                0012      O                  Hc, PIc,
                0013      I                  myThid )
                0014
                0015 C     !DESCRIPTION: \bv
                0016 C     *==========================================================*
                0017 C     | SUBROUTINE LAYERS_DIAPYCNAL
                0018 C     | Calculate the diapycnal velocity in isotracer layers,
                0019 C     | for a chosen tracer.
                0020 C     *==========================================================*
                0021 C     \ev
                0022       IMPLICIT NONE
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
                0027 #include "LAYERS_SIZE.h"
                0028 #include "LAYERS.h"
                0029
                0030 C !INPUT PARAMETERS:
                0031 C     tracer     :: potential temperature, salt or potential density prho
                0032 C     iLa        :: layer coordinate index
                0033 C     TtendSurf  :: Theta tendency due to surface forcing times thickness
                0034 C     TtendDiffh :: Theta tendency due to horiz. diffusion times thickness
                0035 C     TtendDiffr :: Theta tendency due to vert. diffusion times thickness
                0036 C     TtendAdvh  :: Theta tendency due to horiz. advection times thickness
                0037 C     TtendAdvr  :: Theta tendency due to vert. advection times thickness
                0038 C     TtendTot   :: Theta Total tendency
                0039 C     StendSurf  :: Salt  tendency due to surface forcing times thickness
                0040 C     StendDiffh :: Salt  tendency due to horiz. diffusion times thickness
                0041 C     StendDiffr :: Salt  tendency due to vert. diffusion times thickness
                0042 C     StendAdvh  :: Salt  tendency due to horiz. advection times thickness
                0043 C     StendAdvr  :: Salt  tendency due to vert. advection times thickness
                0044 C     StendTot   :: Salt  Total tendency
                0045 C     Hc         :: Layer thickness at the tracer point (m)
                0046 C     PIw        :: 1 if layer exists, 0 otherwise (at tracer point)
                0047 C     myThid     :: my Thread Id number
                0048       INTEGER iLa, myThid
                0049       _RL tracer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,     nSx,nSy)
                0050       _RL TtendSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0051       _RL TtendDiffh(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0052       _RL TtendDiffr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0053       _RL TtendAdvh (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0054       _RL TtendAdvr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0055       _RL TtendTot  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0056       _RL StendSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0057       _RL StendDiffh(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0058       _RL StendDiffr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0059       _RL StendAdvh (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0060       _RL StendAdvr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0061       _RL StendTot  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0062       _RL Hc        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
                0063       _RL PIc       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
                0064 CEOP
                0065
                0066 #ifdef LAYERS_THERMODYNAMICS
                0067
                0068 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0069 C !LOCAL VARIABLES:
                0070 C     bi, bj  :: tile indices
                0071 C     i, j    :: horizontal indices
                0072 C     k       :: vertical index for model grid
                0073 C     kp1     :: vertical index for model grid next cell
                0074 C     kci     :: index from CellIndex
                0075 C     kg      :: index for looping though layers_bounds
                0076 C     kk      :: vertical index for ZZ (fine) grid
                0077 C     kloc    :: local copy of kgu/v to reduce accesses to index arrays
                0078 C     mSteps  :: maximum number of steps for bisection method
                0079 C     TatC    :: temperature at C point
                0080       _RL Hcw       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0081       INTEGER bi, bj
                0082       INTEGER i,j,k,kk,kg,kci,kloc
                0083       INTEGER mSteps
                0084       INTEGER kgc(sNx+1,sNy+1)
                0085       INTEGER kgcw(sNx+1,sNy+1)
                0086       _RL TatC(sNx+1,sNy+1), dzfac, Tfac, Sfac
                0087       LOGICAL errorFlag
                0088       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0089 #ifdef LAYERS_FINEGRID_DIAPYCNAL
                0090       INTEGER kp1
                0091 #endif
                0092
                0093 C --  constants for T and S forcing, gets reset later for rho
                0094       Tfac = 1. _d 0
                0095       Sfac = 1. _d 0
                0096
                0097 C     compute maximum number of steps for bisection method (approx.
                0098 C     log2(Nlayers)) as log2(Nlayers) + 1 for safety
                0099       mSteps = INT(LOG10(DBLE(Nlayers))/LOG10(2. _d 0))+1
                0100
                0101 C      STOP 'DEBUG END: S/R LAYERS_DIAPYCNAL'
                0102
                0103 C --- The tile loops
                0104       DO bj=myByLo(myThid),myByHi(myThid)
                0105       DO bi=myBxLo(myThid),myBxHi(myThid)
                0106
                0107 C     Initialize the search indices
                0108       DO j = 1,sNy+1
                0109         DO i = 1,sNx+1
                0110 C       The temperature index (layer_G) goes from cold to warm.
                0111 C       The water column goes from warm (k=1) to cold (k=Nr).
                0112 C       So initialize the search with the warmest value.
                0113           kgc(i,j) = Nlayers
                0114           kgcw(i,j) = Nlayers-1
                0115         ENDDO
                0116       ENDDO
                0117
                0118 C     Reset the arrays
                0119 C --- These are at the w point
                0120       DO kg=1,Nlayers-1
                0121        DO j = 1-OLy,sNy+OLy
                0122         DO i = 1-OLx,sNx+OLx
                0123          TtendSurf (i,j,kg,bi,bj) = 0. _d 0
                0124          TtendDiffh(i,j,kg,bi,bj) = 0. _d 0
                0125          TtendDiffr(i,j,kg,bi,bj) = 0. _d 0
                0126          TtendAdvh(i,j,kg,bi,bj)  = 0. _d 0
                0127          TtendAdvr(i,j,kg,bi,bj)  = 0. _d 0
                0128          TtendTot(i,j,kg,bi,bj)   = 0. _d 0
                0129          StendSurf (i,j,kg,bi,bj) = 0. _d 0
                0130          StendDiffh(i,j,kg,bi,bj) = 0. _d 0
                0131          StendDiffr(i,j,kg,bi,bj) = 0. _d 0
                0132          StendAdvh(i,j,kg,bi,bj)  = 0. _d 0
                0133          StendAdvr(i,j,kg,bi,bj)  = 0. _d 0
                0134          StendTot(i,j,kg,bi,bj)   = 0. _d 0
                0135          Hcw(i,j,kg,bi,bj) = 0. _d 0
                0136         ENDDO
                0137        ENDDO
                0138       ENDDO
                0139 C --- These are at the c point
                0140       DO kg=1,Nlayers
                0141        DO j = 1-OLy,sNy+OLy
                0142         DO i = 1-OLx,sNx+OLx
                0143          Hc(i,j,kg,bi,bj) = 0. _d 0
                0144          PIc(i,j,kg,bi,bj) = 0. _d 0
                0145         ENDDO
                0146        ENDDO
                0147       ENDDO
                0148
                0149 #ifdef LAYERS_FINEGRID_DIAPYCNAL
                0150       DO kk=1,NZZ
                0151        k = MapIndex(kk)
                0152        kci = CellIndex(kk)
                0153        DO j = 1,sNy+1
                0154         DO i = 1,sNx+1
                0155 C ------ Find theta at the V point (south) on the fine Z grid
                0156          kp1=k+1
                0157          IF (maskC(i,j,kp1,bi,bj).EQ.zeroRS) kp1=k
                0158          TatC(i,j) = MapFact(kk) * tracer(i,j,k,bi,bj) +
                0159      &    (1. _d 0 -MapFact(kk)) * tracer(i,j,kp1,bi,bj)
                0160         ENDDO
                0161        ENDDO
                0162 #else
                0163       DO kk=1,Nr
                0164        k = kk
                0165        kci = kk
                0166        DO j = 1,sNy+1
                0167         DO i = 1,sNx+1
                0168          TatC(i,j) = tracer(i,j,k,bi,bj)
                0169         ENDDO
                0170        ENDDO
                0171 #endif /* LAYERS_FINEGRID_DIAPYCNAL */
                0172
                0173 C ------ debugging stuff
                0174 c         IF (i.EQ.38 .AND. j.EQ.4 .AND. bi.EQ.1 .AND. bj.EQ.1) THEN
                0175 c           i=38
                0176 c           j=4
                0177 c           WRITE(msgBuf,
                0178 c     &       '(A,I3,A,I3,A,I3,A,F7.2,A,F7.2,A,F7.2,A,F7.2,A,F3.1)')
                0179 c     &          'LAYERS_DEBUG: iLa=', iLa,
                0180 c     &          ', kk=', kk,
                0181 c     &          ', k=', k,
                0182 c     &          ', tracer=', tracer(i,j,k,bi,bj),
                0183 c     &          ', TatC=',TatC(i,j),
                0184 c     &          ', hFacC=',hFacC(i,j,k,bi,bj)
                0185 c           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0186 c     &                         SQUEEZE_RIGHT, myThid )
                0187 c         ENDIF
                0188 C ------ Now that we know T everywhere, determine the binning.
                0189 C        find the layer indices kgc for the center point
                0190        CALL LAYERS_LOCATE(
                0191      I      layers_bounds(1,iLa),Nlayers,mSteps,sNx,sNy,TatC,
                0192      O      kgc,
                0193      I      myThid )
                0194 #ifndef TARGET_NEC_SX
                0195 C     check for failures
                0196        IF ( debugLevel .GE. debLevC ) THEN
                0197         errorFlag = .FALSE.
                0198         DO j = 1,sNy+1
                0199          DO i = 1,sNx+1
                0200           IF ( kgc(i,j) .LE. 0 ) THEN
                0201            WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)')
                0202      &          'S/R LAYERS_LOCATE: Could not find a bin in ',
                0203      &          'layers_bounds for TatC(',i,',',j,',)=',TatC(i,j)
                0204            CALL PRINT_ERROR( msgBuf, myThid )
                0205            errorFlag = .TRUE.
                0206           ENDIF
                0207          ENDDO
                0208         ENDDO
                0209         IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_DIAPYCNAL'
                0210        ENDIF
                0211 #endif /* ndef TARGET_NEC_SX */
                0212
                0213 C        find the layer indices kgcw for the w point
                0214        CALL LAYERS_LOCATE(
                0215      I      layers_bounds_w(1,iLa),Nlayers-1,mSteps,sNx,sNy,TatC,
                0216      O      kgcw,
                0217      I      myThid )
                0218 #ifndef TARGET_NEC_SX
                0219 C     check for failures
                0220        IF ( debugLevel .GE. debLevC ) THEN
                0221         errorFlag = .FALSE.
                0222         DO j = 1,sNy+1
                0223          DO i = 1,sNx+1
                0224           IF ( kgcw(i,j) .LE. 0 ) THEN
                0225            WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)')
                0226      &          'S/R LAYERS_LOCATE: Could not find a bin in ',
                0227      &          'layers_bounds for TatC(',i,',',j,',)=',TatC(i,j)
                0228            CALL PRINT_ERROR( msgBuf, myThid )
                0229            errorFlag = .TRUE.
                0230           ENDIF
                0231          ENDDO
                0232         ENDDO
                0233         IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_DIAPYCNAL'
                0234        ENDIF
                0235 #endif /* ndef TARGET_NEC_SX */
                0236
                0237 C ------ Augment the bin values
                0238        DO j = 1,sNy+1
                0239         DO i = 1,sNx+1
                0240 #ifdef LAYERS_FINEGRID_DIAPYCNAL
                0241          dzfac = dZZf(kk) * hFacC(i,j,kci,bi,bj)
                0242 #else
                0243          dzfac = dRf(kci) * hFacC(i,j,kci,bi,bj)
                0244 #endif /* LAYERS_FINEGRID_DIAPYCNAL */
                0245          kloc = kgcw(i,j)
                0246
                0247 C ------- Thickness at w point
                0248          Hcw(i,j,kloc,bi,bj) = Hcw(i,j,kloc,bi,bj)
                0249      &     + dzfac
                0250 C ------- Thickness at c point
                0251          Hc(i,j,kgc(i,j),bi,bj) = Hc(i,j,kgc(i,j),bi,bj)
                0252      &     + dzfac
                0253
                0254 C ------- Now rescale dzfac to include the layer coordinate spacing
                0255          dzfac = dzfac * layers_recip_delta(kloc,iLa)
                0256
                0257 #ifdef LAYERS_PRHO_REF
                0258          IF ( layers_num(iLa) .EQ. 3 ) THEN
                0259            Tfac = layers_alpha(i,j,kci,bi,bj)
                0260            Sfac = layers_beta(i,j,kci,bi,bj)
                0261          ENDIF
                0262 #endif
                0263          IF (kci.EQ.1) THEN
                0264 C ------- We are in the surface layer
                0265           TtendSurf(i,j,kloc,bi,bj) = TtendSurf(i,j,kloc,bi,bj)
                0266      &         + Tfac * dzfac * layers_surfflux(i,j,1,1,bi,bj)
                0267           StendSurf(i,j,kloc,bi,bj) = StendSurf(i,j,kloc,bi,bj)
                0268      &         + Sfac * dzfac * layers_surfflux(i,j,1,2,bi,bj)
                0269          ENDIF
                0270
                0271 #ifdef SHORTWAVE_HEATING
                0272          TtendSurf(i,j,kloc,bi,bj) = TtendSurf(i,j,kloc,bi,bj)
                0273      &        + Tfac * dzfac * layers_sw(i,j,kci,1,bi,bj)
                0274 #endif /* SHORTWAVE_HEATING */
                0275
                0276 C ------- Diffusion
                0277          TtendDiffh(i,j,kloc,bi,bj) = TtendDiffh(i,j,kloc,bi,bj)
                0278      &        + dzfac * Tfac*( layers_dfx(i,j,kci,1,bi,bj)
                0279      &                       + layers_dfy(i,j,kci,1,bi,bj) )
                0280          StendDiffh(i,j,kloc,bi,bj) = StendDiffh(i,j,kloc,bi,bj)
                0281      &        + dzfac * Sfac*( layers_dfx(i,j,kci,2,bi,bj)
                0282      &                       + layers_dfy(i,j,kci,2,bi,bj) )
                0283          TtendDiffr(i,j,kloc,bi,bj) = TtendDiffr(i,j,kloc,bi,bj)
                0284      &        + dzfac * Tfac * layers_dfr(i,j,kci,1,bi,bj)
                0285          StendDiffr(i,j,kloc,bi,bj) = StendDiffr(i,j,kloc,bi,bj)
                0286      &        + dzfac * Sfac * layers_dfr(i,j,kci,2,bi,bj)
                0287 C ------- Advection
                0288          TtendAdvh(i,j,kloc,bi,bj) = TtendAdvh(i,j,kloc,bi,bj)
                0289      &        + dzfac * Tfac*( layers_afx(i,j,kci,1,bi,bj)
                0290      &                       + layers_afy(i,j,kci,1,bi,bj) )
                0291          StendAdvh(i,j,kloc,bi,bj) = StendAdvh(i,j,kloc,bi,bj)
                0292      &        + dzfac * Sfac*( layers_afx(i,j,kci,2,bi,bj)
                0293      &                       + layers_afy(i,j,kci,2,bi,bj) )
                0294          TtendAdvr(i,j,kloc,bi,bj) = TtendAdvr(i,j,kloc,bi,bj)
                0295      &        + dzfac * Tfac * layers_afr(i,j,kci,1,bi,bj)
                0296          StendAdvr(i,j,kloc,bi,bj) = StendAdvr(i,j,kloc,bi,bj)
                0297      &        + dzfac * Sfac * layers_afr(i,j,kci,2,bi,bj)
                0298 C -------- Total Tendency
                0299          TtendTot(i,j,kloc,bi,bj) = TtendTot(i,j,kloc,bi,bj)
                0300      &        + dzfac * Tfac * layers_tottend(i,j,kci,1,bi,bj)
                0301          StendTot(i,j,kloc,bi,bj) = StendTot(i,j,kloc,bi,bj)
                0302      &        + dzfac * Sfac * layers_tottend(i,j,kci,2,bi,bj)
                0303         ENDDO
                0304        ENDDO
                0305 C--   end kk loop
                0306       ENDDO
                0307
                0308 C--   Now that we know the thicknesses, compute the heaviside function
                0309 C--   (Needs another loop through Ng)
                0310       DO kg=1,Nlayers
                0311        DO j = 1,sNy+1
                0312         DO i = 1,sNx+1
                0313          IF (Hc(i,j,kg,bi,bj) .GT. 0.) THEN
                0314           PIc(i,j,kg,bi,bj) = 1. _d 0
                0315          ENDIF
                0316         ENDDO
                0317        ENDDO
                0318       ENDDO
                0319
                0320 C --- End bi,bj loop
                0321       ENDDO
                0322       ENDDO
                0323
                0324 #endif /* LAYERS_THERMODYNAMICS */
                0325
                0326       RETURN
                0327       END