--- MITgcm/pkg/kpp/kpp_routines.F 2003/12/16 20:58:57 1.18 +++ MITgcm/pkg/kpp/kpp_routines.F 2007/10/19 19:11:17 1.35 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/kpp/kpp_routines.F,v 1.18 2003/12/16 20:58:57 heimbach Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/kpp/kpp_routines.F,v 1.35 2007/10/19 19:11:17 jmc Exp $ C $Name: $ #include "KPP_OPTIONS.h" @@ -11,8 +11,7 @@ C-- o WSCALE - Compute turbulent velocity scales. C-- o RI_IWMIX - Compute interior viscosity diffusivity coefficients. C-- o Z121 - Apply 121 vertical smoothing. -C-- o KPP_SMOOTH_HORIZ - Apply horizontal smoothing to KPP array. -C-- o SMOOTH_HORIZ - Apply horizontal smoothing to global array. +C-- o SMOOTH_HORIZ- Apply horizontal smoothing to global array. C-- o BLMIX - Boundary layer mixing coefficients. C-- o ENHANCE - Enhance diffusivity at boundary layer interface. C-- o STATEKPP - Compute buoyancy-related input arrays. @@ -20,13 +19,18 @@ c*********************************************************************** SUBROUTINE KPPMIX ( - I mytime, mythid - I , kmtj, shsq, dvsq, ustar - I , bo, bosol, dbloc, Ritop, coriol + I kmtj, shsq, dvsq, ustar, msk + I , bo, bosol +#ifdef ALLOW_SALT_PLUME + I , boplume,SPDepth +#endif /* ALLOW_SALT_PLUME */ + I , dbloc, Ritop, coriol + I , diffusKzS, diffusKzT I , ikppkey O , diffus U , ghat - O , hbl ) + O , hbl + I , myTime, myIter, myThid ) c----------------------------------------------------------------------- c @@ -47,40 +51,51 @@ #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" -#include "DYNVARS.h" -#include "FFIELDS.h" #include "KPP_PARAMS.h" c input -c myTime - current time in simulation -c myThid - thread number for this instance of the routine -c kmtj (imt) - number of vertical layers on this row -c shsq (imt,Nr) - (local velocity shear)^2 ((m/s)^2) -c dvsq (imt,Nr) - (velocity shear re sfc)^2 ((m/s)^2) -c ustar (imt) - surface friction velocity (m/s) -c bo (imt) - surface turbulent buoy. forcing (m^2/s^3) -c bosol (imt) - radiative buoyancy forcing (m^2/s^3) -c dbloc (imt,Nr) - local delta buoyancy across interfaces (m/s^2) -c dblocSm(imt,Nr) - horizontally smoothed dbloc (m/s^2) -c stored in ghat to save space -c Ritop (imt,Nr) - numerator of bulk Richardson Number -c (zref-z) * delta buoyancy w.r.t. surface ((m/s)^2) -c coriol (imt) - Coriolis parameter (1/s) +c myTime :: Current time in simulation +c myIter :: Current iteration number in simulation +c myThid :: My Thread Id. number +c kmtj (imt) - number of vertical layers on this row +c msk (imt) - surface mask (=1 if water, =0 otherwise) +c shsq (imt,Nr) - (local velocity shear)^2 ((m/s)^2) +c dvsq (imt,Nr) - (velocity shear re sfc)^2 ((m/s)^2) +c ustar (imt) - surface friction velocity (m/s) +c bo (imt) - surface turbulent buoy. forcing (m^2/s^3) +c bosol (imt) - radiative buoyancy forcing (m^2/s^3) +c boplume(imt) - haline buoyancy forcing (m^2/s^3) +c dbloc (imt,Nr) - local delta buoyancy across interfaces (m/s^2) +c dblocSm(imt,Nr) - horizontally smoothed dbloc (m/s^2) +c stored in ghat to save space +c Ritop (imt,Nr) - numerator of bulk Richardson Number +c (zref-z) * delta buoyancy w.r.t. surface ((m/s)^2) +c coriol (imt) - Coriolis parameter (1/s) +c diffusKzS(imt,Nr)- background vertical diffusivity for scalars (m^2/s) +c diffusKzT(imt,Nr)- background vertical diffusivity for theta (m^2/s) c note: there is a conversion from 2-D to 1-D for input output variables, c e.g., hbl(sNx,sNy) -> hbl(imt), c where hbl(i,j) -> hbl((j-1)*sNx+i) - _RL mytime - integer mythid - integer kmtj (imt ) - _KPP_RL shsq (imt,Nr) - _KPP_RL dvsq (imt,Nr) - _KPP_RL ustar (imt ) - _KPP_RL bo (imt ) - _KPP_RL bosol (imt ) - _KPP_RL dbloc (imt,Nr) - _KPP_RL Ritop (imt,Nr) - _KPP_RL coriol(imt ) + _RL myTime + integer myIter + integer myThid + integer kmtj (imt ) + _RL shsq (imt,Nr) + _RL dvsq (imt,Nr) + _RL ustar (imt ) + _RL bo (imt ) + _RL bosol (imt ) +#ifdef ALLOW_SALT_PLUME + _RL boplume (imt ) + _RL SPDepth (imt ) +#endif /* ALLOW_SALT_PLUME */ + _RL dbloc (imt,Nr) + _RL Ritop (imt,Nr) + _RL coriol (imt ) + _RS msk (imt ) + _RL diffusKzS(imt,Nr) + _RL diffusKzT(imt,Nr) integer ikppkey @@ -91,9 +106,9 @@ c ghat (imt) - nonlocal transport coefficient (s/m^2) c hbl (imt) - mixing layer depth (m) - _KPP_RL diffus(imt,0:Nrp1,mdiff) - _KPP_RL ghat (imt,Nr) - _KPP_RL hbl (imt) + _RL diffus(imt,0:Nrp1,mdiff) + _RL ghat (imt,Nr) + _RL hbl (imt) #ifdef ALLOW_KPP @@ -107,14 +122,14 @@ c sigma (imt ) - normalized depth (d / hbl) c Rib (imt,Nr ) - bulk Richardson number - integer kbl (imt ) - _KPP_RL bfsfc (imt ) - _KPP_RL casea (imt ) - _KPP_RL stable(imt ) - _KPP_RL dkm1 (imt, mdiff) - _KPP_RL blmc (imt,Nr,mdiff) - _KPP_RL sigma (imt ) - _KPP_RL Rib (imt,Nr ) + integer kbl(imt ) + _RL bfsfc (imt ) + _RL casea (imt ) + _RL stable (imt ) + _RL dkm1 (imt, mdiff) + _RL blmc (imt,Nr,mdiff) + _RL sigma (imt ) + _RL Rib (imt,Nr ) integer i, k, md @@ -132,8 +147,9 @@ cph) call Ri_iwmix ( I kmtj, shsq, dbloc, ghat + I , diffusKzS, diffusKzT I , ikppkey - O , diffus ) + O , diffus, myThid ) cph( cph these storings avoid recomp. of Ri_iwmix @@ -150,9 +166,9 @@ c----------------------------------------------------------------------- do md = 1, mdiff + do k=1,Nrp1 do i = 1,imt - do k=kmtj(i),Nrp1 - diffus(i,k,md) = 0.0 + if(k.ge.kmtj(i)) diffus(i,k,md) = 0.0 end do end do end do @@ -164,12 +180,15 @@ c----------------------------------------------------------------------- call bldepth ( - I mytime, mythid - I , kmtj - I , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol + I kmtj + I , dvsq, dbloc, Ritop, ustar, bo, bosol +#ifdef ALLOW_SALT_PLUME + I , boplume,SPDepth +#endif /* ALLOW_SALT_PLUME */ + I , coriol I , ikppkey O , hbl, bfsfc, stable, casea, kbl, Rib, sigma - & ) + I , myTime, myIter, myThid ) CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey @@ -180,7 +199,7 @@ call blmix ( I ustar, bfsfc, hbl, stable, casea, diffus, kbl O , dkm1, blmc, ghat, sigma, ikppkey - & ) + I , myThid ) cph( CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey @@ -193,7 +212,8 @@ call enhance ( I dkm1, hbl, kbl, diffus, casea U , ghat - O , blmc ) + O , blmc + I , myThid ) cph( cph avoids recomp. of enhance @@ -209,9 +229,19 @@ do k = 1, Nr do i = 1, imt if (k .lt. kbl(i)) then +#ifdef ALLOW_SHELFICE +C when there is shelfice on top (msk(i)=0), reset the boundary layer +C mixing coefficients blmc to pure Ri-number based mixing + blmc(i,k,1) = max ( blmc(i,k,1)*msk(i), + & diffus(i,k,1) ) + blmc(i,k,2) = max ( blmc(i,k,2)*msk(i), + & diffus(i,k,2) ) + blmc(i,k,3) = max ( blmc(i,k,3)*msk(i), + & diffus(i,k,3) ) +#endif /* not ALLOW_SHELFICE */ diffus(i,k,1) = max ( blmc(i,k,1), viscAr ) - diffus(i,k,2) = max ( blmc(i,k,2), diffKrS ) - diffus(i,k,3) = max ( blmc(i,k,3), diffKrT ) + diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) ) + diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) ) else ghat(i,k) = 0. endif @@ -226,12 +256,15 @@ c************************************************************************* subroutine bldepth ( - I mytime, mythid - I , kmtj - I , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol + I kmtj + I , dvsq, dbloc, Ritop, ustar, bo, bosol +#ifdef ALLOW_SALT_PLUME + I , boplume,SPDepth +#endif /* ALLOW_SALT_PLUME */ + I , coriol I , ikppkey O , hbl, bfsfc, stable, casea, kbl, Rib, sigma - & ) + I , myTime, myIter, myThid ) c the oceanic planetary boundary layer depth, hbl, is determined as c the shallowest depth where the bulk Richardson number is @@ -257,15 +290,13 @@ IMPLICIT NONE #include "SIZE.h" -#include "EEPARAMS.h" -#include "PARAMS.h" #include "KPP_PARAMS.h" -#include "FFIELDS.h" c input c------ -c myTime : current time in simulation -c myThid : thread number for this instance of the routine +c myTime :: Current time in simulation +c myIter :: Current iteration number in simulation +c myThid :: My Thread Id. number c kmtj : number of vertical layers c dvsq : (velocity shear re sfc)^2 ((m/s)^2) c dbloc : local delta buoyancy across interfaces (m/s^2) @@ -275,18 +306,24 @@ c ustar : surface friction velocity (m/s) c bo : surface turbulent buoyancy forcing (m^2/s^3) c bosol : radiative buoyancy forcing (m^2/s^3) +c boplume : haline buoyancy forcing (m^2/s^3) c coriol : Coriolis parameter (1/s) - _RL mytime - integer mythid + _RL myTime + integer myIter + integer myThid integer kmtj(imt) - _KPP_RL dvsq (imt,Nr) - _KPP_RL dbloc (imt,Nr) - _KPP_RL Ritop (imt,Nr) - _KPP_RL ustar (imt) - _KPP_RL bo (imt) - _KPP_RL bosol (imt) - _KPP_RL coriol(imt) - integer ikppkey + _RL dvsq (imt,Nr) + _RL dbloc (imt,Nr) + _RL Ritop (imt,Nr) + _RL ustar (imt) + _RL bo (imt) + _RL bosol (imt) + _RL coriol (imt) + integer ikppkey, kkppkey +#ifdef ALLOW_SALT_PLUME + _RL boplume (imt) + _RL SPDepth (imt) +#endif /* ALLOW_SALT_PLUME */ c output c-------- @@ -297,26 +334,25 @@ c kbl : -1 of first grid level below hbl c Rib : Bulk Richardson number c sigma : normalized depth (d/hbl) - _KPP_RL hbl (imt) - _KPP_RL bfsfc (imt) - _KPP_RL stable(imt) - _KPP_RL casea (imt) - integer kbl (imt) - _KPP_RL Rib (imt,Nr) - _KPP_RL sigma (imt) + _RL hbl (imt) + _RL bfsfc (imt) + _RL stable (imt) + _RL casea (imt) + integer kbl(imt) + _RL Rib (imt,Nr) + _RL sigma (imt) #ifdef ALLOW_KPP c local c------- c wm, ws : turbulent velocity scales (m/s) - _KPP_RL wm(imt), ws(imt) - _RL worka(imt) - - _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2 + _RL wm(imt), ws(imt) + _RL worka(imt) + _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2 integer i, kl - _KPP_RL p5 , eins + _RL p5 , eins parameter ( p5=0.5, eins=1.0 ) _RL minusone parameter ( minusone=-1.0 ) @@ -340,15 +376,21 @@ do kl = 2, Nr +#ifdef ALLOW_AUTODIFF_TAMC + kkppkey = (ikppkey-1)*Nr + kl +#endif + c compute bfsfc = sw fraction at hbf * zgrid do i = 1, imt worka(i) = zgrid(kl) end do +CADJ store worka = comlev1_kpp_k, key = kkppkey call SWFRAC( - I imt, hbf, - I mytime, mythid, - U worka ) + I imt, hbf, + U worka, + I myTime, myIter, myThid ) +CADJ store worka = comlev1_kpp_k, key = kkppkey do i = 1, imt @@ -359,10 +401,26 @@ c compute bfsfc= Bo + radiative contribution down to hbf * hbl bfsfc(i) = bo(i) + bosol(i)*(1. - worka(i)) - stable(i) = p5 + sign(p5,bfsfc(i)) - sigma(i) = stable(i) + (1. - stable(i)) * epsilon end do +#ifdef ALLOW_SALT_PLUME +c compute bfsfc = plume fraction at hbf * zgrid + do i = 1, imt + worka(i) = zgrid(kl) + enddo + call PLUMEFRAC( + I imt, hbf,SPDepth, + U worka, + I myTime, myIter, myThid) + do i = 1, imt + bfsfc(i) = bfsfc(i) + boplume(i)*(1. - worka(i)) + enddo +#endif /* ALLOW_SALT_PLUME */ + + do i = 1, imt + stable(i) = p5 + sign(p5,bfsfc(i)) + sigma(i) = stable(i) + (1. - stable(i)) * epsilon + enddo c----------------------------------------------------------------------- c compute velocity scales at sigma, for hbl= caseA = -zgrid(kl) @@ -370,7 +428,8 @@ call wscale ( I sigma, casea, ustar, bfsfc, - O wm, ws ) + O wm, ws, myThid ) +CADJ store ws = comlev1_kpp_k, key = kkppkey do i = 1, imt @@ -443,14 +502,31 @@ do i = 1, imt worka(i) = hbl(i) end do +CADJ store worka = comlev1_kpp +CADJ & , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /) call SWFRAC( - I imt, minusone, - I mytime, mythid, - U worka ) + I imt, minusone, + U worka, + I myTime, myIter, myThid ) +CADJ store worka = comlev1_kpp +CADJ & , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /) do i = 1, imt bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i)) end do + +#ifdef ALLOW_SALT_PLUME + do i = 1, imt + worka(i) = hbl(i) + enddo + call PLUMEFRAC( + I imt,minusone,SPDepth, + U worka, + I myTime, myIter, myThid ) + do i = 1, imt + bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i)) + enddo +#endif /* ALLOW_SALT_PLUME */ CADJ store bfsfc = comlev1_kpp CADJ & , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /) @@ -511,14 +587,31 @@ do i = 1, imt worka(i) = hbl(i) end do +CADJ store worka = comlev1_kpp +CADJ & , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /) call SWFRAC( I imt, minusone, - I mytime, mythid, - U worka ) + U worka, + I myTime, myIter, myThid ) +CADJ store worka = comlev1_kpp +CADJ & , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /) do i = 1, imt bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i)) end do + +#ifdef ALLOW_SALT_PLUME + do i = 1, imt + worka(i) = hbl(i) + enddo + call PLUMEFRAC( + I imt,minusone,SPDepth, + U worka, + I myTime, myIter, myThid ) + do i = 1, imt + bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i)) + enddo +#endif /* ALLOW_SALT_PLUME */ CADJ store bfsfc = comlev1_kpp CADJ & , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /) @@ -546,7 +639,8 @@ subroutine wscale ( I sigma, hbl, ustar, bfsfc, - O wm, ws ) + O wm, ws, + I myThid ) c compute turbulent velocity scales. c use a 2D-lookup table for wm and ws as functions of ustar and @@ -567,26 +661,28 @@ c hbl : boundary layer depth (m) c ustar : surface friction velocity (m/s) c bfsfc : total surface buoyancy flux (m^2/s^3) - _KPP_RL sigma(imt) - _KPP_RL hbl (imt) - _KPP_RL ustar(imt) - _KPP_RL bfsfc(imt) +c myThid : thread number for this instance of the routine + integer myThid + _RL sigma(imt) + _RL hbl (imt) + _RL ustar(imt) + _RL bfsfc(imt) c output c-------- c wm, ws : turbulent velocity scales at sigma - _KPP_RL wm(imt), ws(imt) + _RL wm(imt), ws(imt) #ifdef ALLOW_KPP c local c------ c zehat : = zeta * ustar**3 - _KPP_RL zehat + _RL zehat integer iz, izp1, ju, i, jup1 - _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam - _KPP_RL wbm, was, wbs, u3, tempVar + _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam + _RL wbm, was, wbs, u3, tempVar c----------------------------------------------------------------------- c use lookup table for zehat < zmax only; otherwise use @@ -641,9 +737,11 @@ c************************************************************************* subroutine Ri_iwmix ( - I kmtj, shsq, dbloc, dblocSm - I , ikppkey - O , diffus ) + I kmtj, shsq, dbloc, dblocSm, + I diffusKzS, diffusKzT, + I ikppkey, + O diffus, + I myThid ) c compute interior viscosity diffusivity coefficients due c to shear instability (dependent on a local Richardson number), @@ -662,28 +760,34 @@ c shsq (imt,Nr) (local velocity shear)^2 ((m/s)^2) c dbloc (imt,Nr) local delta buoyancy (m/s^2) c dblocSm(imt,Nr) horizontally smoothed dbloc (m/s^2) - integer kmtj (imt) - _KPP_RL shsq (imt,Nr) - _KPP_RL dbloc (imt,Nr) - _KPP_RL dblocSm(imt,Nr) +c diffusKzS(imt,Nr)- background vertical diffusivity for scalars (m^2/s) +c diffusKzT(imt,Nr)- background vertical diffusivity for theta (m^2/s) +c myThid :: My Thread Id. number + integer kmtj (imt) + _RL shsq (imt,Nr) + _RL dbloc (imt,Nr) + _RL dblocSm (imt,Nr) + _RL diffusKzS(imt,Nr) + _RL diffusKzT(imt,Nr) integer ikppkey + integer myThid c output c diffus(imt,0:Nrp1,1) vertical viscosivity coefficient (m^2/s) c diffus(imt,0:Nrp1,2) vertical scalar diffusivity (m^2/s) c diffus(imt,0:Nrp1,3) vertical temperature diffusivity (m^2/s) - _KPP_RL diffus(imt,0:Nrp1,3) + _RL diffus(imt,0:Nrp1,3) #ifdef ALLOW_KPP c local variables c Rig local Richardson number c fRi, fcon function of Rig - _KPP_RL Rig - _KPP_RL fRi, fcon - _KPP_RL ratio + _RL Rig + _RL fRi, fcon + _RL ratio integer i, ki - _KPP_RL c1, c0 + _RL c1, c0 #ifdef ALLOW_KPP_VERTICALLY_SMOOTH integer mr @@ -702,7 +806,16 @@ #ifdef ALLOW_AUTODIFF_TAMC C break data flow dependence on diffus diffus(1,1,1) = 0.0 + + do ki = 1, Nr + do i = 1, imt + diffus(i,ki,1) = 0. + diffus(i,ki,2) = 0. + diffus(i,ki,3) = 0. + enddo + enddo #endif + do ki = 1, Nr do i = 1, imt @@ -719,6 +832,7 @@ endif end do end do +CADJ store diffus = comlev1_kpp, key = ikppkey c----------------------------------------------------------------------- c vertically smooth Ri @@ -729,7 +843,8 @@ CADJ & , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /) call z121 ( - U diffus(1,0,1)) + U diffus(1,0,1) + I myThid ) end do #endif @@ -758,13 +873,19 @@ c mixing due to internal waves, and shear and static instability #ifndef EXCLUDE_KPP_SHEAR_MIX - diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0 - diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0 - diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0 + if ( .NOT. inAdMode ) then + diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0 + diffus(i,ki,2) = diffusKzS(i,ki)+fcon*difscon+fRi*difs0 + diffus(i,ki,3) = diffusKzT(i,ki)+fcon*diftcon+fRi*dift0 + else + diffus(i,ki,1) = viscAr + diffus(i,ki,2) = diffusKzS(i,ki) + diffus(i,ki,3) = diffusKzT(i,ki) + endif #else diffus(i,ki,1) = viscAr - diffus(i,ki,2) = diffKrS - diffus(i,ki,3) = diffKrT + diffus(i,ki,2) = diffusKzS(i,ki) + diffus(i,ki,3) = diffusKzT(i,ki) #endif end do @@ -787,7 +908,8 @@ c************************************************************************* subroutine z121 ( - U v ) + U v, + I myThid ) c Apply 121 smoothing in k to 2-d array V(i,k=1,Nr) c top (0) value is used as a dummy @@ -805,16 +927,18 @@ c input/output c------------- c v : 2-D array to be smoothed in Nrp1 direction - _KPP_RL v(imt,0:Nrp1) +c myThid: thread number for this instance of the routine + integer myThid + _RL v(imt,0:Nrp1) #ifdef ALLOW_KPP c local - _KPP_RL zwork, zflag - _KPP_RL KRi_range(1:Nrp1) + _RL zwork, zflag + _RL KRi_range(1:Nrp1) integer i, k, km1, kp1 - _KPP_RL p0 , p25 , p5 , p2 + _RL p0 , p25 , p5 , p2 parameter ( p0 = 0.0, p25 = 0.25, p5 = 0.5, p2 = 2.0 ) KRi_range(Nrp1) = p0 @@ -868,99 +992,23 @@ c************************************************************************* - subroutine kpp_smooth_horiz ( - I k, bi, bj, - U fld ) - -c Apply horizontal smoothing to KPP array - - IMPLICIT NONE -#include "SIZE.h" -#include "KPP_PARAMS.h" - -c input -c bi, bj : array indices -c k : vertical index used for masking - integer k, bi, bj - -c input/output -c fld : 2-D array to be smoothed - _KPP_RL fld( ibot:itop, jbot:jtop ) - -#ifdef ALLOW_KPP - -c local - integer i, j, im1, ip1, jm1, jp1 - _KPP_RL tempVar - _KPP_RL fld_tmp( ibot:itop, jbot:jtop ) - - integer imin , imax , jmin , jmax - parameter( imin=ibot+1, imax=itop-1, jmin=jbot+1, jmax=jtop-1 ) - - _KPP_RL p0 , p5 , p25 , p125 , p0625 - parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 ) - - DO j = jmin, jmax - jm1 = j-1 - jp1 = j+1 - DO i = imin, imax - im1 = i-1 - ip1 = i+1 - tempVar = - & p25 * pMask(i ,j ,k,bi,bj) + - & p125 * ( pMask(im1,j ,k,bi,bj) + - & pMask(ip1,j ,k,bi,bj) + - & pMask(i ,jm1,k,bi,bj) + - & pMask(i ,jp1,k,bi,bj) ) + - & p0625 * ( pMask(im1,jm1,k,bi,bj) + - & pMask(im1,jp1,k,bi,bj) + - & pMask(ip1,jm1,k,bi,bj) + - & pMask(ip1,jp1,k,bi,bj) ) - IF ( tempVar .GE. p25 ) THEN - fld_tmp(i,j) = ( - & p25 * fld(i ,j )*pMask(i ,j ,k,bi,bj) + - & p125 *(fld(im1,j )*pMask(im1,j ,k,bi,bj) + - & fld(ip1,j )*pMask(ip1,j ,k,bi,bj) + - & fld(i ,jm1)*pMask(i ,jm1,k,bi,bj) + - & fld(i ,jp1)*pMask(i ,jp1,k,bi,bj))+ - & p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) + - & fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) + - & fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) + - & fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj))) - & / tempVar - ELSE - fld_tmp(i,j) = fld(i,j) - ENDIF - ENDDO - ENDDO - -c transfer smoothed field to output array - DO j = jmin, jmax - DO i = imin, imax - fld(i,j) = fld_tmp(i,j) - ENDDO - ENDDO - -#endif /* ALLOW_KPP */ - - return - end - -c************************************************************************* - subroutine smooth_horiz ( I k, bi, bj, - U fld ) + U fld, + I myThid ) c Apply horizontal smoothing to global _RL 2-D array IMPLICIT NONE #include "SIZE.h" +#include "GRID.h" #include "KPP_PARAMS.h" c input c bi, bj : array indices c k : vertical index used for masking +c myThid : thread number for this instance of the routine + INTEGER myThid integer k, bi, bj c input/output @@ -987,26 +1035,26 @@ im1 = i-1 ip1 = i+1 tempVar = - & p25 * pMask(i ,j ,k,bi,bj) + - & p125 * ( pMask(im1,j ,k,bi,bj) + - & pMask(ip1,j ,k,bi,bj) + - & pMask(i ,jm1,k,bi,bj) + - & pMask(i ,jp1,k,bi,bj) ) + - & p0625 * ( pMask(im1,jm1,k,bi,bj) + - & pMask(im1,jp1,k,bi,bj) + - & pMask(ip1,jm1,k,bi,bj) + - & pMask(ip1,jp1,k,bi,bj) ) + & p25 * maskC(i ,j ,k,bi,bj) + + & p125 * ( maskC(im1,j ,k,bi,bj) + + & maskC(ip1,j ,k,bi,bj) + + & maskC(i ,jm1,k,bi,bj) + + & maskC(i ,jp1,k,bi,bj) ) + + & p0625 * ( maskC(im1,jm1,k,bi,bj) + + & maskC(im1,jp1,k,bi,bj) + + & maskC(ip1,jm1,k,bi,bj) + + & maskC(ip1,jp1,k,bi,bj) ) IF ( tempVar .GE. p25 ) THEN fld_tmp(i,j) = ( - & p25 * fld(i ,j )*pMask(i ,j ,k,bi,bj) + - & p125 *(fld(im1,j )*pMask(im1,j ,k,bi,bj) + - & fld(ip1,j )*pMask(ip1,j ,k,bi,bj) + - & fld(i ,jm1)*pMask(i ,jm1,k,bi,bj) + - & fld(i ,jp1)*pMask(i ,jp1,k,bi,bj))+ - & p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) + - & fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) + - & fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) + - & fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj))) + & p25 * fld(i ,j )*maskC(i ,j ,k,bi,bj) + + & p125 *(fld(im1,j )*maskC(im1,j ,k,bi,bj) + + & fld(ip1,j )*maskC(ip1,j ,k,bi,bj) + + & fld(i ,jm1)*maskC(i ,jm1,k,bi,bj) + + & fld(i ,jp1)*maskC(i ,jp1,k,bi,bj))+ + & p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) + + & fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) + + & fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) + + & fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj))) & / tempVar ELSE fld_tmp(i,j) = fld(i,j) @@ -1031,7 +1079,7 @@ subroutine blmix ( I ustar, bfsfc, hbl, stable, casea, diffus, kbl O , dkm1, blmc, ghat, sigma, ikppkey - & ) + I , myThid ) c mixing coefficients within boundary layer depend on surface c forcing and the magnitude and gradient of interior mixing below @@ -1053,13 +1101,15 @@ c stable(imt) = 1 in stable forcing c casea (imt) = 1 in case A c diffus(imt,0:Nrp1,mdiff) vertical diffusivities (m^2/s) -c kbl(imt) -1 of first grid level below hbl - _KPP_RL ustar (imt) - _KPP_RL bfsfc (imt) - _KPP_RL hbl (imt) - _KPP_RL stable(imt) - _KPP_RL casea (imt) - _KPP_RL diffus(imt,0:Nrp1,mdiff) +c kbl (imt) -1 of first grid level below hbl +c myThid thread number for this instance of the routine + integer myThid + _RL ustar (imt) + _RL bfsfc (imt) + _RL hbl (imt) + _RL stable(imt) + _RL casea (imt) + _RL diffus(imt,0:Nrp1,mdiff) integer kbl(imt) c output @@ -1067,11 +1117,11 @@ c blmc (imt,Nr,mdiff) boundary layer mixing coefficients (m^2/s) c ghat (imt,Nr) nonlocal scalar transport c sigma(imt) normalized depth (d / hbl) - _KPP_RL dkm1 (imt,mdiff) - _KPP_RL blmc (imt,Nr,mdiff) - _KPP_RL ghat (imt,Nr) - _KPP_RL sigma(imt) - integer ikppkey + _RL dkm1 (imt,mdiff) + _RL blmc (imt,Nr,mdiff) + _RL ghat (imt,Nr) + _RL sigma(imt) + integer ikppkey, kkppkey #ifdef ALLOW_KPP @@ -1079,17 +1129,17 @@ c gat1*(imt) shape function at sigma = 1 c dat1*(imt) derivative of shape function at sigma = 1 c ws(imt), wm(imt) turbulent velocity scales (m/s) - _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt) - _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt) - _KPP_RL ws(imt), wm(imt) + _RL gat1m(imt), gat1s(imt), gat1t(imt) + _RL dat1m(imt), dat1s(imt), dat1t(imt) + _RL ws(imt), wm(imt) integer i, kn, ki - _KPP_RL R, dvdzup, dvdzdn, viscp - _KPP_RL difsp, diftp, visch, difsh, difth - _KPP_RL f1, sig, a1, a2, a3, delhat - _KPP_RL Gm, Gs, Gt - _KPP_RL tempVar + _RL R, dvdzup, dvdzdn, viscp + _RL difsp, diftp, visch, difsh, difth + _RL f1, sig, a1, a2, a3, delhat + _RL Gm, Gs, Gt + _RL tempVar - _KPP_RL p0 , eins + _RL p0 , eins parameter (p0=0.0, eins=1.0) c----------------------------------------------------------------------- @@ -1100,9 +1150,12 @@ sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon end do +CADJ STORE sigma = comlev1_kpp, key = ikppkey call wscale ( I sigma, hbl, ustar, bfsfc, - O wm, ws ) + O wm, ws, myThid ) +CADJ STORE wm = comlev1_kpp, key = ikppkey +CADJ STORE ws = comlev1_kpp, key = ikppkey do i = 1, imt wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i))) @@ -1145,20 +1198,39 @@ & max(ustar(i)**4,phepsi) gat1m(i) = visch / hbl(i) / wm(i) dat1m(i) = -viscp / wm(i) + f1 * visch - dat1m(i) = min(dat1m(i),p0) gat1s(i) = difsh / hbl(i) / ws(i) dat1s(i) = -difsp / ws(i) + f1 * difsh - dat1s(i) = min(dat1s(i),p0) gat1t(i) = difth / hbl(i) / ws(i) dat1t(i) = -diftp / ws(i) + f1 * difth - dat1t(i) = min(dat1t(i),p0) end do +#ifdef KPP_AUTODIFF_EXCESSIVE_STORE +CADJ STORE gat1m = comlev1_kpp, key = ikppkey +CADJ STORE gat1s = comlev1_kpp, key = ikppkey +CADJ STORE gat1t = comlev1_kpp, key = ikppkey +CADJ STORE dat1m = comlev1_kpp, key = ikppkey +CADJ STORE dat1s = comlev1_kpp, key = ikppkey +CADJ STORE dat1t = comlev1_kpp, key = ikppkey +#endif + do i = 1, imt + dat1m(i) = min(dat1m(i),p0) + dat1s(i) = min(dat1s(i),p0) + dat1t(i) = min(dat1t(i),p0) + end do +#ifdef KPP_AUTODIFF_EXCESSIVE_STORE +CADJ STORE dat1m = comlev1_kpp, key = ikppkey +CADJ STORE dat1s = comlev1_kpp, key = ikppkey +CADJ STORE dat1t = comlev1_kpp, key = ikppkey +#endif do ki = 1, Nr +#ifdef ALLOW_AUTODIFF_TAMC + kkppkey = (ikppkey-1)*Nr + ki +#endif + c----------------------------------------------------------------------- c compute turbulent velocity scales on the interfaces c----------------------------------------------------------------------- @@ -1167,9 +1239,16 @@ sig = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i) sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon) end do +#ifdef KPP_AUTODIFF_EXCESSIVE_STORE +CADJ STORE wm = comlev1_kpp_k, key = kkppkey +CADJ STORE ws = comlev1_kpp_k, key = kkppkey +#endif +CADJ STORE sigma = comlev1_kpp_k, key = kkppkey call wscale ( I sigma, hbl, ustar, bfsfc, - O wm, ws ) + O wm, ws, myThid ) +CADJ STORE wm = comlev1_kpp_k, key = kkppkey +CADJ STORE ws = comlev1_kpp_k, key = kkppkey c----------------------------------------------------------------------- c compute the dimensionless shape functions at the interfaces @@ -1213,9 +1292,16 @@ & + (1. - stable(i)) * min(sig,epsilon) end do +#ifdef KPP_AUTODIFF_EXCESSIVE_STORE +CADJ STORE wm = comlev1_kpp, key = ikppkey +CADJ STORE ws = comlev1_kpp, key = ikppkey +#endif +CADJ STORE sigma = comlev1_kpp, key = ikppkey call wscale ( I sigma, hbl, ustar, bfsfc, - O wm, ws ) + O wm, ws, myThid ) +CADJ STORE wm = comlev1_kpp, key = ikppkey +CADJ STORE ws = comlev1_kpp, key = ikppkey do i = 1, imt sig = -zgrid(kbl(i)-1) / hbl(i) @@ -1241,7 +1327,7 @@ I dkm1, hbl, kbl, diffus, casea U , ghat O , blmc - & ) + & , myThid ) c enhance the diffusivity at the kbl-.5 interface @@ -1256,27 +1342,29 @@ c kbl(imt) grid above hbl c diffus(imt,0:Nrp1,mdiff) vertical diffusivities (m^2/s) c casea(imt) = 1 in caseA, = 0 in case B - _KPP_RL dkm1 (imt,mdiff) - _KPP_RL hbl (imt) +c myThid thread number for this instance of the routine + integer myThid + _RL dkm1 (imt,mdiff) + _RL hbl (imt) integer kbl (imt) - _KPP_RL diffus(imt,0:Nrp1,mdiff) - _KPP_RL casea (imt) + _RL diffus(imt,0:Nrp1,mdiff) + _RL casea (imt) c input/output c nonlocal transport, modified ghat at kbl(i)-1 interface (s/m**2) - _KPP_RL ghat (imt,Nr) + _RL ghat (imt,Nr) c output c enhanced bound. layer mixing coeff. - _KPP_RL blmc (imt,Nr,mdiff) + _RL blmc (imt,Nr,mdiff) #ifdef ALLOW_KPP c local c fraction hbl lies beteen zgrid neighbors - _KPP_RL delta + _RL delta integer ki, i, md - _KPP_RL dkmp5, dstar + _RL dkmp5, dstar do i = 1, imt ki = kbl(i)-1 @@ -1302,8 +1390,8 @@ c************************************************************************* SUBROUTINE STATEKPP ( - I bi, bj, myThid, - O RHO1, DBLOC, DBSFC, TTALPHA, SSBETA) + O RHO1, DBLOC, DBSFC, TTALPHA, SSBETA, + I ikppkey, bi, bj, myThid ) c c----------------------------------------------------------------------- c "statekpp" computes all necessary input arrays @@ -1328,6 +1416,9 @@ c written by: jan morzel, feb. 10, 1995 (converted from "sigma" version) c modified by: d. menemenlis, june 1998 : for use with MIT GCM UV c +c 28 april 05: added computation of mixed layer depth KPPmld +c for a deltaRho equivalent to a deltaTheta of 0.8 deg C + c----------------------------------------------------------------------- IMPLICIT NONE @@ -1337,14 +1428,15 @@ #include "PARAMS.h" #include "KPP_PARAMS.h" #include "DYNVARS.h" +#include "GRID.h" c-------------- Routine arguments ----------------------------------------- INTEGER bi, bj, myThid - _KPP_RL RHO1 ( ibot:itop, jbot:jtop ) - _KPP_RL DBLOC ( ibot:itop, jbot:jtop, Nr ) - _KPP_RL DBSFC ( ibot:itop, jbot:jtop, Nr ) - _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 ) - _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 ) + _RL RHO1 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy ) + _RL DBLOC ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr ) + _RL DBSFC ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr ) + _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 ) + _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 ) #ifdef ALLOW_KPP @@ -1355,7 +1447,7 @@ c rhok - density of t(k ) & s(k ) at depth k c rhokm1 - density of t(k-1) & s(k-1) at depth k c rho1k - density of t(1 ) & s(1 ) at depth k -c work1, work2 - work arrays for holding horizontal slabs +c work1,2,3 - work arrays for holding horizontal slabs _RL RHOK (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy) @@ -1363,26 +1455,42 @@ _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) +#ifdef ALLOW_DIAGNOSTICS +c KPPMLD - mixed layer depth based on density criterion + _RL KPPMLD(1 :sNx ,1 :sNy ) +#endif /* ALLOW_DIAGNOSTICS */ + INTEGER I, J, K + INTEGER ikppkey, kkppkey c calculate density, alpha, beta in surface layer, and set dbsfc to zero + kkppkey = (ikppkey-1)*Nr + 1 + +#ifdef KPP_AUTODIFF_EXCESSIVE_STORE +CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey +CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey +#endif /* KPP_AUTODIFF_EXCESSIVE_STORE */ call FIND_RHO( - I bi, bj, ibot, itop, jbot, jtop, 1, 1, + I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, I theta, salt, O WORK1, I myThid ) +#ifdef KPP_AUTODIFF_EXCESSIVE_STORE +CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey +CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey +#endif /* KPP_AUTODIFF_EXCESSIVE_STORE */ call FIND_ALPHA( - I bi, bj, ibot, itop, jbot, jtop, 1, 1, - O WORK2 ) + I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, + O WORK2, myThid ) call FIND_BETA( - I bi, bj, ibot, itop, jbot, jtop, 1, 1, - O WORK3 ) + I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, + O WORK3, myThid ) - DO J = jbot, jtop - DO I = ibot, itop + DO J = 1-OLy, sNy+OLy + DO I = 1-OLx, sNx+OLx RHO1(I,J) = WORK1(I,J) + rhoConst TTALPHA(I,J,1) = WORK2(I,J) SSBETA(I,J,1) = WORK3(I,J) @@ -1390,39 +1498,71 @@ END DO END DO +#ifdef ALLOW_DIAGNOSTICS +c work3 - density of t(1)-.8 & s(1 ) at depth 1 + IF ( useDiagnostics ) THEN + DO J = 1, sNy + DO I = 1, sNx + KPPMLD(I,J) = ABS(R_low(I,J,bi,bj)) + WORK3 (I,J) = WORK1(I,J) - 0.8 _d 0 * WORK2(I,J) + END DO + END DO + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + c calculate alpha, beta, and gradients in interior layers CHPF$ INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2) DO K = 2, Nr + kkppkey = (ikppkey-1)*Nr + k + +#ifdef KPP_AUTODIFF_EXCESSIVE_STORE +CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey +CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey +#endif /* KPP_AUTODIFF_EXCESSIVE_STORE */ call FIND_RHO( - I bi, bj, ibot, itop, jbot, jtop, K, K, + I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K, I theta, salt, O RHOK, I myThid ) +#ifdef KPP_AUTODIFF_EXCESSIVE_STORE +CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey +CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey +#endif /* KPP_AUTODIFF_EXCESSIVE_STORE */ call FIND_RHO( - I bi, bj, ibot, itop, jbot, jtop, K-1, K, + I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K-1, K, I theta, salt, O RHOKM1, I myThid ) +#ifdef KPP_AUTODIFF_EXCESSIVE_STORE +CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey +CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey +#endif /* KPP_AUTODIFF_EXCESSIVE_STORE */ call FIND_RHO( - I bi, bj, ibot, itop, jbot, jtop, 1, K, + I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, K, I theta, salt, O RHO1K, I myThid ) +#ifdef KPP_AUTODIFF_EXCESSIVE_STORE +CADJ STORE rhok (:,:) = comlev1_kpp_k, key=kkppkey +CADJ STORE rhokm1(:,:) = comlev1_kpp_k, key=kkppkey +CADJ STORE rho1k (:,:) = comlev1_kpp_k, key=kkppkey +#endif /* KPP_AUTODIFF_EXCESSIVE_STORE */ + call FIND_ALPHA( - I bi, bj, ibot, itop, jbot, jtop, K, K, - O WORK1 ) + I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K, + O WORK1, myThid ) call FIND_BETA( - I bi, bj, ibot, itop, jbot, jtop, K, K, - O WORK2 ) + I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K, + O WORK2, myThid ) - DO J = jbot, jtop - DO I = ibot, itop + DO J = 1-OLy, sNy+OLy + DO I = 1-OLx, sNx+OLx TTALPHA(I,J,K) = WORK1 (I,J) SSBETA(I,J,K) = WORK2 (I,J) DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) / @@ -1432,17 +1572,51 @@ END DO END DO +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN +c work1 - density of t(k-1) & s(k-1) at depth 1 +c work2 - density of t(k ) & s(k ) at depth 1 +c work3 - density of t(1)-.8 & s(1 ) at depth 1 + call FIND_RHO( + I bi, bj, 1, sNx, 1, sNy, K-1, 1, theta, salt, + O WORK1, + I myThid ) + call FIND_RHO( + I bi, bj, 1, sNx, 1, sNy, K , 1, theta, salt, + O WORK2, + I myThid ) + DO J = 1, sNy + DO I = 1, sNx + IF ( k .LE. klowC(I,J,bi,bj) .AND. + & WORK1(I,J) .LT. WORK3(I,J) .AND. + & WORK2(I,J) .GE. WORK3(I,J) ) THEN + KPPMLD(I,J) = MIN ( KPPMLD(I,J), + & ABS(((WORK3(I,J)-WORK1(I,J))*rC(k)+ + & (WORK2(I,J)-WORK3(I,J))*rC(k-1))/ + & (WORK2(I,J)-WORK1(I,J)))) + ENDIF + END DO + END DO + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + END DO c compute arrays for K = Nrp1 - DO J = jbot, jtop - DO I = ibot, itop + DO J = 1-OLy, sNy+OLy + DO I = 1-OLx, sNx+OLx TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr) SSBETA(I,J,Nrp1) = SSBETA(I,J,Nr) DBLOC(I,J,Nr) = 0. END DO END DO +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + CALL DIAGNOSTICS_FILL(KPPmld,'KPPmld ',0,1,3,bi,bj,myThid) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + #endif /* ALLOW_KPP */ RETURN