pkg/bling/bling_mixedlayer.F

c0d1c06c15 Matt*0001 #include "BLING_OPTIONS.h"
fe1862e69b Mart*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
c0d1c06c15 Matt*0005
                0006 CBOP
4e4ad91a39 Jean*0007       SUBROUTINE BLING_MIXEDLAYER(
e0f9a7ba0b Matt*0008      U           sumMLDepth,
                0009      I           bi, bj, imin, imax, jmin, jmax,
                0010      I           myTime, myIter, myThid)
                0011
c0d1c06c15 Matt*0012 C     =================================================================
                0013 C     | subroutine bling_mixedlayer
e0f9a7ba0b Matt*0014 C     | o Calculate mixed layer depth based on density criterion
                0015 C     |   (default: second derivative criterion; optional: threshold)
c0d1c06c15 Matt*0016 C     =================================================================
                0017
e0f9a7ba0b Matt*0018       IMPLICIT NONE
                0019
c0d1c06c15 Matt*0020 C     === Global variables ===
                0021
                0022 #include "SIZE.h"
                0023 #include "EEPARAMS.h"
                0024 #include "PARAMS.h"
                0025 #include "FFIELDS.h"
                0026 #include "GRID.h"
                0027 #include "DYNVARS.h"
                0028 #include "BLING_VARS.h"
                0029 #include "PTRACERS_SIZE.h"
                0030 #include "PTRACERS_PARAMS.h"
fe1862e69b Mart*0031 #ifdef ALLOW_AUTODIFF_TAMC
                0032 # include "tamc.h"
                0033 #endif
c0d1c06c15 Matt*0034 C     === Routine arguments ===
                0035 C     bi,bj         :: tile indices
                0036 C     iMin,iMax     :: computation domain: 1rst index range
                0037 C     jMin,jMax     :: computation domain: 2nd  index range
                0038 C     myTime        :: current time
                0039 C     myIter        :: current timestep
                0040 C     myThid        :: thread Id. number
                0041       INTEGER bi, bj, imin, imax, jmin, jmax
                0042       INTEGER myThid
                0043       INTEGER myIter
                0044       _RL     myTime
                0045 C     === Output ===
                0046 C      sumMLDepth   :: mixed layer depth
                0047       _RL sumMLDepth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0048
                0049 C     === Local variables ===
