PROGRAM HYCOM_PROFILE_HYBGEN IMPLICIT NONE C C hycom_profile_hybgen - Usage: hycom_profile_hybgen archv.txt blkdat.input archg.txt C C apply the HYCOM hybrid grid generator to a HYCOM profile. C note that the velocities are unchanged on exit. C based on hybgen version 2.3.03 C C archv.txt is assumed to be an HYCOM archive text profile file C blkdat.input is assumed to be an HYCOM parameter file (version 2.3.03) C archg.txt will be the output text profile file C C this version for "serial" Unix systems. C C Alan J. Wallcraft, Naval Research Laboratory, August 2003 C INTEGER IARGC INTEGER NARG CHARACTER*240 CARG C CHARACTER*240 CFILEA,CFILEB,CFILEC CHARACTER*240 CLINE REAL THK,THKIN,FLAG,ONEM REAL U(99),V(99),T(99),S(99),R(99),DP(99),P(99+1) REAL RIN(99),TIN(99),SIN(99),PIN(99+1) REAL*8 SUMI(3),SUMO(3),SUMD(3) INTEGER IOS,K,KDM LOGICAL LDEBUG C LOGICAL ISOPCM INTEGER KDMBLK,NHYBRD,NSIGMA,HYBMAP,HYBFLG REAL HYBISO,HYBRLX,QHYBRLX REAL DP00,DP00X,DP00F,DS00,DS00X,DS00F,DP00I,THETA(999) REAL DP0K(999),DS0K(999),DEPTHS,ISOTOP,DPNS,DSNS REAL THKBOT c integer mnproc,lp common/xxx/ mnproc,lp c integer thflag common/th/thflag c mnproc = 1 lp = 6 C C READ ARGUMENTS. C NARG = IARGC() C CALL GETARG(0,CLINE) K = LEN_TRIM(CLINE) LDEBUG = CLINE(K-5:K) .EQ. "_debug" C IF (NARG.EQ.3) THEN CALL GETARG(1,CFILEA) CALL GETARG(2,CFILEB) CALL GETARG(3,CFILEC) ELSE WRITE(6,*) + 'Usage: hycom_profile_hybgen archv.txt blkdat.input archg.txt' CALL EXIT(1) ENDIF C C OPEN ALL FILES. C OPEN(UNIT=11, FILE=CFILEA, FORM='FORMATTED', STATUS='OLD', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',CFILEA(1:LEN_TRIM(CFILEA)) WRITE(6,*) 'ios = ',ios CALL EXIT(3) ENDIF OPEN(UNIT=99, FILE=CFILEB, FORM='FORMATTED', STATUS='OLD', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',CFILEB(1:LEN_TRIM(CFILEB)) WRITE(6,*) 'ios = ',ios CALL EXIT(3) ENDIF OPEN(UNIT=21, FILE=CFILEC, FORM='FORMATTED', STATUS='NEW', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',CFILEC(1:LEN_TRIM(CFILEC)) WRITE(6,*) 'ios = ',ios CALL EXIT(5) ENDIF C C INPUT BLKDAT PARAMETERS C DO K= 1,10 READ(99,*) ENDDO c c --- 'kdm ' = number of layers c --- 'nhybrd' = number of hybrid levels (0=all isopycnal) c --- 'nsigma' = number of sigma levels (nhybrd-nsigma z-levels) c --- 'dp00' = deep z-level spacing minimum thickness (m) c --- 'dp00x' = deep z-level spacing maximum thickness (m) c --- 'dp00f' = deep z-level spacing stretching factor (1.0=const.z) c --- 'ds00' = shallow z-level spacing minimum thickness (m) c --- 'ds00x' = shallow z-level spacing maximum thickness (m) c --- 'ds00f' = shallow z-level spacing stretching factor (1.0=const.z) c --- 'dp00i' = deep iso-pycnal spacing minimum thickness (m) c --- 'isotop' = shallowest depth for isopycnal layers (m), <0 from file c c --- or, in place of 'dp00','dp00x','dp00f','ds00','ds00x','ds00f' specify: c --- 'dp0k ' = layer k deep z-level spacing minimum thickness (m) c --- k=1,kdm; dp0k must be zero for k>nhybrd c --- 'ds0k ' = layer k shallow z-level spacing minimum thickness (m) c --- k=1,nsigma c call blkini(kdmblk,'kdm ') call blkini(nhybrd,'nhybrd') call blkini(nsigma,'nsigma') call blkinr2(dp00,k, & 'dp00 ','(a6," =",f10.4," m")', & 'dp0k ','(a6," =",f10.4," m")' ) call flush(lp) if (k.eq.1) then !dp00 call blkinr(dp00x, 'dp00x ','(a6," =",f10.4," m")') call blkinr(dp00f, 'dp00f ','(a6," =",f10.4," ")') call blkinr(ds00, 'ds00 ','(a6," =",f10.4," m")') call blkinr(ds00x, 'ds00x ','(a6," =",f10.4," m")') call blkinr(ds00f, 'ds00f ','(a6," =",f10.4," ")') else !dp0k dp0k(1) = dp00 dp00 = -1.0 !signal that dp00 is not input do k=2,kdmblk call blkinr(dp0k(k), 'dp0k ','(a6," =",f10.4," m")') c if (k.gt.nhybrd .and. dp0k(k).ne.0.0) then write(lp,'(/ a,i3 /)') & 'error - dp0k must be zero for k>nhybrd' call flush(lp) stop '(blkdat)' endif !k>nhybrd&dp0k(k)!=0 enddo !k do k=1,nsigma call blkinr(ds0k(k), 'ds0k ','(a6," =",f10.4," m")') enddo !k endif !dp00:dp0k call blkinr(dp00i, 'dp00i ','(a6," =",f10.4," m")') call blkinr(isotop,'isotop','(a6," =",f10.4," m")') call flush(lp) c c --- 'thflag' = reference pressure flag (0=Sigma-0, 2=Sigma-2) c READ(99,*) !'saln0 ' READ(99,*) !'locsig' READ(99,*) !'kapflg' call blkini(thflag,'thflag') READ(99,*) !'thbase' READ(99,*) !'vsigma' c c --- 'sigma ' = isopycnal layer target densities (sigma units) c do k=1,kdmblk call blkinr(THETA(k),'sigma ','(a6," =",f10.4," sigma-0")') enddo call flush(lp) c c --- skip 'iniflg' to 'hybraf'. c DO K= 1,23 READ(99,*) ENDDO c c --- 'hybrlx' = HYBGEN: inverse relaxation coefficient (time steps) c (1.0 for no relaxation) c --- 'hybiso' = HYBGEN: Use PCM if layer is within hybiso of target density c (0.0 for no PCM; large to recover pre-2.2.09 behaviour) c --- 'hybmap' = HYBGEN: remapper flag (0=PCM, 1=PLM, 2=PPM, 3=WENO-like) c --- 'hybflg' = HYBGEN: generator flag (0=T&S, 1=th&S, 2=th&T) c call blkinr(hybrlx,'hybrlx','(a6," =",f10.4," time steps")') call blkinr(hybiso,'hybiso','(a6," =",f10.4," kg/m^3")') call blkini(hybmap,'hybmap') call blkini(hybflg,'hybflg') call flush(lp) c qhybrlx = 1.0/hybrlx isopcm = hybiso.gt.0.0 !use PCM for isopycnal layers? C CLOSE(99) C C INITIALIZE FIXED GRID. C CALL GEOPAR(DP00,DP00X,DP00F,DS00,DS00X,DS00F,DP00I, + NHYBRD,NSIGMA,KDMBLK, + DP0K,DS0K,DPNS,DSNS) * call flush(lp) C C COPY PROFILE HEADER TO OUTPUT. C DO K= 1,99 READ( 11,'(a)') CLINE IF (CLINE(1:5).EQ.'# k ') then EXIT ENDIF WRITE(21,'(a)') TRIM(CLINE) ENDDO C C READ THE PROFILE. C P(1) = 0.0 KDM = -1 DO K= 1,99 READ(11,'(a)',IOSTAT=IOS) CLINE IF (IOS.NE.0) THEN IF (K.NE.KDM+1) THEN WRITE(6,*) 'Error: inconsistent input profile' CALL EXIT(6) ENDIF EXIT ENDIF READ(CLINE,*) KDM,U(K),V(K),T(K),S(K),R(K),THK P(K+1) = P(K) + THK IF (THK.EQ.0.0) THEN U(K) = U(K-1) V(K) = V(K-1) T(K) = T(K-1) S(K) = S(K-1) R(K) = R(K-1) ENDIF ENDDO CLOSE(11) C C TRANSFORM THE PROFILE. C ONEM = 9806.0 !g/thref DO K= 1,KDM DP(K) = (P(K+1) - P(K))*ONEM C RIN(K) = R(K) TIN(K) = T(K) SIN(K) = S(K) PIN(K) = P(K) ENDDO PIN(KDM+1) = P(KDM+1) C DEPTHS = P(KDM+1) ISOTOP = ISOTOP*ONEM THKBOT = 10.0 C CALL HYBGEN(T,S,R,DP,THETA,KDM, + NHYBRD,ISOPCM,HYBMAP,HYBFLG,HYBISO, QHYBRLX, + DP00I,DP0K,DS0K,DPNS,DSNS,DEPTHS,ISOTOP,THKBOT, + LDEBUG) C DO K= 1,KDM P(K+1) = P(K) + DP(K)/ONEM ENDDO C C DIAGNOSTIC OUTPUT C WRITE(6,*) WRITE(6,'(a,a)') & '# k ', & ' sigma dens dens deld thkns thkns' DO K= 1,KDM THK = P( K+1) - P( K) THKIN = PIN(K+1) - PIN(K) IF (ABS(THK-THKIN).LT.0.001) THEN WRITE(6,'(i4,1x,4f8.3,2f9.3)') & K,THETA(K),RIN(K),R(K),R(K)-THETA(K),THKIN,THK ELSE WRITE(6,'(i4,1x,4f8.3,2f9.3,a)') & K,THETA(K),RIN(K),R(K),R(K)-THETA(K),THKIN,THK,' *' ENDIF ENDDO C SUMI(:)=0.0 SUMO(:)=0.0 DO K= 1,KDM THK = P( K+1) - P( K) THKIN = PIN(K+1) - PIN(K) SUMI(1) = SUMI(1) + THKIN*TIN(K) SUMI(2) = SUMI(2) + THKIN*SIN(K) SUMI(3) = SUMI(3) + THKIN*RIN(K) SUMO(1) = SUMO(1) + THK *T( K) SUMO(2) = SUMO(2) + THK *S( K) SUMO(3) = SUMO(3) + THK *R( K) ENDDO DO K=1,3 SUMI(K) = SUMI(K)/PIN(KDM+1) SUMO(K) = SUMO(K)/P( KDM+1) SUMD(K) = SUMO(K) - SUMI(K) ENDDO WRITE(*,'(/A)') ' LABEL: T S R' WRITE(*,'(A,1p3e14.6)') ' INPUT TOTALS: ',SUMI(:) WRITE(*,'(A,1p3e14.6)') 'OUTPUT TOTALS: ',SUMO(:) WRITE(*,'(A,1p3e14.6)') 'OUT-IN TOTALS: ',SUMD(:) WRITE(*,*) C C OUTPUT THE TRANSFORMED PROFILE. C WRITE(21,'(3a)') & '# k', & ' utot vtot p.temp saln p.dens', & ' thkns dpth' DO K= 1,KDM THK = P(K+1) - P(K) WRITE(21,'(i4,2f8.2,3f8.3,1pg15.6,0pf13.6)') & K,U(K),V(K),T(K),S(K),R(K),THK,0.5*(P(K)+P(K+1)) ENDDO END subroutine blkinr(rvar,cvar,cfmt) implicit none c real rvar character cvar*6,cfmt*(*) c c read in one real value c character*6 cvarin c read(99,*) rvar,cvarin write(6,cfmt) cvarin,rvar c if (cvar.ne.cvarin) then write(6,*) write(6,*) 'error in blkinr - input ',cvarin, + ' but should be ',cvar write(6,*) stop endif return end subroutine blkinr2(rvar,nvar,cvar1,cfmt1,cvar2,cfmt2) implicit none c real rvar integer nvar character cvar1*6,cvar2*6,cfmt1*(*),cfmt2*(*) c c read in one real value from stdin, c identified as either cvar1 (return nvar=1) or cvar2 (return nvar=2) c character*6 cvarin c read(99,*) rvar,cvarin c if (cvar1.eq.cvarin) then nvar = 1 write(6,cfmt1) cvarin,rvar elseif (cvar2.eq.cvarin) then nvar = 2 write(6,cfmt2) cvarin,rvar else write(6,*) write(6,*) 'error in blkinr2 - input ',cvarin, + ' but should be ',cvar1,' or ',cvar2 write(6,*) stop endif return end subroutine blkini(ivar,cvar) implicit none c integer ivar character*6 cvar c c read in one integer value c character*6 cvarin c read(99,*) ivar,cvarin write(6,6000) cvarin,ivar c if (cvar.ne.cvarin) then write(6,*) write(6,*) 'error in blkini - input ',cvarin, + ' but should be ',cvar write(6,*) stop endif return 6000 format(a6,' =',i6) end subroutine geopar(dp00,dp00x,dp00f,ds00,ds00x,ds00f,dp00i, & nhybrd,nsigma,kdm, & dp0k,ds0k,dpns,dsns) implicit none c integer nhybrd,nsigma,kdm real dp00,dp00x,dp00f,ds00,ds00x,ds00f,dp00i real dp0k(kdm),ds0k(kdm),dpns,dsns c integer mnproc,lp common/xxx/ mnproc,lp c c --- set up model parameters related to geography c real dp0kf,dpm,dpms,ds0kf,dsm,dsms,onem,qonem integer k c onem = 9806.0 !g/thref qonem = 1.0/onem c c --- calculate dp0k and ds0k? if (dp00.lt.0.0) then c --- dp0k and ds0k already input dp00 =onem*dp0k(1) dp00x=onem*dp0k(kdm-1) dp00i=onem*dp00i dpms = 0.0 do k=1,kdm dpm = dp0k(k) dpms = dpms + dpm dp0k(k) = dp0k(k)*onem if (mnproc.eq.1) then write(lp,135) k,dp0k(k)*qonem,dpm,dpms endif if (mnproc.eq.-99) then ! bugfix that prevents optimization write(6,*) 'geopar: dp0k = ',dp0k(k),k,mnproc endif enddo !k dsms = 0.0 do k=1,nsigma dsm = ds0k(k) dsms = dsms + dsm ds0k(k) = ds0k(k)*onem if (mnproc.eq.1) then write(lp,130) k,ds0k(k)*qonem,dsm,dsms endif if (mnproc.eq.-99) then ! bugfix that prevents optimization write(6,*) 'geopar: ds0k = ',ds0k(k),k,mnproc endif enddo !k if (mnproc.eq.1) then write(lp,*) endif else c --- calculate dp0k and ds0k c c --- logorithmic k-dependence of dp0 (deep z's) dp00 =onem*dp00 dp00x=onem*dp00x dp00i=onem*dp00i if (nhybrd.eq.0) then * dp0k(1)=thkmin*onem dp0k(1)= 20.0*onem else dp0k(1)=dp00 endif dpm = dp0k(1)*qonem dpms = dpm write(lp,*) write(lp,135) 1,dp0k(1)*qonem,dpm,dpms 135 format('dp0k(',i2,') =',f7.2,' m', & ' thkns =',f7.2,' m', & ' depth =',f8.2,' m') c dp0kf=1.0 do k=2,kdm dp0kf=dp0kf*dp00f if (k.le.nhybrd) then dp0k(k)=min(dp00*dp0kf,dp00x) else dp0k(k)=0.0 endif dpm = dp0k(k)*qonem dpms = dpms + dpm write(lp,135) k,dp0k(k)*qonem,dpm,dpms if (mnproc.eq.-99) then ! bugfix that prevents optimization write(6,*) 'geopar: dp0kf = ',dp0kf, mnproc write(6,*) 'geopar: dp0k = ',dp0k(k),k,mnproc endif enddo c c --- logorithmic k-dependence of ds0 (shallow z-s) ds00 =onem*ds00 ds00x=onem*ds00x if (nhybrd.eq.0) then * ds0k(1)=thkmin*onem ds0k(1)= 20.0*onem else ds0k(1)=ds00 endif dsm = ds0k(1)*qonem dsms = dsm if (mnproc.eq.1) then write(lp,*) write(lp,130) 1,ds0k(1)*qonem,dsm,dsms endif 130 format('ds0k(',i2,') =',f7.2,' m', & ' thkns =',f7.2,' m', & ' depth =',f8.2,' m') c ds0kf=1.0 do k=2,nsigma ds0kf=ds0kf*ds00f ds0k(k)=min(ds00*ds0kf,ds00x) dsm = ds0k(k)*qonem dsms = dsms + dsm if (mnproc.eq.1) then write(lp,130) k,ds0k(k)*qonem,dsm,dsms endif if (mnproc.eq.-99) then ! bugfix that prevents optimization write(6,*) 'geopar: ds0kf = ',ds0kf, mnproc write(6,*) 'geopar: ds0k = ',ds0k(k),k,mnproc endif enddo if (mnproc.eq.1) then write(lp,*) endif endif !input:calculate dp0k,ds0k c c --- start and stop depths for terrain following coordinate if (nsigma.eq.0) then dpns = dp0k(1) dsns = 0.0 ds0k(1) = dp0k(1) do k= 2,kdm ds0k(k)=0.0 enddo !k else dpns = 0.0 dsns = 0.0 do k=1,nsigma dpns = dpns + dp0k(k) dsns = dsns + ds0k(k) enddo !k do k= nsigma+1,kdm ds0k(k)=0.0 enddo !k endif !nsigma dpns = dpns*qonem !depths is in m dsns = dsns*qonem !depths is in m c if (mnproc.eq.1) then write(lp,131) nsigma,dpns,dsns endif 131 format('nsigma = ',i2, & ' deep =',f8.2,' m', & ' shallow =',f8.2,' m' ) call flush(lp) c return end real function sig(t,s) real t,s c integer thflag common/th/thflag c c --- equation of state c real sig0,sig2 c if (thflag.eq.0) then sig = sig0(t,s) else sig = sig2(t,s) endif end real function dsigdt(t,s) real t,s c integer thflag common/th/thflag c c --- d(sig)/dt c real dsigdt0,dsigdt2 c if (thflag.eq.0) then dsigdt = dsigdt0(t,s) else dsigdt = dsigdt2(t,s) endif end real function dsigds(t,s) real t,s c integer thflag common/th/thflag c c --- d(sig)/ds c real dsigds0,dsigds2 c if (thflag.eq.0) then dsigds = dsigds0(t,s) else dsigds = dsigds0(t,s) endif end c real function tofsig(r,s) real r,s c integer thflag common/th/thflag c c --- temp (deg c) as a function of sigma and salinity (mil) c real tofsig0,tofsig2 c if (thflag.eq.0) then tofsig = tofsig0(r,s) else tofsig = tofsig2(r,s) endif end real function sofsig(r,t) real r,t c integer thflag common/th/thflag c c --- salinity (mil) as a function of sigma and temperature (deg c) c real sofsig0,sofsig2 c if (thflag.eq.0) then sofsig = sofsig0(r,t) else sofsig = sofsig2(r,t) endif end real function sig0(tx,sx) real tx,sx c c --- equation of state c include 'hycom_profile_hybgen_SIGMA0.h' c sig0 = sig(tx,sx) end real function dsigdt0(tx,sx) real tx,sx c c --- d(sig)/dt c include 'hycom_profile_hybgen_SIGMA0.h' c dsigdt0 = dsigdt(tx,sx) end real function dsigds0(tx,sx) real tx,sx c c --- d(sig)/ds c include 'hycom_profile_hybgen_SIGMA0.h' c dsigds0 = dsigds(tx,sx) end c real function tofsig0(rx,sx) real rx,sx c c --- temp (deg c) as a function of sigma and salinity (mil) c include 'hycom_profile_hybgen_SIGMA0.h' c tofsig0 = tofsig(rx,sx) end real function sofsig0(rx,tx) real rx,tx c c --- salinity (mil) as a function of sigma and temperature (deg c) c include 'hycom_profile_hybgen_SIGMA0.h' c sofsig0 = sofsig(rx,tx) end real function sig2(tx,sx) real tx,sx c c --- equation of state c include 'hycom_profile_hybgen_SIGMA2.h' c sig2 = sig(tx,sx) end real function dsigdt2(tx,sx) real tx,sx c c --- d(sig)/dt c include 'hycom_profile_hybgen_SIGMA2.h' c dsigdt2 = dsigdt(tx,sx) end real function dsigds2(tx,sx) real tx,sx c c --- d(sig)/ds c include 'hycom_profile_hybgen_SIGMA2.h' c dsigds2 = dsigds(tx,sx) end c real function tofsig2(rx,sx) real rx,sx c c --- temp (deg c) as a function of sigma and salinity (mil) c include 'hycom_profile_hybgen_SIGMA2.h' c tofsig2 = tofsig(rx,sx) end real function sofsig2(rx,tx) real rx,tx c c --- salinity (mil) as a function of sigma and temperature (deg c) c include 'hycom_profile_hybgen_SIGMA2.h' c sofsig2 = sofsig(rx,tx) end subroutine hybgen(temp,saln,th3d,dp,theta,kdm, & nhybrd,isopcm,hybmap,hybflg, & hybiso, qhybrlx, & dp00i,dp0k,ds0k,dpns,dsns,depths,topiso, & thkbot, ldebug) implicit none c logical isopcm,ldebug integer kdm, nhybrd,hybmap,hybflg real temp(1,1,kdm,1), & saln(1,1,kdm,1), & th3d(1,1,kdm,1), & dp(1,1,kdm,1), & theta(1,1,kdm) real hybiso,qhybrlx,thkbot real dp00i,dp0k(kdm),ds0k(kdm),dpns,dsns, & depths(1,1),topiso(1,1) c integer mnproc,lp common/xxx/ mnproc,lp c integer klist(1,1) real p(1,1,kdm+1),dpmixl(1,1,1),q2(1,1,1,1),q2l(1,1,1,1) real tracer(1,1,1,1,1),trcflg(1) logical mxlkta,thermo integer j,kk,n, nstep, i0,j0,itest,jtest real epsil,onem,tenm,tencm,onemm,qonem,thbase c real sig,dsigdt,dsigds,tofsig,sofsig external sig,dsigdt,dsigds,tofsig,sofsig c integer, parameter :: mxtrcr=1 integer, parameter :: ntracr=0 logical, parameter :: mxlmy =.false. c c --- --------------------- c --- hybrid grid generator c --- --------------------- c logical, parameter :: lunmix=.true. !unmix a too light deepest layer logical, parameter :: lconserve=.false. !explicitly conserve each column c double precision asum( mxtrcr+4,3) real offset(mxtrcr+4) c logical lcm(kdm) !use PCM for some layers? real s1dsum( mxtrcr+4), !sum of scalar fields & s1d(kdm,mxtrcr+4), !original scalar fields & f1d(kdm,mxtrcr+4), !final scalar fields & c1d(kdm,mxtrcr+4,3), !interpolation coefficients & dpi( kdm), !original layer thicknesses, >= onezm & dprs(kdm), !original layer thicknesses & pres(kdm+1), !original layer interfaces & prsf(kdm+1), !final layer interfaces & qhrlx( kdm+1), !relaxation coefficient, from qhybrlx & dp0ij( kdm), !minimum layer thickness & dp0cum(kdm+1) !minimum interface depth real p_hat,p_hat0,p_hat2,p_hat3,hybrlx, & delt,deltm,dels,delsm,q,qsum,qdep,qtr,qts,thkbop, & zthk,oneum,onezm integer i,k,ka,kp,ktr,l,fixlay,nums1d character*12 cinfo c double precision, parameter :: zp5=0.5 !for sign function c c --- c u s h i o n function (from Bleck & Benjamin, 1992): c --- if delp >= qqmx*dp0 >> dp0, -cushn- returns -delp- c --- if delp <= qqmn*dp0 << -dp0, -cushn- returns -dp0- c real qqmn,qqmx,cusha,cushb parameter (qqmn=-4.0, qqmx=2.0) ! shifted range * parameter (qqmn=-2.0, qqmx=4.0) ! traditional range * parameter (qqmn=-4.0, qqmx=6.0) ! somewhat wider range parameter (cusha=qqmn**2*(qqmx-1.0)/(qqmx-qqmn)**2) parameter (cushb=1.0/qqmn) c real qq,cushn,delp,dp0 * include 'stmt_fns.h' qq( delp,dp0)=max(qqmn, min(qqmx, delp/dp0)) cushn(delp,dp0)=dp0* & (1.0+cusha*(1.0-cushb*qq(delp,dp0))**2)* & max(1.0, delp/(dp0*qqmx)) c mxlkta = .false. thermo = .true. i = 1 j = 1 kk = kdm n = 1 nstep = 1 thbase = 0.0 i0 = 0 j0 = 0 if (ldebug) then itest = 1 jtest = 1 else itest = 0 jtest = 0 endif c epsil = 1.0e-11 onem = 9806.0 !g/thref qonem = 1.0/onem tenm = onem*10.0 tencm = onem/10.0 onemm = onem/1000.0 c onezm = onem*1.e-20 !a zepto-metre, an effectively zero thickness oneum = onem*1.e-6 !a micro-metre, a barely non-zero thickness thkbop = thkbot*onem hybrlx = 1.0/qhybrlx c if (mxlmy) then nums1d = ntracr + 4 else nums1d = ntracr + 2 endif c if (.not.isopcm) then c --- lcm the same for all points do k=1,nhybrd lcm(k) = .false. !use same remapper for all layers enddo !k do k=nhybrd+1,kk lcm(k) = .true. !purely isopycnal layers use PCM enddo !k endif c * do l=1,isp(j) * do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) c c --- terrain following starts at depth dpns and ends at depth dsns qdep = max( 0.0, min( 1.0, & (depths(i,j) - dsns)/ & (dpns - dsns) ) ) c if (qdep.lt.1.0) then c --- terrain following, qhrlx=1 and ignore dp00 p(i,j, 1)=0.0 dp0cum(1)=0.0 qhrlx( 1)=1.0 dp0ij( 1)=qdep*dp0k(1) + (1.0-qdep)*ds0k(1) if (i.eq.itest .and. j.eq.jtest) then k=1 write (lp,*) 'qdep = ',qdep write (lp,'(a/i6,1x,4f9.3/a)') . ' k dp0ij ds0k dp0k p', . k,dp0ij(k)*qonem,ds0k(1)*qonem,dp0k(1)*qonem, . p(i,j,k)*qonem, . ' k dp0ij p-cum p dp0cum' endif !debug dp0cum(2)=dp0cum(1)+dp0ij(1) qhrlx( 2)=1.0 p(i,j, 2)=p(i,j,1)+dp(i,j,1,n) do k=2,kk qhrlx( k+1)=1.0 dp0ij( k) =qdep*dp0k(k) + (1.0-qdep)*ds0k(k) dp0cum(k+1)=dp0cum(k)+dp0ij(k) p(i,j, k+1)=p(i,j,k)+dp(i,j,k,n) if (i.eq.itest .and. j.eq.jtest) then write (lp,'(i6,1x,4f9.3)') . k,dp0ij(k)*qonem,p(i,j,k)*qonem-dp0cum(k)*qonem, . p(i,j,k)*qonem,dp0cum(k)*qonem endif !debug enddo !k else c --- not terrain following p(i,j, 1)=0.0 dp0cum(1)=0.0 qhrlx( 1)=1.0 !no relaxation in top layer dp0ij( 1)=dp0k(1) if (i.eq.itest .and. j.eq.jtest) then k=1 write (lp,*) 'qdep = ',qdep write (lp,'(a/i6,1x,3f8.3,f9.3)') . ' k dp0ij ds0k dp0k p dp0cum', . k,dp0ij(k)*qonem,ds0k(1)*qonem,dp0k(1)*qonem, . p(i,j,k)*qonem write (lp,'(a/i6,1x,f9.3/a)') . ' k dp0ij', . k,dp0ij(k)*qonem, . ' k dp0ij q p-cum p dp0cum' endif !debug dp0cum(2)=dp0cum(1)+dp0ij(1) qhrlx( 2)=1.0 !no relaxation in top layer p(i,j, 2)=p(i,j,1)+dp(i,j,1,n) do k=2,kk c --- q is dp0k(k) when in surface fixed coordinates c --- q is dp00i when much deeper than surface fixed coordinates if (dp0k(k).le.dp00i) then q = dp0k(k) qts= 0.0 !0 at dp0k else q = max( dp00i, & dp0k(k) * dp0k(k)/ & max( dp0k( k), & p(i,j,k)-dp0cum(k) ) ) qts= 1.0 - (q-dp00i)/(dp0k(k)-dp00i) !0 at dp0k, 1 at dp00i endif qhrlx( k+1)=1.0/(1.0 + qts*(hybrlx-1.0)) !1 at dp0k, qhybrlx at dp00i dp0ij( k) =min( q, dp0k(k) ) dp0cum(k+1)=dp0cum(k)+dp0ij(k) p(i,j, k+1)=p(i,j,k)+dp(i,j,k,n) if (i.eq.itest .and. j.eq.jtest) then write (lp,'(i6,1x,5f9.3)') . k,dp0ij(k)*qonem,q*qonem,p(i,j,k)*qonem-dp0cum(k)*qonem, . p(i,j,k)*qonem,dp0cum(k)*qonem endif !debug enddo !k endif !qdep<1:else c c --- identify the always-fixed coordinate layers fixlay = 1 !layer 1 always fixed do k= 2,nhybrd if (dp0cum(k).ge.topiso(i,j)) then exit !layers k to nhybrd can be isopycnal endif qhrlx(k+1)=1.0 !no relaxation in fixed layers fixlay = fixlay+1 enddo !k if (i.eq.itest .and. j.eq.jtest) then write(lp,'(a,i3)') & 'hybgen, always-fixed coordinate layers: 1 to ', & fixlay call flush(lp) endif !debug c if (i.eq.itest .and. j.eq.jtest) then write (lp,'(a/(i6,1x,2f8.3,2f9.3,f9.3))') . 'hybgen: thkns minthk dpth mindpth hybrlx', . (k,dp(i,j,k,n)*qonem, dp0ij(k)*qonem, . p(i,j,k+1)*qonem,dp0cum(k+1)*qonem, . 1.0/qhrlx(k+1), . k=1,kk) endif !debug c c --- identify the deepest layer kp with significant thickness (> oneum) c kp = 2 !minimum allowed value do k=kk,3,-1 if (p(i,j,k+1)-p(i,j,k).ge.oneum) then kp=k exit endif enddo c k=kp !at least 2 c if (i.eq.itest .and. j.eq.jtest) then write(lp,'(a,i3)') & 'hybgen, deepest inflated layer:',k call flush(lp) endif !debug c ka = max(k-2,1) !k might be 2 if (k.gt.fixlay+1 .and. & theta(i,j,k)-epsil.gt.th3d(i,j,k,n) .and. & th3d(i,j,k-1,n) .gt.th3d(i,j,k,n) .and. & th3d(i,j,ka, n) .gt.th3d(i,j,k,n) ) then c c --- water in the deepest inflated layer with significant thickness c --- (kp) is too light, and it is lighter than the two layers above. c --- c --- this should only occur when relaxing or nudging layer thickness c --- and is a bug (bad interaction with tsadvc) even in those cases c --- c --- entrain the entire layer into the one above c--- note the double negative in T=T-q*(T-T'), equiv. to T=T+q*(T'-T) q = (p(i,j,k+1)-p(i,j,k))/(p(i,j,k+1)-p(i,j,k-1)) if (hybflg.eq.0) then !T&S temp(i,j,k-1,n)=temp(i,j,k-1,n)-q*(temp(i,j,k-1,n) - & temp(i,j,k, n) ) saln(i,j,k-1,n)=saln(i,j,k-1,n)-q*(saln(i,j,k-1,n) - & saln(i,j,k, n) ) th3d(i,j,k-1,n)=sig(temp(i,j,k-1,n),saln(i,j,k-1,n))-thbase elseif (hybflg.eq.1) then !th&S th3d(i,j,k-1,n)=th3d(i,j,k-1,n)-q*(th3d(i,j,k-1,n) - & th3d(i,j,k, n) ) saln(i,j,k-1,n)=saln(i,j,k-1,n)-q*(saln(i,j,k-1,n) - & saln(i,j,k, n) ) temp(i,j,k-1,n)=tofsig(th3d(i,j,k-1,n)+thbase,saln(i,j,k-1,n)) elseif (hybflg.eq.2) then !th&T th3d(i,j,k-1,n)=th3d(i,j,k-1,n)-q*(th3d(i,j,k-1,n) - & th3d(i,j,k, n) ) temp(i,j,k-1,n)=temp(i,j,k-1,n)-q*(temp(i,j,k-1,n) - & temp(i,j,k, n) ) saln(i,j,k-1,n)=sofsig(th3d(i,j,k-1,n)+thbase,temp(i,j,k-1,n)) endif do ktr= 1,ntracr tracer(i,j,k-1,n,ktr)=tracer(i,j,k-1,n,ktr)- & q*(tracer(i,j,k-1,n,ktr) - & tracer(i,j,k, n,ktr) ) enddo !ktr if (mxlmy) then q2( i,j,k-1,n)=q2( i,j,k-1,n)- & q*(q2( i,j,k-1,n)-q2( i,j,k,n)) q2l(i,j,k-1,n)=q2l(i,j,k-1,n)- & q*(q2l(i,j,k-1,n)-q2l(i,j,k,n)) endif c --- entrained the entire layer into the one above, so now kp=kp-1 p(i,j,k) = p(i,j,k+1) kp = k-1 if (i.eq.itest .and. j.eq.jtest) then write(lp,'(a,i3,f6.3,5f8.3)') & 'hybgen, 11(+):', & k-1,q,temp(i,j,k-1,n),saln(i,j,k-1,n), & th3d(i,j,k-1,n)+thbase,theta(i,j,k-1)+thbase call flush(lp) endif !debug endif c if (lunmix .and. !usually .true. & k.gt.fixlay+1 .and. & theta(i,j,k)-epsil.gt.th3d(i,j,k, n) .and. & theta(i,j,k-1) .lt.th3d(i,j,k, n) .and. & abs(theta(i,j,k-1)- th3d(i,j,k-1,n)).lt.hybiso .and. & ( th3d(i,j,k,n)- th3d(i,j,k-1,n)).gt. & (theta(i,j,k) - theta(i,j,k-1) )*0.001 ) then c c --- water in the deepest inflated layer with significant thickness c --- (kp) is too light, with the layer above near-isopycnal c --- c --- split layer into 2 sublayers, one near the desired density c --- and one exactly matching the T&S properties of layer k-1. c --- To prevent "runaway" T or S, the result satisfies either c --- abs(T.k - T.k-1) <= abs(T.k-2 - T.k-1) or c --- abs(S.k - S.k-1) <= abs(S.k-2 - S.k-1) c --- It is also limited to a 50% change in layer thickness. c if (i.eq.itest .and. j.eq.jtest) then write(lp,'(a,i3)') & 'hybgen, deepest inflated layer too light (stable):',k call flush(lp) endif !debug c delsm=abs(saln(i,j,k-2,n)-saln(i,j,k-1,n)) dels =abs(saln(i,j,k-1,n)-saln(i,j,k, n)) deltm=abs(temp(i,j,k-2,n)-temp(i,j,k-1,n)) delt =abs(temp(i,j,k-1,n)-temp(i,j,k, n)) c --- sanity check on deltm and delsm q=min(temp(i,j,k-2,n),temp(i,j,k-1,n),temp(i,j,k,n)) if (q.gt. 6.0) then deltm=min( deltm, 6.0*(theta(i,j,k)-theta(i,j,k-1)) ) else !(q.le. 6.0) deltm=min( deltm, 10.0*(theta(i,j,k)-theta(i,j,k-1)) ) endif delsm=min( delsm, 1.3*(theta(i,j,k)-theta(i,j,k-1)) ) qts=0.0 if (delt.gt.epsil) then qts=max(qts, (min(deltm, 2.0*delt)-delt)/delt) ! qts<=1.0 endif if (dels.gt.epsil) then qts=max(qts, (min(delsm, 2.0*dels)-dels)/dels) ! qts<=1.0 endif q=(theta(i,j,k)-th3d(i,j,k, n))/ & (theta(i,j,k)-th3d(i,j,k-1,n)) q=min(q,qts/(1.0+qts)) ! upper sublayer <= 50% of total q=qhrlx(k)*q c --- qhrlx is relaxation coefficient (inverse baroclinic time steps) p_hat=q*(p(i,j,k+1)-p(i,j,k)) p(i,j,k)=p(i,j,k)+p_hat if (hybflg.eq.0) then !T&S temp(i,j,k,n)=temp(i,j,k,n)+(q/(1.0-q))*(temp(i,j,k,n) - & temp(i,j,k-1,n) ) saln(i,j,k,n)=saln(i,j,k,n)+(q/(1.0-q))*(saln(i,j,k,n) - & saln(i,j,k-1,n) ) th3d(i,j,k,n)=sig(temp(i,j,k,n),saln(i,j,k,n))-thbase elseif (hybflg.eq.1) then !th&S th3d(i,j,k,n)=th3d(i,j,k,n)+(q/(1.0-q))*(th3d(i,j,k,n) - & th3d(i,j,k-1,n) ) saln(i,j,k,n)=saln(i,j,k,n)+(q/(1.0-q))*(saln(i,j,k,n) - & saln(i,j,k-1,n) ) temp(i,j,k,n)=tofsig(th3d(i,j,k,n)+thbase,saln(i,j,k,n)) elseif (hybflg.eq.2) then !th&T th3d(i,j,k,n)=th3d(i,j,k,n)+(q/(1.0-q))*(th3d(i,j,k,n) - & th3d(i,j,k-1,n) ) temp(i,j,k,n)=temp(i,j,k,n)+(q/(1.0-q))*(temp(i,j,k,n) - & temp(i,j,k-1,n) ) saln(i,j,k,n)=sofsig(th3d(i,j,k,n)+thbase,temp(i,j,k,n)) endif if (ntracr.gt.0 .and. p_hat.ne.0.0) then qtr=p_hat/(p(i,j,k)-p(i,j,k-1)) !ok because <1.0 and >0.0 do ktr= 1,ntracr if (trcflg(ktr).eq.2) then !temperature tracer tracer(i,j,k,n,ktr)=tracer(i,j,k,n,ktr)+ & (q/(1.0-q))*(tracer(i,j,k, n,ktr)- & tracer(i,j,k-1,n,ktr)) else !standard tracer - not split into two sub-layers tracer(i,j,k-1,n,ktr)=tracer(i,j,k-1,n,ktr)+ & qtr*(tracer(i,j,k, n,ktr)- & tracer(i,j,k-1,n,ktr)) if (i.eq.itest .and. j.eq.jtest) then write(lp,'(a,i4,i3,5e12.3)') & 'hybgen, 10(+):', & k,ktr,p_hat,p(i,j,k),p(i,j,k-1), & qtr,tracer(i,j,k-1,n,ktr) call flush(lp) endif !debug endif !trcflg enddo !ktr endif !tracers if (mxlmy .and. p_hat.ne.0.0) then qtr=p_hat/(p(i,j,k)-p(i,j,k-1)) !ok because <1.0 and >0.0 q2( i,j,k-1,n)=q2( i,j,k-1,n)+ & qtr*(q2( i,j,k,n)-q2( i,j,k-1,n)) q2l(i,j,k-1,n)=q2l(i,j,k-1,n)+ & qtr*(q2l(i,j,k,n)-q2l(i,j,k-1,n)) if (i.eq.itest .and. j.eq.jtest) then write(lp,'(a,i4,i3,6e12.3)') & 'hybgen, 10(+):', & k,0,p_hat,p(i,j,k)-p(i,j,k-1),p(i,j,k+1)-p(i,j,k), & qtr,q2(i,j,k-1,n),q2l(i,j,k-1,n) call flush(lp) endif !debug endif if (i.eq.itest .and. j.eq.jtest) then write(lp,'(a,i3,f6.3,5f8.3)') & 'hybgen, 10(+):', & k,q,temp(i,j,k,n),saln(i,j,k,n), & th3d(i,j,k,n)+thbase,theta(i,j,k)+thbase call flush(lp) endif !debug if (i.eq.itest .and. j.eq.jtest) then write(lp,'(a,i3,f6.3,5f8.3)') & 'hybgen, 10(-):', & k,0.0,temp(i,j,k,n),saln(i,j,k,n), & th3d(i,j,k,n)+thbase,theta(i,j,k)+thbase call flush(lp) endif !debug endif !too light c c --- massless or near-massless (thickness < oneum) layers c --- conserve mass, but make them all the same as deepest significant layer c qsum = 0.0 s1dsum(:) = 0.0 do k=kp,kk q = max( p(i,j,k+1)-p(i,j,k), 0.0 ) !enforce interface order qsum = qsum + q if (hybflg.eq.0) then !T&S s1dsum(1) = s1dsum(1) + temp(i,j,k,n)*q s1dsum(2) = s1dsum(2) + saln(i,j,k,n)*q elseif (hybflg.eq.1) then !th&S s1dsum(1) = s1dsum(1) + th3d(i,j,k,n)*q s1dsum(2) = s1dsum(2) + saln(i,j,k,n)*q elseif (hybflg.eq.2) then !th&T s1dsum(1) = s1dsum(1) + th3d(i,j,k,n)*q s1dsum(2) = s1dsum(2) + temp(i,j,k,n)*q endif do ktr= 1,ntracr s1dsum(2+ktr) = s1dsum(2+ktr) + tracer(i,j,k,n,ktr)*q enddo if (mxlmy) then s1dsum(ntracr+3) = s1dsum(ntracr+3) + q2( i,j,k,n)*q s1dsum(ntracr+4) = s1dsum(ntracr+4) + q2l(i,j,k,n)*q endif enddo !k q = 1.0/qsum !qsum must be > onemu and so > 0 s1dsum(:) = s1dsum(:)*q do k=kp,kk if (k.le.nhybrd) then c --- homogenize near-bottom layers if (hybflg.eq.0) then !T&S temp(i,j,k,n)=s1dsum(1) saln(i,j,k,n)=s1dsum(2) th3d(i,j,k,n)=sig(s1dsum(1),s1dsum(2))-thbase elseif (hybflg.eq.1) then !th&S th3d(i,j,k,n)=s1dsum(1) saln(i,j,k,n)=s1dsum(2) temp(i,j,k,n)=tofsig(s1dsum(1)+thbase, & s1dsum(2)) elseif (hybflg.eq.2) then !th&T th3d(i,j,k,n)=s1dsum(1) temp(i,j,k,n)=s1dsum(2) saln(i,j,k,n)=sofsig(s1dsum(1)+thbase, & s1dsum(2)) endif elseif (th3d(i,j,k,n).ne.theta(i,j,k)) then if (hybflg.ne.2) then c --- homogenize saln th3d(i,j,k,n)=max(theta(i,j,k), th3d(i,j,k-1,n)) saln(i,j,k,n)=s1dsum(2) temp(i,j,k,n)=tofsig(th3d(i,j,k,n)+thbase,saln(i,j,k,n)) saln(i,j,k,n)=sofsig(th3d(i,j,k,n)+thbase,temp(i,j,k,n)) else c --- homogenize temp th3d(i,j,k,n)=max(theta(i,j,k), th3d(i,j,k-1,n)) temp(i,j,k,n)=s1dsum(2) saln(i,j,k,n)=sofsig(th3d(i,j,k,n)+thbase,temp(i,j,k,n)) endif endif do ktr= 1,ntracr tracer(i,j,k,n,ktr)=s1dsum(2+ktr) enddo if (mxlmy) then q2 (i,j,k,n)=s1dsum(ntracr+3) q2l(i,j,k,n)=s1dsum(ntracr+4) endif enddo !k c c --- store one-dimensional arrays of -temp-, -saln-, and -p- for the 'old' c --- vertical grid before attempting to restore isopycnic conditions pres(1)=p(i,j,1) do k=1,kk if (hybflg.eq.0) then !T&S s1d(k,1) = temp(i,j,k,n) s1d(k,2) = saln(i,j,k,n) elseif (hybflg.eq.1) then !th&S s1d(k,1) = th3d(i,j,k,n) s1d(k,2) = saln(i,j,k,n) elseif (hybflg.eq.2) then !th&T s1d(k,1) = th3d(i,j,k,n) s1d(k,2) = temp(i,j,k,n) endif do ktr= 1,ntracr s1d(k,2+ktr) = tracer(i,j,k,n,ktr) enddo if (mxlmy) then s1d(k,ntracr+3) = q2( i,j,k,n) s1d(k,ntracr+4) = q2l(i,j,k,n) endif pres(k+1)=p(i,j,k+1) dprs(k) =pres(k+1)-pres(k) dpi( k) =max(dprs(k),onezm) c if (isopcm) then if (k.le.fixlay) then lcm(k) = .false. !fixed layers are never PCM else c --- thin and isopycnal layers remapped with PCM. lcm(k) = k.gt.nhybrd & .or. dprs(k).le.oneum & .or. abs(th3d(i,j,k,n)-theta(i,j,k)).lt.hybiso endif !k<=fixlay:else endif !isopcm enddo !k c c --- try to restore isopycnic conditions by moving layer interfaces c --- qhrlx(k) are relaxation coefficients (inverse baroclinic time steps) c if (fixlay.ge.1) then c c --- maintain constant thickness, layer k = 1 k=1 p_hat=p(i,j,k)+dp0ij(k) p(i,j,k+1)=min(p_hat,p(i,j,k+2)) endif c do k=2,nhybrd c if (i.eq.itest .and. j.eq.jtest) then write(cinfo,'(a9,i2.2,1x)') ' do 88 k=',k 109 format (i9,2i5,a,a/(i9,8x,a,a,i3,f9.2,f8.2,f9.2,f8.2)) write(lp,109) nstep,itest+i0,jtest+j0, . cinfo,': othkns odpth nthkns ndpth', . (nstep,cinfo,':',ka, . (pres(ka+1)- . pres(ka) )*qonem, . pres(ka+1) *qonem, . (p(itest,jtest,ka+1)- . p(itest,jtest,ka) )*qonem, . p(itest,jtest,ka+1) *qonem,ka=1,kk) call flush(lp) endif !debug c if (k.le.fixlay) then c c --- maintain constant thickness, k <= fixlay if (k.lt.kk) then !p.kk+1 not changed p(i,j,k+1)=min(dp0cum(k+1),p(i,j,kk+1)) if (k.eq.fixlay) then c --- enforce interface order (may not be necessary). do ka= k+2,kk if (p(i,j,ka).ge.p(i,j,k+1)) then exit ! usually get here quickly else p(i,j,ka) = p(i,j,k+1) endif enddo !ka endif !k.eq.fixlay endif !k.lt.kk c if (i.eq.itest .and. j.eq.jtest) then write(lp,'(a,i3.2,f8.2)') 'hybgen, fixlay :', & k+1,p(i,j,k+1)*qonem call flush(lp) endif !debug else c c --- do not maintain constant thickness, k > fixlay c if (th3d(i,j,k,n).gt.theta(i,j,k)+epsil .and. & k.gt.fixlay+1) then c c --- water in layer k is too dense c --- try to dilute with water from layer k-1 c --- do not move interface if k = fixlay + 1 c if (th3d(i,j,k-1,n).ge.theta(i,j,k-1) .or. & p(i,j,k).le.dp0cum(k)+onem .or. & p(i,j,k+1)-p(i,j,k).le.p(i,j,k)-p(i,j,k-1)) then c c --- if layer k-1 is too light, thicken the thinner of the two, c --- i.e. skip this layer if it is thicker. c if (i.eq.itest .and. j.eq.jtest) then write(lp,'(a,3x,i2.2,1pe13.5)') & 'hybgen, too dense:',k,th3d(i,j,k,n)-theta(i,j,k) call flush(lp) endif !debug c if ((theta(i,j,k)-th3d(i,j,k-1,n)).le.epsil) then c layer k-1 too dense, take entire layer p_hat=p(i,j,k-1)+dp0ij(k-1) else q=(theta(i,j,k)-th3d(i,j,k, n))/ & (theta(i,j,k)-th3d(i,j,k-1,n)) ! -1 <= q < 0 p_hat0=p(i,j,k)+q*(p(i,j,k+1)-p(i,j,k)) !
p(i,j,k+1) endif c c --- if layer k+1 does not touch the bottom then maintain minimum c --- thicknesses of layers k and k+1 as much as possible, c --- but permit layers to collapse to zero thickness at the bottom c if (p(i,j,k+2).lt.p(i,j,kk+1)) then if (p(i,j,kk+1)-p(i,j,k).gt. & dp0ij(k)+dp0ij(k+1) ) then p_hat=p(i,j,k+2)-cushn(p(i,j,k+2)-p_hat,dp0ij(k+1)) endif p_hat=p(i,j,k)+max(p_hat-p(i,j,k),dp0ij(k)) p_hat=min(p_hat, & max(0.5*(p(i,j,k+1)+p(i,j,k+2)), & p(i,j,k+2)-dp0ij(k+1))) else p_hat=min(p(i,j,k+2),p_hat) endif !p.k+2