C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/kpp/kpp_calc.F,v 1.18 2003/10/02 14:38:37 adcroft Exp $ C $Name: $ #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_SEAICE #include "SEAICE_EXTERNAL.h" #endif /* ALLOW_SEAICE */ #ifdef ALLOW_AUTODIFF_TAMC #include "tamc.h" #include "tamc_keys.h" #else /* ALLOW_AUTODIFF_TAMC */ integer ikppkey #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 constants c minusone, p0, p5, p25, p125, p0625 c imin, imax, jmin, jmax - array computation indices _RL minusone parameter( minusone=-1.0) _KPP_RL p0 , p5 , p25 , p125 , p0625 parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 ) integer imin ,imax ,jmin ,jmax #ifdef FRUGAL_KPP parameter(imin=1 ,imax=sNx ,jmin=1 ,jmax=sNy ) #else parameter(imin=2-OLx,imax=sNx+OLx-1,jmin=2-OLy,jmax=sNy+OLy-1) #endif 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) _RL worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy ) integer work1 ( ibot:itop , jbot:jtop ) _KPP_RL work2 ( ibot:itop , jbot:jtop ) _KPP_RL work3 ( ibot:itop , jbot:jtop ) _KPP_RL ustar ( ibot:itop , jbot:jtop ) _KPP_RL bo ( ibot:itop , jbot:jtop ) _KPP_RL bosol ( ibot:itop , jbot:jtop ) _KPP_RL shsq ( ibot:itop , jbot:jtop , Nr ) _KPP_RL dVsq ( ibot:itop , jbot:jtop , Nr ) _KPP_RL dbloc ( ibot:itop , jbot:jtop , Nr ) _KPP_RL Ritop ( ibot:itop , jbot:jtop , Nr ) _KPP_RL vddiff( ibot:itop , jbot:jtop , 0:Nrp1, mdiff ) _KPP_RL ghat ( ibot:itop , jbot:jtop , Nr ) _KPP_RL hbl ( ibot:itop , jbot:jtop ) #ifdef KPP_ESTIMATE_UREF _KPP_RL z0 ( ibot:itop , jbot:jtop ) _KPP_RL zRef ( ibot:itop , jbot:jtop ) _KPP_RL uRef ( ibot:itop , jbot:jtop ) _KPP_RL vRef ( ibot:itop , jbot:jtop ) #endif /* KPP_ESTIMATE_UREF */ _KPP_RL tempvar2 integer i, j, k, kp1, im1, ip1, jm1, jp1 #ifdef KPP_ESTIMATE_UREF _KPP_RL tempvar1, dBdz1, dBdz2, ustarX, ustarY #endif #ifdef ALLOW_AUTODIFF_TAMC act1 = bi - myBxLo(myThid) max1 = myBxHi(myThid) - myBxLo(myThid) + 1 act2 = bj - myByLo(myThid) max2 = myByHi(myThid) - myByLo(myThid) + 1 act3 = myThid - 1 max3 = nTx*nTy act4 = ikey_dynamics - 1 ikppkey = (act1 + 1) + act2*max1 & + act3*max1*max2 & + act4*max1*max2*max3 #endif /* ALLOW_AUTODIFF_TAMC */ 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 O , vddiff(ibot,jbot,1,1), vddiff(ibot,jbot,1,2) & ) CALL TIMER_STOP ('STATEKPP [KPP_CALC]', myThid) DO k = 1, Nr DO j = jbot, jtop DO i = ibot, itop ghat(i,j,k) = dbloc(i,j,k) ENDDO ENDDO ENDDO #ifdef KPP_SMOOTH_DBLOC c horizontally smooth dbloc with a 121 filter c smooth dbloc stored in ghat to save space c dbloc(k) is buoyancy gradientnote between k and k+1 c levels therefore k+1 mask must be used DO k = 1, Nr-1 CALL KPP_SMOOTH_HORIZ ( I k+1, bi, bj, U ghat (ibot,jbot,k) ) ENDDO #endif /* KPP_SMOOTH_DBLOC */ #ifdef KPP_SMOOTH_DENS c horizontally smooth density related quantities with 121 filters CALL KPP_SMOOTH_HORIZ ( I 1, bi, bj, U work2 ) DO k = 1, Nr CALL KPP_SMOOTH_HORIZ ( I k+1, bi, bj, U dbloc (ibot,jbot,k) ) CALL KPP_SMOOTH_HORIZ ( I k, bi, bj, U Ritop (ibot,jbot,k) ) CALL KPP_SMOOTH_HORIZ ( I k, bi, bj, U vddiff(ibot,jbot,k,1) ) CALL KPP_SMOOTH_HORIZ ( I k, bi, bj, U vddiff(ibot,jbot,k,2) ) ENDDO #endif /* KPP_SMOOTH_DENS */ DO k = 1, Nr DO j = jbot, jtop DO i = ibot, itop 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 cph( cph this avoids a single or double recomp./call of statekpp CADJ store work2 = comlev1_kpp, key = ikppkey #ifdef ALLOW_AUTODIFF_KPP_EXTENSIVE_STORE CADJ store dbloc, Ritop, ghat = comlev1_kpp, key = ikppkey CADJ store vddiff = comlev1_kpp, key = ikppkey #endif cph) c------------------------------------------------------------------------ c friction velocity, turbulent and radiative surface buoyancy forcing c ------------------------------------------------------------------- c taux / rho = SurfaceTendencyU * drF(1) (N/m^2) c tauy / rho = SurfaceTendencyV * drF(1) (N/m^2) c ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho ) (m/s) c bo = - g * ( alpha*SurfaceTendencyT + c beta *SurfaceTendencyS ) * drF(1) / rho (m^2/s^3) c bosol = - g * alpha * Qsw * drF(1) / rho (m^2/s^3) c------------------------------------------------------------------------ c initialize arrays to zero DO j = jbot, jtop DO i = ibot, itop ustar(i,j) = p0 bo (I,J) = p0 bosol(I,J) = p0 END DO END DO DO j = jmin, jmax jp1 = j + 1 DO i = imin, imax ip1 = i+1 work3(i,j) = & (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)) END DO END DO cph( CADJ store work3 = comlev1_kpp, key = ikppkey cph) DO j = jmin, jmax jp1 = j + 1 DO i = imin, imax ip1 = i+1 if ( work3(i,j) .lt. (phepsi*phepsi) ) then ustar(i,j) = SQRT( phepsi * p5 * drF(1) ) else tempVar2 = SQRT( work3(i,j) ) * p5 * drF(1) ustar(i,j) = SQRT( tempVar2 ) endif if ( .NOT. useSEAICE ) & bo(I,J) = - gravity * & ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) + & vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj) & ) * & drF(1) / work2(I,J) #ifdef ALLOW_SEAICE if ( useSEAICE ) & bo(I,J) = - gravity * & ( vddiff(I,J,1,1) * (SurfaceTendencyT(i,j,bi,bj)+ & SurfaceTendencyTice(i,j,bi,bj)) + & vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj) & ) * & drF(1) / work2(I,J) #endif /* ALLOW_SEAICE */ bosol(I,J) = gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) * & recip_Cp*recip_rhoConst*recip_dRf(1) * & drF(1) / work2(I,J) END DO END DO cph( CADJ store ustar = comlev1_kpp, key = ikppkey cph) 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 initialize arrays to zero DO k = 1, Nr DO j = jbot, jtop DO i = ibot, itop shsq(i,j,k) = p0 dVsq(i,j,k) = p0 END DO END DO END DO 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. DO j = jmin, jmax jp1 = j + 1 DO i = imin, imax 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 tempVar1 = 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 ( tempVar1 .lt. (epsln*epsln) ) then tempVar2 = epsln else tempVar2 = SQRT ( tempVar1 ) endif z0(i,j) = rF(2) * & ( rF(3) * LOG ( rF(3) / rF(2) ) / & ( rF(3) - rF(2) ) - & tempVar2 * 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 tempVar1 = ustarX * ustarX + ustarY * ustarY if ( tempVar1 .lt. (epsln*epsln) ) then tempVar2 = epsln else tempVar2 = SQRT ( tempVar1 ) endif tempVar2 = ustar(i,j) * & ( LOG ( zRef(i,j) / rF(2) ) + & z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) / & vonK / tempVar2 uRef(i,j) = uRef(i,j) + ustarX * tempVar2 vRef(i,j) = vRef(i,j) + ustarY * tempVar2 ENDIF END DO END DO DO k = 1, Nr DO j = jmin, jmax jm1 = j - 1 jp1 = j + 1 DO i = imin, imax 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 DO j = jmin, jmax jm1 = j - 1 jp1 = j + 1 DO i = imin, imax 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 DO j = jmin, jmax jm1 = j - 1 jp1 = j + 1 DO i = imin, imax 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 cph( #ifdef ALLOW_AUTODIFF_KPP_EXTENSIVE_STORE CADJ store dvsq, shsq = comlev1_kpp, key = ikppkey #endif cph) c----------------------------------------------------------------------- c solve for viscosity, diffusivity, ghat, and hbl on "t-grid" c----------------------------------------------------------------------- DO j = jbot, jtop DO i = ibot, itop 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 mytime, mythid I , work1, shsq, dVsq, ustar I , bo, bosol, dbloc, Ritop, work2 I , ikppkey O , vddiff U , ghat O , hbl ) CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid) c----------------------------------------------------------------------- c zero out land values and transfer to global variables c----------------------------------------------------------------------- DO j = jmin, jmax DO i = imin, imax DO k = 1, Nr KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * pMask(i,j,k,bi,bj) KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * pMask(i,j,k,bi,bj) KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * pMask(i,j,k,bi,bj) KPPghat(i,j,k,bi,bj) = ghat(i,j,k) * pMask(i,j,k,bi,bj) END DO KPPhbl(i,j,bi,bj) = hbl(i,j) * pMask(i,j,1,bi,bj) 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 DO k = 1, Nr CALL SMOOTH_HORIZ ( I k, bi, bj, U KPPviscAz(1-OLx,1-OLy,k,bi,bj) ) END DO _EXCH_XYZ_R8(KPPviscAz , myThid ) #endif /* KPP_SMOOTH_VISC */ #ifdef KPP_SMOOTH_DIFF c horizontal smoothing of vertical diffusivity DO k = 1, Nr CALL SMOOTH_HORIZ ( I k, bi, bj, U KPPdiffKzS(1-OLx,1-OLy,k,bi,bj) ) CALL SMOOTH_HORIZ ( I k, bi, bj, U KPPdiffKzT(1-OLx,1-OLy,k,bi,bj) ) END DO _EXCH_XYZ_R8(KPPdiffKzS , myThid ) _EXCH_XYZ_R8(KPPdiffKzT , myThid ) #endif /* KPP_SMOOTH_DIFF */ cph( cph crucial: this avoids full recomp./call of kppmix CADJ store KPPhbl = comlev1_kpp, key = ikppkey cph) 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), minusone, I mytime, mythid, U worka ) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx KPPfrac(i,j,bi,bj) = worka(i,j) ENDDO ENDDO ENDIF #endif /* ALLOW_KPP */ RETURN END subroutine KPP_CALC_DUMMY( I bi, bj, myTime, myThid ) C /==========================================================\ C | SUBROUTINE KPP_CALC_DUMMY | C | o Compute all KPP fields defined in KPP.h | C | o Dummy routine for TAMC C |==========================================================| C | This subroutine serves as an interface between MITGCMUV | C | code and NCOM 1-D routines in kpp_routines.F | C \==========================================================/ IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "KPP.h" #include "KPP_PARAMS.h" #include "GRID.h" 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 constants integer i, j, k DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx KPPhbl (i,j,bi,bj) = 1.0 KPPfrac(i,j,bi,bj) = 0.0 DO k = 1,Nr KPPghat (i,j,k,bi,bj) = 0.0 KPPviscAz (i,j,k,bi,bj) = viscAz KPPdiffKzT(i,j,k,bi,bj) = diffKzT KPPdiffKzS(i,j,k,bi,bj) = diffKzS ENDDO ENDDO ENDDO #endif RETURN END