PROGRAM PROFILEI IMPLICIT NONE C C hycom_profile2zi - Usage: hycom_profile2zi archv.txt zi.txt archz.txt [itype] C C converts an HYCOM isopycnal text profile file to C a finite volume z-level text profile file. C C archv.txt is assumed to be an HYCOM archive text profile file C zi.txt is assumed to contain a list of interface depths, one per line C the first layer is assumed to be from the surface. C the last layer need not reach the bottom, so set the C last interface depth very deep to guarentee full coverage. C z-levels are nominally half way between interface depths. C archz.txt will be the output, new cell, text profile file C itype is the input interpolation type (default 1) C =0; piecewise constant method (PCM) or donor cell C =1; piecewise linear method (PLM) or VanLeer C =2; piecewise parabolic method (PPM) C =3; weighted essentially non-oscillatory (WENO) C C all input interpolation schemes conserve the original cell average, C and the output is the average of the interpolation profile across C each new cell. So the total depth average is conserved (provided C the last interface is at or below the total depth). C C If some or all of the target coordinate system is in density space, C use hycom_profile_remap to both choose the interface depths and C perform the averaging over the new cells. C C this version for "serial" Unix systems. C C Tim Campbell, Mississippi State University C Alan J. Wallcraft, Naval Research Laboratory, Sep. 2002 (July 2008). C INTEGER IARGC INTEGER NARG CHARACTER*240 CARG C CHARACTER*240 CFILEA,CFILEB,CFILEC,CFILED CHARACTER*240 CLINE REAL HYBISO,THK,DUM5(5),FLAG REAL SIGISO(9999),PO(9999) REAL SIT(5),SOT(5) INTEGER IOS,ITYPE,K,KK,KDM,KO INTEGER I,KMAX C REAL, ALLOCATABLE :: SI(:,:),PI(:),SO(:,:),ZO(:) C INTEGER NSAMP REAL PMAX PARAMETER(NSAMP=1,PMAX=500.0) C C READ ARGUMENTS. C NARG = IARGC() C IF (NARG.EQ.4) THEN CALL GETARG(1,CFILEA) CALL GETARG(2,CFILEB) CALL GETARG(3,CFILEC) CALL GETARG(4,CLINE) READ(CLINE,*) ITYPE CFILED = " " HYBISO = -1.0 ELSEIF (NARG.EQ.3) THEN CALL GETARG(1,CFILEA) CALL GETARG(2,CFILEB) CALL GETARG(3,CFILEC) ITYPE = 1 CFILED = " " HYBISO = -1.0 ELSEIF (NARG.EQ.6) THEN !undocumented option CALL GETARG(1,CFILEA) CALL GETARG(2,CFILEB) CALL GETARG(3,CFILEC) CALL GETARG(4,CLINE) READ(CLINE,*) ITYPE CALL GETARG(5,CFILED) CALL GETARG(6,CLINE) READ(CLINE,*) HYBISO ELSE WRITE(6,*) + 'Usage: hycom_profile2zi archv.txt zi.txt archz.txt [itype]' 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 ',TRIM(CFILEA) WRITE(6,*) 'ios = ',ios CALL EXIT(3) ENDIF OPEN(UNIT=12, FILE=CFILEB, FORM='FORMATTED', STATUS='OLD', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',TRIM(CFILEB) WRITE(6,*) 'ios = ',ios CALL EXIT(4) ENDIF IF (CFILED.NE." ") THEN OPEN(UNIT=13, FILE=CFILED, FORM='FORMATTED', STATUS='OLD', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',TRIM(CFILED) WRITE(6,*) 'ios = ',ios CALL EXIT(5) ENDIF ENDIF !cfiled OPEN(UNIT=21, FILE=CFILEC, FORM='FORMATTED', STATUS='NEW', + IOSTAT=IOS) IF (IOS.NE.0) THEN WRITE(6,*) 'Error: can''t open ',TRIM(CFILEC) WRITE(6,*) 'ios = ',ios CALL EXIT(5) ENDIF C C COPY PROFILE HEADER TO OUTPUT. C DO K= 1,4 READ( 11,'(a)') CLINE WRITE(21,'(a)') TRIM(CLINE) ENDDO READ( 11,'(a)') CLINE C C READ THE ISOPYCNAL PROFILE, TO GET KDM. C 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 ENDDO C C RE-READ THE ISOPYCNAL PROFILE. C ALLOCATE( PI(KDM+1), SI(KDM,5) ) C REWIND(11) DO K= 1,5 READ(11,*) ENDDO PI(1) = 0.0 DO K= 1,KDM READ(11,'(a)',IOSTAT=IOS) CLINE IF (IOS.NE.0) THEN WRITE(6,*) 'Error: inconsistent input profile' CALL EXIT(6) ENDIF READ(CLINE,*) KDM,SI(K,1),SI(K,2),SI(K,3),SI(K,4),SI(K,5),THK PI(K+1) = PI(K) + THK IF (THK.EQ.0.0) THEN DO KK= 1,5 SI(K,KK)=SI(K-1,KK) ENDDO !kk ENDIF ENDDO CLOSE(11) C C READ SIGISO C IF (CFILED.NE." ") THEN DO K= 1,KDM READ(13,*) SIGISO(K) ENDDO CLOSE(13) ENDIF C C STORE COPY OF INPUT AVERAGES AND C COMPUTE INPUT TOTAL COLUMN AVERAGES C KMAX=0 DO KK= 1,5 SIT(KK)=0.0 ENDDO !kk DO K=1,KDM IF (PI(K).LE.PMAX) KMAX=KMAX+1 THK=PI(K+1)-PI(K) DO KK= 1,5 SIT(KK)=SIT(KK)+THK*SI(K,KK) ENDDO !kk ENDDO DO KK= 1,5 SIT(KK)=SIT(KK)/PI(KDM+1) ENDDO !kk C C READ THE ZI FILE. C PO(1) = 0.0 DO K= 2,99999 READ(12,'(a)',IOSTAT=IOS) CLINE IF (IOS.NE.0) THEN KO = K-2 EXIT ENDIF READ(CLINE,*) PO(K) IF (K.EQ.2 .AND. PO(K).LE.0.0) THEN WRITE(6,*) 'Error: inconsistent output profile' WRITE(6,*) '1st value must be > 0.0' WRITE(6,*) PO(2) CALL EXIT(7) ELSEIF (PO(K).LT.PO(K-1)) THEN WRITE(6,*) 'Error: inconsistent output profile' WRITE(6,*) k-1,'-th value < previous value' WRITE(6,*) PO(2:K) CALL EXIT(8) ENDIF PO(K) =MIN( PO(K), PI(KDM+1) ) ENDDO CLOSE(12) C ALLOCATE( ZO(KO), SO(KO,5) ) C DO K= 1,KO ZO(K) = 0.5*(PO(K)+PO(K+1)) ENDDO C C FORM NEW CELL AVERAGES C FLAG = 9999.99 IF (ITYPE.EQ.0) THEN CALL LAYER2ZI_PCM(SI,PI,KDM,5, & SO,PO,KO, FLAG) ELSEIF (ITYPE.EQ.1) THEN CALL LAYER2ZI_PLM(SI,PI,KDM,5, & SO,PO,KO, FLAG, SIGISO,HYBISO) ELSEIF (ITYPE.EQ.2) THEN CALL LAYER2ZI_PPM(SI,PI,KDM,5, & SO,PO,KO, FLAG, SIGISO,HYBISO) ELSEIF (ITYPE.EQ.3) THEN CALL LAYER2ZI_WENO(SI,PI,KDM,5, & SO,PO,KO, FLAG, SIGISO,HYBISO) ELSE WRITE(6,*) + 'Usage: hycom_profile2zi archv.txt zi.txt archz.txt [itype]' WRITE(6,*) 'unsupported itype value' CALL EXIT(1) ENDIF C C COMPUTE AND PRINT NEW TOTAL COLUMN AVERAGE AND DEVIATION METRIC C DO KK= 1,5 SOT(KK)=0.0 DO K=1,KO THK=PO(K+1)-PO(K) SOT(KK)=SOT(KK)+THK*SO(K,KK) ENDDO !k SOT(KK)=SOT(KK)/PO(KO+1) ENDDO !kk WRITE(*,'(/A)') ' LABEL: U V T S R' WRITE(*,'(A,5e14.6)') ' INPUT TOTALS: ',SIT(:) WRITE(*,'(A,5e14.6)') 'OUTPUT TOTALS: ',SOT(:) C C OUTPUT. C WRITE(21,'(a,a)') & '# k', & ' utot vtot temp saln dens thkns dpth' DO K= 1,KO WRITE(21,'(i4,1x,2f8.2,3f8.3,f9.3,f10.3)') & K, & SO(K,1), !cm/s & SO(K,2), !cm/s & SO(K,3), !degC & SO(K,4), !psu & SO(K,5), !SigmaT & PO(K+1)-PO(K), !m & ZO(K) !m ENDDO END subroutine layer2zi_pcm(si,pi,ki,ks, & so,po,ko, flag) implicit none c integer ki,ks,ko real si(ki,ks),pi(ki+1), & so(ko,ks),po(ko+1),flag c c********** c* c 1) remap from one set of vertical cells to another. c method: piecewise constant across each input cell c the output is the average of the interpolation c profile across each output cell. c c 2) input arguments: c si - scalar fields in pi-layer space c pi - layer interface depths (non-negative m) c pi( 1) is the surface c pi(ki+1) is the bathymetry c ki - 1st dimension of si (number of input layers) c ks - 2nd dimension of si,so (number of fields) c po - target interface depths (non-negative m) c po(k+1) >= po(k) c ko - 1st dimension of so (number of output layers) c flag - data void (land) marker c c 3) output arguments: c so - scalar fields in po-layer space c c 4) except at data voids, must have: c pi( 1) == zero (surface) c pi( l+1) >= pi(l) c pi(ki+1) == bathymetry c 0 <= po(k) <= po(k+1) c output layers completely below the bathymetry inherit values c from the layer above. c c 5) Alan J. Wallcraft, Naval Research Laboratory, Sep. 2002 (Aug. 2005). c* c********** c real thin parameter (thin=1.e-6) ! minimum layer thickness c integer i,k real dpi(ki) c if (si(1,1).eq.flag) then do k= 1,ko do i= 1,ks so(k,i) = flag !land enddo !i enddo !k else do k= 1,ki dpi(k) = max( pi(k+1) - pi(k), thin ) enddo !k call hybgen_pcm_remap(si,pi,dpi, & so,po,ki,ko,ks,thin) endif return end subroutine layer2zi_pcm subroutine layer2zi_plm(si,pi,ki,ks, & so,po,ko, flag, sigiso,hybiso) implicit none c integer ki,ks,ko real si(ki,ks),pi(ki+1), & so(ko,ks),po(ko+1),flag, & sigiso(ki),hybiso c c********** c* c 1) remap from one set of vertical cells to another. c method: piecewise linear across each input cell c the output is the average of the interpolation c profile across each output cell. c c 2) input arguments: c si - scalar fields in pi-layer space c pi - layer interface depths (non-negative m) c pi( 1) is the surface c pi(ki+1) is the bathymetry c ki - 1st dimension of si (number of input layers) c ks - 2nd dimension of si,so (number of fields) c po - target interface depths (non-negative m) c po(k+1) >= po(k) c ko - 1st dimension of so (number of output layers) c flag - data void (land) marker c sigiso - isopycnal target densities c hybiso - Use PCM if layer is within hybiso of target density c <0.0 to never use PCM c c 3) output arguments: c so - scalar fields in po-layer space c c 4) except at data voids, must have: c pi( 1) == zero (surface) c pi( l+1) >= pi(l) c pi(ki+1) == bathymetry c 0 <= po(k) <= po(k+1) c output layers completely below the bathymetry inherit values c from the layer above. c c 5) Tim Campbell, Mississippi State University, October 2002. C Alan J. Wallcraft, Naval Research Laboratory, Aug. 2005. c* c********** c real thin parameter (thin=1.e-6) ! minimum layer thickness c integer i,k logical lcm(ki) real dpi(ki),c1d(ki,ks,3) c if (si(1,1).eq.flag) then do k= 1,ko do i= 1,ks so(k,i) = flag !land enddo !i enddo !k else do k= 1,ki dpi(k) = max( pi(k+1) - pi(k), thin ) * lcm(k) = .false. !no PCM lcm(k) = hybiso.gt.0.0 .and. & abs(si(k,5)-sigiso(k)).lt.hybiso enddo !k call hybgen_plm_coefs(si, dpi,lcm,c1d, & ki, ks,thin) call hybgen_plm_remap(si,pi,dpi, c1d, & so,po,ki,ko,ks,thin) endif return end subroutine layer2zi_plm subroutine layer2zi_ppm(si,pi,ki,ks, & so,po,ko, flag, sigiso,hybiso) implicit none c integer ki,ks,ko real si(ki,ks),pi(ki+1), & so(ko,ks),po(ko+1),flag, & sigiso(ki),hybiso c c********** c* c 1) remap from one set of vertical cells to another. c method: piecewise parabolic method across each input cell c the output is the average of the interpolation c profile across each output cell. c c 2) input arguments: c si - scalar fields in pi-layer space c pi - layer interface depths (non-negative m) c pi( 1) is the surface c pi(ki+1) is the bathymetry c ki - 1st dimension of si (number of input layers) c ks - 2nd dimension of si,so (number of fields) c po - target interface depths (non-negative m) c po(k+1) >= po(k) c ko - 1st dimension of so (number of output layers) c flag - data void (land) marker c sigiso - isopycnal target densities c hybiso - Use PCM if layer is within hybiso of target density c <0.0 to never use PCM c c 3) output arguments: c so - scalar fields in po-layer space c c 4) except at data voids, must have: c pi( 1) == zero (surface) c pi( l+1) >= pi(l) c pi(ki+1) == bathymetry c 0 <= po(k) <= po(k+1) c output layers completely below the bathymetry inherit values c from the layer above. c c 5) Tim Campbell, Mississippi State University, October 2002. C Alan J. Wallcraft, Naval Research Laboratory, Aug. 2005. c* c********** c real thin parameter (thin=1.e-6) ! minimum layer thickness c integer i,k logical lcm(ki) real dpi(ki),c1d(ki,ks,3) c if (si(1,1).eq.flag) then do k= 1,ko do i= 1,ks so(k,i) = flag !land enddo !i enddo !k else do k= 1,ki dpi(k) = max( pi(k+1) - pi(k), thin ) * lcm(k) = .false. !no PCM lcm(k) = hybiso.gt.0.0 .and. & abs(si(k,5)-sigiso(k)).lt.hybiso enddo !k call hybgen_ppm_coefs(si, dpi,lcm,c1d, & ki, ks,thin) call hybgen_ppm_remap(si,pi,dpi, c1d, & so,po,ki,ko,ks,thin) endif return end subroutine layer2zi_ppm subroutine layer2zi_weno(si,pi,ki,ks, & so,po,ko, flag, sigiso,hybiso) implicit none c integer ki,ks,ko real si(ki,ks),pi(ki+1), & so(ko,ks),po(ko+1),flag, & sigiso(ki),hybiso c c********** c* c 1) remap from one set of vertical cells to another. c method: WENO c c 2) input arguments: c si - scalar fields in pi-layer space c pi - layer interface depths (non-negative m) c pi( 1) is the surface c pi(ki+1) is the bathymetry c ki - 1st dimension of si (number of input layers) c ks - 2nd dimension of si,so (number of fields) c po - target interface depths (non-negative m) c po(k+1) >= po(k) c ko - 1st dimension of so (number of output layers) c flag - data void (land) marker c sigiso - isopycnal target densities c hybiso - Use PCM if layer is within hybiso of target density c <0.0 to never use PCM c c 3) output arguments: c so - scalar fields in po-layer space c c 4) except at data voids, must have: c pi( 1) == zero (surface) c pi( l+1) >= pi(l) c pi(ki+1) == bathymetry c 0 <= po(k) <= po(k+1) c output layers completely below the bathymetry inherit values c from the layer above. c c 5) Tim Campbell, Mississippi State University, October 2002. C Alan J. Wallcraft, Naval Research Laboratory, Aug. 2005. c* c********** c real thin parameter (thin=1.e-6) ! minimum layer thickness c integer i,k logical lcm(ki) real dpi(ki),c1d(ki,ks,2) c if (si(1,1).eq.flag) then do k= 1,ko do i= 1,ks so(k,i) = flag !land enddo !i enddo !k else do k= 1,ki dpi(k) = max( pi(k+1) - pi(k), thin ) * lcm(k) = .false. !no PCM lcm(k) = hybiso.gt.0.0 .and. & abs(si(k,5)-sigiso(k)).lt.hybiso enddo !k call hybgen_weno_coefs(si, dpi,lcm,c1d, & ki, ks,thin) call hybgen_weno_remap(si,pi,dpi, c1d, & so,po,ki,ko,ks,thin) endif return end subroutine layer2zi_weno subroutine hybgen_pcm_remap(si,pi,dpi, & so,po,ki,ko,ks,thin) implicit none c integer ki,ko,ks real si(ki,ks),pi(ki+1),dpi(ki), & so(ko,ks),po(ko+1),thin c c----------------------------------------------------------------------- c 1) remap from one set of vertical cells to another. c method: piecewise constant across each input cell c the output is the average of the interpolation c profile across each output cell. c c 2) input arguments: c si - initial scalar fields in pi-layer space c pi - initial layer interface depths (non-negative) c pi( 1) is the surface c pi(ki+1) is the bathymetry c pi(k+1) >= pi(k) c dpi - initial layer thicknesses (dpi(k)=pi(k+1)-pi(k)) c ki - number input of layers c ko - number output of layers c ks - number of fields c po - target interface depths (non-negative) c po( 1) is the surface c po(ko+1) is the bathymetry (== pi(ki+1)) c po(k+1) >= po(k) c thin - layer thickness (>0) that can be ignored c c 3) output arguments: c so - scalar fields in po-layer space c c 4) Tim Campbell, Mississippi State University, October 2002. c Alan J. Wallcraft, Naval Research Laboratory, Aug. 2007. c----------------------------------------------------------------------- c integer i,k,l,lb,lt real dpb,dpt,xb,xt,zb,zt,zx,o real*8 sz c zx=pi(ki+1) !maximum depth zb=max(po(1),pi(1)) lb=1 do while (pi(lb+1).lt.zb .and. lb.lt.ki) lb=lb+1 enddo do k= 1,ko !output layers zt = zb zb = min(po(k+1),zx) * write(lp,*) 'k,zt,zb = ',k,zt,zb lt=lb !top will always correspond to bottom of previous !find input layer containing bottom output interface do while (pi(lb+1).lt.zb .and. lb.lt.ki) lb=lb+1 enddo if (zb-zt.le.thin .or. zt.ge.zx) then if (k.ne.1) then c c --- thin or bottomed layer, values taken from layer above c do i= 1,ks so(k,i) = so(k-1,i) enddo !i else !thin surface layer do i= 1,ks so(k,i) = si(k,i) enddo !i endif else c c form layer averages. c use PPM-like logic (may not have minimum operation count) c * if (pi(lb).gt.zt) then * write(lp,*) 'bad lb = ',lb * stop * endif if (lt.ne.lb) then !multiple layers xt=(zt-pi(lt))/max(dpi(lt),thin) xb=(zb-pi(lb))/max(dpi(lb),thin) dpt=pi(lt+1)-zt dpb=zb-pi(lb) do i= 1,ks o = si((lt+lb)/2,i) !offset to reduce round-off sz = dpt*(si(lt,i)-o) do l=lt+1,lb-1 sz = sz+dpi(l)*(si(l,i)-o) enddo !l sz = sz+dpb*(si(lb,i)-o) so(k,i) = o + sz/(zb-zt) !zb-zt>=thin enddo !i else !single layer do i= 1,ks so(k,i) = si(lt,i) enddo !i endif endif !thin:std layer enddo !k return end subroutine hybgen_pcm_remap subroutine hybgen_plm_coefs(si,dpi,lc,ci,kk,ks,thin) implicit none c integer kk,ks logical lc(kk) real si(kk,ks),dpi(kk),ci(kk,ks),thin c c----------------------------------------------------------------------- c 1) coefficents for remaping from one set of vertical cells to another. c method: piecewise linear across each input cell with c monotonized central-difference limiter. c the output is the average of the interpolation c profile across each output cell. c c 2) input arguments: c si - initial scalar fields in pi-layer space c dpi - initial layer thicknesses (dpi(k)=pi(k+1)-pi(k)) c lc - use PCM for selected layers c kk - number of layers c ks - number of fields c thin - layer thickness (>0) that can be ignored c c 3) output arguments: c ci - coefficents (slopes) for hybgen_plm_remap c profile(y)=si+ci*(y-1), 0<=y<=1 c c 4) Tim Campbell, Mississippi State University, October 2002. c Alan J. Wallcraft, Naval Research Laboratory, Aug. 2007. c----------------------------------------------------------------------- c integer k,i real qcen,zbot,zcen,ztop c do i= 1,ks ci(1, i) = 0.0 ci(kk,i) = 0.0 enddo !i do k= 2,kk-1 if (lc(k) .or. dpi(k).le.thin) then !use PCM do i= 1,ks ci(k,i) = 0.0 enddo !i else c --- use qcen in place of 0.5 to allow for non-uniform grid qcen = dpi(k)/(dpi(k)+0.5*(dpi(k-1)+dpi(k+1))) !dpi(k)>thin do i= 1,ks c --- PLM (non-zero slope, but no new extrema) c --- layer value is si-0.5*ci at top interface, c --- and si+0.5*ci at bottom interface. c c --- monotonized central-difference limiter (van Leer, 1977, c --- JCP 23 pp 276-299). For a discussion of PLM limiters, see c --- Finite Volume Methods for Hyperbolic Problems by R.J. Leveque. ztop = 2.0*(si(k, i)-si(k-1,i)) zbot = 2.0*(si(k+1,i)-si(k, i)) zcen =qcen*(si(k+1,i)-si(k-1,i)) if (ztop*zbot.gt.0.0) then !ztop,zbot are the same sign ci(k,i)=sign(min(abs(zcen),abs(zbot),abs(ztop)),zbot) else ci(k,i)=0.0 !local extrema, so no slope endif enddo !i endif !PCM:PLM enddo !k return end subroutine hybgen_plm_coefs subroutine hybgen_plm_remap(si,pi,dpi,ci, & so,po,ki,ko,ks,thin) implicit none c integer ki,ko,ks real si(ki,ks),pi(ki+1),dpi(ki),ci(ki,ks), & so(ko,ks),po(ko+1),thin c c----------------------------------------------------------------------- c 1) remap from one set of vertical cells to another. c method: piecewise linear across each input cell c the output is the average of the interpolation c profile across each output cell. c c 2) input arguments: c si - initial scalar fields in pi-layer space c pi - initial layer interface depths (non-negative) c pi( 1) is the surface c pi(ki+1) is the bathymetry c pi(k+1) >= pi(k) c dpi - initial layer thicknesses (dpi(k)=pi(k+1)-pi(k)) c ci - coefficents (slopes) from hybgen_plm_coefs c profile(y)=si+ci*(y-1), 0<=y<=1 c ki - number input of layers c ko - number output of layers c ks - number of fields c po - target interface depths (non-negative) c po( 1) is the surface c po(ko+1) is the bathymetry (== pi(ki+1)) c po(k+1) >= po(k) c thin - layer thickness (>0) that can be ignored c c 3) output arguments: c so - scalar fields in po-layer space c c 4) Tim Campbell, Mississippi State University, October 2002. c Alan J. Wallcraft, Naval Research Laboratory, Aug. 2007. c----------------------------------------------------------------------- c integer i,k,l,lb,lt real c0,qb0,qb1,qt0,qt1,xb,xt,zb,zt,zx,o real*8 sz c zx=pi(ki+1) !maximum depth zb=max(po(1),pi(1)) lb=1 do while (pi(lb+1).lt.zb .and. lb.lt.ki) lb=lb+1 enddo do k= 1,ko !output layers zt = zb zb = min(po(k+1),zx) * write(lp,*) 'k,zt,zb = ',k,zt,zb lt=lb !top will always correspond to bottom of previous !find input layer containing bottom output interface do while (pi(lb+1).lt.zb .and. lb.lt.ki) lb=lb+1 enddo if (zb-zt.le.thin .or. zt.ge.zx) then if (k.ne.1) then c c --- thin or bottomed layer, values taken from layer above c do i= 1,ks so(k,i) = so(k-1,i) enddo !i else !thin surface layer do i= 1,ks so(k,i) = si(k,i) enddo !i endif else c c form layer averages. c use PPM-like logic (may not have minimum operation count) c * if (pi(lb).gt.zt) then * write(lp,*) 'bad lb = ',lb * stop * endif xt=(zt-pi(lt))/max(dpi(lt),thin) xb=(zb-pi(lb))/max(dpi(lb),thin) if (lt.ne.lb) then !multiple layers qt0 = (1.0-xt) qt1 = (1.0-xt**2)*0.5 qb0 = xb qb1 = xb**2 *0.5 do i= 1,ks o = si((lt+lb)/2,i) !offset to reduce round-off c0 = si(lt,i) - o - 0.5*ci(lt,i) sz= dpi(lt)*(c0*qt0 + ci(lt,i)*qt1) do l=lt+1,lb-1 sz = sz+dpi(l)*(si(l,i) - o) enddo !l c0 = si(lb,i) - o - 0.5*ci(lb,i) sz = sz+dpi(lb)*(c0*qb0 + ci(lb,i)*qb1) so(k,i) = o + sz/(zb-zt) !zb-zt>=thin enddo !i else !single layer qt1 = (xb**2-xt**2 - (xb-xt))*0.5 do i= 1,ks sz = dpi(lt)*(ci(lt,i)*qt1) so(k,i) = si(lt,i) + sz/(zb-zt) !zb-zt>=thin enddo !i endif endif !thin:std layer enddo !k return end subroutine hybgen_plm_remap subroutine hybgen_ppm_coefs(s,dp,lc,ci,kk,ks,thin) implicit none c integer kk,ks logical lc(kk) real s(kk,ks),dp(kk),ci(kk,ks,3),thin c c----------------------------------------------------------------------- c 1) coefficents for remaping from one set of vertical cells to another. c method: monotonic piecewise parabolic across each input cell c the output is the average of the interpolation c profile across each output cell. c c Colella, P. & P.R. Woodward, 1984, J. Comp. Phys., 54, 174-201. c c 2) input arguments: c s - initial scalar fields in pi-layer space c dp - initial layer thicknesses (>=thin) c lc - use PCM for selected layers c kk - number of layers c ks - number of fields c thin - layer thickness (>0) that can be ignored c c 3) output arguments: c ci - coefficents for hybgen_ppm_remap c profile(y)=ci.1+ci.2*y+ci.3*y^2, 0<=y<=1 c c 4) Tim Campbell, Mississippi State University, October 2002. c Alan J. Wallcraft, Naval Research Laboratory, Aug. 2007. c----------------------------------------------------------------------- c integer j,i real da,a6,slj,scj,srj real as(kk),al(kk),ar(kk) real dpjp(kk), dp2jp(kk), dpj2p(kk), & qdpjp(kk),qdp2jp(kk),qdpj2p(kk),dpq3(kk),qdp4(kk) c !compute grid metrics do j=1,kk-1 dpjp( j) = dp(j) + dp(j+1) dp2jp(j) = dp(j) + dpjp(j) dpj2p(j) = dpjp(j) + dp(j+1) qdpjp( j) = 1.0/dpjp( j) qdp2jp(j) = 1.0/dp2jp(j) qdpj2p(j) = 1.0/dpj2p(j) enddo !j dpq3(2) = dp(2)/(dp(1)+dpjp(2)) do j=3,kk-1 dpq3(j) = dp(j)/(dp(j-1)+dpjp(j)) !dp(j)/ (dp(j-1)+dp(j)+dp(j+1)) qdp4(j) = 1.0/(dpjp(j-2)+dpjp(j)) !1.0/(dp(j-2)+dp(j-1)+dp(j)+dp(j+1)) enddo !j c do i= 1,ks !Compute average slopes: Colella, Eq. (1.8) as(1)=0. do j=2,kk-1 if (lc(j) .or. dp(j).le.thin) then !use PCM as(j) = 0.0 else slj=s(j, i)-s(j-1,i) srj=s(j+1,i)-s(j, i) if (slj*srj.gt.0.) then scj=dpq3(j)*( dp2jp(j-1)*srj*qdpjp(j) & +dpj2p(j) *slj*qdpjp(j-1) ) as(j)=sign(min(abs(2.0*slj),abs(scj),abs(2.0*srj)),scj) else as(j)=0. endif endif !PCM:PPM enddo !j as(kk)=0. !Compute "first guess" edge values: Colella, Eq. (1.6) al(1)=s(1,i) !1st layer PCM ar(1)=s(1,i) !1st layer PCM al(2)=s(1,i) !1st layer PCM do j=3,kk-1 al(j)=s(j-1,i)+dp(j-1)*(s(j,i)-s(j-1,i))*qdpjp(j-1) & +qdp4(j)*( & 2.*dp(j)*dp(j-1)*qdpjp(j-1)*(s(j,i)-s(j-1,i))* & ( dpjp(j-2)*qdp2jp(j-1) & -dpjp(j) *qdpj2p(j-1) ) & -dp(j-1)*as(j) *dpjp(j-2)*qdp2jp(j-1) & +dp(j) *as(j-1)*dpjp(j) *qdpj2p(j-1) & ) ar(j-1)=al(j) enddo !j ar(kk-1)=s(kk,i) !last layer PCM al(kk) =s(kk,i) !last layer PCM ar(kk) =s(kk,i) !last layer PCM !Impose monotonicity: Colella, Eq. (1.10) do j=2,kk-1 if (lc(j) .or. dp(j).le.thin) then !use PCM al(j)=s(j,i) ar(j)=s(j,i) elseif ((s(j+1,i)-s(j,i))*(s(j,i)-s(j-1,i)).le.0.) then !local extremum al(j)=s(j,i) ar(j)=s(j,i) else da=ar(j)-al(j) a6=6.0*s(j,i)-3.0*(al(j)+ar(j)) if (da*a6 .gt. da*da) then !peak in right half of zone al(j)=3.0*s(j,i)-2.0*ar(j) elseif (da*a6 .lt. -da*da) then !peak in left half of zone ar(j)=3.0*s(j,i)-2.0*al(j) endif endif enddo !j !Set coefficients do j=1,kk if (al(j).ne.ar(j)) then ci(j,i,1)=al(j) ci(j,i,2)=ar(j)-al(j) ci(j,i,3)=6.0*s(j,i)-3.0*(al(j)+ar(j)) else !PCM ci(j,i,1)=al(j) ci(j,i,2)=0.0 ci(j,i,3)=0.0 endif enddo !j enddo !i return end subroutine hybgen_ppm_coefs subroutine hybgen_ppm_remap(si,pi,dpi,ci, & so,po,ki,ko,ks,thin) implicit none c integer ki,ko,ks real si(ki,ks),pi(ki+1),dpi(ki),ci(ki,ks,3), & so(ko,ks),po(ko+1),thin c c----------------------------------------------------------------------- c 1) remap from one set of vertical cells to another. c method: monotonic piecewise parabolic across each input cell c the output is the average of the interpolation c profile across each output cell. c Colella, P. & P.R. Woodward, 1984, J. Comp. Phys., 54, 174-201. c c 2) input arguments: c si - initial scalar fields in pi-layer space c pi - initial layer interface depths (non-negative) c pi( 1) is the surface c pi(ki+1) is the bathymetry c pi(k+1) >= pi(k) c dpi - initial layer thicknesses (dpi(k)=pi(k+1)-pi(k)) c ci - coefficents from hybgen_ppm_coefs c profile(y)=ci.1+ci.2*y+ci.3*y^2, 0<=y<=1 c ki - number input of layers c ko - number output of layers c ks - number of fields c po - target interface depths (non-negative) c po( 1) is the surface c po(ko+1) is the bathymetry (== pi(ki+1)) c po(k+1) >= po(k) c thin - layer thickness (>0) that can be ignored c c 3) output arguments: c so - scalar fields in po-layer space c c 4) Tim Campbell, Mississippi State University, October 2002. c Alan J. Wallcraft, Naval Research Laboratory, Aug. 2007. c----------------------------------------------------------------------- c integer i,k,l,lb,lt real qb0,qb1,qb2,qt0,qt1,qt2,xb,xt,zb,zt,zx,o real*8 sz c zx=pi(ki+1) !maximum depth zb=max(po(1),pi(1)) lb=1 do while (pi(lb+1).lt.zb .and. lb.lt.ki) lb=lb+1 enddo do k= 1,ko !output layers zt = zb zb = min(po(k+1),zx) * write(lp,*) 'k,zt,zb = ',k,zt,zb lt=lb !top will always correspond to bottom of previous !find input layer containing bottom output interface do while (pi(lb+1).lt.zb .and. lb.lt.ki) lb=lb+1 enddo if (zb-zt.le.thin .or. zt.ge.zx) then if (k.ne.1) then c c --- thin or bottomed layer, values taken from layer above c do i= 1,ks so(k,i) = so(k-1,i) enddo !i else !thin surface layer do i= 1,ks so(k,i) = si(k,i) enddo !i endif else c c form layer averages. c * if (pi(lb).gt.zt) then * write(lp,*) 'bad lb = ',lb * stop * endif xt=(zt-pi(lt))/max(dpi(lt),thin) xb=(zb-pi(lb))/max(dpi(lb),thin) if (lt.ne.lb) then !multiple layers qt0 = (1.0-xt) qt1 = (1.-xt**2)*0.5 qt2 = (1.-xt**3)/3.0 qb0 = xb qb1 = xb**2 *0.5 qb2 = xb**3 /3.0 do i= 1,ks o = si((lt+lb)/2,i) !offset to reduce round-off sz = dpi(lt)*( (ci(lt,i,1)-o)*qt0 & +(ci(lt,i,2)+ & ci(lt,i,3) ) *qt1 & -ci(lt,i,3) *qt2 ) do l=lt+1,lb-1 sz = sz+dpi(l)*(si(l,i)-o) enddo !l sz = sz+dpi(lb)*( (ci(lb,i,1)-o)*qb0 & +(ci(lb,i,2)+ & ci(lb,i,3) ) *qb1 & -ci(lb,i,3) *qb2 ) so(k,i) = o + sz/(zb-zt) !zb-zt>=thin enddo !i else !single layer qt0 = (xb-xt) qt1 = (xb**2-xt**2)*0.5 qt2 = (xb**3-xt**3)/3.0 do i= 1,ks sz = dpi(lt)*( (ci(lt,i,1)-o)*qt0 & +(ci(lt,i,2)+ & ci(lt,i,3) ) *qt1 & -ci(lt,i,3) *qt2 ) so(k,i) = o + sz/(zb-zt) !zb-zt>=thin enddo !i endif endif !thin:std layer enddo !k return end subroutine hybgen_ppm_remap subroutine hybgen_weno_coefs(s,dp,lc,ci,kk,ks,thin) implicit none c integer kk,ks logical lc(kk) real s(kk,ks),dp(kk),ci(kk,ks,2),thin c c----------------------------------------------------------------------- c 1) coefficents for remaping from one set of vertical cells to another. c method: WENO c c REFERENCE? c c 2) input arguments: c s - initial scalar fields in pi-layer space c dp - initial layer thicknesses (>=thin) c lc - use PCM for selected layers c kk - number of layers c ks - number of fields c thin - layer thickness (>0) that can be ignored c c 3) output arguments: c ci - coefficents for hybgen_weno_remap c ci.1 is value at interface above c ci.2 is value at interface below c profile(y)=ci.1+ci.2*y+ci.3*y^2, 0<=y<=1 c c 4) Laurent Debreu, Grenoble. c Alan J. Wallcraft, Naval Research Laboratory, July 2008. c----------------------------------------------------------------------- c real, parameter :: dsmll=1.0e-8 c integer j,i real q,q01,q02,q001,q002 real qdpjm(kk),qdpjmjp(kk),dpjm2jp(kk) real zw(kk+1,3) !compute grid metrics do j=2,kk-1 qdpjm( j) = 1.0/(dp(j-1) + dp(j)) qdpjmjp(j) = 1.0/(dp(j-1) + dp(j) + dp(j+1)) dpjm2jp(j) = dp(j-1) + 2.0*dp(j) + dp(j+1) enddo !j j=kk qdpjm( j) = 1.0/(dp(j-1) + dp(j)) c do i= 1,ks do j=2,kk zw(j,3) = qdpjm(j)*(s(j,i)-s(j-1,i)) enddo !j j = 1 !PCM first layer ci(j,i,1) = s(j,i) ci(j,i,2) = s(j,i) zw(j, 1) = 0.0 zw(j, 2) = 0.0 do j=2,kk-1 if (lc(j) .or. dp(j).le.thin) then !use PCM ci(j,i,1) = s(j,i) ci(j,i,2) = s(j,i) zw(j, 1) = 0.0 zw(j, 2) = 0.0 else q001 = dp(j)*zw(j+1,3) q002 = dp(j)*zw(j, 3) if (q001*q002 < 0.0) then q001 = 0.0 q002 = 0.0 endif q01 = dpjm2jp(j)*zw(j+1,3) q02 = dpjm2jp(j)*zw(j, 3) if (abs(q001) > abs(q02)) then q001 = q02 endif if (abs(q002) > abs(q01)) then q002 = q01 endif q = (q001-q002)*qdpjmjp(j) q001 = q001-q*dp(j+1) q002 = q002+q*dp(j-1) ci(j,i,2) = s(j,i)+q001 ci(j,i,1) = s(j,i)-q002 zw( j,1) = (2.0*q001-q002)**2 zw( j,2) = (2.0*q002-q001)**2 endif !PCM:WEND enddo !j j = kk !PCM last layer ci(j,i,1) = s(j,i) ci(j,i,2) = s(j,i) zw(j, 1) = 0.0 zw(j, 2) = 0.0 do j=2,kk q002 = max(zw(j-1,2),dsmll) q001 = max(zw(j, 1),dsmll) zw(j,3) = (q001*ci(j-1,i,2)+q002*ci(j,i,1))/(q001+q002) enddo !j zw( 1,3) = 2.0*s( 1,i)-zw( 2,3) !not used? zw(kk+1,3) = 2.0*s(kk,i)-zw(kk,3) !not used? do j=2,kk-1 if (.not.(lc(j) .or. dp(j).le.thin)) then !don't use PCM q01 = zw(j+1,3)-s(j,i) q02 = s(j,i)-zw(j,3) q001 = 2.0*q01 q002 = 2.0*q02 if (q01*q02 < 0.0) then q01 = 0.0 q02 = 0.0 elseif (abs(q01) > abs(q002)) then q01 = q002 elseif (abs(q02) > abs(q001)) then q02 = q001 endif ci(j,i,1) = s(j,i)-q02 ci(j,i,2) = s(j,i)+q01 endif !PCM:WEND enddo !j enddo !i return end subroutine hybgen_weno_coefs subroutine hybgen_weno_remap(si,pi,dpi,ci, & so,po,ki,ko,ks,thin) implicit none c integer ki,ko,ks real si(ki,ks),pi(ki+1),dpi(ki),ci(ki,ks,2), & so(ko,ks),po(ko+1),thin c c----------------------------------------------------------------------- c 1) remap from one set of vertical cells to another. c method: WENO c c REFERENCE? c c 2) input arguments: c si - initial scalar fields in pi-layer space c pi - initial layer interface depths (non-negative) c pi( 1) is the surface c pi(ki+1) is the bathymetry c pi(k+1) >= pi(k) c dpi - initial layer thicknesses (dpi(k)=pi(k+1)-pi(k)) c ci - coefficents from hybgen_weno_coefs c ci.1 is value at interface above c ci.2 is value at interface below c profile(y)=ci.1+ci.2*y+ci.3*y^2, 0<=y<=1 c ki - number input of layers c ko - number output of layers c ks - number of fields c po - target interface depths (non-negative) c po( 1) is the surface c po(ko+1) is the bathymetry (== pi(ki+1)) c po(k+1) >= po(k) c thin - layer thickness (>0) that can be ignored c c 3) output arguments: c so - scalar fields in po-layer space c c 4) Laurent Debreu, Grenoble. c Alan J. Wallcraft, Naval Research Laboratory, Aug. 2007. c----------------------------------------------------------------------- c integer i,k,l,lb,lt real dpb,dpt,qb0,qb1,qb2,qt0,qt1,qt2,xb,xt,zb,zt,zx,o real*8 sz c zx=pi(ki+1) !maximum depth zb=max(po(1),pi(1)) lb=1 do while (pi(lb+1).lt.zb .and. lb.lt.ki) lb=lb+1 enddo do k= 1,ko !output layers zt = zb zb = min(po(k+1),zx) * write(lp,*) 'k,zt,zb = ',k,zt,zb lt=lb !top will always correspond to bottom of previous !find input layer containing bottom output interface do while (pi(lb+1).lt.zb .and. lb.lt.ki) lb=lb+1 enddo if (zb-zt.le.thin .or. zt.ge.zx) then if (k.ne.1) then c c --- thin or bottomed layer, values taken from layer above c do i= 1,ks so(k,i) = so(k-1,i) enddo !i else !thin surface layer do i= 1,ks so(k,i) = si(k,i) enddo !i endif else c c form layer averages. c * if (pi(lb).gt.zt) then * write(lp,*) 'bad lb = ',lb * stop * endif xt=(zt-pi(lt))/max(dpi(lt),thin) xb=(zb-pi(lb))/max(dpi(lb),thin) if (lt.ne.lb) then !multiple layers dpt = pi(lt+1)-zt dpb = zb-pi(lb) qt1 = xt*(xt-1.0) qt2 = qt1+xt qt0 = 1.0-qt1-qt2 qb1 = (xb-1.0)**2 qb2 = qb1-1.0+xb qb0 = 1.0-qb1-qb2 do i= 1,ks o = si((lt+lb)/2,i) !offset to reduce round-off sz = dpt*(qt0*(si(lt,i) -o) + & qt1*(ci(lt,i,1)-o) + & qt2*(ci(lt,i,2)-o) ) do l=lt+1,lb-1 sz = sz+dpi(l)*(si(l,i) - o) enddo !l sz = sz + dpb*(qb0*(si(lb,i) -o) + & qb1*(ci(lb,i,1)-o) + & qb2*(ci(lb,i,2)-o) ) so(k,i) = o + sz/(zb-zt) !zb-zt>=thin enddo !i else !single layer qt1 = xb**2 + xt**2 + xb*xt + 1.0 - 2.0*(xb+xt) qt2 = qt1 - 1.0 + (xb+xt) qt0 = 1.0 - qt1 - qt2 do i= 1,ks sz=qt0*(si(lt,i) -o) + & qt1*(ci(lt,i,1)-o) + & qt2*(ci(lt,i,2)-o) so(k,i) = o + sz enddo !i endif !layers endif !thin:std layer enddo !k return end subroutine hybgen_weno_remap