e0f9a7ba0b Matt*0050 #ifdef BLING_USE_THRESHOLD_MLD
c0d1c06c15 Matt*0051       _RL dens_surf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0052       _RL dens_z    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0053       _RL delta_dens(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
e0f9a7ba0b Matt*0054 #else
                0055       _RL blg_stra   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0056       _RL rhoKm1     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0057       _RL rhoKp1     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0058       _RL blg_minstra(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0059       _RL blg_str2   (Nr)
                0060       _RL blg_str2max
fe1862e69b Mart*0061       INTEGER blgI,blgJ
                0062 #endif
                0063       INTEGER i,j,k
                0064 #ifdef ALLOW_AUTODIFF_TAMC
                0065       INTEGER ijkkey, ikey, kkey
e0f9a7ba0b Matt*0066 #endif
c0d1c06c15 Matt*0067 CEOP
                0068
fe1862e69b Mart*0069 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0070       ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
fe1862e69b Mart*0071       ikey = (ikey-1)*Nr
                0072 #endif /* ALLOW_AUTODIFF_TAMC */
e0f9a7ba0b Matt*0073 # ifndef BLING_USE_THRESHOLD_MLD
                0074 c ---------------------------------------------------------------------
                0075 c  Mixed layer depth initialization
                0076
                0077       DO j=jmin,jmax
                0078         DO i=imin,imax
                0079           sumMLDepth(i,j) = drF(1)
                0080           rhoKm1 (i,j)   = 0. _d 0
                0081           rhoKp1 (i,j)   = 0. _d 0
                0082           blg_minstra(i,j) = 0. _d 0
                0083         ENDDO
                0084       ENDDO
                0085
                0086       DO k=1,Nr
                0087         DO j=jmin,jmax
                0088           DO i=imin,imax
                0089             blg_stra(i,j,k) = 0. _d 0
                0090           ENDDO
                0091         ENDDO
                0092         blg_str2(k) = 0
                0093       ENDDO
                0094
                0095 c  get drhdr
                0096       DO k=1,Nr-1
                0097        DO j=jmin,jmax
                0098         DO i=imin,imax
                0099           rhoKm1(i,j) = rhoKp1(i,j)
                0100         ENDDO
                0101        ENDDO
                0102        CALL FIND_RHO_2D(
                0103      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
                0104      I        theta(1-OLx,1-OLy,k+1,bi,bj), salt(1-OLx,1-OLy,k+1,bi,bj),
                0105      O        rhoKp1,
                0106      I        k+1, bi, bj, myThid )
                0107
                0108        IF (k.EQ.1) THEN
                0109         DO j=jmin,jmax
                0110          DO i=imin,imax
                0111           blg_stra(i,j,k)= 0. _d 0
                0112          ENDDO
                0113         ENDDO
                0114        ELSE
                0115         DO j=jmin,jmax
                0116          DO i=imin,imax
                0117           blg_stra(i,j,k)= maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
                0118      &                *recip_drC(k)*rkSign
                0119      &                *(rhoKp1(i,j)-rhoKm1(i,j))
                0120          ENDDO
                0121         ENDDO
                0122        ENDIF
                0123       ENDDO
                0124 CMM Strat at k is on bottom of cell
                0125
                0126 C  % 1. find min(drho/dz)
fe1862e69b Mart*0127       DO k = 1,Nr
                0128 #ifdef ALLOW_AUTODIFF_TAMC
                0129        kkey = ikey + k
                0130 CADJ STORE blg_minstra = comlev1_bibj_k, key = kkey
                0131 #endif
                0132        DO j=jmin,jmax
                0133         DO i=imin,imax
                0134          IF(blg_stra(i,j,k).LT.blg_minstra(i,j))
e0f9a7ba0b Matt*0135      &        blg_minstra(i,j) = blg_stra(i,j,k)
                0136         ENDDO
fe1862e69b Mart*0137        ENDDO
                0138       ENDDO
e0f9a7ba0b Matt*0139
                0140 CMM NOW LOOP TO GET MLD
                0141         DO j=jmin,jmax
                0142          DO i=imin,imax
                0143 C  % 2. check that we have at least 3 grid cells of water
                0144           IF (hFacC(i,j,3,bi,bj).GT. 0. _d 0)  THEN
                0145 C  %  check that we have stable strat
                0146            IF (blg_minstra(i,j).LT. 0. _d 0) THEN
                0147 C  % 3. find the index of minimum stratification
                0148             blgI = 3
fe1862e69b Mart*0149             DO k = Nr,1,-1
e0f9a7ba0b Matt*0150              IF (blg_stra(i,j,k).EQ.blg_minstra(i,j))  blgI = k
                0151             ENDDO
                0152 C        %if deep enough see strat starts before pynocline
                0153             IF (blgI.GT.3) THEN
                0154 C % 4. compute the second derivative (d2rho/dz2)..
                0155 CMM Only look up
                0156             blg_str2max = -999
                0157             DO k = 3,blgI
fe1862e69b Mart*0158 #ifdef ALLOW_AUTODIFF_TAMC
                0159              ijkkey = ikey * (2*OLx+sNx) * (2*OLy+sNy)
                0160      &            + ( (j-1)*(2*OLx+sNx) + (i-1) ) + k
                0161 CADJ STORE blg_str2max = comlev1_bibj_ijk, key = ijkkey
                0162 #endif
e0f9a7ba0b Matt*0163 CMM blg_str2 will be in cell center
                0164              blg_str2(k) = (blg_stra(i,j,k-1) -blg_stra(i,j,k))
                0165      &                                             *recip_drF(k)
                0166 C     &                                             *blg_invdz(k)
                0167 Cblg_invdz should be DRF
                0168              IF (blg_str2(k).GT.blg_str2max)
                0169      &            blg_str2max = blg_str2(k)
                0170             ENDDO
fe1862e69b Mart*0171             DO k = Nr,1,-1
e0f9a7ba0b Matt*0172              IF (blg_str2(k).EQ.blg_str2max)  blgJ = k
                0173             ENDDO
                0174
                0175 CMM %TAKE WHERE STRAT IS FASTEST INCREASING ABOVE PYCNOCLINE
                0176             IF (blgJ.LT.blgI)  blgI = blgJ
                0177
                0178             IF ( (blg_str2(blgI).LT. 0. _d 0) .AND. (blgI.EQ.blgJ) )
                0179 CMM; %error: strat should increase or be flat
                0180      &       print*,'error: strat should increase or be flat'
                0181
                0182             sumMLDepth(i,j) = -RF(blgI)
                0183            ENDIF
                0184           ENDIF
                0185          ENDIF
                0186         ENDDO
                0187        ENDDO
                0188
                0189 #else /* BLING_USE_THRESHOLD_MLD */
c0d1c06c15 Matt*0190 c ---------------------------------------------------------------------
e0f9a7ba0b Matt*0191 c  Mixed layer depth
c0d1c06c15 Matt*0192
                0193       DO j=jmin,jmax
                0194         DO i=imin,imax
                0195           SumMLDepth(i,j) = drf(1)
                0196         ENDDO
                0197       ENDDO
                0198
                0199 c  Surface density
                0200       CALL FIND_RHO_2D(
                0201      I     1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
                0202      I     theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
                0203      O     dens_surf,
                0204      I     1, bi, bj, myThid )
                0205
                0206       DO k=1,Nr
                0207         DO j=jmin,jmax
                0208           DO i=imin,imax
                0209              if (k.eq.1) then
                0210               delta_dens(i,j,1) = 0. _d 0
                0211              else
                0212               delta_dens(i,j,k) = 9999. _d 0
                0213              endif
                0214           ENDDO
                0215         ENDDO
                0216       ENDDO
                0217
                0218       DO k = 2,Nr
                0219
e0f9a7ba0b Matt*0220 c  Potential density
c0d1c06c15 Matt*0221          CALL FIND_RHO_2D(
                0222      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
                0223      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
                0224      O        dens_z,
                0225      I        k, bi, bj, myThid )
                0226
                0227         DO j=jmin,jmax
                0228           DO i=imin,imax
e0f9a7ba0b Matt*0229
c0d1c06c15 Matt*0230 c           SumMLDepth(i,j) = 0. _d 0
                0231
                0232            IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
                0233             delta_dens(i,j,k) = dens_z(i,j)-dens_surf(i,j)
                0234             IF (delta_dens(i,j,k) .LT. 0.03 _d 0) THEN
                0235              SumMLDepth(i,j) = SumMLDepth(i,j)+drF(k)
e0f9a7ba0b Matt*0236             ENDIF
c0d1c06c15 Matt*0237            ENDIF
e0f9a7ba0b Matt*0238
c0d1c06c15 Matt*0239           ENDDO
                0240         ENDDO
                0241       ENDDO
                0242
e0f9a7ba0b Matt*0243 #endif /* BLING_USE_THRESHOLD_MLD */
                0244
c0d1c06c15 Matt*0245       RETURN
                0246       END