C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/kpp/kpp_calc.F,v 1.3 2000/09/13 17:07:11 heimbach Exp $ #include "KPP_OPTIONS.h" subroutine KPP_CALC( I bi, bj, myTime, myThid ) C /==========================================================\ C | SUBROUTINE KPP_CALC | C | o Compute all KPP fields defined in KPP.h | C |==========================================================| C | This subroutine serves as an interface between MITGCMUV | C | code and NCOM 1-D routines in kpp_routines.F | C \==========================================================/ IMPLICIT NONE c======================================================================= c c written by : jan morzel, august 11, 1994 c modified by : jan morzel, january 25, 1995 : "dVsq" and 1d code c detlef stammer, august, 1997 : for MIT GCM Classic c d. menemenlis, july, 1998 : for MIT GCM UV c c compute vertical mixing coefficients based on the k-profile c and oceanic planetary boundary layer scheme by large & mcwilliams. c c summary: c - compute interior mixing everywhere: c interior mixing gets computed at all interfaces due to constant c internal wave background activity ("fkpm" and "fkph"), which c is enhanced in places of static instability (local richardson c number < 0). c Additionally, mixing can be enhanced by adding contribution due c to shear instability which is a function of the local richardson c number c - double diffusivity: c interior mixing can be enhanced by double diffusion due to salt c fingering and diffusive convection (ifdef "kmixdd"). c - kpp scheme in the boundary layer: c c a.boundary layer depth: c at every gridpoint the depth of the oceanic boundary layer c ("hbl") gets computed by evaluating bulk richardson numbers. c b.boundary layer mixing: c within the boundary layer, above hbl, vertical mixing is c determined by turbulent surface fluxes, and interior mixing at c the lower boundary, i.e. at hbl. c c this subroutine provides the interface between the MIT GCM UV and the c subroutine "kppmix", where boundary layer depth, vertical c viscosity, vertical diffusivity, and counter gradient term (ghat) c are computed slabwise. c note: subroutine "kppmix" uses m-k-s units. c c time level: c input tracer and velocity profiles are evaluated at time level c tau, surface fluxes come from tau or tau-1. c c grid option: c in this "1-grid" implementation, diffusivity and viscosity c profiles are computed on the "t-grid" (by using velocity shear c profiles averaged from the "u,v-grid" onto the "t-grid"; note, that c the averaging includes zero values on coastal and seafloor grid c points). viscosity on the "u,v-grid" is computed by averaging the c "t-grid" viscosity values onto the "u,v-grid". c c vertical grid: c mixing coefficients get evaluated at the bottom of the lowest c layer, i.e., at depth zw(Nr). these values are only useful when c the model ocean domain does not include the entire ocean down to c the seafloor ("upperocean" setup) and allows flux through the c bottom of the domain. for full-depth runs, these mixing c coefficients are being zeroed out before leaving this subroutine. c c------------------------------------------------------------------------- c global parameters updated by kpp_calc c KPPviscAz - KPP eddy viscosity coefficient (m^2/s) c KPPdiffKzT - KPP diffusion coefficient for temperature (m^2/s) c KPPdiffKzS - KPP diffusion coefficient for salt and tracers (m^2/s) c KPPghat - Nonlocal transport coefficient (s/m^2) c KPPhbl - Boundary layer depth on "t-grid" (m) c KPPfrac - Fraction of short-wave flux penetrating mixing layer c-- KPP_CALC computes vertical viscosity and diffusivity for region c (-2:sNx+3,-2:sNy+3) as required by CALC_DIFFUSIVITY and requires c values of uVel, vVel, SurfaceTendencyU, SurfaceTendencyV in the c region (-2:sNx+4,-2:sNy+4). c Hence overlap region needs to be set OLx=4, OLy=4. c When option FRUGAL_KPP is used, computation in overlap regions c is replaced with exchange calls hence reducing overlap requirements c to OLx=1, OLy=1. #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "DYNVARS.h" #include "KPP.h" #include "KPP_PARAMS.h" #include "FFIELDS.h" #include "GRID.h" #ifdef ALLOW_AUTODIFF_TAMC #include "tamc.h" #include "tamc_keys.h" INTEGER isbyte PARAMETER( isbyte = 4 ) #else /* ALLOW_AUTODIFF_TAMC */ integer ikey #endif /* ALLOW_AUTODIFF_TAMC */ EXTERNAL DIFFERENT_MULTIPLE LOGICAL DIFFERENT_MULTIPLE c Routine arguments c bi, bj - array indices on which to apply calculations c myTime - Current time in simulation INTEGER bi, bj INTEGER myThid _RL myTime #ifdef ALLOW_KPP c Local arrays and variables c work? (nx,ny) - horizontal working arrays c ustar (nx,ny) - surface friction velocity (m/s) c bo (nx,ny) - surface turbulent buoyancy forcing (m^2/s^3) c bosol (nx,ny) - surface radiative buoyancy forcing (m^2/s^3) c shsq (nx,ny,Nr) - local velocity shear squared c at interfaces for ri_iwmix (m^2/s^2) c dVsq (nx,ny,Nr) - velocity shear re surface squared c at grid levels for bldepth (m^2/s^2) c dbloc (nx,ny,Nr) - local delta buoyancy at interfaces c for ri_iwmix and bldepth (m/s^2) c Ritop (nx,ny,Nr) - numerator of bulk richardson number c at grid levels for bldepth c vddiff (nx,ny,Nrp2,1)- vertical viscosity on "t-grid" (m^2/s) c vddiff (nx,ny,Nrp2,2)- vert. diff. on next row for temperature (m^2/s) c vddiff (nx,ny,Nrp2,3)- vert. diff. on next row for salt&tracers (m^2/s) c ghat (nx,ny,Nr) - nonlocal transport coefficient (s/m^2) c hbl (nx,ny) - mixing layer depth (m) c kmtj (nx,ny) - maximum number of wet levels in each column c z0 (nx,ny) - Roughness length (m) c zRef (nx,ny) - Reference depth: Hmix * epsilon (m) c uRef (nx,ny) - Reference zonal velocity (m/s) c vRef (nx,ny) - Reference meridional velocity (m/s) _RS worka (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS workb (1-OLx:sNx+OLx,1-OLy:sNy+OLy) #ifdef FRUGAL_KPP integer work1(sNx,sNy) _RS work2 (sNx,sNy) _RS ustar (sNx,sNy) _RS bo (sNx,sNy) _RS bosol (sNx,sNy) _RS shsq (sNx,sNy,Nr) _RS dVsq (sNx,sNy,Nr) _RS dbloc (sNx,sNy,Nr) _RS Ritop (sNx,sNy,Nr) _RS vddiff (sNx,sNy,0:Nrp1,mdiff) _RS ghat (sNx,sNy,Nr) _RS hbl (sNx,sNy) #ifdef KPP_ESTIMATE_UREF _RS z0 (sNx,sNy) _RS zRef (sNx,sNy) _RS uRef (sNx,sNy) _RS vRef (sNx,sNy) #endif /* KPP_ESTIMATE_UREF */ #else /* FRUGAL_KPP */ integer work1(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS work2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS ustar (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS bo (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS bosol (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS shsq (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RS dVsq (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RS dbloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RS Ritop (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RS vddiff (1-OLx:sNx+OLx,1-OLy:sNy+OLy,0:Nrp1,mdiff) _RS ghat (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RS hbl (1-OLx:sNx+OLx,1-OLy:sNy+OLy) #ifdef KPP_ESTIMATE_UREF _RS z0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS zRef (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS uRef (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS vRef (1-OLx:sNx+OLx,1-OLy:sNy+OLy) #endif /* KPP_ESTIMATE_UREF */ #endif /* FRUGAL_KPP */ c imin,imax,jmin,jmax - array indices integer imin , imax , jmin , jmax parameter( imin=-2, imax=sNx+3, jmin=-2, jmax=sNy+3 ) c mixing process switches logical lri parameter( lri = .true. ) _RS m1 parameter( m1=-1.0) _RS p0 , p5 , p25 , p125 , p0625 parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 ) _RL tempVar integer i, j, k, kp1, im1, ip1, jm1, jp1 #ifdef KPP_ESTIMATE_UREF _RS dBdz1, dBdz2, ustarX, ustarY #endif c Check to see if new vertical mixing coefficient should be computed now? IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,myTime-deltaTClock) .OR. 1 myTime .EQ. startTime ) THEN c----------------------------------------------------------------------- c prepare input arrays for subroutine "kppmix" to compute c viscosity and diffusivity and ghat. c All input arrays need to be in m-k-s units. c c note: for the computation of the bulk richardson number in the c "bldepth" subroutine, gradients of velocity and buoyancy are c required at every depth. in the case of very fine vertical grids c (thickness of top layer < 2m), the surface reference depth must c be set to zref=epsilon/2*zgrid(k), and the reference value c of velocity and buoyancy must be computed as vertical average c between the surface and 2*zref. in the case of coarse vertical c grids zref is zgrid(1)/2., and the surface reference value is c simply the surface value at zgrid(1). c----------------------------------------------------------------------- c------------------------------------------------------------------------ c density related quantities c -------------------------- c c work2 - density of surface layer (kg/m^3) c dbloc - local buoyancy gradient at Nr interfaces c g/rho{k+1,k+1} * [ drho{k,k+1}-drho{k+1,k+1} ] (m/s^2) c dbsfc (stored in Ritop to conserve stack memory) c - buoyancy difference with respect to the surface c g * [ drho{1,k}/rho{1,k} - drho{k,k}/rho{k,k} ] (m/s^2) c ttalpha (stored in vddiff(:,:,:,1) to conserve stack memory) c - thermal expansion coefficient without 1/rho factor c d(rho{k,k})/d(T(k)) (kg/m^3/C) c ssbeta (stored in vddiff(:,:,:,2) to conserve stack memory) c - salt expansion coefficient without 1/rho factor c d(rho{k,k})/d(S(k)) (kg/m^3/PSU) c------------------------------------------------------------------------ CALL TIMER_START('STATEKPP [KPP_CALC]', myThid) CALL STATEKPP( I bi, bj, myThid O , work2, dbloc, Ritop #ifdef FRUGAL_KPP O , vddiff(1 ,1 ,1,1), vddiff(1 ,1 ,1,2) #else O , vddiff(1-OLx,1-OLy,1,1), vddiff(1-OLx,1-OLy,1,2) #endif & ) CALL TIMER_STOP ('STATEKPP [KPP_CALC]', myThid) #ifdef KPP_SMOOTH_DBLOC c horizontally smooth dbloc with a 121 filter c (stored in ghat to save space) DO k = 1, Nr CALL SMOOTH_HORIZ_RS ( I k, bi, bj, I dbloc(1-OLx,1-OLy,k), O ghat (1-OLx,1-OLy,k) ) ENDDO #else /* KPP_SMOOTH_DBLOC */ DO k = 1, Nr #ifdef FRUGAL_KPP DO j = 1, sNy DO i = 1, sNx #else DO j = 1-OLy, sNy+OLy DO i = imin, imax #endif ghat(i,j,k) = dbloc(i,j,k) ENDDO ENDDO ENDDO #endif /* KPP_SMOOTH_DBLOC */ #ifdef KPP_SMOOTH_DENS c horizontally smooth density related quantities with 121 filters CALL SMOOTH_HORIZ_RS ( I k, bi, bj, I work2, O work2 ) DO k = 1, Nr CALL SMOOTH_HORIZ_RS ( I k, bi, bj, I dbloc (1-OLx,1-OLy,k) , O dbloc (1-OLx,1-OLy,k) ) CALL SMOOTH_HORIZ_RS ( I k, bi, bj, I Ritop (1-OLx,1-OLy,k) , O Ritop (1-OLx,1-OLy,k) ) CALL SMOOTH_HORIZ_RS ( I k, bi, bj, I vddiff(1-OLx,1-OLy,k,1), O vddiff(1-OLx,1-OLy,k,1) ) CALL SMOOTH_HORIZ_RS ( I k, bi, bj, I vddiff(1-OLx,1-OLy,k,2), O vddiff(1-OLx,1-OLy,k,2) ) ENDDO #endif /* KPP_SMOOTH_DENS */ DO k = 1, Nr #ifdef FRUGAL_KPP DO j = 1, sNy DO i = 1, sNx #else DO j = 1-OLy, sNy+OLy DO i = 1-OLx, sNx+OLx #endif c zero out dbloc over land points (so that the convective c part of the interior mixing can be diagnosed) dbloc(i,j,k) = dbloc(i,j,k) * pMask(i,j,k,bi,bj) ghat(i,j,k) = ghat(i,j,k) * pMask(i,j,k,bi,bj) Ritop(i,j,k) = Ritop(i,j,k) * pMask(i,j,k,bi,bj) if(k.eq.nzmax(i,j,bi,bj)) then dbloc(i,j,k) = p0 ghat(i,j,k) = p0 Ritop(i,j,k) = p0 endif c numerator of bulk richardson number on grid levels c note: land and ocean bottom values need to be set to zero c so that the subroutine "bldepth" works correctly Ritop(i,j,k) = (zgrid(1)-zgrid(k)) * Ritop(i,j,k) END DO END DO END DO c------------------------------------------------------------------------ c friction velocity, turbulent and radiative surface buoyancy forcing c ------------------------------------------------------------------- c taux / rho = SurfaceTendencyU * delZ(1) (N/m^2) c tauy / rho = SurfaceTendencyV * delZ(1) (N/m^2) c ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho ) (m/s) c bo = - g * ( alpha*SurfaceTendencyT + c beta *SurfaceTendencyS ) * delZ(1) / rho (m^2/s^3) c bosol = - g * alpha * Qsw * delZ(1) / rho (m^2/s^3) c------------------------------------------------------------------------ #ifdef FRUGAL_KPP DO j = 1, sNy jp1 = j + 1 DO i = 1, sNx #else DO j = jmin, jmax jp1 = j + 1 DO i = imin, imax #endif ip1 = i+1 tempVar = & (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) * & (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) + & (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) * & (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) if ( tempVar .lt. (epsln*epsln) ) then ustar(i,j) = SQRT( epsln * p5 * delZ(1) ) else ustar(i,j) = SQRT( SQRT( tempVar ) * p5 * delZ(1) ) endif bo(I,J) = - gravity * & ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) + & vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj) & ) * & delZ(1) / work2(I,J) bosol(I,J) = - gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) * & delZ(1) / work2(I,J) END DO END DO #ifndef FRUGAL_KPP c set array edges to zero DO j = jmin, jmax DO i = 1-OLx, imin-1 ustar(i,j) = p0 bo (I,J) = p0 bosol(I,J) = p0 END DO DO i = imax+1, sNx+OLx ustar(i,j) = p0 bo (I,J) = p0 bosol(I,J) = p0 END DO END DO DO i = 1-OLx, sNx+OLx DO j = 1-OLy, jmin-1 ustar(i,j) = p0 bo (I,J) = p0 bosol(I,J) = p0 END DO DO j = jmax+1, sNy+OLy ustar(i,j) = p0 bo (I,J) = p0 bosol(I,J) = p0 END DO END DO #endif c------------------------------------------------------------------------ c velocity shear c -------------- c Get velocity shear squared, averaged from "u,v-grid" c onto "t-grid" (in (m/s)**2): c dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2 at grid levels c shsq(k)=(U(k)-U(k+1))**2+(V(k)-V(k+1))**2 at interfaces c------------------------------------------------------------------------ c dVsq computation #ifdef KPP_ESTIMATE_UREF c Get rid of vertical resolution dependence of dVsq term by c estimating a surface velocity that is independent of first level c thickness in the model. First determine mixed layer depth hMix. c Second zRef = espilon * hMix. Third determine roughness length c scale z0. Third estimate reference velocity. #ifdef FRUGAL_KPP DO j = 1, sNy jp1 = j + 1 DO i = 1, sNx #else DO j = jmin, jmax jp1 = j + 1 DO i = imin, imax #endif /* FRUGAL_KPP */ ip1 = i + 1 c Determine mixed layer depth hMix as the shallowest depth at which c dB/dz exceeds 5.2e-5 s^-2. work1(i,j) = nzmax(i,j,bi,bj) DO k = 1, Nr IF ( k .LT. nzmax(i,j,bi,bj) .AND. & dbloc(i,j,k) / drC(k+1) .GT. dB_dz ) & work1(i,j) = k END DO c Linearly interpolate to find hMix. k = work1(i,j) IF ( k .EQ. 0 .OR. nzmax(i,j,bi,bj) .EQ. 1 ) THEN zRef(i,j) = p0 ELSEIF ( k .EQ. 1) THEN dBdz2 = dbloc(i,j,1) / drC(2) zRef(i,j) = drF(1) * dB_dz / dBdz2 ELSEIF ( k .LT. nzmax(i,j,bi,bj) ) THEN dBdz1 = dbloc(i,j,k-1) / drC(k ) dBdz2 = dbloc(i,j,k ) / drC(k+1) zRef(i,j) = rF(k) + drF(k) * (dB_dz - dBdz1) / & MAX ( phepsi, dBdz2 - dBdz1 ) ELSE zRef(i,j) = rF(k+1) ENDIF c Compute roughness length scale z0 subject to 0 < z0 tempVar = p5 * ( & (uVel(i, j, 1,bi,bj)-uVel(i, j, 2,bi,bj)) * & (uVel(i, j, 1,bi,bj)-uVel(i, j, 2,bi,bj)) + & (uVel(ip1,j, 1,bi,bj)-uVel(ip1,j, 2,bi,bj)) * & (uVel(ip1,j, 1,bi,bj)-uVel(ip1,j, 2,bi,bj)) + & (vVel(i, j, 1,bi,bj)-vVel(i, j, 2,bi,bj)) * & (vVel(i, j, 1,bi,bj)-vVel(i, j, 2,bi,bj)) + & (vVel(i, jp1,1,bi,bj)-vVel(i, jp1,2,bi,bj)) * & (vVel(i, jp1,1,bi,bj)-vVel(i, jp1,2,bi,bj)) ) if ( tempVar .lt. (epsln*epsln) ) then tempVar = epsln else tempVar = SQRT ( tempVar ) endif z0(i,j) = rF(2) * & ( rF(3) * LOG ( rF(3) / rF(2) ) / & ( rF(3) - rF(2) ) - & tempVar * vonK / & MAX ( ustar(i,j), phepsi ) ) z0(i,j) = MAX ( z0(i,j), phepsi ) c zRef is set to 0.1 * hMix subject to z0 <= zRef <= drF(1) zRef(i,j) = MAX ( epsilon * zRef(i,j), z0(i,j) ) zRef(i,j) = MIN ( zRef(i,j), drF(1) ) c Estimate reference velocity uRef and vRef. uRef(i,j) = p5 * & ( uVel(i,j,1,bi,bj) + uVel(ip1,j,1,bi,bj) ) vRef(i,j) = p5 * & ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) ) IF ( zRef(i,j) .LT. drF(1) ) THEN ustarX = ( SurfaceTendencyU(i, j,bi,bj) + & SurfaceTendencyU(ip1,j,bi,bj) ) * p5 ustarY = ( SurfaceTendencyV(i,j, bi,bj) + & SurfaceTendencyU(i,jp1,bi,bj) ) * p5 tempVar = ustarX * ustarX + ustarY * ustarY if ( tempVar .lt. (epsln*epsln) ) then tempVar = epsln else tempVar = SQRT ( tempVar ) endif tempVar = ustar(i,j) * & ( LOG ( zRef(i,j) / rF(2) ) + & z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) / & vonK / tempVar uRef(i,j) = uRef(i,j) + ustarX * tempVar vRef(i,j) = vRef(i,j) + ustarY * tempVar ENDIF END DO END DO IF (KPPmixingMaps) THEN #ifdef FRUGAL_KPP CALL PRINT_MAPRS( I zRef, 'zRef', PRINT_MAP_XY, I 1, sNx, 1, sNy, 1, 1, 1, 1, I 1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) CALL PRINT_MAPRS( I z0, 'z0', PRINT_MAP_XY, I 1, sNx, 1, sNy, 1, 1, 1, 1, I 1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) CALL PRINT_MAPRS( I uRef, 'uRef', PRINT_MAP_XY, I 1, sNx, 1, sNy, 1, 1, 1, 1, I 1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) CALL PRINT_MAPRS( I vRef, 'vRef', PRINT_MAP_XY, I 1, sNx, 1, sNy, 1, 1, 1, 1, I 1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) #else CALL PRINT_MAPRS( I zRef, 'zRef', PRINT_MAP_XY, I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1, I 1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) CALL PRINT_MAPRS( I z0, 'z0', PRINT_MAP_XY, I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1, I 1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) CALL PRINT_MAPRS( I uRef, 'uRef', PRINT_MAP_XY, I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1, I 1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) CALL PRINT_MAPRS( I vRef, 'vRef', PRINT_MAP_XY, I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1, I 1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) #endif ENDIF DO k = 1, Nr #ifdef FRUGAL_KPP DO j = 1, sNy jm1 = j - 1 jp1 = j + 1 DO i = 1, sNx #else DO j = jmin, jmax jm1 = j - 1 jp1 = j + 1 DO i = imin, imax #endif /* FRUGAL_KPP */ im1 = i - 1 ip1 = i + 1 dVsq(i,j,k) = p5 * ( $ (uRef(i,j) - uVel(i, j, k,bi,bj)) * $ (uRef(i,j) - uVel(i, j, k,bi,bj)) + $ (uRef(i,j) - uVel(ip1,j, k,bi,bj)) * $ (uRef(i,j) - uVel(ip1,j, k,bi,bj)) + $ (vRef(i,j) - vVel(i, j, k,bi,bj)) * $ (vRef(i,j) - vVel(i, j, k,bi,bj)) + $ (vRef(i,j) - vVel(i, jp1,k,bi,bj)) * $ (vRef(i,j) - vVel(i, jp1,k,bi,bj)) ) #ifdef KPP_SMOOTH_DVSQ dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * ( $ (uRef(i,j) - uVel(i, jm1,k,bi,bj)) * $ (uRef(i,j) - uVel(i, jm1,k,bi,bj)) + $ (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) * $ (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) + $ (uRef(i,j) - uVel(i, jp1,k,bi,bj)) * $ (uRef(i,j) - uVel(i, jp1,k,bi,bj)) + $ (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) * $ (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) + $ (vRef(i,j) - vVel(im1,j, k,bi,bj)) * $ (vRef(i,j) - vVel(im1,j, k,bi,bj)) + $ (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) * $ (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) + $ (vRef(i,j) - vVel(ip1,j, k,bi,bj)) * $ (vRef(i,j) - vVel(ip1,j, k,bi,bj)) + $ (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) * $ (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) ) #endif /* KPP_SMOOTH_DVSQ */ END DO END DO END DO #else /* KPP_ESTIMATE_UREF */ DO k = 1, Nr #ifdef FRUGAL_KPP DO j = 1, sNy jm1 = j - 1 jp1 = j + 1 DO i = 1, sNx #else DO j = jmin, jmax jm1 = j - 1 jp1 = j + 1 DO i = imin, imax #endif /* FRUGAL_KPP */ im1 = i - 1 ip1 = i + 1 dVsq(i,j,k) = p5 * ( $ (uVel(i, j, 1,bi,bj)-uVel(i, j, k,bi,bj)) * $ (uVel(i, j, 1,bi,bj)-uVel(i, j, k,bi,bj)) + $ (uVel(ip1,j, 1,bi,bj)-uVel(ip1,j, k,bi,bj)) * $ (uVel(ip1,j, 1,bi,bj)-uVel(ip1,j, k,bi,bj)) + $ (vVel(i, j, 1,bi,bj)-vVel(i, j, k,bi,bj)) * $ (vVel(i, j, 1,bi,bj)-vVel(i, j, k,bi,bj)) + $ (vVel(i, jp1,1,bi,bj)-vVel(i, jp1,k,bi,bj)) * $ (vVel(i, jp1,1,bi,bj)-vVel(i, jp1,k,bi,bj)) ) #ifdef KPP_SMOOTH_DVSQ dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * ( $ (uVel(i, jm1,1,bi,bj)-uVel(i, jm1,k,bi,bj)) * $ (uVel(i, jm1,1,bi,bj)-uVel(i, jm1,k,bi,bj)) + $ (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) * $ (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) + $ (uVel(i, jp1,1,bi,bj)-uVel(i, jp1,k,bi,bj)) * $ (uVel(i, jp1,1,bi,bj)-uVel(i, jp1,k,bi,bj)) + $ (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) * $ (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) + $ (vVel(im1,j, 1,bi,bj)-vVel(im1,j, k,bi,bj)) * $ (vVel(im1,j, 1,bi,bj)-vVel(im1,j, k,bi,bj)) + $ (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) * $ (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) + $ (vVel(ip1,j, 1,bi,bj)-vVel(ip1,j, k,bi,bj)) * $ (vVel(ip1,j, 1,bi,bj)-vVel(ip1,j, k,bi,bj)) + $ (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) * $ (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) ) #endif /* KPP_SMOOTH_DVSQ */ END DO END DO END DO #endif /* KPP_ESTIMATE_UREF */ c shsq computation DO k = 1, Nrm1 kp1 = k + 1 #ifdef FRUGAL_KPP DO j = 1, sNy jm1 = j - 1 jp1 = j + 1 DO i = 1, sNx #else DO j = jmin, jmax jm1 = j - 1 jp1 = j + 1 DO i = imin, imax #endif /* FRUGAL_KPP */ im1 = i - 1 ip1 = i + 1 shsq(i,j,k) = p5 * ( $ (uVel(i, j, k,bi,bj)-uVel(i, j, kp1,bi,bj)) * $ (uVel(i, j, k,bi,bj)-uVel(i, j, kp1,bi,bj)) + $ (uVel(ip1,j, k,bi,bj)-uVel(ip1,j, kp1,bi,bj)) * $ (uVel(ip1,j, k,bi,bj)-uVel(ip1,j, kp1,bi,bj)) + $ (vVel(i, j, k,bi,bj)-vVel(i, j, kp1,bi,bj)) * $ (vVel(i, j, k,bi,bj)-vVel(i, j, kp1,bi,bj)) + $ (vVel(i, jp1,k,bi,bj)-vVel(i, jp1,kp1,bi,bj)) * $ (vVel(i, jp1,k,bi,bj)-vVel(i, jp1,kp1,bi,bj)) ) #ifdef KPP_SMOOTH_SHSQ shsq(i,j,k) = p5 * shsq(i,j,k) + p125 * ( $ (uVel(i, jm1,k,bi,bj)-uVel(i, jm1,kp1,bi,bj)) * $ (uVel(i, jm1,k,bi,bj)-uVel(i, jm1,kp1,bi,bj)) + $ (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) * $ (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) + $ (uVel(i, jp1,k,bi,bj)-uVel(i, jp1,kp1,bi,bj)) * $ (uVel(i, jp1,k,bi,bj)-uVel(i, jp1,kp1,bi,bj)) + $ (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) * $ (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) + $ (vVel(im1,j, k,bi,bj)-vVel(im1,j, kp1,bi,bj)) * $ (vVel(im1,j, k,bi,bj)-vVel(im1,j, kp1,bi,bj)) + $ (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) * $ (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) + $ (vVel(ip1,j, k,bi,bj)-vVel(ip1,j, kp1,bi,bj)) * $ (vVel(ip1,j, k,bi,bj)-vVel(ip1,j, kp1,bi,bj)) + $ (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) * $ (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) ) #endif END DO END DO END DO c shsq @ Nr computation #ifdef FRUGAL_KPP DO j = 1, sNy DO i = 1, sNx #else DO j = jmin, jmax DO i = imin, imax #endif shsq(i,j,Nr) = p0 END DO END DO #ifndef FRUGAL_KPP c set array edges to zero DO k = 1, Nr DO j = jmin, jmax DO i = 1-OLx, imin-1 shsq(i,j,k) = p0 dVsq(i,j,k) = p0 END DO DO i = imax+1, sNx+OLx shsq(i,j,k) = p0 dVsq(i,j,k) = p0 END DO END DO DO i = 1-OLx, sNx+OLx DO j = 1-OLy, jmin-1 shsq(i,j,k) = p0 dVsq(i,j,k) = p0 END DO DO j = jmax+1, sNy+OLy shsq(i,j,k) = p0 dVsq(i,j,k) = p0 END DO END DO END DO #endif c----------------------------------------------------------------------- c solve for viscosity, diffusivity, ghat, and hbl on "t-grid" c----------------------------------------------------------------------- #ifdef FRUGAL_KPP DO j = 1, sNy DO i = 1, sNx #else DO j = 1-OLy, sNy+OLy DO i = 1-OLx, sNx+OLx #endif work1(i,j) = nzmax(i,j,bi,bj) work2(i,j) = Fcori(i,j,bi,bj) END DO END DO CALL TIMER_START('KPPMIX [KPP_CALC]', myThid) CALL KPPMIX ( I lri, work1, shsq, dVsq, ustar I , bo, bosol, dbloc, Ritop, work2 I , ikey O , vddiff U , ghat O , hbl & ) CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid) IF (KPPmixingMaps) THEN #ifdef FRUGAL_KPP CALL PRINT_MAPRS( I hbl, 'hbl', PRINT_MAP_XY, I 1, sNx, 1, sNy, 1, 1, 1, 1, I 1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) #else CALL PRINT_MAPRS( I hbl, 'hbl', PRINT_MAP_XY, I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1, I 1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ) #endif ENDIF #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE vddiff, ghat = comlev1_kpp, key = ikey #endif /* ALLOW_AUTODIFF_TAMC */ c----------------------------------------------------------------------- c zero out land values, c make sure coefficients are within reasonable bounds, c and transfer to global variables c----------------------------------------------------------------------- #ifdef FRUGAL_KPP DO j = 1, sNy DO i = 1, sNx #else DO j = jmin, jmax DO i = imin, imax #endif DO k = 1, Nr c KPPviscAz tempVar = min( maxKPPviscAz(k), vddiff(i,j,k-1,1) ) tempVar = max( minKPPviscAz, tempVar ) KPPviscAz(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj) c KPPdiffKzS tempVar = min( maxKPPdiffKzS, vddiff(i,j,k-1,2) ) tempVar = max( minKPPdiffKzS, tempVar ) KPPdiffKzS(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj) c KPPdiffKzT tempVar = min( maxKPPdiffKzT, vddiff(i,j,k-1,3) ) tempVar = max( minKPPdiffKzT, tempVar ) KPPdiffKzT(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj) c KPPghat tempVar = min( maxKPPghat, ghat(i,j,k) ) tempVar = max( minKPPghat, tempVar ) KPPghat(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj) END DO c KPPhbl: set to -zgrid(1) over land KPPhbl(i,j,bi,bj) = (hbl(i,j) + zgrid(1)) & * pMask(i,j,1,bi,bj) - & zgrid(1) END DO END DO #ifdef FRUGAL_KPP _EXCH_XYZ_R8(KPPviscAz , myThid ) _EXCH_XYZ_R8(KPPdiffKzS , myThid ) _EXCH_XYZ_R8(KPPdiffKzT , myThid ) _EXCH_XYZ_R8(KPPghat , myThid ) _EXCH_XY_R8 (KPPhbl , myThid ) #endif #ifdef KPP_SMOOTH_VISC c horizontal smoothing of vertical viscosity c as coded requires FRUGAL_KPP and OLx=4, OLy=4 c alternatively could recode with OLx=5, OLy=5 DO k = 1, Nr CALL SMOOTH_HORIZ_RL ( I k, bi, bj, I KPPviscAz(1-OLx,1-OLy,k,bi,bj), O KPPviscAz(1-OLx,1-OLy,k,bi,bj) ) END DO #endif /* KPP_SMOOTH_VISC */ #ifdef KPP_SMOOTH_DIFF c horizontal smoothing of vertical diffusivity c as coded requires FRUGAL_KPP and OLx=4, OLy=4 c alternatively could recode with OLx=5, OLy=5 DO k = 1, Nr CALL SMOOTH_HORIZ_RL ( I k, bi, bj, I KPPdiffKzS(1-OLx,1-OLy,k,bi,bj), O KPPdiffKzS(1-OLx,1-OLy,k,bi,bj) ) CALL SMOOTH_HORIZ_RL ( I k, bi, bj, I KPPdiffKzT(1-OLx,1-OLy,k,bi,bj), O KPPdiffKzT(1-OLx,1-OLy,k,bi,bj) ) END DO #endif /* KPP_SMOOTH_DIFF */ C Compute fraction of solar short-wave flux penetrating to C the bottom of the mixing layer. DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx worka(i,j) = KPPhbl(i,j,bi,bj) ENDDO ENDDO CALL SWFRAC( I (sNx+2*OLx)*(sNy+2*OLy), m1, worka, O workb ) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx KPPfrac(i,j,bi,bj) = workb(i,j) ENDDO ENDDO ENDIF #endif ALLOW_KPP RETURN END