diff --git a/ROMS/Nonlinear/Sediment/sed_bed2.F b/ROMS/Nonlinear/Sediment/sed_bed2.F
index e68fbd483fd2bd340583f575cc2d1852ba83d618..ef0f1ed18a84ec3e272bc7e892d26b150ca36b66 100755
--- a/ROMS/Nonlinear/Sediment/sed_bed2.F
+++ b/ROMS/Nonlinear/Sediment/sed_bed2.F
@@ -63,9 +63,9 @@
      &                   GRID(ng) % rmask_wet,                          &
 # endif
 # ifdef SEDBIO_COUP
-     &                   GRID(ng) % Hz,                                 &
+     &                  GRID(ng) % Hz,                                  &
 #  ifdef DIAGNOSTICS_BIO
-     &                   DIAGS(ng) % DiaBio2d,                          &
+     &                  DIAGS(ng) % DiaBio2d,                           &
 #   ifdef WET_DRY
      &                   GRID(ng) % rmask_full,                         &
 #   endif
@@ -216,7 +216,7 @@
       real(r8), intent(in) :: bvstr(LBi:UBi,LBj:UBj)
 #  endif
 #  if defined SED_MORPH
-      real(r8), intent(inout):: bed_thick(LBi:UBi,LBj:UBj,3)
+      real(r8), intent(inout):: bed_thick(LBi:UBi,LBj:UBj,2)
 #  endif
       real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
 #  ifdef SUSPLOAD
@@ -330,6 +330,9 @@ minlayer_thick(ng) = newlayer_thick(ng)
       DO i=Istr,Iend
 #  ifdef SEDBIO_COUP
          ndep_thck_tot(i,j)=0.0_r8
+!IF ((i.eq.3).and.(j.eq.3)) THEN
+!write(6,*),"sd conc 0",bed(i,j,1,iboxy)/bed(i,j,1,ithck)/bed(i,j,1,iporo)
+!ENDIF  
 #  endif
          SED_LOOP: DO ised=1,NST  
             dep_mass(i,ised)=0.0_r8
@@ -874,6 +877,12 @@ minlayer_thick(ng) = newlayer_thick(ng)
              cff2=MIN(1.0_r8,MAX(0.0_r8,Hz(i,j,1)/                      &
      &            (Hz(i,j,1)+MAX(bed(i,j,1,ithck)*bed(i,j,1,iporo),     &
      &            0.0_r8))))
+!IF ((i.eq.3).and.(j.eq.3).and.(ised.eq.iboxy)) THEN
+!write(6,*),"w.c.conc 1",t(i,j,1,nnew,iwc)/Hz(i,j,1)
+!write(6,*),"sd conc 1",bed(i,j,1,ised)/bed(i,j,1,ithck)/bed(i,j,1,iporo)
+!write(6,*),"flux",cff1*cff2-t(i,j,1,nnew,iwc)
+!write(6,*),"DiaBio2d 0",DiaBio2d(i,j,idf)
+!ENDIF 
 #   if defined DIAGNOSTICS_BIO
              DiaBio2d(i,j,idf)=DiaBio2d(i,j,idf) +                      &
 #    ifdef WET_DRY
@@ -883,6 +892,12 @@ minlayer_thick(ng) = newlayer_thick(ng)
 #   endif
              t(i,j,1,nnew,iwc)=cff1*cff2
              bed(i,j,1,ised)=cff1*(1.0_r8-cff2)
+!IF ((i.eq.3).and.(j.eq.3).and.(ised.eq.iboxy)) THEN
+!write(6,*),"w.c.conc 2",t(i,j,1,nnew,iwc)/Hz(i,j,1)
+!write(6,*),"sd conc 2",bed(i,j,1,ised)/bed(i,j,1,ithck)/bed(i,j,1,iporo)
+!write(6,*),"flux",cff1*cff2-t(i,j,1,nnew,iwc)
+!write(6,*),"DiaBio2d 2",DiaBio2d(i,j,idf)
+!ENDIF 
 #  else
 !            cff1=dC/dz in mmol/m3/m
 !            cff2=diffusive flux out of seabed in mmol/m2
@@ -902,6 +917,9 @@ minlayer_thick(ng) = newlayer_thick(ng)
              cff4=(t(i,j,1,nnew,iwc)+bed(i,j,1,ised))                   &
      &            *(bed(i,j,1,ithck)*bed(i,j,1,iporo)                   &
      &            /MAX(eps,Hz(i,j,1)+bed(i,j,1,ithck)*bed(i,j,1,iporo)))
+!IF ((i.eq.3).and.(j.eq.3).and.(ised.eq.iboxy)) THEN
+!write(6,*),"DiaBio2d 0",DiaBio2d(i,j,idf)
+!ENDIF 
 #   if defined DIAGNOSTICS_BIO
              DiaBio2d(i,j,idf)=DiaBio2d(i,j,idf) +                      &
 #    ifdef WET_DRY
@@ -909,12 +927,25 @@ minlayer_thick(ng) = newlayer_thick(ng)
 #    endif
      &            cff2
 #   endif
+!IF ((i.eq.3).and.(j.eq.3).and.(ised.eq.iboxy)) THEN
+!write(6,*),"w.c.conc 1",t(i,j,1,nnew,iwc)/Hz(i,j,1)
+!write(6,*),"sd conc 1",bed(i,j,1,ised)/bed(i,j,1,ithck)/bed(i,j,1,iporo)
+!write(6,*) "dC/dz",cff1
+!write(6,*),"flux",cff2
+!write(6,*),"DiaBio2d 1",DiaBio2d(i,j,idf)
+!ENDIF 
              t(i,j,1,nnew,iwc)=t(i,j,1,nnew,iwc) + cff2
              bed(i,j,1,ised)=bed(i,j,1,ised) - cff2
              IF ( ((cff2.gt.eps).and.(t(i,j,1,nnew,iwc).gt.cff3)).or.   &
      &       ((cff2.lt.(eps*-1.0_r8)).and.(bed(i,j,1,ised).gt.cff4)))   &
      &       THEN
 !            diffusion exceeded equilibrium
+!IF ((i.eq.3).and.(j.eq.3).and.(ised.eq.iboxy)) THEN
+!write(6,*) "diffusion exceeded equilibrium"
+!write(6,*),"w.c. conc",cff3/Hz(i,j,1)
+!write(6,*),"seab conc",cff4/(bed(i,j,1,ithck)*bed(i,j,1,iporo))
+!write(6,*),"flux",cff3-t(i,j,1,nnew,iwc)
+!ENDIF
 #   if defined DIAGNOSTICS_BIO
                 DiaBio2d(i,j,idf)=DiaBio2d(i,j,idf) +                   &
 #    ifdef WET_DRY
@@ -925,6 +956,12 @@ minlayer_thick(ng) = newlayer_thick(ng)
                 t(i,j,1,nnew,iwc)=cff3
                 bed(i,j,1,ised)=cff4
              ENDIF
+!IF ((i.eq.3).and.(j.eq.3).and.(ised.eq.iboxy)) THEN
+!write(6,*),"DiaBio2d 2",DiaBio2d(i,j,idf)
+!write(6,*),"w.c.conc 2",t(i,j,1,nnew,iOxyg)/Hz(i,j,1)
+!write(6,*),"sd conc 2",bed(i,j,1,iboxy)/bed(i,j,1,ithck)/bed(i,j,1,iporo)
+!ENDIF
+
 #  endif
           END DO
         END DO
diff --git a/ROMS/Nonlinear/Sediment/sed_biodiff.F b/ROMS/Nonlinear/Sediment/sed_biodiff.F
index 8c6c7f0a0666fda6c986aa79be331d6eff07db6e..991a39150d9150f06f8388d761ce5181f48b8828 100755
--- a/ROMS/Nonlinear/Sediment/sed_biodiff.F
+++ b/ROMS/Nonlinear/Sediment/sed_biodiff.F
@@ -118,7 +118,7 @@
       real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
 #  endif
 #  if defined SED_MORPH
-      real(r8), intent(inout):: bed_thick(LBi:UBi,LBj:UBj,3)
+      real(r8), intent(inout):: bed_thick(LBi:UBi,LBj:UBj,2)
 #  endif
       real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
       real(r8), intent(inout) :: bed(LBi:UBi,LBj:UBj,Nbed,MBEDP)
@@ -321,6 +321,11 @@
 !
 !   ...make working copies 
               DO k=1,Nbed
+!IF ((i.eq.3).and.(j.eq.3).and.(k.eq.1).and.(ised.eq.iboxy)) THEN
+!write(6,*) "Db", Db(k)
+!write(6,*) "Biodiff sd conc 5",bed(i,j,k,ised)/bed(i,j,k,ithck)/        &
+!     &                         bed(i,j,k,iporo)
+!ENDIF
                 cc(k) = bed(i,j,k,ised)/bed(i,j,k,ithck)/               &
      &                  bed(i,j,k,iporo)
                 dd(k)= d(k)
@@ -345,12 +350,18 @@
 !   ...solution stored in cc; copy out
               DO k = 1,Nbed
                 bed(i,j,k,ised)=cc(k)*bed(i,j,k,ithck)*bed(i,j,k,iporo)
+!IF ((i.eq.3).and.(j.eq.3).and.(k.eq.1).and.(ised.eq.iboxy)) THEN
+!write(6,*) "Db", Db(k)
+!write(6,*) "Biodiff sd conc 6",bed(i,j,k,ised)/bed(i,j,k,ithck)/        &
+!     &                         bed(i,j,k,iporo)
+!ENDIF
               ENDDO
             ENDDO
 
 # endif
 ! Recompute bed masses
             DO k=1,Nbed
+!             write(*,*) i,j,k,(bed_frac(i,j,k,ised),ised=1,NST)
 !            debugging: ensure fracs add up to 1
               cff3 = 0.0_r8
               DO ised=1,NST
diff --git a/ROMS/Nonlinear/Sediment/sedbed_mod.h b/ROMS/Nonlinear/Sediment/sedbed_mod.h
index 8ecf59cc5b8ae33ecd678434d4b7c087d347fda3..e54db75fd7c604b4fba67c9b01d9b8ceff7f242b 100755
--- a/ROMS/Nonlinear/Sediment/sedbed_mod.h
+++ b/ROMS/Nonlinear/Sediment/sedbed_mod.h
@@ -45,13 +45,22 @@
 !  ursell_no      Ursell number of the asymmetric wave.                !
 !  RR_asymwave    Velocity skewness parameter of the asymmetric wave.  !
 !  beta_asymwave  Accleration assymetry parameter.                     !
+!  Zr_wbl         Reference height to get near bottom current vel.(m/s)! 
+!  ksd_wbl        Bed roughness from the currents (m).                 !
+!  ustrc_wbl      Current friction vel. (m/s).                         ! 
+!  thck_wbl       Thickness at WBL edge (m).                           ! 
+!  udelta_wbl     Current vel. at wave boundary layer(WBL) edge (m/s). ! 
+!  phic_sgwbl     angle between waves/currents. at user input elevation!
+!                 to get near-bottom current velocity.                 ! 
+!  phi_wc         angle between waves/currents.                        !
+!  fd_wbl         Friction factor at WBL edge (m).                     ! 
 !  ucrest_r       Crest velocity of the asymmetric wave form (m/s).    !
 !  utrough_r      Trough velocity of the asymmetric wave form (m/s).   !
 !  T_crest        Crest time period of the asymmetric wave form (s).   !
 !  T_trough       Trough time period of the asymmetric wave form (s).  !
-# endif
+# endif 
 #endif
-!
+! 
 !  bottom         Exposed sediment layer properties:                   !
 !                   bottom(:,:,isd50) => mean grain diameter           !
 !                   bottom(:,:,idens) => mean grain density            !
@@ -60,25 +69,15 @@
 !                   bottom(:,:,irlen) => ripple length                 !
 !                   bottom(:,:,irhgt) => ripple height                 !
 !                   bottom(:,:,ibwav) => bed wave excursion amplitude  !
-!                   bottom(:,:,izdef) => default bottom roughness      !
-!                   bottom(:,:,izapp) => apparent bottom roughness     !
 !                   bottom(:,:,izNik) => Nikuradse bottom roughness    !
 !                   bottom(:,:,izbio) => biological bottom roughness   !
 !                   bottom(:,:,izbfm) => bed form bottom roughness     !
 !                   bottom(:,:,izbld) => bed load bottom roughness     !
+!                   bottom(:,:,izapp) => apparent bottom roughness     !
 !                   bottom(:,:,izwbl) => wave bottom roughness         !
+!                   bottom(:,:,izdef) => default bottom roughness      !
 !                   bottom(:,:,iactv) => active layer thickness        !
 !                   bottom(:,:,ishgt) => saltation height              !
-!                   bottom(:,:,imaxD) => maximum inundation depth      !
-!                   bottom(:,:,idnet) => Erosion or deposition         !
-!                   bottom(:,:,idtbl) => Thickness of wbl              !
-!                   bottom(:,:,idubl) => Current velocity at wbl       !
-!                   bottom(:,:,idfdw) => Friction factor from currents !
-!                   bottom(:,:,idzrw) => Ref height for near bottom vel!
-!                   bottom(:,:,idksd) => Bed roughness for wbl         !
-!                   bottom(:,:,idusc) => Current friction velocity wbl !
-!                   bottom(:,:,idpcx) => Angle between currents and xi !
-!                   bottom(:,:,idpwc) => Angle between waves / currents!
 #if defined COHESIVE_BED || defined SED_BIODIFF || defined MIXED_BED
 !                   bottom(:,:,idoff) => tau critical offset           !
 !                   bottom(:,:,idslp) => tau critical slope            !
@@ -142,10 +141,18 @@
 #ifdef BEDLOAD
         real(r8), pointer :: bedldu(:,:,:)
         real(r8), pointer :: bedldv(:,:,:)
-# ifdef BEDLOAD_VANDERA
+# ifdef BEDLOAD_VANDERA 
         real(r8), pointer :: ursell_no(:,:)
         real(r8), pointer :: RR_asymwave(:,:)
         real(r8), pointer :: beta_asymwave(:,:)
+        real(r8), pointer :: Zr_wbl(:,:)
+        real(r8), pointer :: ksd_wbl(:,:)
+        real(r8), pointer :: ustrc_wbl(:,:)
+        real(r8), pointer :: thck_wbl(:,:)
+        real(r8), pointer :: udelta_wbl(:,:)
+        real(r8), pointer :: phic_sgwbl(:,:)
+        real(r8), pointer :: phi_wc(:,:)
+        real(r8), pointer :: fd_wbl(:,:)
         real(r8), pointer :: ucrest_r(:,:)
         real(r8), pointer :: utrough_r(:,:)
         real(r8), pointer :: T_crest(:,:)
@@ -259,20 +266,28 @@
 #endif
 #if defined SEDIMENT && defined SED_MORPH
       allocate ( SEDBED(ng) % bed_thick0(LBi:UBi,LBj:UBj) )
-      allocate ( SEDBED(ng) % bed_thick(LBi:UBi,LBj:UBj,1:3) )
+      allocate ( SEDBED(ng) % bed_thick(LBi:UBi,LBj:UBj,1:2) )
 #endif
 #ifdef BEDLOAD
       allocate ( SEDBED(ng) % bedldu(LBi:UBi,LBj:UBj,NST) )
       allocate ( SEDBED(ng) % bedldv(LBi:UBi,LBj:UBj,NST) )
-# ifdef BEDLOAD_VANDERA
+# ifdef BEDLOAD_VANDERA 
       allocate ( SEDBED(ng) % ursell_no(LBi:UBi,LBj:UBj) )
       allocate ( SEDBED(ng) % RR_asymwave(LBi:UBi,LBj:UBj) )
       allocate ( SEDBED(ng) % beta_asymwave(LBi:UBi,LBj:UBj) )
+      allocate ( SEDBED(ng) % Zr_wbl(LBi:UBi,LBj:UBj) )
+      allocate ( SEDBED(ng) % ksd_wbl(LBi:UBi,LBj:UBj) )
+      allocate ( SEDBED(ng) % ustrc_wbl(LBi:UBi,LBj:UBj) )
+      allocate ( SEDBED(ng) % thck_wbl(LBi:UBi,LBj:UBj) )
+      allocate ( SEDBED(ng) % udelta_wbl(LBi:UBi,LBj:UBj) )
+      allocate ( SEDBED(ng) % phic_sgwbl(LBi:UBi,LBj:UBj) )
+      allocate ( SEDBED(ng) % phi_wc(LBi:UBi,LBj:UBj) )
+      allocate ( SEDBED(ng) % fd_wbl(LBi:UBi,LBj:UBj) )
       allocate ( SEDBED(ng) % ucrest_r(LBi:UBi,LBj:UBj) )
       allocate ( SEDBED(ng) % utrough_r(LBi:UBi,LBj:UBj) )
       allocate ( SEDBED(ng) % T_crest(LBi:UBi,LBj:UBj) )
       allocate ( SEDBED(ng) % T_trough(LBi:UBi,LBj:UBj) )
-# endif
+# endif 
 #endif
       allocate ( SEDBED(ng) % bottom(LBi:UBi,LBj:UBj,MBOTP) )
 #if defined SEDIMENT && defined SUSPLOAD
@@ -294,7 +309,7 @@
 # endif
 # if defined SEDIMENT && defined SED_MORPH
       allocate ( SEDBED(ng) % tl_bed_thick0(LBi:UBi,LBj:UBj) )
-      allocate ( SEDBED(ng) % tl_bed_thick(LBi:UBi,LBj:UBj,1:3) )
+      allocate ( SEDBED(ng) % tl_bed_thick(LBi:UBi,LBj:UBj,1:2) )
 # endif
 # ifdef BEDLOAD
       allocate ( SEDBED(ng) % tl_bedldu(LBi:UBi,LBj:UBj,NST) )
@@ -318,7 +333,7 @@
 # endif
 # if defined SEDIMENT && defined SED_MORPH
       allocate ( SEDBED(ng) % ad_bed_thick0(LBi:UBi,LBj:UBj) )
-      allocate ( SEDBED(ng) % ad_bed_thick(LBi:UBi,LBj:UBj,1:3) )
+      allocate ( SEDBED(ng) % ad_bed_thick(LBi:UBi,LBj:UBj,1:2) )
 # endif
 # ifdef BEDLOAD
       allocate ( SEDBED(ng) % ad_bedldu(LBi:UBi,LBj:UBj,NST) )
@@ -449,7 +464,6 @@
             SEDBED(ng) % bed_thick0(i,j) = IniVal
             SEDBED(ng) % bed_thick(i,j,1) = IniVal
             SEDBED(ng) % bed_thick(i,j,2) = IniVal
-            SEDBED(ng) % bed_thick(i,j,3) = IniVal
           END DO
 #endif
 #ifdef BEDLOAD
@@ -459,11 +473,19 @@
               SEDBED(ng) % bedldv(i,j,itrc) = IniVal
             END DO
           END DO
-# ifdef BEDLOAD_VANDERA
+# ifdef BEDLOAD_VANDERA 
           DO i=Imin,Imax
             SEDBED(ng) % ursell_no(i,j)    = IniVal
             SEDBED(ng) % RR_asymwave(i,j)  = IniVal
             SEDBED(ng) % beta_asymwave(i,j)= IniVal
+            SEDBED(ng) % Zr_wbl(i,j)       = IniVal
+            SEDBED(ng) % ksd_wbl(i,j)      = IniVal
+            SEDBED(ng) % ustrc_wbl(i,j)    = IniVal
+            SEDBED(ng) % thck_wbl(i,j)     = IniVal
+            SEDBED(ng) % udelta_wbl(i,j)   = IniVal
+            SEDBED(ng) % phic_sgwbl(i,j)   = IniVal
+            SEDBED(ng) % phi_wc(i,j)       = IniVal
+            SEDBED(ng) % fd_wbl(i,j)       = IniVal
             SEDBED(ng) % ucrest_r(i,j)     = IniVal
             SEDBED(ng) % utrough_r(i,j)    = IniVal
             SEDBED(ng) % T_crest(i,j)      = IniVal
@@ -523,7 +545,6 @@
             SEDBED(ng) % tl_bed_thick0(i,j) = IniVal
             SEDBED(ng) % tl_bed_thick(i,j,1) = IniVal
             SEDBED(ng) % tl_bed_thick(i,j,2) = IniVal
-            SEDBED(ng) % tl_bed_thick(i,j,3) = IniVal
           END DO
 # endif
 # ifdef BEDLOAD
@@ -580,7 +601,6 @@
             SEDBED(ng) % ad_bed_thick0(i,j) = IniVal
             SEDBED(ng) % ad_bed_thick(i,j,1) = IniVal
             SEDBED(ng) % ad_bed_thick(i,j,2) = IniVal
-            SEDBED(ng) % ad_bed_thick(i,j,3) = IniVal
           END DO
 # endif
 # ifdef BEDLOAD
diff --git a/ROMS/Nonlinear/Sediment/sediment_inp.h b/ROMS/Nonlinear/Sediment/sediment_inp.h
index 84de2e8908d57d5bd18012abf13ceeb279ad43d7..75b9fb89e308032003dc1425b3efc7379cb5dd2c 100755
--- a/ROMS/Nonlinear/Sediment/sediment_inp.h
+++ b/ROMS/Nonlinear/Sediment/sediment_inp.h
@@ -65,6 +65,60 @@
           SELECT CASE (TRIM(KeyWord))
             CASE ('Lsediment')
               Npts=load_l(Nval, Cval, Ngrids, Lsediment)
+            CASE ('NEWLAYER_THICK')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+              DO ng=1,Ngrids
+                newlayer_thick(ng)=Rbed(ng)
+              END DO
+            CASE ('MINLAYER_THICK')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+              DO ng=1,Ngrids
+                minlayer_thick(ng)=Rbed(ng)
+              END DO
+#ifdef MIXED_BED
+            CASE ('TRANSC')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+              DO ng=1,Ngrids
+                transC(ng)=Rbed(ng)
+              END DO
+            CASE ('TRANSN')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+              DO ng=1,Ngrids
+                transN(ng)=Rbed(ng)
+              END DO
+#endif
+            CASE ('BEDLOAD_COEFF')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+              DO ng=1,Ngrids
+                bedload_coeff(ng)=Rbed(ng)
+              END DO
+#ifdef SED_BIODIFF
+	        CASE ('DBMAX')
+	          Npts=load_r(Nval, Rval, Ngrids, Rbed)
+              DO ng=1,Ngrids
+	            Dbmx(ng)=Rbed(ng)
+              END DO
+            CASE ('DBMIN')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+              DO ng=1,Ngrids
+	            Dbmm(ng)=Rbed(ng)
+              END DO
+            CASE ('DBZS')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+              DO ng=1,Ngrids
+                Dbzs(ng)=Rbed(ng)
+              END DO
+            CASE ('DBZM')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+              DO ng=1,Ngrids
+                Dbzm(ng)=Rbed(ng)
+              END DO
+            CASE ('DBZP')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+              DO ng=1,Ngrids
+                Dbzp(ng)=Rbed(ng)
+              END DO
+#endif
             CASE ('Hadvection')
               IF (itracer.lt.NST) THEN
                 itracer=itracer+1
@@ -133,379 +187,471 @@
      &                      idsed(iTrcStr), idsed(iTrcEnd),             &
      &                      Vname(1,idTvar(idsed(itracer))), ad_LBC)
 #endif
-            CASE ('NEWLAYER_THICK')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('MUD_SD50')
+              IF (.not.allocated(Sd50)) allocate (Sd50(NST,Ngrids))
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                newlayer_thick(ng)=Rbed(ng)
+                DO itrc=1,NCS
+                  Sd50(itrc,ng)=Rmud(itrc,ng)
+                END DO
               END DO
-            CASE ('MINLAYER_THICK')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('MUD_CSED')
+              IF (.not.allocated(Csed)) allocate (Csed(NST,Ngrids))
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud )
               DO ng=1,Ngrids
-                minlayer_thick(ng)=Rbed(ng)
+                DO itrc=1,NCS
+                  Csed(itrc,ng)=Rmud(itrc,ng)
+                END DO
               END DO
-            CASE ('BEDLOAD_COEFF')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('MUD_SRHO')
+              IF (.not.allocated(Srho)) allocate (Srho(NST,Ngrids))
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                bedload_coeff(ng)=Rbed(ng)
+                DO itrc=1,NCS
+                  Srho(itrc,ng)=Rmud(itrc,ng)
+                END DO
               END DO
-#ifdef SED_BIODIFF
-	        CASE ('DBMAX')
-	          Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('MUD_WSED')
+              IF (.not.allocated(Wsed)) allocate (Wsed(NST,Ngrids))
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-	            Dbmx(ng)=Rbed(ng)
+                DO itrc=1,NCS
+                  Wsed(itrc,ng)=Rmud(itrc,ng)
+                END DO
               END DO
-            CASE ('DBMIN')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('MUD_ERATE')
+              IF (.not.allocated(Erate)) allocate (Erate(NST,Ngrids))
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-	            Dbmm(ng)=Rbed(ng)
+                DO itrc=1,NCS
+                  Erate(itrc,ng)=Rmud(itrc,ng)
+                END DO
               END DO
-            CASE ('DBZS')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('MUD_TAU_CE')
+              IF (.not.allocated(tau_ce)) allocate (tau_ce(NST,Ngrids))
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                Dbzs(ng)=Rbed(ng)
+                DO itrc=1,NCS
+                  tau_ce(itrc,ng)=Rmud(itrc,ng)
+                END DO
               END DO
-            CASE ('DBZM')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('MUD_TAU_CD')
+              IF (.not.allocated(tau_cd)) allocate (tau_cd(NST,Ngrids))
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                Dbzm(ng)=Rbed(ng)
+                DO itrc=1,NCS
+                  tau_cd(itrc,ng)=Rmud(itrc,ng)
+                END DO
               END DO
-            CASE ('DBZP')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('MUD_POROS')
+              IF (.not.allocated(poros)) allocate (poros(NST,Ngrids))
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                Dbzp(ng)=Rbed(ng)
-              END DO
-#endif
-            CASE ('SG_ZWBL')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
-                DO ng=1,Ngrids
-                  sg_zwbl(ng)=Rbed(ng)
-                END DO
-#ifdef BEDLOAD
-            CASE ('SEDSLOPE_CRIT_WET')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
-                DO ng=1,Ngrids
-                  sedslope_crit_wet(ng)=Rbed(ng)
-                END DO
-            CASE ('SEDSLOPE_CRIT_DRY')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
-                DO ng=1,Ngrids
-                  sedslope_crit_dry(ng)=Rbed(ng)
-                END DO
-            CASE ('SLOPEFAC_WET')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
-                DO ng=1,Ngrids
-                  slopefac_wet(ng)=Rbed(ng)
-                END DO
-            CASE ('SLOPEFAC_DRY')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
-                DO ng=1,Ngrids
-                  slopefac_dry(ng)=Rbed(ng)
-                END DO
-            CASE ('BEDLOAD_VANDERA_ALPHAW')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
-                DO ng=1,Ngrids
-                  bedload_vandera_alphaw(ng)=Rbed(ng)
-                END DO
-            CASE ('BEDLOAD_VANDERA_ALPHAC')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
-                DO ng=1,Ngrids
-                  bedload_vandera_alphac(ng)=Rbed(ng)
+                DO itrc=1,NCS
+                  poros(itrc,ng)=Rmud(itrc,ng)
                 END DO
-#endif
-            CASE ('Hout(ithck)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbed)
-              i=idSbed(ithck)
-              DO ng=1,Ngrids
-                Hout(i,ng)=Lbed(ng)
               END DO
-            CASE ('Hout(iaged)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbed)
-              i=idSbed(iaged)
+            CASE ('MUD_RXN')
+              IF (.not.allocated(sed_rxn)) allocate (sed_rxn(NST,Ngrids))
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbed(ng)
+                DO itrc=1,NCS
+                  sed_rxn(itrc,ng)=Rmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(iporo)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbed)
-              i=idSbed(iporo)
+            CASE ('MUD_TNU2')
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbed(ng)
+                DO itrc=1,NCS
+                  i=idsed(itrc)
+                  nl_tnu2(i,ng)=Rmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(idiff)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbed)
-              i=idSbed(idiff)
+            CASE ('MUD_TNU4')
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbed(ng)
+                DO itrc=1,NCS
+                  i=idsed(itrc)
+                  nl_tnu4(i,ng)=Rmud(itrc,ng)
+                END DO
               END DO
-# ifdef BEDLOAD_VANDERA
-            CASE ('Hout(idsurs)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+            CASE ('ad_MUD_TNU2')
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                DO itrc=1,NNS
-                  i=idsurs
-                  Hout(i,ng)=Lsand(itrc,ng)
+                DO itrc=1,NCS
+                  i=idsed(itrc)
+                  ad_tnu2(i,ng)=Rmud(itrc,ng)
+                  tl_tnu2(i,ng)=Rmud(itrc,ng)
                 END DO
               END DO
-            CASE ('Hout(idsrrw)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+            CASE ('ad_MUD_TNU4')
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                DO itrc=1,NNS
-                  i=idsrrw
-                  Hout(i,ng)=Lsand(itrc,ng)
+                DO itrc=1,NCS
+                  i=idsed(itrc)
+                  ad_tnu4(i,ng)=Rmud(itrc,ng)
+                  tl_tnu4(i,ng)=Rmud(itrc,ng)
                 END DO
               END DO
-            CASE ('Hout(idsbtw)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+            CASE ('MUD_Sponge')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                DO itrc=1,NNS
-                  i=idsbtw
-                  Hout(i,ng)=Lsand(itrc,ng)
+                DO itrc=1,NCS
+                  i=idsed(itrc)
+                  LtracerSponge(i,ng)=Lmud(itrc,ng)
                 END DO
               END DO
-            CASE ('Hout(idsucr)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+            CASE ('MUD_AKT_BAK')
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                DO itrc=1,NNS
-                  i=idsucr
-                  Hout(i,ng)=Lsand(itrc,ng)
+                DO itrc=1,NCS
+                  i=idsed(itrc)
+                  Akt_bak(i,ng)=Rmud(itrc,ng)
                 END DO
               END DO
-            CASE ('Hout(idsutr)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+            CASE ('MUD_AKT_fac')
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                DO itrc=1,NNS
-                  i=idsutr
-                  Hout(i,ng)=Lsand(itrc,ng)
+                DO itrc=1,NCS
+                  i=idsed(itrc)
+                  ad_Akt_fac(i,ng)=Rmud(itrc,ng)
+                  tl_Akt_fac(i,ng)=Rmud(itrc,ng)
                 END DO
               END DO
-            CASE ('Hout(idstcr)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+            CASE ('MUD_TNUDG')
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                DO itrc=1,NNS
-                  i=idstcr
-                  Hout(i,ng)=Lsand(itrc,ng)
+                DO itrc=1,NCS
+                  i=idsed(itrc)
+                  Tnudg(i,ng)=Rmud(itrc,ng)
                 END DO
               END DO
-            CASE ('Hout(idsttr)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+            CASE ('MUD_MORPH_FAC')
+              IF (.not.allocated(morph_fac)) THEN
+                allocate (morph_fac(NST,Ngrids))
+              END IF
+              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
               DO ng=1,Ngrids
-                DO itrc=1,NNS
-                  i=idsttr
-                  Hout(i,ng)=Lsand(itrc,ng)
+                DO itrc=1,NCS
+                  morph_fac(itrc,ng)=Rmud(itrc,ng)
                 END DO
               END DO
-# endif
-            CASE ('Hout(isd50)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(isd50)
+#if defined COHESIVE_BED || defined MIXED_BED
+            CASE ('MUD_TAUCR_MIN')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                tcr_min(ng)=Rbed(ng)
               END DO
-            CASE ('Hout(idens)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idens)
+            CASE ('MUD_TAUCR_MAX')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
-             END DO
-            CASE ('Hout(iwsed)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(iwsed)
+                tcr_max(ng)=Rbed(ng)
+              END DO
+            CASE ('MUD_TAUCR_SLOPE')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                tcr_slp(ng)=Rbed(ng)
               END DO
-            CASE ('Hout(itauc)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(itauc)
+            CASE ('MUD_TAUCR_OFF')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                tcr_off(ng)=Rbed(ng)
               END DO
-            CASE ('Hout(irlen)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(irlen)
+            CASE ('MUD_TAUCR_TIME')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                tcr_tim(ng)=Rbed(ng)
               END DO
-            CASE ('Hout(irhgt)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(irhgt)
+#endif
+            CASE ('MUD_Ltsrc', 'MUD_Ltracer')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idsed(itrc)
+                  LtracerSrc(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(ibwav)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(ibwav)
+            CASE ('MUD_Ltclm')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idsed(itrc)
+                  LtracerCLM(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(izdef)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(izdef)
+            CASE ('MUD_Tnudge')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idsed(itrc)
+                  LnudgeTCLM(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(izapp)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(izapp)
+            CASE ('Hout(idmud)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idTvar(idsed(itrc))
+                  Hout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(izNik)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(izNik)
-             DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+            CASE ('Hout(iMfrac)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+              DO ng=1,Ngrids
+                DO itrc=1,NCS
+                  i=idfrac(itrc)
+                  Hout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(izbio)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(izbio)
+            CASE ('Hout(iMmass)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idBmas(itrc)
+                  Hout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(izbfm)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(izbfm)
+#ifdef BEDLOAD
+            CASE ('Hout(iMUbld)')
               DO ng=1,Ngrids
-               Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  IF (idUbld(itrc).eq.0) THEN
+                    IF (Master) WRITE (out,30) 'idUbld'
+                    exit_flag=5
+                    RETURN
+                  END IF
+                END DO
               END DO
-            CASE ('Hout(izbld)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(izbld)
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idUbld(itrc)
+                  Hout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(izwbl)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(izwbl)
+            CASE ('Hout(iMVbld)')
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  IF (idVbld(itrc).eq.0) THEN
+                    IF (Master) WRITE (out,30) 'idVbld'
+                    exit_flag=5
+                    RETURN
+                  END IF
+                END DO
               END DO
-            CASE ('Hout(iactv)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(iactv)
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idVbld(itrc)
+                  Hout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(ishgt)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(ishgt)
+#endif
+            CASE ('Qout(idmud)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idTvar(idsed(itrc))
+                  Qout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(imaxD)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(imaxD)
+            CASE ('Qout(iSmud)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idsurT(idsed(itrc))
+                  Qout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(idnet)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idnet)
+            CASE ('Qout(iMfrac)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idfrac(itrc)
+                  Qout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(idtbl)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idtbl)
+            CASE ('Qout(iMmass)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idBmas(itrc)
+                  Qout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(idubl)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idubl)
+#ifdef BEDLOAD
+            CASE ('Qout(iMUbld)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idUbld(itrc)
+                  Qout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(idfdw)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idfdw)
+            CASE ('Qout(iMVbld)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idVbld(itrc)
+                  Qout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(idzrw)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idzrw)
+#endif
+#if defined AVERAGES    || \
+   (defined AD_AVERAGES && defined ADJOINT) || \
+   (defined RP_AVERAGES && defined TL_IOMS) || \
+   (defined TL_AVERAGES && defined TANGENT)
+            CASE ('Aout(idmud)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idTvar(idsed(itrc))
+                  Aout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(idksd)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idksd)
+            CASE ('Aout(iMTTav)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idTTav(idsed(itrc))
+                  Aout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(idusc)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idusc)
+            CASE ('Aout(iMUTav)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idUTav(idsed(itrc))
+                  Aout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(idpcx)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idpcx)
+            CASE ('Aout(iMVTav)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO itrc=1,NCS
+                  i=idVTav(idsed(itrc))
+                  Aout(i,ng)=Lmud(itrc,ng)
+                END DO
+              END DO
+            CASE ('Aout(MHUTav)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+              DO ng=1,Ngrids
+                DO itrc=1,NCS
+                  i=iHUTav(idsed(itrc))
+                  Aout(i,ng)=Lmud(itrc,ng)
+                END DO
+              END DO
+            CASE ('Aout(MHVTav)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+              DO ng=1,Ngrids
+                DO itrc=1,NCS
+                  i=iHVTav(idsed(itrc))
+                  Aout(i,ng)=Lmud(itrc,ng)
+                END DO
+              END DO
+# ifdef BEDLOAD
+            CASE ('Aout(iMUbld)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+              DO ng=1,Ngrids
+                DO itrc=1,NCS
+                  i=idUbld(itrc)
+                  Aout(i,ng)=Lmud(itrc,ng)
+                END DO
+              END DO
+            CASE ('Aout(iMVbld)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+              DO ng=1,Ngrids
+                DO itrc=1,NCS
+                  i=idVbld(itrc)
+                  Aout(i,ng)=Lmud(itrc,ng)
+                END DO
               END DO
-            CASE ('Hout(idpwc)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idpwc)
+# endif
+#endif
+#ifdef DIAGNOSTICS_TS
+            CASE ('Dout(MTrate)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO i=1,NCS
+                  itrc=idsed(i)
+                  Dout(idDtrc(itrc,iTrate),ng)=Lmud(i,ng)
+                END DO
               END DO
-#if defined COHESIVE_BED || defined SED_BIODIFF || defined MIXED_BED
-            CASE ('Hout(idoff)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idoff)
+            CASE ('Dout(MThadv)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO i=1,NCS
+                  itrc=idsed(i)
+                  Dout(idDtrc(itrc,iThadv),ng)=Lmud(i,ng)
+                END DO
               END DO
-            CASE ('Hout(idslp)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idslp)
+            CASE ('Dout(MTxadv)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO i=1,NCS
+                  itrc=idsed(i)
+                  Dout(idDtrc(itrc,iTxadv),ng)=Lmud(i,ng)
+                END DO
               END DO
-            CASE ('Hout(idtim)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idtim)
+            CASE ('Dout(MTyadv)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO i=1,NCS
+                  itrc=idsed(i)
+                  Dout(idDtrc(itrc,iTyadv),ng)=Lmud(i,ng)
+                END DO
               END DO
-            CASE ('Hout(idbmx)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idbmx)
+            CASE ('Dout(MTvadv)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO i=1,NCS
+                  itrc=idsed(i)
+                  Dout(idDtrc(itrc,iTvadv),ng)=Lmud(i,ng)
+                END DO
               END DO
-            CASE ('Hout(idbmm)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idbmm)
+# if defined TS_DIF2 || defined TS_DIF4
+            CASE ('Dout(MThdif)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO i=1,NCS
+                  itrc=idsed(i)
+                  Dout(idDtrc(itrc,iThdif),ng)=Lmud(i,ng)
+                END DO
               END DO
-            CASE ('Hout(idbzs)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idbzs)
+            CASE ('Dout(MTxdif)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO i=1,NCS
+                  itrc=idsed(i)
+                  Dout(idDtrc(itrc,iTxdif),ng)=Lmud(i,ng)
+                END DO
               END DO
-            CASE ('Hout(idbzm)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idbzm)
+            CASE ('Dout(MTydif)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO i=1,NCS
+                  itrc=idsed(i)
+                  Dout(idDtrc(itrc,iTydif),ng)=Lmud(i,ng)
+                END DO
               END DO
-            CASE ('Hout(idbzp)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idbzp)
+#  if defined MIX_GEO_TS || defined MIX_ISO_TS
+            CASE ('Dout(MTsdif)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO i=1,NCS
+                  itrc=idsed(i)
+                  Dout(idDtrc(itrc,iTsdif),ng)=Lmud(i,ng)
+                END DO
               END DO
-#endif
-#if defined MIXED_BED
-            CASE ('Hout(idprp)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
-              i=idBott(idprp)
+#  endif
+# endif
+            CASE ('Dout(MTvdif)')
+              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
               DO ng=1,Ngrids
-                Hout(i,ng)=Lbottom(ng)
+                DO i=1,NCS
+                  itrc=idsed(i)
+                  Dout(idDtrc(itrc,iTvdif),ng)=Lmud(i,ng)
+                END DO
               END DO
 #endif
             CASE ('SAND_SD50')
@@ -682,74 +828,223 @@
                   i=idsed(NCS+itrc)
                   LtracerCLM(i,ng)=Lsand(itrc,ng)
                 END DO
-              END DO
-            CASE ('SAND_Tnudge')
+              END DO
+            CASE ('SAND_Tnudge')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idsed(NCS+itrc)
+                  LnudgeTCLM(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+            CASE ('Hout(idsand)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idTvar(idsed(NCS+itrc))
+                  Hout(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+            CASE ('Hout(iSfrac)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idfrac(NCS+itrc)
+                  Hout(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+            CASE ('Hout(iSmass)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idBmas(NCS+itrc)
+                  Hout(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+#ifdef BEDLOAD
+            CASE ('Hout(iSUbld)')
+              DO ng=1,Ngrids
+                DO itrc=NCS+1,NST
+                  IF (idUbld(itrc).eq.0) THEN
+                    IF (Master) WRITE (out,30) 'idUbld'
+                    exit_flag=5
+                    RETURN
+                  END IF
+                END DO
+              END DO
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idUbld(NCS+itrc)
+                  Hout(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+            CASE ('Hout(iSVbld)')
+              DO ng=1,Ngrids
+                DO itrc=NCS+1,NST
+                  IF (idVbld(itrc).eq.0) THEN
+                    IF (Master) WRITE (out,30) 'idVbld'
+                    exit_flag=5
+                    RETURN
+                  END IF
+                END DO
+              END DO
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idVbld(NCS+itrc)
+                  Hout(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+!# ifdef BEDLOAD_VANDERA
+            CASE ('Hout(idsurs)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idsurs
+                  Hout(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+            CASE ('Hout(idsrrw)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idsrrw
+                  Hout(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+            CASE ('Hout(idsbtw)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idsbtw
+                  Hout(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+            CASE ('Hout(idsksd)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idsksd
+                  Hout(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+            CASE ('Hout(idszrw)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idszrw
+                  Hout(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+            CASE ('Hout(idsusc)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
+              DO ng=1,Ngrids
+                DO itrc=1,NNS
+                  i=idsusc
+                  Hout(i,ng)=Lsand(itrc,ng)
+                END DO
+              END DO
+            CASE ('SG_ZWBL')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+                DO ng=1,Ngrids
+                  sg_zwbl(ng)=Rbed(ng)
+                END DO
+            CASE ('SEDSLOPE_CRIT_WET')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+                DO ng=1,Ngrids
+                  sedslope_crit_wet(ng)=Rbed(ng)
+                END DO
+            CASE ('SEDSLOPE_CRIT_DRY')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+                DO ng=1,Ngrids
+                  sedslope_crit_dry(ng)=Rbed(ng)
+                END DO
+            CASE ('SLOPEFAC_WET')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+                DO ng=1,Ngrids
+                  slopefac_wet(ng)=Rbed(ng)
+                END DO
+            CASE ('SLOPEFAC_DRY')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+                DO ng=1,Ngrids
+                  slopefac_dry(ng)=Rbed(ng)
+                END DO
+            CASE ('BEDLOAD_VANDERA_ALPHAC')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+                DO ng=1,Ngrids
+                  bedload_vandera_alphac(ng)=Rbed(ng)
+                END DO
+            CASE ('BEDLOAD_VANDERA_ALPHAW')
+              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+                DO ng=1,Ngrids
+                  bedload_vandera_alphaw(ng)=Rbed(ng)
+                END DO
+            CASE ('Hout(idstbl)')
               Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
                 DO itrc=1,NNS
-                  i=idsed(NCS+itrc)
-                  LnudgeTCLM(i,ng)=Lsand(itrc,ng)
+                  i=idstbl
+                  Hout(i,ng)=Lsand(itrc,ng)
                 END DO
               END DO
-            CASE ('Hout(idsand)')
+            CASE ('Hout(idsubl)')
               Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
                 DO itrc=1,NNS
-                  i=idTvar(idsed(NCS+itrc))
+                  i=idsubl
                   Hout(i,ng)=Lsand(itrc,ng)
                 END DO
               END DO
-            CASE ('Hout(iSfrac)')
+            CASE ('Hout(idspwc)')
               Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
                 DO itrc=1,NNS
-                  i=idfrac(NCS+itrc)
+                  i=idspwc
                   Hout(i,ng)=Lsand(itrc,ng)
                 END DO
               END DO
-            CASE ('Hout(iSmass)')
+            CASE ('Hout(idsfdw)')
               Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
                 DO itrc=1,NNS
-                  i=idBmas(NCS+itrc)
+                  i=idsfdw
                   Hout(i,ng)=Lsand(itrc,ng)
                 END DO
               END DO
-#ifdef BEDLOAD
-            CASE ('Hout(iSUbld)')
+            CASE ('Hout(idsucr)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=NCS+1,NST
-                  IF (idUbld(itrc).eq.0) THEN
-                    IF (Master) WRITE (out,30) 'idUbld'
-                    exit_flag=5
-                    RETURN
-                  END IF
+                DO itrc=1,NNS
+                  i=idsucr
+                  Hout(i,ng)=Lsand(itrc,ng)
                 END DO
               END DO
+            CASE ('Hout(idsutr)')
               Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
                 DO itrc=1,NNS
-                  i=idUbld(NCS+itrc)
+                  i=idsutr
                   Hout(i,ng)=Lsand(itrc,ng)
                 END DO
               END DO
-            CASE ('Hout(iSVbld)')
+            CASE ('Hout(idstcr)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=NCS+1,NST
-                  IF (idVbld(itrc).eq.0) THEN
-                    IF (Master) WRITE (out,30) 'idVbld'
-                    exit_flag=5
-                    RETURN
-                  END IF
+                DO itrc=1,NNS
+                  i=idstcr
+                  Hout(i,ng)=Lsand(itrc,ng)
                 END DO
               END DO
+            CASE ('Hout(idsttr)')
               Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
                 DO itrc=1,NNS
-                  i=idVbld(NCS+itrc)
+                  i=idsttr
                   Hout(i,ng)=Lsand(itrc,ng)
                 END DO
               END DO
+!# endif
 #endif
             CASE ('Qout(idsand)')
               Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
@@ -865,576 +1160,327 @@
             CASE ('Aout(iSVbld)')
               Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=1,NNS
-                  i=idVbld(NCS+itrc)
-                  Aout(i,ng)=Lsand(itrc,ng)
-                END DO
-              END DO
-# endif
-#endif
-#ifdef DIAGNOSTICS_TS
-            CASE ('Dout(STrate)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
-              DO ng=1,Ngrids
-                DO i=1,NNS
-                  itrc=idsed(NCS+i)
-                  Dout(idDtrc(itrc,iTrate),ng)=Lsand(i,ng)
-                END DO
-              END DO
-            CASE ('Dout(SThadv)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
-              DO ng=1,Ngrids
-                DO i=1,NNS
-                  itrc=idsed(NCS+i)
-                  Dout(idDtrc(itrc,iThadv),ng)=Lsand(i,ng)
-                END DO
-              END DO
-            CASE ('Dout(STxadv)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
-              DO ng=1,Ngrids
-                DO i=1,NNS
-                  itrc=idsed(NCS+i)
-                  Dout(idDtrc(itrc,iTxadv),ng)=Lsand(i,ng)
-                END DO
-              END DO
-            CASE ('Dout(STyadv)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
-              DO ng=1,Ngrids
-                DO i=1,NNS
-                  itrc=idsed(NCS+i)
-                  Dout(idDtrc(itrc,iTyadv),ng)=Lsand(i,ng)
-                END DO
-              END DO
-            CASE ('Dout(STvadv)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
-              DO ng=1,Ngrids
-                DO i=1,NNS
-                  itrc=idsed(NCS+i)
-                  Dout(idDtrc(itrc,iTvadv),ng)=Lsand(i,ng)
-                END DO
-              END DO
-# if defined TS_DIF2 || defined TS_DIF4
-            CASE ('Dout(SThdif)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
-              DO ng=1,Ngrids
-                DO i=1,NNS
-                  itrc=idsed(NCS+i)
-                  Dout(idDtrc(itrc,iThdif),ng)=Lsand(i,ng)
-                END DO
-              END DO
-            CASE ('Dout(STxdif)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
-              DO ng=1,Ngrids
-                DO i=1,NNS
-                  itrc=idsed(NCS+i)
-                  Dout(idDtrc(itrc,iTxdif),ng)=Lsand(i,ng)
-                END DO
-              END DO
-            CASE ('Dout(STydif)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
-              DO ng=1,Ngrids
-                DO i=1,NNS
-                  itrc=idsed(NCS+i)
-                  Dout(idDtrc(itrc,iTydif),ng)=Lsand(i,ng)
-                END DO
-              END DO
-#  if defined MIX_GEO_TS || defined MIX_ISO_TS
-            CASE ('Dout(STsdif)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
-              DO ng=1,Ngrids
-                DO i=1,NNS
-                  itrc=idsed(NCS+i)
-                  Dout(idDtrc(itrc,iTsdif),ng)=Lsand(i,ng)
-                END DO
-              END DO
-#  endif
-# endif
-            CASE ('Dout(STvdif)')
-              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
-              DO ng=1,Ngrids
-                DO i=1,NNS
-                  itrc=idsed(NCS+i)
-                  Dout(idDtrc(itrc,iTvdif),ng)=Lsand(i,ng)
-                END DO
-              END DO
-#endif
-            CASE ('MUD_SD50')
-              IF (.not.allocated(Sd50)) allocate (Sd50(NST,Ngrids))
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  Sd50(itrc,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_CSED')
-              IF (.not.allocated(Csed)) allocate (Csed(NST,Ngrids))
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud )
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  Csed(itrc,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_SRHO')
-              IF (.not.allocated(Srho)) allocate (Srho(NST,Ngrids))
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  Srho(itrc,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_WSED')
-              IF (.not.allocated(Wsed)) allocate (Wsed(NST,Ngrids))
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  Wsed(itrc,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_ERATE')
-              IF (.not.allocated(Erate)) allocate (Erate(NST,Ngrids))
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  Erate(itrc,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_TAU_CE')
-              IF (.not.allocated(tau_ce)) allocate (tau_ce(NST,Ngrids))
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  tau_ce(itrc,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_TAU_CD')
-              IF (.not.allocated(tau_cd)) allocate (tau_cd(NST,Ngrids))
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  tau_cd(itrc,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_POROS')
-              IF (.not.allocated(poros)) allocate (poros(NST,Ngrids))
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  poros(itrc,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_RXN')
-              IF (.not.allocated(sed_rxn)) allocate (sed_rxn(NST,Ngrids))
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  sed_rxn(itrc,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_TNU2')
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsed(itrc)
-                  nl_tnu2(i,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_TNU4')
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsed(itrc)
-                  nl_tnu4(i,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('ad_MUD_TNU2')
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsed(itrc)
-                  ad_tnu2(i,ng)=Rmud(itrc,ng)
-                  tl_tnu2(i,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('ad_MUD_TNU4')
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsed(itrc)
-                  ad_tnu4(i,ng)=Rmud(itrc,ng)
-                  tl_tnu4(i,ng)=Rmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_Sponge')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsed(itrc)
-                  LtracerSponge(i,ng)=Lmud(itrc,ng)
-                END DO
-              END DO
-            CASE ('MUD_AKT_BAK')
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
-              DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsed(itrc)
-                  Akt_bak(i,ng)=Rmud(itrc,ng)
+                DO itrc=1,NNS
+                  i=idVbld(NCS+itrc)
+                  Aout(i,ng)=Lsand(itrc,ng)
                 END DO
               END DO
-            CASE ('MUD_AKT_fac')
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
+# endif
+#endif
+#ifdef DIAGNOSTICS_TS
+            CASE ('Dout(STrate)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsed(itrc)
-                  ad_Akt_fac(i,ng)=Rmud(itrc,ng)
-                  tl_Akt_fac(i,ng)=Rmud(itrc,ng)
+                DO i=1,NNS
+                  itrc=idsed(NCS+i)
+                  Dout(idDtrc(itrc,iTrate),ng)=Lsand(i,ng)
                 END DO
               END DO
-            CASE ('MUD_TNUDG')
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
+            CASE ('Dout(SThadv)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsed(itrc)
-                  Tnudg(i,ng)=Rmud(itrc,ng)
+                DO i=1,NNS
+                  itrc=idsed(NCS+i)
+                  Dout(idDtrc(itrc,iThadv),ng)=Lsand(i,ng)
                 END DO
               END DO
-            CASE ('MUD_MORPH_FAC')
-              IF (.not.allocated(morph_fac)) THEN
-                allocate (morph_fac(NST,Ngrids))
-              END IF
-              Npts=load_r(Nval, Rval, NCS, Ngrids, Rmud)
+            CASE ('Dout(STxadv)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  morph_fac(itrc,ng)=Rmud(itrc,ng)
+                DO i=1,NNS
+                  itrc=idsed(NCS+i)
+                  Dout(idDtrc(itrc,iTxadv),ng)=Lsand(i,ng)
                 END DO
               END DO
-            CASE ('MUD_Ltsrc', 'MUD_Ltracer')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Dout(STyadv)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsed(itrc)
-                  LtracerSrc(i,ng)=Lmud(itrc,ng)
+                DO i=1,NNS
+                  itrc=idsed(NCS+i)
+                  Dout(idDtrc(itrc,iTyadv),ng)=Lsand(i,ng)
                 END DO
               END DO
-            CASE ('MUD_Ltclm')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Dout(STvadv)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsed(itrc)
-                  LtracerCLM(i,ng)=Lmud(itrc,ng)
+                DO i=1,NNS
+                  itrc=idsed(NCS+i)
+                  Dout(idDtrc(itrc,iTvadv),ng)=Lsand(i,ng)
                 END DO
               END DO
-            CASE ('MUD_Tnudge')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+# if defined TS_DIF2 || defined TS_DIF4
+            CASE ('Dout(SThdif)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsed(itrc)
-                  LnudgeTCLM(i,ng)=Lmud(itrc,ng)
+                DO i=1,NNS
+                  itrc=idsed(NCS+i)
+                  Dout(idDtrc(itrc,iThdif),ng)=Lsand(i,ng)
                 END DO
               END DO
-            CASE ('Hout(idmud)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Dout(STxdif)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idTvar(idsed(itrc))
-                  Hout(i,ng)=Lmud(itrc,ng)
+                DO i=1,NNS
+                  itrc=idsed(NCS+i)
+                  Dout(idDtrc(itrc,iTxdif),ng)=Lsand(i,ng)
                 END DO
               END DO
-            CASE ('Hout(iMfrac)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Dout(STydif)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idfrac(itrc)
-                  Hout(i,ng)=Lmud(itrc,ng)
+                DO i=1,NNS
+                  itrc=idsed(NCS+i)
+                  Dout(idDtrc(itrc,iTydif),ng)=Lsand(i,ng)
                 END DO
               END DO
-            CASE ('Hout(iMmass)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+#  if defined MIX_GEO_TS || defined MIX_ISO_TS
+            CASE ('Dout(STsdif)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idBmas(itrc)
-                  Hout(i,ng)=Lmud(itrc,ng)
+                DO i=1,NNS
+                  itrc=idsed(NCS+i)
+                  Dout(idDtrc(itrc,iTsdif),ng)=Lsand(i,ng)
                 END DO
               END DO
-#ifdef BEDLOAD
-            CASE ('Hout(iMUbld)')
+#  endif
+# endif
+            CASE ('Dout(STvdif)')
+              Npts=load_l(Nval, Cval, NNS, Ngrids, Lsand)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  IF (idUbld(itrc).eq.0) THEN
-                    IF (Master) WRITE (out,30) 'idUbld'
-                    exit_flag=5
-                    RETURN
-                  END IF
+                DO i=1,NNS
+                  itrc=idsed(NCS+i)
+                  Dout(idDtrc(itrc,iTvdif),ng)=Lsand(i,ng)
                 END DO
               END DO
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+#endif
+            CASE ('Hout(ithck)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbed)
+              i=idSbed(ithck)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idUbld(itrc)
-                  Hout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbed(ng)
               END DO
-            CASE ('Hout(iMVbld)')
+            CASE ('Hout(iaged)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbed)
+              i=idSbed(iaged)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  IF (idVbld(itrc).eq.0) THEN
-                    IF (Master) WRITE (out,30) 'idVbld'
-                    exit_flag=5
-                    RETURN
-                  END IF
-                END DO
+                Hout(i,ng)=Lbed(ng)
               END DO
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(iporo)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbed)
+              i=idSbed(iporo)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idVbld(itrc)
-                  Hout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbed(ng)
+              END DO
+#if defined COHESIVE_BED || defined SED_BIODIFF || defined MIXED_BED
+            CASE ('Hout(ibtcr)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbed)
+              i=idSbed(ibtcr)
+              DO ng=1,Ngrids
+                Hout(i,ng)=Lbed(ng)
               END DO
 #endif
-            CASE ('Qout(idmud)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+#if defined SEDBIO_COUP
+            CASE ('Hout(iboxy)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbed)
+              i=idSbed(iboxy)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idTvar(idsed(itrc))
-                  Qout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbed(ng)
               END DO
-            CASE ('Qout(iSmud)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(ibno3)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbed)
+              i=idSbed(ibno3)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idsurT(idsed(itrc))
-                  Qout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbed(ng)
               END DO
-            CASE ('Qout(iMfrac)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(ibnh4)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbed)
+              i=idSbed(ibnh4)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idfrac(itrc)
-                  Qout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbed(ng)
               END DO
-            CASE ('Qout(iMmass)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(ibodu)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbed)
+              i=idSbed(ibodu)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idBmas(itrc)
-                  Qout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbed(ng)
               END DO
-#ifdef BEDLOAD
-            CASE ('Qout(iMUbld)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+#endif
+            CASE ('Hout(idiff)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbed)
+              i=idSbed(idiff)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idUbld(itrc)
-                  Qout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbed(ng)
               END DO
-            CASE ('Qout(iMVbld)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(isd50)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(isd50)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idVbld(itrc)
-                  Qout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-#endif
-#if defined AVERAGES    || \
-   (defined AD_AVERAGES && defined ADJOINT) || \
-   (defined RP_AVERAGES && defined TL_IOMS) || \
-   (defined TL_AVERAGES && defined TANGENT)
-            CASE ('Aout(idmud)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(itpor)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(itpor)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idTvar(idsed(itrc))
-                  Aout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Aout(iMTTav)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(idens)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(idens)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idTTav(idsed(itrc))
-                  Aout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
+             END DO
+            CASE ('Hout(iwsed)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(iwsed)
+              DO ng=1,Ngrids
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Aout(iMUTav)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(itauc)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(itauc)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idUTav(idsed(itrc))
-                  Aout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Aout(iMVTav)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(irlen)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(irlen)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idVTav(idsed(itrc))
-                  Aout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Aout(MHUTav)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(irhgt)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(irhgt)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=iHUTav(idsed(itrc))
-                  Aout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Aout(MHVTav)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(ibwav)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(ibwav)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=iHVTav(idsed(itrc))
-                  Aout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-# ifdef BEDLOAD
-            CASE ('Aout(iMUbld)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(izdef)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(izdef)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idUbld(itrc)
-                  Aout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Aout(iMVbld)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(izapp)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(izapp)
               DO ng=1,Ngrids
-                DO itrc=1,NCS
-                  i=idVbld(itrc)
-                  Aout(i,ng)=Lmud(itrc,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-# endif
-#endif
-#ifdef DIAGNOSTICS_TS
-            CASE ('Dout(MTrate)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(izNik)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(izNik)
+             DO ng=1,Ngrids
+                Hout(i,ng)=Lbottom(ng)
+              END DO
+            CASE ('Hout(izbio)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(izbio)
               DO ng=1,Ngrids
-                DO i=1,NCS
-                  itrc=idsed(i)
-                  Dout(idDtrc(itrc,iTrate),ng)=Lmud(i,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Dout(MThadv)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(izbfm)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(izbfm)
               DO ng=1,Ngrids
-                DO i=1,NCS
-                  itrc=idsed(i)
-                  Dout(idDtrc(itrc,iThadv),ng)=Lmud(i,ng)
-                END DO
+               Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Dout(MTxadv)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(izbld)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(izbld)
               DO ng=1,Ngrids
-                DO i=1,NCS
-                  itrc=idsed(i)
-                  Dout(idDtrc(itrc,iTxadv),ng)=Lmud(i,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Dout(MTyadv)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(izwbl)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(izwbl)
               DO ng=1,Ngrids
-                DO i=1,NCS
-                  itrc=idsed(i)
-                  Dout(idDtrc(itrc,iTyadv),ng)=Lmud(i,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Dout(MTvadv)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(iactv)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(iactv)
               DO ng=1,Ngrids
-                DO i=1,NCS
-                  itrc=idsed(i)
-                  Dout(idDtrc(itrc,iTvadv),ng)=Lmud(i,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-# if defined TS_DIF2 || defined TS_DIF4
-            CASE ('Dout(MThdif)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(ishgt)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(ishgt)
               DO ng=1,Ngrids
-                DO i=1,NCS
-                  itrc=idsed(i)
-                  Dout(idDtrc(itrc,iThdif),ng)=Lmud(i,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Dout(MTxdif)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(imaxD)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(imaxD)
               DO ng=1,Ngrids
-                DO i=1,NCS
-                  itrc=idsed(i)
-                  Dout(idDtrc(itrc,iTxdif),ng)=Lmud(i,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('Dout(MTydif)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(idnet)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(idnet)
               DO ng=1,Ngrids
-                DO i=1,NCS
-                  itrc=idsed(i)
-                  Dout(idDtrc(itrc,iTydif),ng)=Lmud(i,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-#  if defined MIX_GEO_TS || defined MIX_ISO_TS
-            CASE ('Dout(MTsdif)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+#if defined COHESIVE_BED || defined SED_BIODIFF || defined MIXED_BED
+            CASE ('Hout(idoff)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(idoff)
               DO ng=1,Ngrids
-                DO i=1,NCS
-                  itrc=idsed(i)
-                  Dout(idDtrc(itrc,iTsdif),ng)=Lmud(i,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-#  endif
-# endif
-            CASE ('Dout(MTvdif)')
-              Npts=load_l(Nval, Cval, NCS, Ngrids, Lmud)
+            CASE ('Hout(idslp)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(idslp)
               DO ng=1,Ngrids
-                DO i=1,NCS
-                  itrc=idsed(i)
-                  Dout(idDtrc(itrc,iTvdif),ng)=Lmud(i,ng)
-                END DO
+                Hout(i,ng)=Lbottom(ng)
               END DO
-#endif
-#ifdef MIXED_BED
-            CASE ('TRANSC')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('Hout(idtim)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(idtim)
               DO ng=1,Ngrids
-                transC(ng)=Rbed(ng)
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('TRANSN')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('Hout(idbmx)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(idbmx)
               DO ng=1,Ngrids
-                transN(ng)=Rbed(ng)
+                Hout(i,ng)=Lbottom(ng)
               END DO
-#endif
-#if defined COHESIVE_BED || defined MIXED_BED
-            CASE ('MUD_TAUCR_MIN')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('Hout(idbmm)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(idbmm)
               DO ng=1,Ngrids
-                tcr_min(ng)=Rbed(ng)
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('MUD_TAUCR_MAX')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('Hout(idbzs)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(idbzs)
               DO ng=1,Ngrids
-                tcr_max(ng)=Rbed(ng)
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('MUD_TAUCR_SLOPE')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('Hout(idbzm)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(idbzm)
               DO ng=1,Ngrids
-                tcr_slp(ng)=Rbed(ng)
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('MUD_TAUCR_OFF')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+            CASE ('Hout(idbzp)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(idbzp)
               DO ng=1,Ngrids
-                tcr_off(ng)=Rbed(ng)
+                Hout(i,ng)=Lbottom(ng)
               END DO
-            CASE ('MUD_TAUCR_TIME')
-              Npts=load_r(Nval, Rval, Ngrids, Rbed)
+#endif
+#if defined MIXED_BED
+            CASE ('Hout(idprp)')
+              Npts=load_l(Nval, Cval, Ngrids, Lbottom)
+              i=idBott(idprp)
               DO ng=1,Ngrids
-                tcr_tim(ng)=Rbed(ng)
+                Hout(i,ng)=Lbottom(ng)
               END DO
 #endif
 #if defined SED_FLOCS
@@ -1643,41 +1689,6 @@
               DO ng=1,Ngrids
                 Qout(i,ng)=Lbottom(ng)
               END DO
-#if defined COHESIVE_BED || defined SED_BIODIFF || defined MIXED_BED
-            CASE ('Hout(ibtcr)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbed)
-              i=idSbed(ibtcr)
-              DO ng=1,Ngrids
-                Hout(i,ng)=Lbed(ng)
-              END DO
-#endif
-#if defined SEDBIO_COUP
-            CASE ('Hout(iboxy)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbed)
-              i=idSbed(iboxy)
-              DO ng=1,Ngrids
-                Hout(i,ng)=Lbed(ng)
-              END DO
-            CASE ('Hout(ibno3)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbed)
-              i=idSbed(ibno3)
-              DO ng=1,Ngrids
-                Hout(i,ng)=Lbed(ng)
-              END DO
-            CASE ('Hout(ibnh4)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbed)
-              i=idSbed(ibnh4)
-              DO ng=1,Ngrids
-                Hout(i,ng)=Lbed(ng)
-              END DO
-            CASE ('Hout(ibodu)')
-              Npts=load_l(Nval, Cval, Ngrids, Lbed)
-              i=idSbed(ibodu)
-              DO ng=1,Ngrids
-                Hout(i,ng)=Lbed(ng)
-              END DO
-#endif
-
           END SELECT
         END IF
       END DO
@@ -1790,7 +1801,7 @@
               DO itrc=1,NST
                 i=idBmas(itrc)
                 IF (Hout(i,ng)) WRITE (out,160) Hout(i,ng),             &
-     &              'Hout(idmass)',                                     &
+     &              'Hout(idfrac)',                                     &
      &              'Write out mass, sediment ', itrc,                  &
      &              TRIM(Vname(1,i))
               END DO
diff --git a/ROMS/Nonlinear/Sediment/sediment_mod.h b/ROMS/Nonlinear/Sediment/sediment_mod.h
index 36b1101f90bb9adcf21c3091a03cf7befbc4f851..257b5afe93409e68763441c12857ba2aa2641cb1 100755
--- a/ROMS/Nonlinear/Sediment/sediment_mod.h
+++ b/ROMS/Nonlinear/Sediment/sediment_mod.h
@@ -48,10 +48,10 @@
 !   idiff           Sediment layer bio-diffusivity (m2/s).             !
 !   ibtcr           Sediment critical stress for erosion (Pa).         !
 #if defined SEDBIO_COUP
-!   iboxy           Sediment porewater oxygen (mmol O2/m2)             !
-!   ibno3           Sediment porewater nitrate (mmol NO3/m2)           !
-!   ibnh4           Sediment porewater ammonium (mmol NH4/m2)          !
-!   ibodu           Sediment porewater oxygen demand units (mmol O2/m2)!
+!   iboxy           Sediment porewater oxygen (mmol O2/m2)             ! 
+!   ibno3           Sediment porewater nitrate (mmol NO3/m2)           ! 
+!   ibnh4           Sediment porewater ammonium (mmol NH4/m2)          ! 
+!   ibodu           Sediment porewater oxygen demand units (mmol O2/m2)! 
 #endif
 !                                                                      !
 !  BOTTOM properties indices:                                          !
@@ -76,18 +76,13 @@
 !   iactv           Active layer thickness for erosive potential (m).  !
 !   ishgt           Sediment saltation height (m).                     !
 !   imaxD           Maximum inundation depth.                          !
-!   idnet           Erosion/deposition                                 !
-!   idtbl           Thickness at wave boundary layer                   !
-!   idubl           Current velocity at wave boundary layer            !
-!   idfdw           Friction factor from the currents                  !
-!   idzrw           Reference height to get near bottom current vel    !
-!   idksd           Bed roughness (Zo) for the wave boundary layer calc!
-!   idusc           Current friction velocity the wave boundary layer  !
-!   idpcx           Anlge between currents and xi axis                 !
-!   idpwc           Angle between waves/currents                       !
-!
-!   idoff           Offset for calc of dmix erodibility profile (m).   !
-!   idslp           Slope  for calc of dmix or erodibility profile.    !
+!   isgrH           Seagrass height.                                   !
+!   isgrD           Seagrass shoot density.                            !
+!   idnet           Erosion or deposition.                             !
+!   idoff           Offset for calculation of dmix erodibility         !
+!                     profile (m).                                     !
+!   idslp           Slope  for calculation of dmix or erodibility      !
+!                     profile.                                         !
 !   idtim           Time scale for restoring erodibility profile (s).  !
 !   idbmx           Bed biodifusivity maximum.                         !
 !   idbmm           Bed biodifusivity minimum.                         !
@@ -96,8 +91,6 @@
 !   idbzp           Bed biodifusivity phi.                             !
 !   idprp           Cohesive behavior.                                 !
 !                                                                      !
-!   isgrH           Seagrass height.                                   !
-!   isgrD           Seagrass shoot density.                            !
 !   nTbiom          Number of hours for depth integration              !
 !=======================================================================
 !
@@ -123,7 +116,7 @@
       integer :: MBEDP                     ! Number of bed properties
       integer :: ithck, iaged, iporo, idiff
 #if defined COHESIVE_BED || defined SED_BIODIFF || defined MIXED_BED
-      integer :: ibtcr
+      integer :: ibtcr 
 #endif
 #if defined SEDBIO_COUP
       integer :: iboxy, ibno3, ibnh4, ibodu
@@ -138,9 +131,8 @@
       integer :: irlen, irhgt, ibwav, izdef
       integer :: izapp, izNik, izbio, izbfm
       integer :: izbld, izwbl, iactv, ishgt
-      integer :: imaxD, idnet
-      integer :: idtbl, idubl, idfdw, idzrw
-      integer :: idksd, idusc, idpcx, idpwc
+      integer :: imaxD, idnet  
+      integer :: itpor
 #if defined COHESIVE_BED || defined SED_BIODIFF || defined MIXED_BED
       integer :: idoff, idslp, idtim, idbmx
       integer :: idbmm, idbzs, idbzm, idbzp
@@ -161,6 +153,13 @@
       integer :: idsurs                    ! Ursell number of the asymmetric wave
       integer :: idsrrw                    ! velocity skewness of the asymmetric wave
       integer :: idsbtw                    ! acceleration asymmetry parameter
+      integer :: idszrw                    ! Reference height to get near bottom current velocity
+      integer :: idsksd                    ! Bed roughness (zo) to calc. wave boundary layer
+      integer :: idsusc                    ! Current friction velocity at wave boundary layer
+      integer :: idstbl                    ! Thickness at wave boundary layer
+      integer :: idsubl                    ! Current velocity at wave boundary layer
+      integer :: idspwc                    ! Angle between waves/currents
+      integer :: idsfdw                    ! Friction factor from the current cycle
       integer :: idsucr                    ! Crest velocity of the asymmetric wave
       integer :: idsutr                    ! Trough velocity of the asymmetric wave
       integer :: idstcr                    ! Crest time period of the asymmetric wave
@@ -175,14 +174,17 @@
       real(r8), allocatable :: newlayer_thick(:)   ! deposit thickness criteria
       real(r8), allocatable :: minlayer_thick(:)   ! 2nd layer thickness criteria
       real(r8), allocatable :: bedload_coeff(:)    ! bedload rate coefficient
-      real(r8), allocatable :: sg_zwbl(:)         ! input elevation to get near-bottom current vel
+!
 #if defined BEDLOAD
+!# if defined BEDLOAD_VANDERA
+      real(r8), allocatable :: sg_zwbl(:)         ! input elevation to get near-bottom current vel
       real(r8), allocatable :: sedslope_crit_wet(:) ! critical wet bed slope for slumping
       real(r8), allocatable :: sedslope_crit_dry(:) ! critical dry bed slope for slumping
       real(r8), allocatable :: slopefac_wet(:)    ! bedload wet bed slumping factor
       real(r8), allocatable :: slopefac_dry(:)    ! bedload dry bed slumping factor
       real(r8), allocatable :: bedload_vandera_alphaw(:)    ! bedload scale factor for waves contribution
       real(r8), allocatable :: bedload_vandera_alphac(:)    ! bedload scale factor for currs contribution
+!# endif
 #endif
 !
       real(r8), allocatable :: Csed(:,:)       ! initial concentration
@@ -254,14 +256,14 @@
       ibtcr    = counter1    ! layer critical stress
 #endif
 #if defined SEDBIO_COUP
-      counter1 = counter1+1
-      iboxy    = counter1    ! porewater oxygen
-      counter1 = counter1+1
+      counter1 = counter1+1 
+      iboxy    = counter1    ! porewater oxygen 
+      counter1 = counter1+1 
       ibno3    = counter1    ! porewater nitrate
-      counter1 = counter1+1
+      counter1 = counter1+1 
       ibnh4    = counter1    ! porewater ammonium
-      counter1 = counter1+1
-      ibodu    = counter1    ! porewater oxygen demand units
+      counter1 = counter1+1 
+      ibodu    = counter1    ! porewater oxygen demand units 
 #endif
 !
 !-----------------------------------------------------------------------
@@ -304,22 +306,6 @@
       imaxD    = counter2    ! Maximum inundation depth.
       counter2 = counter2+1
       idnet    = counter2    ! Erosion/deposition
-      counter2 = counter2+1
-      idtbl    = counter2    ! Thickness at wave boundary layer
-      counter2 = counter2+1
-      idubl    = counter2    ! Current velocity at wave boundary layer
-      counter2 = counter2+1
-      idfdw    = counter2    ! Friction factor from the currents
-      counter2 = counter2+1
-      idzrw    = counter2    ! Reference height to get near bottom current velocity
-      counter2 = counter2+1
-      idksd    = counter2    ! Bed roughness (zo) to calc. wave boundary layer
-      counter2 = counter2+1
-      idusc    = counter2    ! Current friction velocity at wave boundary layer
-      counter2 = counter2+1
-      idpcx    = counter2    ! Anlge between currents and xi axis.
-      counter2 = counter2+1
-      idpwc    = counter2    ! Angle between waves/currents
 #if defined COHESIVE_BED || defined SED_BIODIFF || defined MIXED_BED
       counter2 = counter2+1
       idoff    = counter2    ! Offset for calculation of dmix erodibility profile (m).
@@ -406,13 +392,13 @@
         Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
       END IF
 #endif
+#if defined BEDLOAD
+!# if defined BEDLOAD_VANDERA
       IF (.not.allocated(sg_zwbl)) THEN
         allocate ( sg_zwbl(Ngrids) )
         sg_zwbl = 0.1_r8
         Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
       END IF
-#if defined BEDLOAD
-!# if defined BEDLOAD_VANDERA
       IF (.not.allocated(sedslope_crit_wet)) THEN
         allocate ( sedslope_crit_wet(Ngrids) )
         sedslope_crit_wet = IniVal
@@ -467,7 +453,7 @@
         Dbzp = IniVal
         Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
       END IF
-#endif
+#endif 
 #if defined COHESIVE_BED || defined MIXED_BED
       IF (.not.allocated(tcr_min)) THEN
         allocate ( tcr_min(Ngrids) )
diff --git a/ROMS/Nonlinear/Sediment/sediment_var.h b/ROMS/Nonlinear/Sediment/sediment_var.h
index 480180a70972d171a9cc4e615505a1396cc50f21..1da8f0e38eccad6516fa0fbcb575a6a37db7078b 100755
--- a/ROMS/Nonlinear/Sediment/sediment_var.h
+++ b/ROMS/Nonlinear/Sediment/sediment_var.h
@@ -28,13 +28,27 @@
                 idSbed(iporo)=varid
               CASE ('idSbed(idiff)')
                 idSbed(idiff)=varid
-#if defined BEDLOAD_VANDERA
+#if defined BEDLOAD && defined BEDLOAD_VANDERA
               CASE ('idsurs')
                 idsurs=varid
               CASE ('idsrrw')
                 idsrrw=varid
               CASE ('idsbtw')
                 idsbtw=varid
+              CASE ('idszrw')
+                idszrw=varid
+              CASE ('idsksd')
+                idsksd=varid
+              CASE ('idsusc')
+                idsusc=varid
+              CASE ('idstbl')
+                idstbl=varid
+              CASE ('idsubl')
+                idsubl=varid
+              CASE ('idspwc')
+                idspwc=varid
+              CASE ('idsfdw')
+                idsfdw=varid
               CASE ('idsucr')
                 idsucr=varid
               CASE ('idsutr')
@@ -43,7 +57,7 @@
                 idstcr=varid
               CASE ('idsttr')
                 idsttr=varid
-#endif
+#endif 
 #if defined COHESIVE_BED || defined SED_BIODIFF || defined MIXED_BED
               CASE ('idSbed(ibtcr)')
                 idSbed(ibtcr)=varid
@@ -94,22 +108,6 @@
                 idBott(imaxD)=varid
               CASE ('idBott(idnet)')
                 idBott(idnet)=varid
-              CASE ('idBott(idtbl)')
-                idBott(idtbl)=varid
-              CASE ('idBott(idubl)')
-                idBott(idubl)=varid
-              CASE ('idBott(idfdw)')
-                idBott(idfdw)=varid
-              CASE ('idBott(idzrw)')
-                idBott(idzrw)=varid
-              CASE ('idBott(idksd)')
-                idBott(idksd)=varid
-              CASE ('idBott(idusc)')
-                idBott(idusc)=varid
-              CASE ('idBott(idpcx)')
-                idBott(idpcx)=varid
-              CASE ('idBott(idpwc)')
-                idBott(idpwc)=varid
 #if defined COHESIVE_BED || defined SED_BIODIFF || defined MIXED_BED
               CASE ('idBott(idoff)')
                 idBott(idoff)=varid
diff --git a/ROMS/Nonlinear/Vegetation/vegarr_mod.h b/ROMS/Nonlinear/Vegetation/vegarr_mod.h
index 32bd03db408088e1517158cc251b9a9d4b9e0200..7f29e35f93a0b36c1b3297c6333f6cf11282f262 100755
--- a/ROMS/Nonlinear/Vegetation/vegarr_mod.h
+++ b/ROMS/Nonlinear/Vegetation/vegarr_mod.h
@@ -27,6 +27,7 @@
 !  step2d_vveg    Momentum term for 2d y direction                     !
 #ifdef VEG_FLEX
 !  bend           Bending for each vegetation                          !
+!  plant_hght_flex Effective plant height modified                     !
 !  Lveg           Effective blade length                               ! 
 # endif
 # ifdef VEG_FLEX
@@ -96,6 +97,7 @@
 # endif 
 # ifdef VEG_FLEX 
         real(r8), pointer :: bend(:,:,:)
+        real(r8), pointer :: plant_hght_flex(:,:)
 # endif         
 # ifdef VEG_TURB
         real(r8), pointer :: tke_veg(:,:,:)
@@ -195,6 +197,7 @@
       allocate ( VEG(ng) % Lveg(LBi:UBi,LBj:UBj,N(ng)) )
 #  ifdef VEG_FLEX
       allocate ( VEG(ng) % bend(LBi:UBi,LBj:UBj,NVEG) )
+      allocate ( VEG(ng) % plant_hght_flex(LBi:UBi,LBj:UBj) )
 #  endif
 #  ifdef VEG_HMIXING
       allocate ( VEG(ng) % visc2d_r_veg(LBi:UBi,LBj:UBj) )
@@ -376,6 +379,11 @@
             END DO 
           END DO 
         END DO 
+        DO j=Jmin,Jmax
+          DO i=Imin,Imax
+            VEG(ng) % plant_hght_flex(i,j) = IniVal
+          END DO 
+        END DO 
 #  endif 
 !  
 #  ifdef VEG_TURB 
diff --git a/ROMS/Nonlinear/Vegetation/vegetation_def_his.h b/ROMS/Nonlinear/Vegetation/vegetation_def_his.h
index b443ab285e7b7b70de0daba58cb986f89baeb76e..f86e0932c3ebac4bb2586c01b5982ed99e8a5fc6 100755
--- a/ROMS/Nonlinear/Vegetation/vegetation_def_his.h
+++ b/ROMS/Nonlinear/Vegetation/vegetation_def_his.h
@@ -37,6 +37,70 @@
         END IF
       END DO
 #endif
+#ifdef VEG_FLEX
+!
+!  Define flexible height variable.
+!
+        IF (Hout(idhgtf,ng)) THEN 
+          Vinfo( 1)=Vname(1,idhgtf)
+          Vinfo( 2)=Vname(2,idhgtf)
+          Vinfo( 3)=Vname(3,idhgtf)
+          Vinfo(14)=Vname(4,idhgtf)
+!         Vinfo(16)=Vname(1,idhgtf)
+          Vinfo(16)=Vname(1,idtime)
+# if defined WRITE_WATER && defined MASKING
+          Vinfo(20)='mask_rho'
+# endif
+          Vinfo(22)='coordinates'
+          Aval(5)=REAL(Iinfo(1,idhgtf,ng),r8)
+          status=def_var(ng, iNLM, HIS(ng)%ncid, HIS(ng)%Vid(idhgtf),   &
+     &                   NF_FOUT, nvd3, t2dgrd, Aval, Vinfo, ncname,    &
+     &                   SetFillVal = .FALSE.)
+          IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN
+
+        END IF
+#endif 
+#if defined VEG_TURB && defined VEG_TURB_WRITEHIS
+!
+!  Define turbulent kinetic energy due to veg.
+!
+        IF (Hout(idvtke,ng)) THEN
+          Vinfo( 1)=Vname(1,idvtke)
+          Vinfo( 2)=Vname(2,idvtke)
+          Vinfo( 3)=Vname(3,idvtke)
+          Vinfo(14)=Vname(4,idvtke)
+          Vinfo(16)=Vname(1,idtime)
+#  if defined WRITE_WATER && defined MASKING
+          Vinfo(20)='mask_rho'
+#  endif
+          Vinfo(22)='coordinates'
+          Aval(5)=REAL(Iinfo(1,idvtke,ng),r8)
+          status=def_var(ng, iNLM, HIS(ng)%ncid, HIS(ng)%Vid(idvtke),   &
+     &                   NF_FOUT, nvd4, w3dgrd, Aval, Vinfo, ncname,    &
+     &                   SetFillVal = .FALSE.)
+          IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN
+        END IF
+!
+!  Define length scale change due to veg.
+!
+        IF (Hout(idvgls,ng)) THEN
+          Vinfo( 1)=Vname(1,idvgls)
+          Vinfo( 2)=Vname(2,idvgls)
+          Vinfo( 3)=Vname(3,idvgls)
+          Vinfo(14)=Vname(4,idvgls)
+          Vinfo(16)=Vname(1,idtime)
+#  if defined WRITE_WATER && defined MASKING
+          Vinfo(20)='mask_rho'
+#  endif
+          Vinfo(22)='coordinates'
+          Aval(5)=REAL(Iinfo(1,idvgls,ng),r8)
+          status=def_var(ng, iNLM, HIS(ng)%ncid, HIS(ng)%Vid(idvgls),   &
+     &                   NF_FOUT, nvd4, w3dgrd, Aval, Vinfo, ncname,    &
+     &                   SetFillVal = .FALSE.)
+          IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN
+       ENDIF 
+#endif
+!
 #if defined VEG_STREAMING
 !
 !  Define wave dissipation due to vegetation.
diff --git a/ROMS/Nonlinear/Vegetation/vegetation_drag.F b/ROMS/Nonlinear/Vegetation/vegetation_drag.F
index 7d6355ab8757654fc56c4ec52ad9f07a0c8e7dd3..e281742054374e1402b2c37de3d4b6880e813992 100755
--- a/ROMS/Nonlinear/Vegetation/vegetation_drag.F
+++ b/ROMS/Nonlinear/Vegetation/vegetation_drag.F
@@ -61,6 +61,7 @@
      &                        VEG(ng) % plant,                          &
 # ifdef VEG_FLEX
      &                        VEG(ng) % bend,                           &
+     &                        VEG(ng) % plant_hght_flex,                &
 # endif
 # if defined MARSH_DYNAMICS && defined MARSH_RETREAT
      &                        VEG(ng) % marsh_mask,                     &
@@ -89,6 +90,7 @@
      &                              plant,                              &
 # ifdef VEG_FLEX
      &                              bend,                               &
+     &                              plant_hght_flex,                    &
 # endif
 # if defined MARSH_DYNAMICS && defined MARSH_RETREAT
      &                              marsh_mask,                         &
@@ -103,6 +105,7 @@
       USE mod_scalars
       USE mod_vegetation
       USE mod_vegarr
+      USE bc_2d_mod, ONLY: bc_r2d_tile
       USE bc_3d_mod, ONLY: bc_r3d_tile
       USE exchange_2d_mod
       USE exchange_3d_mod
@@ -125,6 +128,7 @@
       real(r8), intent(inout) :: plant(LBi:,LBj:,:,:)
 #  ifdef VEG_FLEX
       real(r8), intent(inout) :: bend(LBi:,LBj:,:)
+      real(r8), intent(inout) :: plant_hght_flex(LBi:,LBj:)
 #  endif
 # if defined MARSH_DYNAMICS && defined MARSH_RETREAT
       real(r8), intent(inout) :: marsh_mask(LBi:,LBj:)
@@ -143,6 +147,7 @@
       real(r8), intent(inout) :: plant(LBi:UBi,LBj:UBj,NVEG,NVEGP)
 #  ifdef VEG_FLEX
       real(r8), intent(inout) :: bend(LBi:UBi,LBj:UBj,N(ng),NVEG)
+      real(r8), intent(inout) :: plant_hght_flex(LBi:UBi,LBj:UBj)
 #  endif
 # if defined MARSH_DYNAMICS && defined MARSH_RETREAT 
       real(r8), intent(inout) :: marsh_mask(LBi:UBi,LBj:UBj)
@@ -191,7 +196,7 @@
 !  to zero. The value of 0.75 is arbitrary limitation assigment
 !  (same as for bottom stress).
 !
-      cff=0.75_r8/dt(ng)
+      cff=0.25_r8/dt(ng)
 # endif
 !
 # if defined MARSH_DYNAMICS && defined MARSH_RETREAT 
@@ -267,9 +272,10 @@
               cflex=1.0_r8-((1.0_r8-0.9_r8*cff6)/(cff5))
               cflex=MIN(cflex, 1.0_r8)
 !
-! Effective blade length
+! Effective blade length and save it
 !
               plant_height_eff=cflex*plant(i,j,iveg,phght)
+              plant_hght_flex(i,j)=plant_height_eff
 !
 ! Blade bending angle
 !
@@ -391,12 +397,22 @@
      &                          rv_veg(:,:,:))
       END IF
 !
-#  if defined MARSH_DYNAMICS && defined MARSH_RETREAT 
+#  if defined MARSH_RETREAT
       DO ivpr=1,NVEGP
         CALL bc_r3d_tile (ng, tile,                                     &
      &                          LBi, UBi, LBj, UBj, 1, NVEG,            &
      &                          plant(:,:,:,ivpr))
       END DO
+#  endif 
+#  if defined VEG_FLEX 
+      CALL bc_r2d_tile (ng, tile,                                       &
+     &                          LBi, UBi, LBj, UBj,                     &
+                                plant_hght_flex)
+      CALL mp_exchange2d (ng, tile, iNLM,  1,                           &
+     &                    LBi, UBi, LBj, UBj,                           &
+     &                    NghostPoints,                                 &
+     &                    EWperiodic(ng), NSperiodic(ng),               &
+     &                    plant_hght_flex)
 #  endif
 !
 # ifdef DISTRIBUTE
@@ -414,7 +430,7 @@
      &                    EWperiodic(ng), NSperiodic(ng),               &
      &                    ru_veg, rv_veg)
 ! 
-#  if defined MARSH_DYNAMICS && defined MARSH_RETREAT 
+#  if defined MARSH_RETREAT || defined VEG_FLEX
       CALL mp_exchange4d (ng, tile, iNLM, 1,                            &
      &                    LBi, UBi, LBj, UBj, 1, NVEG, 1, NVEGP,        &
      &                    NghostPoints,                                 &
diff --git a/ROMS/Nonlinear/Vegetation/vegetation_inp.h b/ROMS/Nonlinear/Vegetation/vegetation_inp.h
index e71c98e0268b8851d98fd424aa2633038479638a..d7a230d19ba8fd842de71060c405a65bd5150420 100755
--- a/ROMS/Nonlinear/Vegetation/vegetation_inp.h
+++ b/ROMS/Nonlinear/Vegetation/vegetation_inp.h
@@ -266,6 +266,31 @@
               END IF
               Npts=load_l(Nval, Cval, Ngrids, Hout(idvprp(pthck),:))
 #endif                                         
+#if defined VEG_TURB && defined VEG_TURB_WRITEHIS 
+            CASE ('Hout(idvtke)')
+              IF (idvtke.eq.0) THEN
+                IF (Master) WRITE (out,30) 'idvtke'
+                exit_flag=5
+                RETURN
+              END IF
+              Npts=load_l(Nval, Cval, Ngrids, Hout(idvtke,:))
+            CASE ('Hout(idvgls)')
+              IF (idvgls.eq.0) THEN
+                IF (Master) WRITE (out,30) 'idvgls'
+                exit_flag=5
+                RETURN
+              END IF
+              Npts=load_l(Nval, Cval, Ngrids, Hout(idvgls,:))
+#endif
+#if defined VEG_FLEX 
+            CASE ('Hout(idhgtf)')
+              IF (idhgtf.eq.0) THEN
+                IF (Master) WRITE (out,30) 'idhgtf'
+                exit_flag=5
+                RETURN
+              END IF
+              Npts=load_l(Nval, Cval, Ngrids, Hout(idhgtf,:))
+#endif
 #ifdef VEG_STREAMING
             CASE ('Hout(idWdvg)')
               IF ((idWdvg).eq.0) THEN
diff --git a/ROMS/Nonlinear/Vegetation/vegetation_mod.h b/ROMS/Nonlinear/Vegetation/vegetation_mod.h
index c0cb13ccd5a5dd544f3293d726d569e09f9e370b..4f0bb31866310d2d2ddacd184dbe4c6020a74dbf 100755
--- a/ROMS/Nonlinear/Vegetation/vegetation_mod.h
+++ b/ROMS/Nonlinear/Vegetation/vegetation_mod.h
@@ -62,6 +62,9 @@
 !  ipdwbm         Id to output below ground biomass                    !
 !  idWdvg         Id to output wave dissipation from vegetation        !
 !  idCdvg         Id to output spectral Cd from waves vegetation       !
+!  idvgls         Id to output length scale change due to veg.         !
+!  idvtke         Id to output tke change due to veg.                  !
+!  idhgtf         Id to output height change due to flex.              !
 !                                                                      !
 !  Marsh wave induced erosion Output:                                  !
 !  ==========================                                          !
@@ -90,10 +93,16 @@
       integer :: phght, pdens, pdiam, pthck
       integer :: ipdens, iphght, ipdiam, ipthck
 #endif
+#if defined VEG_TURB  && defined VEG_TURB_WRITEHIS  
+      integer :: idvgls, idvtke
+#endif
 !
 #ifdef VEG_STREAMING 
       integer :: idWdvg, idCdvg
 #endif 
+#ifdef VEG_FLEX
+      integer :: idhgtf
+#endif 
 #if defined VEG_DRAG || defined VEG_BIOMASS  
       integer, allocatable :: idvprp(:)
 #endif
diff --git a/ROMS/Nonlinear/Vegetation/vegetation_turb_cal.F b/ROMS/Nonlinear/Vegetation/vegetation_turb_cal.F
index 0fd89e29076c341a638cba80462a0a09654b722c..2973d482b0808f3b0ea88df6445a4096b8dc559e 100755
--- a/ROMS/Nonlinear/Vegetation/vegetation_turb_cal.F
+++ b/ROMS/Nonlinear/Vegetation/vegetation_turb_cal.F
@@ -99,6 +99,13 @@
       USE mod_vegetation 
       USE mod_vegarr
       USE vegetation_drag_mod
+      USE exchange_3d_mod, ONLY : exchange_w3d_tile
+# ifdef DISTRIBUTE
+      USE mp_exchange_mod, ONLY : mp_exchange3d
+# endif
+
+
+!      USE tkebc_mod, ONLY : tkebc_tile
 !
 !  Imported variable declarations.
 !
@@ -251,7 +258,36 @@
           END DO
         END DO 
       END DO VEG_LOOP
+
+# if defined VEG_TURB && defined VEG_TURB_WRITEHIS 
 !
+!-----------------------------------------------------------------------
+!  Set lateral boundary conditions.
+!-----------------------------------------------------------------------
+!
+!      CALL tkebc_tile (ng, tile,                                        &
+!     &                 LBi, UBi, LBj, UBj, N(ng),                       &
+!     &                 IminS, ImaxS, JminS, JmaxS,                      &
+!     &                 nnew, nstp,                                      &
+!     &                 gls_veg, tke_veg)
+
+      IF (EWperiodic(ng).or.NSperiodic(ng)) THEN
+        CALL exchange_w3d_tile (ng, tile,                               &
+     &                          LBi, UBi, LBj, UBj, 0, N(ng),           &
+     &                          tke_veg(:,:,:))
+        CALL exchange_w3d_tile (ng, tile,                               &
+     &                          LBi, UBi, LBj, UBj, 0, N(ng),           &
+     &                          gls_veg(:,:,:))
+      END IF
+
+#  ifdef DISTRIBUTE
+      CALL mp_exchange3d (ng, tile, iNLM, 2,                            &
+     &                    LBi, UBi, LBj, UBj, 0, N(ng),                 &
+     &                    NghostPoints,                                 &
+     &                    EWperiodic(ng), NSperiodic(ng),               &
+     &                    tke_veg(:,:,:),    gls_veg(:,:,:))
+#  endif 
+# endif               
       RETURN
       END SUBROUTINE vegetation_turb_tile
 #endif
diff --git a/ROMS/Nonlinear/Vegetation/vegetation_var.h b/ROMS/Nonlinear/Vegetation/vegetation_var.h
index 433cc4c5927001332ec9c83da40ad5bef4c97b67..d7e7d2c73ad63b53d79c99b7b02f396edcf8803f 100755
--- a/ROMS/Nonlinear/Vegetation/vegetation_var.h
+++ b/ROMS/Nonlinear/Vegetation/vegetation_var.h
@@ -34,12 +34,23 @@
 !                idvprp(pbgbm)=varid
 !#endif
 #endif 
+#if defined VEG_TURB && defined VEG_TURB_WRITEHIS
+            CASE ('idvtke')
+              idvtke=varid
+            CASE ('idvgls')
+              idvgls=varid
+#endif
 #if defined VEG_STREAMING 
             CASE ('idWdvg')
               idWdvg=varid
             CASE ('idCdvg')
               idCdvg=varid
 #endif 
+#if defined VEG_FLEX
+            CASE ('idhgtf')
+              idhgtf=varid
+#endif 
+	       
 ! 
 #if defined MARSH_DYNAMICS
             CASE ('idTims')
diff --git a/ROMS/Nonlinear/Vegetation/vegetation_wrt_his.h b/ROMS/Nonlinear/Vegetation/vegetation_wrt_his.h
index 867242d5a50f60477fed56b21d2494f552ca45fa..17de63c24b41c84a150e8d4761a44123c6d6271a 100755
--- a/ROMS/Nonlinear/Vegetation/vegetation_wrt_his.h
+++ b/ROMS/Nonlinear/Vegetation/vegetation_wrt_his.h
@@ -42,13 +42,37 @@
       END DO
 # endif
 !
+# ifdef VEG_FLEX 
+!
+!  Write out height modification due to flexibility.
+!
+      IF (Hout(idhgtf,ng)) THEN
+        scale=1.0_r8
+        gtype=gfactor*r2dvar
+        status=nf_fwrite2d(ng, iNLM, HIS(ng)%ncid, HIS(ng)%Vid(idhgtf), &
+     &                     HIS(ng)%Rindex, gtype,                       &
+     &                     LBi, UBi, LBj, UBj, scale,                   &
+# ifdef MASKING
+     &                     GRID(ng) % rmask,                            &
+# endif
+     &                     VEG(ng)%plant_hght_flex)
+        IF (FoundError(status, nf90_noerr, __LINE__, MyFile)) THEN
+          IF (Master) THEN
+            WRITE (stdout,10) TRIM(Vname(1,idhgtf)), HIS(ng)%Rindex
+          END IF
+          exit_flag=3
+          ioerror=status
+          RETURN
+        END IF
+      END IF
+# endif 
+!
 # ifdef VEG_STREAMING
 !
 !  Write out wave dissipation due to vegetation
 !
       IF (Hout(idWdvg,ng)) THEN
-        scale=rho0                          ! W m /kg to W/m2
-!       scale=1.0_r8
+        scale=1.0_r8
         gtype=gfactor*r2dvar
         status=nf_fwrite2d(ng, iNLM, HIS(ng)%ncid, HIS(ng)%Vid(idWdvg), &
      &                     HIS(ng)%Rindex, gtype,                       &
@@ -69,7 +93,7 @@
 !
 !  Write out spectral Cd due to vegetation.
 !
-      IF (Hout(idCdvg,ng)) THEN
+      IF (Hout(idCdvg,ng)) THEN 
         scale=1.0_r8
         gtype=gfactor*r2dvar
         status=nf_fwrite2d(ng, iNLM, HIS(ng)%ncid, HIS(ng)%Vid(idCdvg), &
@@ -80,7 +104,7 @@
 # endif
      &                     VEG(ng)%Cdwave_veg)
         IF (FoundError(status, nf90_noerr, __LINE__, MyFile)) THEN
-          IF (Master) THEN
+          IF (Master) THEN 
             WRITE (stdout,10) TRIM(Vname(1,idCdvg)), HIS(ng)%Rindex
           END IF
           exit_flag=3
@@ -89,7 +113,56 @@
         END IF
       END IF
 # endif
+# if defined VEG_TURB && defined VEG_TURB_WRITEHIS 
+!
+!  Write out turbulent kinetic energy modification due to veg.
 !
+      IF (Hout(idvtke,ng)) THEN
+        scale=1.0_dp
+        gtype=gfactor*w3dvar
+        status=nf_fwrite3d(ng, iNLM, HIS(ng)%ncid, HIS(ng)%Vid(idvtke), &
+     &                     HIS(ng)%Rindex, gtype,                       &
+     &                     LBi, UBi, LBj, UBj, 0, N(ng), scale,         &
+#  ifdef MASKING
+     &                     GRID(ng) % rmask_full,                       &
+#  endif
+     &                     VEG(ng) % tke_veg(:,:,:),                    &
+     &                     SetFillVal = .FALSE.)
+        IF (FoundError(status, nf90_noerr, __LINE__, MyFile)) THEN
+          IF (Master) THEN
+            WRITE (stdout,10) TRIM(Vname(1,idvtke)), HIS(ng)%Rindex
+          END IF
+          exit_flag=3
+          ioerror=status
+          RETURN
+        END IF
+      END IF
+!
+!  Write out turbulent length scale field due to veg. 
+!
+      IF (Hout(idvgls,ng)) THEN
+        scale=1.0_dp
+        gtype=gfactor*w3dvar
+        status=nf_fwrite3d(ng, iNLM, HIS(ng)%ncid, HIS(ng)%Vid(idvgls), &
+     &                     HIS(ng)%Rindex, gtype,                       &
+     &                     LBi, UBi, LBj, UBj, 0, N(ng), scale,         &
+#  ifdef MASKING
+     &                     GRID(ng) % rmask_full,                       &
+#  endif
+     &                     VEG(ng) % gls_veg(:,:,:),                    &
+     &                     SetFillVal = .FALSE.)
+        IF (FoundError(status, nf90_noerr, __LINE__, MyFile)) THEN
+          IF (Master) THEN
+            WRITE (stdout,10) TRIM(Vname(1,idvgls)), HIS(ng)%Rindex
+          END IF
+          exit_flag=3
+          ioerror=status
+          RETURN
+        END IF
+      END IF
+!
+# endif 
+! 
 # ifdef MARSH_DYNAMICS
 !
 !  Write out masking for marsh cells.
@@ -216,9 +289,9 @@
           RETURN
         END IF
       END IF
-#  endif
+#  endif 
 !
-#  if  defined MARSH_VERT_GROWTH
+#  if  defined MARSH_VERT_GROWTH 
 !
 !  Write mean high high water over a given frequency.
 !
@@ -332,5 +405,5 @@
           RETURN
         END IF
       END IF
-#  endif
-# endif
+#  endif  
+# endif  
diff --git a/ROMS/Nonlinear/bulk_flux.F b/ROMS/Nonlinear/bulk_flux.F
index 280ec74d8a68e1e06faf9cadedb3cca98e17eef1..61d56c9dc256e9f56cfbb4cd2df769b59b3df2a3 100755
--- a/ROMS/Nonlinear/bulk_flux.F
+++ b/ROMS/Nonlinear/bulk_flux.F
@@ -906,14 +906,12 @@
 #  ifdef WET_DRY
           evap(i,j)=evap(i,j)*rmask_wet(i,j)
 #  endif
-#  ifdef SALINITY
           stflux(i,j,isalt)=cff*(evap(i,j)-rain(i,j))
-#   ifdef MASKING
+#  ifdef MASKING
           stflux(i,j,isalt)=stflux(i,j,isalt)*rmask(i,j)
-#   endif
-#   ifdef WET_DRY
+#  endif
+#  ifdef WET_DRY
           stflux(i,j,isalt)=stflux(i,j,isalt)*rmask_wet(i,j)
-#   endif
 #  endif
 # endif
         END DO
@@ -966,11 +964,9 @@
         CALL exchange_r2d_tile (ng, tile,                               &
      &                          LBi, UBi, LBj, UBj,                     &
      &                          evap)
-#  ifdef SALINITY
         CALL exchange_r2d_tile (ng, tile,                               &
      &                          LBi, UBi, LBj, UBj,                     &
      &                          stflux(:,:,isalt))
-#  endif
 # endif
         CALL exchange_u2d_tile (ng, tile,                               &
      &                          LBi, UBi, LBj, UBj,                     &
@@ -988,18 +984,12 @@
      &                    lrflx, lhflx, shflx,                          &
      &                    stflux(:,:,itemp))
 #  ifdef EMINUSP
-      CALL mp_exchange2d (ng, tile, iNLM, 1,                            &
-     &                    LBi, UBi, LBj, UBj,                           &
-     &                    NghostPoints,                                 &
-     &                    EWperiodic(ng), NSperiodic(ng),               &
-     &                    evap)
-#   ifdef SALINITY
-      CALL mp_exchange2d (ng, tile, iNLM, 1,                            &
+      CALL mp_exchange2d (ng, tile, iNLM, 2,                            &
      &                    LBi, UBi, LBj, UBj,                           &
      &                    NghostPoints,                                 &
      &                    EWperiodic(ng), NSperiodic(ng),               &
+     &                    evap,                                         &
      &                    stflux(:,:,isalt))
-#   endif
 #  endif
       CALL mp_exchange2d (ng, tile, iNLM, 2,                            &
      &                    LBi, UBi, LBj, UBj,                           &
diff --git a/SWAN/Src/swancom1.F b/SWAN/Src/swancom1.F
index f6ea5297ec515f597b9f5f079517f8e951a78916..f636181ba1a0577d99869c7b611ba5f6041f48a1 100755
--- a/SWAN/Src/swancom1.F
+++ b/SWAN/Src/swancom1.F
@@ -3219,7 +3219,7 @@
      &  XIS                 ,COMPDA(1,JFRC2)     ,IT                  ,   40.00
      &  COMPDA(1,JNPLA2)    ,COMPDA(1,JNPHT2)                         ,   40.59 40.35
      &  COMPDA(1,JNPDA2)    ,COMPDA(1,JNPTK2)                         ,   40.59 40.35
-     &  COMPDA(1,JDSXC)     ,                                             40.59 
+     &  COMPDA(1,JDSXC)     ,COMPDA(1,JDSKC)     ,COMPDA(1,JDQQQ)     ,
      &  COMPDA(1,JTURB2)    ,COMPDA(1,JMUDL2)                         ,   40.59 40.35 40.55
      &  COMPDA(1,JAICE2)    ,COMPDA(1,JHICE2)    ,                        41.75
      &  COMPDA(1,JURSEL)    ,LSWMAT(1,1,JABIN)   ,REFLSO              ,   40.41 40.03
@@ -3366,6 +3366,7 @@
      &                  COMPDA(1,JDSXS)    ,                              40.67 40.61
      &                  COMPDA(1,JDSXW)    ,                              40.67 40.61
      &                  COMPDA(1,JDSXV)                        ,          40.35 
+!     &                  COMPDA(1,JDSXV)    ,COMPDA(1,JDSXC)    ,          40.35 40.67 40.61
      &                  COMPDA(1,JDSXT)    ,COMPDA(1,JDSXM)    ,          40.35 40.67 40.61
      &                  COMPDA(1,JDSXI)    ,                              41.75
      &                  COMPDA(1,JDSXL)    ,                              40.88
@@ -3672,6 +3673,15 @@
       IF (IVEG.EQ.1) THEN                                                 40.55
          WRITE(PRINTF,7006)                                               40.55
  7006    FORMAT(' Vegetation due to Dalrymple (1984)')                    40.55
+      ELSEIF (IVEG.EQ.2) THEN                                             41.77
+         WRITE(PRINTF,7009)                                               41.77
+ 7009    FORMAT(' Vegetation due to Jacobsen et al. (2019)')              41.77
+      ELSEIF (IVEG.EQ.3) THEN                                             41.77
+         WRITE(PRINTF,7010)                                               41.77
+ 7010    FORMAT(' Vegetation-COAWST due to Dalrymple (1984)')             41.77
+      ELSEIF (IVEG.EQ.4) THEN                                             41.77
+         WRITE(PRINTF,7011)                                               41.77
+ 7011    FORMAT(' Vegetation-COAWST due to Jacobsen et al.(2019)')        41.77
       ELSE                                                                40.55
          WRITE (PRINTF, *) 'Vegetation is off'                            40.55
       ENDIF                                                               40.55
@@ -6612,7 +6622,7 @@
      &                     XIS        ,FRCOEF     ,IT         ,           30.00
      &                     NPLA2      ,NPHT2                  ,           40.59 40.35 40.55
      &                     NPDA2      ,NPTK2                  ,           40.59 40.35 40.55
-     &                     CDVEG_COAWST                       ,           40.59  
+     &                     CDVEG_COAWST, KC_COAWST, Q_COAWST  ,           40.59  
      &                     TURBV2     ,MUDL2                  ,           40.59 40.35 40.55
      &                     AICE2      ,HICE2                  ,           41.75
      &                     URSELL     ,ANYBIN     ,REFLSO     ,           40.41 40.03
@@ -6895,7 +6905,9 @@
 !   If SVEG is on (IVEG == 1 ) then,
 !     Call SVEG  to compute the source term due to vegetation
 !      dissipation according to Dalrymple (1984)
-!   If SVEG_COAWST is on (IVEG == 2 ) then,
+!   If SVEG is on (IVEG == 2 ) then,
+!      dissipation according to Jacobsen et al. (2019)
+!   If SVEG_COAWST is on (IVEG == 3 ) then,
 !     Call SVEG_COAWST spatially varying veg to compute the source 
 !      term due to vegetation
 !      dissipation according to Dalrymple (1984)
@@ -7026,7 +7038,7 @@
       REAL  :: MUDL2(MCGRD)                                               40.59
       REAL, INTENT(IN)  :: AICE2(MCGRD)                                   41.75
       REAL, INTENT(IN)  :: HICE2(MCGRD)                                   41.75
-      REAL  :: CDVEG_COAWST(MCGRD)                                        41.75
+      REAL  :: CDVEG_COAWST(MCGRD), KC_COAWST(MCGRD), Q_COAWST(MCGRD)     41.75
       REAL  :: REFLSO(MDC,MSC)                                            40.41
       REAL  :: URMSTOP(MCGRD)                                             40.88
 !
@@ -7104,29 +7116,34 @@
 !TIMG      CALL SWTSTO(138)                                                    40.59
 !
 !TIMG      CALL SWTSTA(139)                                                    40.55
-      IF ( IVEG.EQ.1 ) THEN                                               40.55
-!
-!     *** wave-vegetation interactions ***                                40.55
+      IF ( IVEG.EQ.1.OR.IVEG.EQ.2 ) THEN                                  40.55
 !
+!        *** wave-vegetation interactions ***                             40.55
+!        iveg=1
 !        *** energy dissipation according to Dalrymple (1984)             40.55
+!        iveg=2
+!        *** energy dissipation according to JACOBSEN (2019)              40.55
 !
          CALL SVEG (DEP2   ,IMATDA   ,ETOT   ,SMEBRK    ,                 40.58 40.55
-     &              KMESPC ,PLVEGT                      ,                 40.58 40.55
+     &              KWAVE  ,KMESPC   ,PLVEGT ,                            41.77 40.58 40.55
      &              IDCMIN ,IDCMAX   ,ISSTOP ,DISSC1    ,                 40.55
      &              NPLA2  )                                              40.67 40.61 40.55
       END IF                                                              40.55
-      IF ( IVEG.EQ.2 ) THEN                                               40.55
+      IF ( IVEG.EQ.3.OR.IVEG.EQ.4 ) THEN                                  40.55
 !
 !     *** wave-vegetation interactions ***                                40.55
-!
+!         with spatially varying veg and spectral Cd
+!        iveg=3
 !        *** energy dissipation according to Dalrymple (1984)             40.55
+!        iveg=4
+!        *** energy dissipation according to JACOBSEN (2019)              40.55
 !
          CALL SVEG_COAWST (DEP2   ,IMATDA   ,ETOT   ,SMEBRK    ,          40.58 40.55
-     &              KMESPC ,PLVEGT                             ,          40.58 40.55
+     &              KWAVE,  KMESPC   ,PLVEGT                   ,          40.58 40.55
      &              IDCMIN ,IDCMAX   ,ISSTOP ,DISSC1           ,          40.55
-     &              ABRBOT ,KWAVE                              ,          40.55
+     &              ABRBOT                                     ,          40.55
      &              NPLA2  ,NPHT2    ,NPDA2  ,NPTK2            ,          40.67 40.61 40.55
-     &              CDVEG_COAWST)
+     &              CDVEG_COAWST, KC_COAWST, Q_COAWST)
       END IF                                                              40.55
 !TIMG      CALL SWTSTO(139)                                                    40.55
 !
diff --git a/SWAN/Src/swancom2.F b/SWAN/Src/swancom2.F
index 02346e03129784236ca39df90663bccd27f39189..537d3fbba90d71494802d9166f24ad8855525f8a 100755
--- a/SWAN/Src/swancom2.F
+++ b/SWAN/Src/swancom2.F
@@ -553,7 +553,7 @@
 !****************************************************************
 !
       SUBROUTINE SVEG ( DEP2   ,IMATDA   ,ETOT   ,SMEBRK    ,
-     &                  KMESPC ,PLVEGT   ,
+     &                  KWAVE  ,KMESPC   ,PLVEGT ,
      &                  IDCMIN ,IDCMAX   ,ISSTOP ,DISSC1    ,
      &                  NPLA2  )
 !
@@ -668,6 +668,69 @@
 !
 !     and is linearized by means of the Picard iteration
 !
+!     IVEG = 2
+!     Jacobsen et al. (2019) published a frequency-dependent
+!     dissipation model for waves propagating over a canopy
+!
+!     The dissipation per frequency is
+!
+!                              ----------
+!                              | 2 mu,0
+!     dv_spectral =   2 F Su   | -------
+!                              |   pi
+!                             \|
+!
+!     where
+!
+!     F = 1/2 * rho * alpha^3 * Cd * bv * Nv
+!
+!     The velocity spectrum Su is associated with the
+!     water surface spectral distribution Sn, as given by
+!
+!                                 2
+!               w * cosh k(h+z)
+!     Su =  (  ----------------- )   * Sn
+!                   sinh kh
+!
+!     where
+!
+!     dv_spectral =   dissipation per frequency and layer [kg^2/ms]
+!     F           =   canopy factor                        [kg/m^4]
+!     Su          =   velocity spectrum                     [m^2/s]
+!     mu,0        =   first moment velocity spectrum      [m^2/s^2]
+!     rho         =   water density                        [kg/m^3]
+!     Nv          =   vegetation density                [stems/m^2]
+!     bv          =   stem thickness                            [m]
+!     Cd          =   drag coefficient vegetation               [-]
+!     alpha = 1   =   velocity reduction factor by canopies     [-]
+!     w           =   radian frequency (per component)        [1/s]
+!     k           =   wave number (per component)             [1/m]
+!     z           =   elevation of still water level            [m]
+!                     (z = 0 at water surface)
+!     h           =   water depth                               [m]
+!     Sn          =   variance density spectrum             [m^2*s]
+!
+!     The depth-integrated and frequency-dependent dissipation is then
+!
+!
+!                   -h + ah
+!                   |\
+!                   |
+!     Sveg  =   -   |      dv_spectral dz / rho / g
+!                   |
+!                  \|
+!                   -h
+!
+!     where
+!
+!     ah =  canopy height (under water)                         [m]
+!
+!     The vertical integration is approximated using the Simpson's rule.
+!     A minimum number of integration points would be needed to reduce
+!     the error of the approximation; 21 points appeared to be sufficient
+!
+!     Note: current effects are not included in Jacobsen et al. (2019)
+
 !
 !  4. Argument variables
 !
@@ -680,6 +743,7 @@
 !     ISSTOP      maximum counter of wave component in frequency
 !                 space that is propagated
 !     KMESPC      mean average wavenumber according to the WAM-formulation
+!     KWAVE       wave number
 !     NPLA2       number of plants per square meter (depth-averaged)
 !     PLVEGT      array containing the vegetation source term for test-output
 !     SMEBRK      mean frequency according to first order moment
@@ -687,38 +751,57 @@
       INTEGER ISSTOP, IDCMIN(MSC), IDCMAX(MSC)
       REAL    DEP2(MCGRD)          ,
      &        IMATDA(MDC,MSC)      ,
+     &        KWAVE(MSC,MICMAX)    ,
      &        DISSC1(MDC,MSC,MDISP),
      &        PLVEGT(MDC,MSC,NPTST),
      &        NPLA2 (MCGRD)
       REAL    ETOT, SMEBRK, KMESPC
 !
+!  5. Parameter variables
+!
+!     ALFU        (orbital) velocity reduction factor by canopies
+!     NIP         total number of subintervals used by Simpson's rule
+!
+      INTEGER, PARAMETER :: NIP  = 20
+      REAL   , PARAMETER :: ALFU = 1.
+!
 !  6. Local variables
 !
 !     A     :     auxiliary variable
 !     B     :     auxiliary variable
 !     C     :     auxiliary variable
 !     D     :     auxiliary variable
+!     DCIP  :     frequency-dependent dissipation for each integration point
+!     DZ    :     interval for vertical integration
+!     EKZ   :     exponential of k(h+z)
+!     FDD   :     factor with orbital velocity to determine Su from Sn
 !     ID    :     counter of the spectral direction
 !     IDDUM :     counter
 !     IENT  :     number of entries
 !     IK    :     counter
 !     IL    :     counter
 !     IS    :     counter of relative frequency band
+!     KC    :     wavenumber times canopy layer height
 !     KD    :     wavenumber times water depth
 !     KVEGH :     wavenumber times plant height
+!     KZ    :     wavenumber times z
 !     LAYPRT:     part of layer below water level
+!     MU    :     first order moment of velocity spectrum
 !     SINHK :     sinh(kh)
 !     SLAYH :     total sum of layer thicknesses
 !     SLAYH1:     sum of layer thicknesses below water level
-!     SLAYH2:     sum of layer thicknesses below water level
-!     SVEG1 :     layer-independent dissipation factor
 !     SVEG2 :     total sum of dissipation factor over layers
 !     SVEGET:     source term containing dissipation due to vegetation
 !                 to be stored in the array IMATDA
+!     ZDH   :     vertical point z between -depth to 0 or - value water clearance in water column
+!     ZH    :     cumulative layer thickness for velocities, bottom up
+!                 (z+d between 0 and vegetation height or depth value in water column)
 !
       INTEGER ID, IDDUM, IENT, IK, IL, IS
       REAL    A, B, C, D, KD, KVEGH, LAYPRT, SINHK, SLAYH,
-     &        SLAYH1, SLAYH2, SVEG1, SVEG2, SVEGET
+     &        SLAYH1, SLAYH2, SVEG1, SVEG2
+      REAL    DZ, EKZ, KC, KZ, MU, ZDH, ZH
+      REAL    DCIP(0:NIP,MSC), FDD(MSC), SVEGET(MSC)
 !
 !  9. Subroutines calling
 !
@@ -756,61 +839,69 @@
       DATA IENT/0/
       IF (LTRACE) CALL STRACE (IENT,'SVEG')
 
-!     --- compute layer-independent vegetation dissipation factor
-
-      KD    = KMESPC * DEP2(KCGRD(1))
-      IF ( KD.GT.10. ) RETURN
-      C     = 3.*KMESPC*(COSH(KD))**3
-      SVEG1 = SQRT(2./PI) * GRAV**2 * (KMESPC/SMEBRK)**3 * SQRT(ETOT)/ C
-      IF ( VARNPL ) SVEG1 = SVEG1 * NPLA2(KCGRD(1))
-
 !     --- compute dissipation factor for each layer and summed up
 
       SLAYH = 0.
       DO IL = 1, ILMAX
          SLAYH = SLAYH + LAYH(IL)
       ENDDO
+      IF ( IVEG.EQ.1 ) THEN
+!     --- Suzuki et al. (2011)
+!     --- compute layer-independent vegetation dissipation factor
+         KD    = KMESPC * DEP2(KCGRD(1))
+         IF ( KD.GT.10. ) RETURN
+         C     = 3.*KMESPC*(COSH(KD))**3
+         SVEG1 = SQRT(2./PI)*GRAV**2 * (KMESPC/SMEBRK)**3 * SQRT(ETOT)/C
+         IF ( VARNPL ) SVEG1 = SVEG1 * NPLA2(KCGRD(1))
 
-      KVEGH = 0.
-      C     = 0.
-      D     = 0.
-      SVEG2 = 0.
+         KVEGH = 0.
+         C     = 0.
+         D     = 0.
+         SVEG2 = 0.
 
-      IF ( DEP2(KCGRD(1)).GT.SLAYH ) THEN
+         IF ( DEP2(KCGRD(1)).GT.SLAYH ) THEN
 
-         DO IL = 1, ILMAX
-            KVEGH = KVEGH + KMESPC * LAYH(IL)
-            SINHK = SINH(KVEGH)
-            A     = C
-            B     = D
-            C     = SINHK**3
-            D     = 3.*SINHK
-            A     = C - A
-            B     = D - B
-            SVEG2 = SVEG2 + VEGDRL(IL)*VEGDIL(IL)*VEGNSL(IL)*(A + B)
-         END DO
+            DO IL = 1, ILMAX
+               KVEGH = KVEGH + KMESPC * LAYH(IL)
+               SINHK = SINH(KVEGH)
+               A     = C
+               B     = D
+               A     = C - A
+               B     = D - B
+               SVEG2 = SVEG2 + VEGDRL(IL)*VEGDIL(IL)*VEGNSL(IL)*(A + B)
+            END DO
 
-      ELSE IF ( DEP2(KCGRD(1)).LT.LAYH(1) ) THEN
+         ELSE IF ( DEP2(KCGRD(1)).LT.LAYH(1) ) THEN
 
-         SINHK = SINH(KD)
-         A     = SINHK**3
-         B     = 3.*SINHK
-         SVEG2 = VEGDRL(1)*VEGDIL(1)*VEGNSL(1)*(A + B)
+            SINHK = SINH(KD)
+            A     = SINHK**3
+            B     = 3.*SINHK
+            SVEG2 = VEGDRL(1)*VEGDIL(1)*VEGNSL(1)*(A + B)
 
-      ELSE
+         ELSE
 
-         SLAYH1 = 0.
-         SLAYH2 = 0.
-         LAYPRT = 0.
-         VGLOOP : DO IL = 1, ILMAX
-            SLAYH1 = SLAYH1 + LAYH(IL)
-            IF (DEP2(KCGRD(1)).LE.SLAYH1) THEN
-               DO IK = 1, IL-1
-                  SLAYH2 = SLAYH2 + LAYH(IK)
-               END DO
-               LAYPRT = DEP2(KCGRD(1)) - SLAYH2
-               DO IK = 1, IL-1
-                  KVEGH = KVEGH + KMESPC * LAYH(IK)
+            SLAYH1 = 0.
+            SLAYH2 = 0.
+            LAYPRT = 0.
+            VGLOOP : DO IL = 1, ILMAX
+               SLAYH1 = SLAYH1 + LAYH(IL)
+               IF (DEP2(KCGRD(1)).LE.SLAYH1) THEN
+                  DO IK = 1, IL-1
+                     SLAYH2 = SLAYH2 + LAYH(IK)
+                  END DO
+                  LAYPRT = DEP2(KCGRD(1)) - SLAYH2
+                  DO IK = 1, IL-1
+                    KVEGH = KVEGH + KMESPC * LAYH(IK)
+                    SINHK = SINH(KVEGH)
+                    A     = C
+                    B     = D
+                    C     = SINHK**3
+                    D     = 3.*SINHK
+                    A     = C - A
+                    B     = D - B
+                    SVEG2 = SVEG2+VEGDRL(IK)*VEGDIL(IK)*VEGNSL(IK)*(A+B)
+                  END DO
+                  KVEGH = KVEGH + KMESPC * LAYPRT
                   SINHK = SINH(KVEGH)
                   A     = C
                   B     = D
@@ -818,44 +909,106 @@
                   D     = 3.*SINHK
                   A     = C - A
                   B     = D - B
-                  SVEG2 = SVEG2+VEGDRL(IK)*VEGDIL(IK)*VEGNSL(IK)*(A + B)
+                  SVEG2 = SVEG2 + VEGDRL(IL)*VEGDIL(IL)*VEGNSL(IL)*(A+B)
+                  EXIT VGLOOP
+               END IF
+            END DO VGLOOP
+
+         END IF
+
+!        --- compute total dissipation
+
+         SVEGET(1:MSC) = SVEG1 * SVEG2
+
+      ELSE IF ( IVEG.EQ.2 ) THEN
+
+!     --- Jacobsen et al. (2019)
+!     --- compute layer- and frequency-independent canopy dissipation factor
+
+         SVEG1 = SQRT(2./PI)*(1/GRAV) * ALFU**3 *
+     &                                 VEGDRL(1) * VEGDIL(1) * VEGNSL(1)
+         IF ( VARNPL ) SVEG1 = SVEG1 * NPLA2(KCGRD(1))
+
+!        --- determine integration interval (submerged vegetation is assumed)
+
+         DZ = MIN( SLAYH, DEP2(KCGRD(1)) ) / REAL(NIP)
+
+!        --- integration from bottom to surface using Simpson's rule
+
+         DO IK = 0, NIP
+
+            ZH  = DZ * REAL(IK)
+            ZDH = ZH - DEP2(KCGRD(1))
+
+            MU = 0.
+
+            DO IS = 1, ISSTOP
+
+               KD = KWAVE(IS,1) * DEP2(KCGRD(1))
+               KC = KWAVE(IS,1) * ZH
+               KZ = KWAVE(IS,1) * ZDH
+
+               IF ( KD.LT.20. ) THEN
+                  ! for all wave-water regimes
+                  FDD(IS) = ( SPCSIG(IS) * COSH(KC)/SINH(KD) )**2      
+                 ! coshk/sinhk should be smaller than 1 but for large kd numbers, almost 1
+               ELSE
+                  !option A: deep water orbital velocity
+                  EKZ     = EXP(KZ)                         ! argument of exp is negative
+                  FDD(IS) = ( SPCSIG(IS) * EKZ )**2
+                  !option B: very small energy in very high frequencies, neglect its dissipation
+                  !FDD(IS) = 0.
+               END IF
+
+!              --- compute first order moment
+
+               DO ID = 1, MDC
+                 MU = MU + FDD(IS) * SPCSIG(IS)**2 * AC2(ID,IS,KCGRD(1))
+                 ! based on velocity spectrum Su
                END DO
-               KVEGH = KVEGH + KMESPC * LAYPRT
-               SINHK = SINH(KVEGH)
-               A     = C
-               B     = D
-               C     = SINHK**3
-               D     = 3.*SINHK
-               A     = C - A
-               B     = D - B
-               SVEG2 = SVEG2 + VEGDRL(IL)*VEGDIL(IL)*VEGNSL(IL)*(A + B)
-               EXIT VGLOOP
+
+            END DO
+
+!           --- integrate Su in frequencies and directions
+
+            MU = MU * DDIR * FRINTF
+
+!           --- determine weight coefficient for integration based on Simpson's rule
+
+            IF ( IK.EQ.0 .OR. IK.EQ.NIP ) THEN
+               C = 1. / 3.
+            ELSE IF ( MOD(IK,2).EQ.0 ) THEN
+               C = 2. / 3.
+            ELSE
+               C = 4. / 3.
             END IF
-         END DO VGLOOP
+!           --- compute frequency-distributed dissipation per integration point
+            DO IS = 1, ISSTOP
+               DCIP(IK,IS) = C * FDD(IS) * SQRT(MU)
+            END DO
+         END DO
 
-      END IF
+         ! --- compute dissipation per frequency
 
-!     --- compute total dissipation
+         SVEGET = SVEG1 * SUM(DCIP,1) * DZ
 
-      SVEGET = SVEG1 * SVEG2
+      END IF
 !
 !     *** test output ***
 !
       IF (TESTFL .AND. ITEST.GE.60) THEN
-         WRITE (PRTEST, 110) IVEG, KCGRD(1), DEP2(KCGRD(1)), SVEGET
+         WRITE (PRTEST, 110) IVEG, KCGRD(1), DEP2(KCGRD(1)), SVEGET(1)
  110     FORMAT (' SVEG :IVEG INDX DEP VEGFAC:', 2I5, 2E12.4)
       END IF
 
       DO IS = 1, ISSTOP
          DO IDDUM = IDCMIN(IS), IDCMAX(IS)
             ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1
-
 !           *** store the results in the array IMATDA ***
 !           *** if testfl store results in array for isoline plot ***
-
-            IMATDA(ID,IS) = IMATDA(ID,IS) + SVEGET
-            IF (TESTFL) PLVEGT(ID,IS,IPTST) = -1.* SVEGET
-            DISSC1(ID,IS,5) = DISSC1(ID,IS,5) + SVEGET
+            IMATDA(ID,IS) = IMATDA(ID,IS) + SVEGET(IS)
+            IF (TESTFL) PLVEGT(ID,IS,IPTST) = -1.* SVEGET(IS)
+            DISSC1(ID,IS,5) = DISSC1(ID,IS,5) + SVEGET(IS)
          END DO
       END DO
 
@@ -865,11 +1018,11 @@
 !****************************************************************
 !
       SUBROUTINE SVEG_COAWST ( DEP2   ,IMATDA   ,ETOT ,SMEBRK,
-     &                  KMESPC ,PLVEGT                       ,
+     &                  KWAVE  ,KMESPC   ,PLVEGT             ,
      &                  IDCMIN ,IDCMAX   ,ISSTOP ,DISSC1     ,
-     &                  ABRBOT ,KWAVE                        ,
+     &                  ABRBOT                               ,
      &                  NPLA2  ,NPHT2    ,NPDA2  ,NPTK2      ,
-     &                  CDVEG_COAWST )
+     &                  CDVEG_COAWST, KC_COAWST, Q_COAWST )
 !
 !****************************************************************
 !
@@ -923,10 +1076,17 @@
 !     40.55, May. 05: implementation of vegetation dissipation formula
 !     40.61, Sep. 06: introduce DISVEG variable for output purposes
 !     40.58, Nov. 08: some modifications and corrections
+!     41.77, Feb. 20: adding Jacobsen vegetation formulation
+!
+!     This routine has been modified for COAWST model 
+!     by Tarandeep Kalra, Nic Jeffress, Jim Hench, Neil Ganju and
+!     John Warner (2021). Cite Nic's future paper. 
 !
 !  2. Purpose
 !
 !     Computation of the source term due to vegetation dissipation
+!     that can account for spatially varying veg. with one vertical
+!     layer and spectral Cd (drag coefficient). 
 !
 !  3. Method
 !
@@ -936,6 +1096,7 @@
 !     waves. Vegetation characteristics that are used as input are
 !     drag coefficient, vegetation height, plant density and diameter.
 !
+!     IVEG=3, COAWST modification to account for spatially varying veg.
 !     The formula used in SWAN is due to Dalrymple (1984)
 !     (see Mendez and Losada, 2004):
 !
@@ -983,6 +1144,68 @@
 !
 !     and is linearized by means of the Picard iteration
 !
+!     IVEG=4, COAWST modification to account for spatially varying veg.
+!     Jacobsen et al. (2019) published a frequency-dependent
+!     dissipation model for waves propagating over a canopy
+!
+!     The dissipation per frequency is
+!
+!                              ----------
+!                              | 2 mu,0
+!     dv_spectral =   2 F Su   | -------
+!                              |   pi
+!                             \|
+!
+!     where
+!
+!     F = 1/2 * rho * alpha^3 * Cd * bv * Nv
+!
+!     The velocity spectrum Su is associated with the
+!     water surface spectral distribution Sn, as given by
+!
+!                                 2
+!               w * cosh k(h+z)
+!     Su =  (  ----------------- )   * Sn
+!                   sinh kh
+!
+!     where
+!
+!     dv_spectral =   dissipation per frequency and layer [kg^2/ms]
+!     F           =   canopy factor                        [kg/m^4]
+!     Su          =   velocity spectrum                     [m^2/s]
+!     mu,0        =   first moment velocity spectrum      [m^2/s^2]
+!     rho         =   water density                        [kg/m^3]
+!     Nv          =   vegetation density                [stems/m^2]
+!     bv          =   stem thickness                            [m]
+!     Cd          =   drag coefficient vegetation               [-]
+!     alpha = 1   =   velocity reduction factor by canopies     [-]
+!     w           =   radian frequency (per component)        [1/s]
+!     k           =   wave number (per component)             [1/m]
+!     z           =   elevation of still water level            [m]
+!                     (z = 0 at water surface)
+!     h           =   water depth                               [m]
+!     Sn          =   variance density spectrum             [m^2*s]
+!
+!     The depth-integrated and frequency-dependent dissipation is then
+!
+!
+!                   -h + ah
+!                   |\
+!                   |
+!     Sveg  =   -   |      dv_spectral dz / rho / g
+!                   |
+!                  \|
+!                   -h
+!
+!     where
+!
+!     ah =  canopy height (under water)                         [m]
+!
+!     The vertical integration is approximated using the Simpson's rule.
+!     A minimum number of integration points would be needed to reduce
+!     the error of the approximation; 21 points appeared to be sufficient
+!
+!     Note: current effects are not included in Jacobsen et al. (2019)
 !
 !  4. Argument variables
 !
@@ -1012,43 +1235,62 @@
      &        KWAVE(MSC,MICMAX)    ,
      &        NPLA2 (MCGRD), NPHT2 (MCGRD), 
      &        NPDA2 (MCGRD), NPTK2 (MCGRD), 
-     &        CDVEG_COAWST(MCGRD)
+     &        CDVEG_COAWST(MCGRD),
+     &        KC_COAWST(MCGRD), 
+     &        Q_COAWST(MCGRD)
       REAL    ETOT, SMEBRK, KMESPC, ABRBOT
 !
+!  5. Parameter variables
+!
+!     ALFU        (orbital) velocity reduction factor by canopies
+!     NIP         total number of subintervals used by Simpson's rule
+!
+      INTEGER, PARAMETER :: NIP  = 20
+      REAL   , PARAMETER :: ALFU = 1.
+!
 !  6. Local variables
 !
 !     A     :     auxiliary variable
 !     B     :     auxiliary variable
 !     C     :     auxiliary variable
 !     D     :     auxiliary variable
+!     DCIP  :     frequency-dependent dissipation for each integration point
+!     DZ    :     interval for vertical integration
+!     EKZ   :     exponential of k(h+z)
+!     FDD   :     factor with orbital velocity to determine Su from Sn
 !     ID    :     counter of the spectral direction
 !     IDDUM :     counter
 !     IENT  :     number of entries
 !     IK    :     counter
 !     IL    :     counter
 !     IS    :     counter of relative frequency band
+!     KC    :     wavenumber times canopy layer height
 !     KD    :     wavenumber times water depth
 !     KVEGH :     wavenumber times plant height
+!     KZ    :     wavenumber times z
 !     LAYPRT:     part of layer below water level
+!     MU    :     first order moment of velocity spectrum
 !     SINHK :     sinh(kh)
 !     SLAYH :     total sum of layer thicknesses
 !     SLAYH1:     sum of layer thicknesses below water level
-!     SLAYH2:     sum of layer thicknesses below water level
-!     SVEG1 :     layer-independent dissipation factor
 !     SVEG2 :     total sum of dissipation factor over layers
 !     SVEGET:     source term containing dissipation due to vegetation
 !                 to be stored in the array IMATDA
+!     ZDH   :     vertical point z between -depth to 0 or - value water clearance in water column
+!     ZH    :     cumulative layer thickness for velocities, bottom up
+!                 (z+d between 0 and vegetation height or depth value in water column)
 !
       INTEGER ID, IDDUM, IENT, IK, IL, IS
       REAL    A, B, C, D, KD, KVEGH, LAYPRT, SINHK, SLAYH,
-     &        SLAYH1, SLAYH2, SVEG1, SVEG2, SVEGET
+     &        SLAYH1, SLAYH2, SVEG1, SVEG2
       REAL    VEGDIL_LOC, VEGTHK_LOC, VEGHGT_LOC 
       REAL    LBYWE 
+      REAL    DZ, EKZ, KC, KZ, MU, ZDH, ZH
       REAL    CAUCHY_NUM, KC_NUM, REY_NUM
       REAL    CFF1, CFF2, CFF3, LAYH_RAT 
       REAL    ALP, ORBVEL, SMA, Q, VEGDRL_CAL  
       REAL    SVEGDR
-!      REAL    CDVEG_COAWST(MCGRD)
+      REAL    DCIP(0:NIP,MSC), FDD(MSC), SVEGET(MSC)
 !
 ! Viscosity of water
 !       
@@ -1061,29 +1303,8 @@
 !
 ! 12. Structure
 !
-!     Vegetation parameters are given per layer thickness, so for
-!     each layer the contribution to wave damping is calculated
-!
-!     This routine checks in which layer the water level is present
-!
-!            -----------
-!                        ILMAX
-!            -----------
-!              _ d       2
-!            --|--------
-!              |         1
-!            --|--------
-!
-!     d     = waterdepth
-!     ILMAX = number of layers in grid point
-!
-!     Subsequently, the vegetation parameters up to the layer where the
-!     water level is in, are used to calculate dissipation for each layer
-!
-!     Thereafter, the contributions to disspation are summed up
-!
-!     With this summation the total dissipation due to vertical varying
-!     vegetation is calculated
+!     Vegetation properties are assumed to be in 1 vertical layer
+!     for COAWST-Veg routine. 
 !
 ! 13. Source text
 !
@@ -1091,87 +1312,100 @@
       DATA IENT/0/
       IF (LTRACE) CALL STRACE (IENT,'SVEG')
 
-!     --- compute layer-independent vegetation dissipation factor
+! --- compute layer-independent vegetation dissipation factor
 
 !   The current logic only works if there is 1 vertical layer. 
-! Tsk KD is KDMEAN 
-      KD     = KMESPC * DEP2(KCGRD(1))
+!   Use local arrays to save spatially varying vegetation properties.
+!   First calculate drag coeff. based on KC and Reynolds number 
+!
       VEGDIL_LOC = NPDA2(KCGRD(1))
       VEGTHK_LOC = NPTK2(KCGRD(1))
       VEGHGT_LOC = NPHT2(KCGRD(1))
- 
-      IF ( KD.GT.10. ) RETURN
-      C     = 3.*KMESPC*(COSH(KD))**3
-      SVEG1 = SQRT(2./PI) * GRAV**2 * (KMESPC/SMEBRK)**3 * SQRT(ETOT)/ C
-
-      IF ( VARNPL ) SVEG1 = SVEG1 * NPLA2(KCGRD(1))
-!
-!     --- compute dissipation factor for each layer and summed up
-!
-      SLAYH = 0.
-      DO IL = 1, ILMAX
-!
-!     --- Modify the veg height that is spatially varying. 
-!
-        SLAYH   = SLAYH + LAYH(IL)
-      ENDDO
-!
-      SLAYH  =VEGHGT_LOC
+      SLAYH      = VEGHGT_LOC
 
-      KVEGH = 0.
-      C     = 0.
-      D     = 0.
-      SVEG2 = 0.
-!
-! If total veg height < depth of water, all the vertical layers 
-! of veg are fully submerged. 
-! 
       IF ( DEP2(KCGRD(1)).GT.SLAYH ) THEN
+#     ifdef VEG_FLEX_SWAN 
 !
-!           iffective blade length from Luhar and Nepf 2016.
+!       Effective blade length from Luhar and Nepf 2016.
 !          
-            ALP        = SLAYH/DEP2(KCGRD(1))
-            ORBVEL     = ( SQRT(2*ETOT)*SMEBRK*COSH(KD*ALP)/SINH(KD) )
+        ALP        = SLAYH/DEP2(KCGRD(1))
+        ORBVEL     = ( SQRT(2*ETOT)*SMEBRK*COSH(KD*ALP)/SINH(KD) )
 !
 ! Second moment of area. 
 !           
-            SMA        = VEGDIL_LOC*VEGTHK_LOC**(3./12.) 
-            CFF1       = RHO*VEGDIL_LOC*(ORBVEL**2.)*(SLAYH**3.)
+        SMA        = VEGDIL_LOC*VEGTHK_LOC**(3./12.) 
+        CFF1       = RHO*VEGDIL_LOC*(ORBVEL**2.)*(SLAYH**3.)
 !            
 !  Avoid division by zero later. 
 !    
-            CFF2       = MAX( (VEGYML(1)*SMA) , MIN_EPS)
+        CFF2       = MAX( (VEGYML(1)*SMA) , MIN_EPS)
 !
 ! Cauchy number. 
 !
-            CAUCHY_NUM = MIN(CFF1/CFF2, MAX_EPS)
+        CAUCHY_NUM = MIN(CFF1/CFF2, MAX_EPS)
 ! 
-            LAYH_RAT   = SLAYH/MIN(ABRBOT, MAX_EPS)
-            CFF3       = MIN( ( 0.7*(CAUCHY_NUM*LAYH_RAT)**(-0.21) ),  
+        LAYH_RAT   = SLAYH/MIN(ABRBOT, MAX_EPS)
+        CFF3       = MIN( ( 0.7*(CAUCHY_NUM*LAYH_RAT)**(-0.21) ),  
      &                               MAX_EPS  )
-            CFF3       = MIN(CFF3,1.0) ! to avoid NaN when Cauchy is zero 
+        CFF3       = MIN(CFF3,1.0) ! to avoid NaN when Cauchy is zero 
+#       else
+        CFF3       = 1.0
+#       endif 
 !
 ! Keulegan Carpenter number. 
 !
-            KC_NUM   = ( SQRT(2.*ETOT)*SMEBRK*COSH(KD*ALP)/ 
+        KC_NUM   = ( SQRT(2.*ETOT)*SMEBRK*COSH(KD*ALP)/ 
      &                   SINH(KD) )*(2.*pi/SMEBRK)/VEGDIL_LOC
-            KC_NUM   = MIN(KC_NUM, MAX_EPS)
+        KC_NUM   = MIN(KC_NUM, MAX_EPS)
 !
 ! Compute Reynolds number.
 !
-            REY_NUM  = ( SQRT(2.*ETOT)*SMEBRK*COSH(KD*ALP)/
+        REY_NUM  = ( SQRT(2.*ETOT)*SMEBRK*COSH(KD*ALP)/
      &                   SINH(KD) )*VEGDIL_LOC/NU  
-            REY_NUM   = MIN(REY_NUM, MAX_EPS)
+        REY_NUM   = MIN(REY_NUM, MAX_EPS)
 
-            IF    ( VEGDRL(1).eq.1. ) THEN 
-                    VEGDRL_CAL   = 1.
-            ELSEIF( VEGDRL(1).eq.2. ) THEN        
+        IF    ( VEGDRL(1).eq.1. ) THEN 
+           VEGDRL_CAL   = 1.
+        ELSEIF( VEGDRL(1).eq.2. ) THEN        
+           Q            = MIN(KC_NUM/(ALP**0.76), MAX_EPS)
+           VEGDRL_CAL   = MIN(EXP(-0.0138*Q)/(Q**0.3),MAX_EPS)
+        ENDIF 
+      ELSE IF ( DEP2(KCGRD(1)).LT.SLAYH ) THEN
+        IF    ( VEGDRL(1).eq.1. ) THEN 
+           VEGDRL_CAL   = 1.
+        ELSEIF( VEGDRL(1).eq.2. ) THEN        
+           VEGDRL_CAL   = 1.95
+        ENDIF 
+      END IF 
+!
+!     --- store values for output, compute drag coefficient in vegetation
+!
+        SVEGDR = VEGDRL_CAL
+        CDVEG_COAWST(KCGRD(1))=MIN(SVEGDR, MAX_EPS)
+        Q_COAWST(KCGRD(1))=Q
+        KC_COAWST(KCGRD(1))=KC_NUM
 !            
 ! Suzuki et al 2012 (actually Mendez and Losada 2004)- almost identical to B&H2009 below
 !  
-                    Q            = MIN(KC_NUM/(ALP**0.76), MAX_EPS)
-                    VEGDRL_CAL   = MIN(EXP(-0.0138*Q)/(Q**0.3),MAX_EPS)
-            ENDIF 
+      IF ( IVEG.EQ.3 ) THEN
+        KD     = KMESPC * DEP2(KCGRD(1))
+ 
+        IF ( KD.GT.10. ) RETURN
+        C     = 3.*KMESPC*(COSH(KD))**3
+        SVEG1 = SQRT(2./PI) * GRAV**2 * (KMESPC/SMEBRK)**3*SQRT(ETOT)/C
+
+        IF ( VARNPL ) SVEG1 = SVEG1 * NPLA2(KCGRD(1))
+
+        SLAYH  =VEGHGT_LOC
+
+        KVEGH = 0.
+        C     = 0.
+        D     = 0.
+        SVEG2 = 0.
+!
+! If total veg height < depth of water, all of the veg
+! are fully submerged. 
+        IF ( DEP2(KCGRD(1)).GT.SLAYH ) THEN
              
             KVEGH = KVEGH + KMESPC * SLAYH * CFF3 
             SINHK = SINH(KVEGH)
@@ -1185,27 +1419,96 @@
 !
 ! If veg height > depth. 
 !
-      ELSE IF ( DEP2(KCGRD(1)).LT.SLAYH ) THEN
-         VEGDRL_CAL=1.95
+        ELSE IF ( DEP2(KCGRD(1)).LT.SLAYH ) THEN
 !
-         SINHK = SINH(KD)
-         A     = SINHK**3
-         B     = 3.*SINHK
-         SVEG2 = VEGDRL_CAL*VEGDIL_LOC*(A + B)
+            SINHK = SINH(KD)
+            A     = SINHK**3
+            B     = 3.*SINHK
+            SVEG2 = VEGDRL_CAL*VEGDIL_LOC*(A + B)
+        END IF
+!        
+!        --- compute total dissipation
+!
+        SVEGET(1:MSC) = SVEG1 * SVEG2
 
-      END IF
+      ELSE IF ( IVEG.EQ.4 ) THEN
+!
+!     --- Jacobsen et al. (2019)
+!
+         SVEG1 = SQRT(2./PI)*(1/GRAV) * ALFU**3 *
+     &                        VEGDRL_CAL * VEGDIL_LOC 
+         IF ( VARNPL ) SVEG1 = SVEG1 * NPLA2(KCGRD(1))
+
+!        --- determine integration interval (submerged vegetation is assumed)
 
-!     --- compute total dissipation
+         DZ = MIN( SLAYH, DEP2(KCGRD(1)) ) / REAL(NIP)
 
-      SVEGET = SVEG1 * SVEG2
-!     --- compute drag coefficient in vegetation
-      SVEGDR = VEGDRL_CAL
-      CDVEG_COAWST(KCGRD(1))=MIN(SVEGDR, MAX_EPS)
+!        --- integration from bottom to surface using Simpson's rule
+
+         DO IK = 0, NIP
+
+            ZH  = DZ * REAL(IK)
+            ZDH = ZH - DEP2(KCGRD(1))
+
+            MU = 0.
+
+            DO IS = 1, ISSTOP
+
+               KD = KWAVE(IS,1) * DEP2(KCGRD(1))
+               KC = KWAVE(IS,1) * ZH
+               KZ = KWAVE(IS,1) * ZDH
+
+               IF ( KD.LT.20. ) THEN
+                  ! for all wave-water regimes
+                  FDD(IS) = ( SPCSIG(IS) * COSH(KC)/SINH(KD) )**2
+                 ! coshk/sinhk should be smaller than 1 but for large kd numbers, almost 1
+               ELSE
+                  !option A: deep water orbital velocity
+                  EKZ     = EXP(KZ)                         ! argument of exp is negative
+                  FDD(IS) = ( SPCSIG(IS) * EKZ )**2
+                  !option B: very small energy in very high frequencies, neglect its dissipation
+                  !FDD(IS) = 0.
+               END IF
+
+!              --- compute first order moment
+
+               DO ID = 1, MDC
+                 MU = MU + FDD(IS) * SPCSIG(IS)**2 * AC2(ID,IS,KCGRD(1))
+                 ! based on velocity spectrum Su
+               END DO
+
+            END DO
+
+!           --- integrate Su in frequencies and directions
+
+            MU = MU * DDIR * FRINTF
+
+!           --- determine weight coefficient for integration based on Simpson's rule
+
+            IF ( IK.EQ.0 .OR. IK.EQ.NIP ) THEN
+               C = 1. / 3.
+            ELSE IF ( MOD(IK,2).EQ.0 ) THEN
+               C = 2. / 3.
+            ELSE
+               C = 4. / 3.
+            END IF
+
+!           --- compute frequency-distributed dissipation per integration point
+
+            DO IS = 1, ISSTOP
+               DCIP(IK,IS) = C * FDD(IS) * SQRT(MU)
+            END DO
+         END DO
+!
+! --- compute dissipation per frequency
+!
+         SVEGET = SVEG1 * SUM(DCIP,1) * DZ
+      END IF
 !
 !     *** test output ***
 !
       IF (TESTFL .AND. ITEST.GE.60) THEN
-         WRITE (PRTEST, 110) IVEG, KCGRD(1), DEP2(KCGRD(1)), SVEGET
+         WRITE (PRTEST, 110) IVEG, KCGRD(1), DEP2(KCGRD(1)), SVEGET(1)
  110     FORMAT (' SVEG :IVEG INDX DEP VEGFAC:', 2I5, 2E12.4)
       END IF
 
@@ -1216,10 +1519,10 @@
 !           *** store the results in the array IMATDA ***
 !           *** if testfl store results in array for isoline plot ***
 
-            IMATDA(ID,IS) = IMATDA(ID,IS) + SVEGET
-            IF (TESTFL) PLVEGT(ID,IS,IPTST) = -1.* SVEGET
-            DISSC1(ID,IS,5) = DISSC1(ID,IS,5) + SVEGET
-!           DISSC1(ID,IS,9) = DISSC1(ID,IS,9) + SVEGDR 
+            IMATDA(ID,IS) = IMATDA(ID,IS) + SVEGET(IS) 
+            IF (TESTFL) PLVEGT(ID,IS,IPTST) = -1.* SVEGET(IS)
+            DISSC1(ID,IS,5) = DISSC1(ID,IS,5) + SVEGET(IS)
+!            DISSC1(ID,IS,9) = DISSC1(ID,IS,9) + SVEGDR 
          END DO
       END DO
 
diff --git a/SWAN/Src/swanmain.F b/SWAN/Src/swanmain.F
index 4cab1a50a1992918c55bdd5126817846e6e3445c..8f46cea3e16efc7624fcc9841582b8f219faad48 100755
--- a/SWAN/Src/swanmain.F
+++ b/SWAN/Src/swanmain.F
@@ -2880,8 +2880,10 @@
       JNPDA3= 41                                                          40.55
       JNPTK2= 42                                                          40.55
       JNPTK3= 43                                                          40.55
-      JDSXC = 44                                                          40.55
-      MCMVAR = 44
+      JDSXC=  44                                                          40.55
+      JDSKC=  45                                                          40.55
+      JDQQQ=  46                                                          40.55
+      MCMVAR = 46
       LADDS = .TRUE.
 #else
       MCMVAR = 28                                                         41.38 40.65 40.61 40.51 40.13
@@ -2940,6 +2942,8 @@
       JDSXW  = 1                                                          40.65 40.61
       JDSXV  = 1                                                          40.65 40.61
       JDSXC  = 1                                                          40.65 40.61
+      JDSKC  = 1                                                          40.65 40.61
+      JDQQQ  = 1                                                          40.65 40.61
 #endif
       JDSXM  = 1                                                          40.65 40.61
       JDSXT  = 1                                                          40.35
@@ -3709,6 +3713,30 @@
       OVLEXP(IVTYPE) = 0.
       OVHEXP(IVTYPE) = 1.
       OVEXCV(IVTYPE) = -9.
+!
+      IVTYPE = 81
+      OVKEYW(IVTYPE) = 'DIKC'                                             40.61
+      OVSNAM(IVTYPE) = 'dikc'                                             40.85
+      OVLNAM(IVTYPE) = 'Kuelgan Carpenter number (veg)'
+      OVUNIT(IVTYPE) = ''
+      OVSVTY(IVTYPE) = 1
+      OVLLIM(IVTYPE) = 0.
+      OVULIM(IVTYPE) = 1.
+      OVLEXP(IVTYPE) = 0.
+      OVHEXP(IVTYPE) = 1.
+      OVEXCV(IVTYPE) = -9.
+!
+      IVTYPE = 82
+      OVKEYW(IVTYPE) = 'DICQ '                                             40.61
+      OVSNAM(IVTYPE) = 'dicq'                                              40.85
+      OVLNAM(IVTYPE) = 'Q number'
+      OVUNIT(IVTYPE) = ''
+      OVSVTY(IVTYPE) = 1
+      OVLLIM(IVTYPE) = 0.
+      OVULIM(IVTYPE) = 1.
+      OVLEXP(IVTYPE) = 0.
+      OVHEXP(IVTYPE) = 1.
+      OVEXCV(IVTYPE) = -9.
 !
       IVTYPE = 58                                                         40.64
       OVKEYW(IVTYPE) = 'QP'
@@ -3974,24 +4002,6 @@
       OVLEXP(IVTYPE) = 0.
       OVHEXP(IVTYPE) = 360.
       OVEXCV(IVTYPE) = -999.
-
-!      As a note, 80 is listed earlier with Veg.
-!      IVTYPE = 80
-!      OVKEYW(IVTYPE) = 'DICD'                                             40.61
-!      ......
-#endif
-#if defined SWAN_MODEL && defined SPECTRUM_STOKES
-      IVTYPE = 81
-      OVKEYW(IVTYPE) = 'SSPEC'
-      OVSNAM(IVTYPE) = 'SSpec'
-      OVLNAM(IVTYPE) = 'Spectrum of Stokes drift'
-      OVUNIT(IVTYPE) = 'm2/s'
-      OVSVTY(IVTYPE) = 3
-      OVLLIM(IVTYPE) = -100.
-      OVULIM(IVTYPE) = 100.
-      OVLEXP(IVTYPE) = -2.
-      OVHEXP(IVTYPE) = 2.
-      OVEXCV(IVTYPE) = 0.
 #endif
 !
       IVTYPE = 100                                                        41.62
@@ -6422,6 +6432,8 @@
 !     initialize CDVEG_COAWST
       DO INDX = 1, MCGRD                                                  32.02
         COMPDA(INDX,JDSXC) =  0.                                          32.02
+        COMPDA(INDX,JDSKC) =  0.                                          32.02
+        COMPDA(INDX,JDQQQ) =  0.                                          32.02
       END DO 
 !
 !     --- initialize GAMBR and RESPL                                      41.38
diff --git a/SWAN/Src/swanout1.F b/SWAN/Src/swanout1.F
index c8c4b1c6a06191cb9f330bf825f523b6267ed0ae..37bf6f5df27692900ddf5f7d139843e54f4e03ed 100755
--- a/SWAN/Src/swanout1.F
+++ b/SWAN/Src/swanout1.F
@@ -261,9 +261,6 @@
       INTEGER, ALLOCATABLE :: IONOD(:)                                    40.51
       REAL, ALLOCATABLE :: ACLOC(:), AUX1(:), VOQ(:)                      40.31
       REAL, ALLOCATABLE :: FORCE(:,:)                                     40.80
-#ifdef SPECTRUM_STOKES
-      REAL, ALLOCATABLE :: KSS(:), USS(:), VSS(:)                         
-#endif
       REAL              :: KNUM(MSC), CG(MSC), NE(MSC), NED(MSC)          40.31
       TYPE(OPSDAT), POINTER :: CUOPS                                      40.31
       TYPE(ORQDAT), POINTER :: CORQ                                       40.31
@@ -404,11 +401,6 @@
 !         assign memory to array ACLOC (contains spectrum for one output point)
 !
           ALLOCATE(ACLOC(MDC*MSC))                                        40.31
-#ifdef SPECTRUM_STOKES
-          ALLOCATE(KSS(MIP*MSC))                                          
-          ALLOCATE(USS(MIP*MSC))
-          ALLOCATE(VSS(MIP*MSC))
-#endif
 !
 !         call SWOEXA to compute quantities for which spectrum is needed
 !         (except wave-induced force)
@@ -421,14 +413,7 @@
      &                 KNUM                ,CG                  ,         40.31 30.90
      &                 SPCDIR              ,NE                  ,         40.31 30.90
      &                 NED                 ,KGRPNT              ,         40.31 30.90
-     &                 COMPDA(1,JDP2)      ,CROSS                         40.86
-#ifdef SPECTRUM_STOKES
-     &                                                          ,
-     &                 KSS(1)                                   ,
-     &                 USS(1)                                   ,
-     &                 VSS(1)                                   
-#endif
-     &                                                          )
+     &                 COMPDA(1,JDP2)      ,CROSS               )         40.86
 !
 !         call SWOEXF to compute wave-driven force on regular grid
 !
@@ -445,9 +430,6 @@
      &                                                             )
 !
           DEALLOCATE(ACLOC)                                               40.31
-#ifdef SPECTRUM_STOKES
-          DEALLOCATE(KSS,USS,VSS)
-#endif
         ENDIF
 !
         IF (ITEST.GE.100 ) THEN
@@ -935,7 +917,7 @@
                                          ! output point and computational 40.86
                                          ! grid point                     40.86
 !
-      INTEGER :: Numcouple, IVTYP, ic
+      INTEGER :: Numcouple, IVTYP
       INTEGER, ALLOCATABLE :: CplType(:)
       INTEGER, ALLOCATABLE :: OQI(:)
       REAL*8, ALLOCATABLE :: OQR(:)
@@ -944,9 +926,6 @@
       INTEGER, ALLOCATABLE :: IONOD(:)                                    40.51
       REAL, ALLOCATABLE :: ACLOC(:), AUX1(:), VOQ(:)                      40.31
       REAL, ALLOCATABLE :: FORCE(:,:)                                     40.80
-#ifdef SPECTRUM_STOKES
-      REAL, ALLOCATABLE :: KSS(:), USS(:), VSS(:)
-#endif
       REAL              :: KNUM(MSC), CG(MSC), NE(MSC), NED(MSC)          40.31
       TYPE(OPSDAT), POINTER :: CUOPS                                      40.31
       TYPE(ORQDAT), POINTER :: CORQ                                       40.31
@@ -960,49 +939,27 @@
 
       Numcouple=14
 #if defined VEGETATION && defined VEG_SWAN_COUPLING 
-      Numcouple=Numcouple+2
-#endif
-#if defined SPECTRUM_STOKES
-      Numcouple=Numcouple+1
+!      Numcouple=15
+      Numcouple=16
 #endif 
       ALLOCATE (Cpltype(Numcouple))
-      ic=1
-      Cpltype(ic)=54                      ! DISBOT  SWOEXD
-      ic=ic+1
-      Cpltype(ic)=55                      ! DISSURF SWOEXD
-      ic=ic+1
-      Cpltype(ic)=56                      ! DISWCAP SWOEXD
-      ic=ic+1
-      Cpltype(ic)=50                      ! TMBOT   SWOEXD
-      ic=ic+1
-      Cpltype(ic)=6                       ! UBOT    SWOEXD
-      ic=ic+1
-      Cpltype(ic)=8                       ! QB      SWOEXD
-      ic=ic+1
-      Cpltype(ic)=10                      ! HSIGN   SWOEXA
-      ic=ic+1
-      Cpltype(ic)=12                      ! RTP     SWOEXA
-      ic=ic+1
-      Cpltype(ic)=13                      ! DIR     SWOEXA
-      ic=ic+1
-      Cpltype(ic)=14                     ! PDI     SWOEXA
-      ic=ic+1
-      Cpltype(ic)=17                     ! WLEN    SWOEXA
-      ic=ic+1
-      Cpltype(ic)=71                     ! WLENP   SWOEXA
-      ic=ic+1
-      Cpltype(ic)=16                     ! WDSPR   SWOEXA
-      ic=ic+1
-      Cpltype(ic)=58                     ! WQP     SWOEXA
+      Cpltype(1)=54                      ! DISBOT  SWOEXD
+      Cpltype(2)=55                      ! DISSURF SWOEXD
+      Cpltype(3)=56                      ! DISWCAP SWOEXD
+      Cpltype(4)=50                      ! TMBOT   SWOEXD
+      Cpltype(5)=6                       ! UBOT    SWOEXD
+      Cpltype(6)=8                       ! QB      SWOEXD
+      Cpltype(7)=10                      ! HSIGN   SWOEXA
+      Cpltype(8)=12                      ! RTP     SWOEXA
+      Cpltype(9)=13                      ! DIR     SWOEXA
+      Cpltype(10)=14                     ! PDI     SWOEXA
+      Cpltype(11)=17                     ! WLEN    SWOEXA
+      Cpltype(12)=71                     ! WLENP   SWOEXA
+      Cpltype(13)=16                     ! WDSPR   SWOEXA
+      Cpltype(14)=58                     ! WQP     SWOEXA
 #if defined VEGETATION && defined VEG_SWAN_COUPLING 
-      ic=ic+1
-      Cpltype(ic)=57                     ! DISVEG  SWOEXD         
-      ic=ic+1
-      Cpltype(ic)=80                     ! DICDVG  SWOEXD         
-#endif 
-#if defined SPECTRUM_STOKES
-      ic=ic+1
-      Cpltype(ic)=81                     ! SSPEC SWOEXA
+      Cpltype(15)=57                     ! DISVEG  SWOEXD         
+      Cpltype(16)=80                     ! DICDVG  SWOEXD         
 #endif 
 !
       ALLOCATE (OQI(4))
@@ -1053,12 +1010,7 @@
 !       assign memory to array VOQ (contains output quantities for all
 !                                   output points)
         ALLOCATE(VOQ(MIP*NVOQP))                                          40.31
-#ifdef SPECTRUM_STOKES
-!       assign memory to array KSS, USS, VSS
-        ALLOCATE(KSS(MIP*MSC))
-        ALLOCATE(USS(MIP*MSC))
-        ALLOCATE(VSS(MIP*MSC))
-#endif
+!       VOQ = 0.
 !
 !       assign memory to array CROSS (indicates crossing of obstacles     40.86
 !                                     in between output and grid points)  40.86
@@ -1108,14 +1060,7 @@
      &                 KNUM                ,CG                  ,         40.31 30.90
      &                 SPCDIR              ,NE                  ,         40.31 30.90
      &                 NED                 ,KGRPNT              ,         40.31 30.90
-     &                 COMPDA(1,JDP2)      ,CROSS                         40.86
-#ifdef SPECTRUM_STOKES
-     &                                                          ,
-     &                 KSS(1)                                   ,
-     &                 USS(1)                                   ,
-     &                 VSS(1)                                   
-#endif
-     &                                                          )
+     &                 COMPDA(1,JDP2)      ,CROSS               )         40.86
 !
           DEALLOCATE(ACLOC)                                               40.31
         ENDIF
@@ -1124,15 +1069,9 @@
 
 !       Couple to ocean model
         CALL WAV2OCN_COUPLING (MIP, NVOQP, VOQR, VOQ(1), IRQ,
-#ifdef SPECTRUM_STOKES
-     &                         KSS(1), USS(1), VSS(1),
-#endif
      &                         IVTYPE, COMPDA, Numcouple, ng, io)
         IF (STPNOW()) RETURN
         DEALLOCATE(VOQ,CROSS,IONOD)
-#ifdef SPECTRUM_STOKES
-        DEALLOCATE(KSS,USS,VSS)
-#endif
   70  CONTINUE
 !
 !     Termination of output
@@ -1388,9 +1327,6 @@
       INTEGER, ALLOCATABLE :: IONOD(:)                                    40.51
       REAL, ALLOCATABLE :: ACLOC(:), AUX1(:), VOQ(:)                      40.31
       REAL, ALLOCATABLE :: FORCE(:,:)                                     40.80
-#ifdef SPECTRUM_STOKES
-      REAL, ALLOCATABLE :: KSS(:), USS(:), VSS(:)
-#endif
       REAL              :: KNUM(MSC), CG(MSC), NE(MSC), NED(MSC)          40.31
       TYPE(OPSDAT), POINTER :: CUOPS                                      40.31
       TYPE(ORQDAT), POINTER :: CORQ                                       40.31
@@ -1471,11 +1407,7 @@
 !       assign memory to array VOQ (contains output quantities for all
 !                                   output points)
         ALLOCATE(VOQ(MIP*NVOQP))                                          40.31
-#ifdef SPECTRUM_STOKES
-       ALLOCATE(KSS(MIP*MSC))
-       ALLOCATE(USS(MIP*MSC))
-       ALLOCATE(VSS(MIP*MSC))
-#endif
+!       VOQ = 0.
 !
 !       assign memory to array CROSS (indicates crossing of obstacles     40.86
 !                                     in between output and grid points)  40.86
@@ -1525,19 +1457,10 @@
      &                 KNUM                ,CG                  ,         40.31 30.90
      &                 SPCDIR              ,NE                  ,         40.31 30.90
      &                 NED                 ,KGRPNT              ,         40.31 30.90
-     &                 COMPDA(1,JDP2)      ,CROSS                         40.86
-#ifdef SPECTRUM_STOKES
-     &                 KSS(1)                                   ,
-     &                 USS(1)                                   ,
-     &                 VSS(1)                                   
-#endif
-     &                                                          )
+     &                 COMPDA(1,JDP2)      ,CROSS               )         40.86
 !
           DEALLOCATE(ACLOC)                                               40.31
         ENDIF
-#ifdef SPECTRUM_STOKES
-        DEALLOCATE(KSS,USS,VSS)
-#endif
 !
   68    CONTINUE                                                          40.31
 
@@ -1834,18 +1757,6 @@
            ENDIF                                                          41.62
          ENDIF                                                            41.62
 !
-!        for spectrum based Stokes drift estimate
-!
-#if defined SPECTRUM_STOKES
-         IF (IVTYPE.EQ.81) THEN
-           IF (.NOT.OQPROC(81)) THEN
-             NVOQP = NVOQP + 1
-             VOQR(81) = NVOQP
-             OQPROC(81) = .TRUE.
-           ENDIF
-         ENDIF
-#endif
-!
 !        for some quantities compute action densities
 !
          IF (IVTYPE.EQ.10 .OR. IVTYPE.EQ.11 .OR. IVTYPE.EQ.12 .OR.
@@ -1862,10 +1773,6 @@
      &       IVTYPE.EQ.32 .OR. IVTYPE.EQ.33 .OR. IVTYPE.EQ.42 .OR.
      &       IVTYPE.EQ.47 .OR. IVTYPE.EQ.59 .OR. IVTYPE.EQ.71      )      41.15 40.64 40.41 40.00
      &   BKC = 2
-
-#ifdef SPECTRUM_STOKES
-         IF (IVTYPE.EQ.81) BKC =2 
-#endif
 !
          IF (IVTYPE.EQ.11 .AND. ICUR.GT.0) BKC = 2                        20.36
          IF (BKC.GT.0) THEN
@@ -3036,6 +2943,48 @@
         ENDIF
       ENDIF
 !
+!     Keulegan Carpenter number 
+!
+      IF (OQPROC(81)) THEN                                               40.61
+        IF (ITEST.GE.50 .OR. IOUTES .GE. 10) WRITE (PRTEST, 121) 81,
+     &  VOQR(81), JDSKC
+        IF (JDSKC.GT.1) THEN                                              40.65
+           IF (OPTG.NE.5) THEN                                            40.80
+             CALL SWIPOL(COMPDA(1,JDSKC),OVEXCV(81), XC, YC, MIP, CROSS,  40.86
+     &                    VOQ(1,VOQR(81)) ,KGRPNT, COMPDA(1,JDP2))
+           ELSE                                                           40.80
+              CALL SwanInterpolateOutput ( VOQ(1, VOQR(81)), VOQ(1,1),    40.80
+     &                                     VOQ(1,2), COMPDA(1,JDSKC),     40.80
+     &                                     MIP, KVERT, OVEXCV(81) )       40.80
+          ENDIF                                                           40.80
+        ELSE
+          DO IP = 1, MIP                                                  31.02
+            VOQ(IP,VOQR(81))  = OVEXCV(81)                                31.02
+          ENDDO                                                           31.02
+        ENDIF
+      ENDIF
+!
+!      Q  
+!
+      IF (OQPROC(82)) THEN                                               40.61
+        IF (ITEST.GE.50 .OR. IOUTES .GE. 10) WRITE (PRTEST, 121) 82,
+     &  VOQR(82), JDQQQ
+        IF (JDQQQ.GT.1) THEN                                              40.65
+           IF (OPTG.NE.5) THEN                                            40.80
+             CALL SWIPOL(COMPDA(1,JDQQQ),OVEXCV(82), XC, YC, MIP, CROSS,  40.86
+     &                    VOQ(1,VOQR(82)) ,KGRPNT, COMPDA(1,JDP2))
+           ELSE                                                           40.80
+              CALL SwanInterpolateOutput ( VOQ(1, VOQR(82)), VOQ(1,1),    40.80
+     &                                     VOQ(1,2), COMPDA(1,JDQQQ),     40.80
+     &                                     MIP, KVERT, OVEXCV(82) )       40.80
+          ENDIF                                                           40.80
+        ELSE
+          DO IP = 1, MIP                                                  31.02
+            VOQ(IP,VOQR(82))  = OVEXCV(82)                                31.02
+          ENDDO                                                           31.02
+        ENDIF
+      ENDIF
+!
 !     turbulent dissipation
 !
       IF (OQPROC(72)) THEN                                                40.35
@@ -4223,15 +4172,7 @@
      &                   WK         ,CG         ,
      &                   SPCDIR     ,NE         ,
      &                   NED        ,KGRPNT     ,
-     &                   DEPXY      ,CROSS                                40.86 30.50
-#ifdef SPECTRUM_STOKES
-     &                                          ,
-     &                   KSS                    ,
-     &                   USS                    ,
-     &                   VSS
-#endif
-     &                                          )
-     
+     &                   DEPXY      ,CROSS      )                         40.86 30.50
 !                                                                      *
 !***********************************************************************
 !
@@ -4387,11 +4328,7 @@
 !
 ! 13. Source text
 !
-#ifdef SPECTRUM_STOKES
-      PARAMETER  (NVOTP=94)
-#else
-      PARAMETER  (NVOTP=93)                                               41.15 40.64 40.51 40.41 40.00
-#endif
+      PARAMETER  (NVOTP=93)                                               41.62 41.15 40.64 40.51 40.41 40.00
       REAL       XC(MIP)        ,YC(MIP)       ,AC2(MDC,MSC,MCGRD),
      &           VOQ(MIP,*)     ,
      &           WK(*)          ,
@@ -4404,12 +4341,6 @@
 !
       INTEGER    NP             ,DIMXPT                                   41.62
       REAL       UABS           ,UDIR                                     41.62
-#ifdef SPECTRUM_STOKES
-      REAL       ETOTX, ETOTY, ETOTXM1, ETOTYM1
-      REAL       USS(MIP,MSC)
-      REAL       VSS(MIP,MSC)
-      REAL       KSS(MIP,MSC)
-#endif
       REAL, ALLOCATABLE :: XPT(:,:)                                       41.62
 !
       REAL, ALLOCATABLE :: FLUX(:,:,:), FLOC(:,:)
@@ -4418,17 +4349,6 @@
       LOGICAL :: EXCPT     ! if true value in point is undefined          40.86
       SAVE IENT, IVOTP
       DATA IENT /0/
-#ifdef SPECTRUM_STOKES
-      DATA IVOTP /10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 28, 32, 33,     20.61
-     &            100, 101, 102, 103, 104, 105, 106, 107, 108, 109,       41.62
-     &            110, 111, 112, 113, 114, 115, 116, 117, 118, 119,       41.62
-     &            120, 121, 122, 123, 124, 125, 126, 127, 128, 129,       41.62
-     &            130, 131, 132, 133, 134, 135, 136, 137, 138, 139,       41.62
-     &            140, 141, 142, 143, 144, 145, 146, 147, 148, 149,       41.62
-     &            150, 151, 152, 153, 154, 155, 156, 157, 158, 159,       41.62
-     &            160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 171,  41.62
-     &            42, 43, 44, 47, 48, 53, 58, 59, 71, 81/                     41.15 40.64 40.51 40.41 40.00
-#else
       DATA IVOTP /10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 28, 32, 33,     20.61
      &            100, 101, 102, 103, 104, 105, 106, 107, 108, 109,       41.62
      &            110, 111, 112, 113, 114, 115, 116, 117, 118, 119,       41.62
@@ -4437,8 +4357,7 @@
      &            140, 141, 142, 143, 144, 145, 146, 147, 148, 149,       41.62
      &            150, 151, 152, 153, 154, 155, 156, 157, 158, 159,       41.62
      &            160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 171,  41.62
-     &            42, 43, 44, 47, 48, 53, 58, 59, 71/                     41.15 40.64 40.51 40.41 40.00
-#endif
+     &            42, 43, 44, 47, 48, 53, 58, 59, 71 /                    41.15 40.64 40.51 40.41 40.00
       CALL STRACE (IENT, 'SWOEXA')
 !
 !     in case of transport of energy, compute the energy flux in all grid points
@@ -5798,38 +5717,6 @@
              VOQ(IP,VOQR(IVTYPE)) = OVEXCV(IVTYPE)
            ENDIF
         ENDIF
-#ifdef SPECTRUM_STOKES
-!
-!       Stokes drift spectrum in shallow water
-!
-        IVTYPE = 81
-        IF (OQPROC(IVTYPE)) THEN
-           ETOTXM1 = 0.
-           ETOTYM1 = 0.
-           DO IS = 2, MSC
-              ETOTX = 0.
-              ETOTY = 0.
-              DO ID = 1, MDC
-                 EADX  = ACLOC(ID,IS)*SPCDIR(ID,2)*DDIR
-                 EADY  = ACLOC(ID,IS)*SPCDIR(ID,3)*DDIR
-                 ETOTX = ETOTX + EADX
-                 ETOTY = ETOTY + EADY   
-              ENDDO
-              DS = SPCSIG(IS) - SPCSIG(IS-1)
-              USS(IP,IS-1) = DS*(WK(IS)*ETOTX*SPCSIG(IS)**2.0
-     &                     + WK(IS-1)*ETOTXM1*SPCSIG(IS-1)**2.0)
-              VSS(IP,IS-1) = DS*(WK(IS)*ETOTY*SPCSIG(IS)**2.0
-     &                     + WK(IS-1)*ETOTYM1*SPCSIG(IS-1)**2.0)
-              ETOTXM1 = ETOTX
-              ETOTYM1 = ETOTY
-              KSS(IP,IS-1) = 0.5*(WK(IS)+WK(IS-1))
-           END DO
-
-           USS(IP,MSC) = 2.0*WK(MSC)*ETOTX*SPCSIG(MSC)**3.0
-           VSS(IP,MSC) = 2.0*WK(MSC)*ETOTY*SPCSIG(MSC)**3.0
-           KSS(IP,MSC) = WK(MSC)
-        ENDIF
-#endif
 !
 !       partitioning output according to Hanson and Phillips (2001)
 !
@@ -5921,16 +5808,6 @@
           ENDIF
  730    CONTINUE
 !
-#ifdef SPECTRUM_STK
-        IF (IVTYPE.EQ.80) THEN
-          DO IS = 1, MSC
-             USS(IP,IS) = OVEXCV(IVTYPE)
-             VSS(IP,IS) = OVEXCV(IVTYPE)
-             KSS(IP,IS) = OVEXCV(IVTYPE)
-          END DO
-        ENDIF
-#endif
-! 
  800  CONTINUE
 !
       IF (ALLOCATED(FLUX)) DEALLOCATE(FLUX)
diff --git a/SWAN/Src/swanpre1.F b/SWAN/Src/swanpre1.F
index 7ce1b6fc94bfbcf6ba0b91ae851ddd4b84562338..f952432273424b33b97e02e68291dca0d816304c 100755
--- a/SWAN/Src/swanpre1.F
+++ b/SWAN/Src/swanpre1.F
@@ -356,6 +356,7 @@
       TYPE VEGPT                                                          40.55
         INTEGER              :: N                                         40.55
         REAL                 :: H, D, C, Y                                40.55
+!        REAL                 :: C, Y                                      40.55
         TYPE(VEGPT), POINTER :: NEXTV                                     40.55
       END TYPE VEGPT                                                      40.55
 
@@ -1530,7 +1531,7 @@
           PSNAME = 'NPHTGRID'                                             40.55
 #ifdef COAWST_COUPLING
           CALL SINPGR (IGRD, IGR2, PSNAME, ng)
-          VARNPL = .TRUE.                                                   40.55
+          VARNPH = .TRUE.                                                   40.55
           IF (JNPHT2.LE.1) THEN                                             40.55
             MCMVAR = MCMVAR + 2                                             40.55
             JNPHT2 = MCMVAR - 1                                             40.55
@@ -2774,12 +2775,14 @@
 !
 ! ==========================================
 !
-!   VEGEtation   < [height]  [diamtr]  [nstems]  [drag] >
+!   VEGEtation   < [iveg] [height]  [diamtr]  [nstems]  [drag] >
 !
 ! ==========================================
 !
       IF (KEYWIS ('VEGE')) THEN
-        IVEG  = 1
+!        IVEG  = 1
+!        -        IVEG  = 1
+        CALL ININTG ('IVEG', IVEG, 'STA', 1)                        
         ILMAX = 0
         FRSTV%H = 0.
         FRSTV%D = 0.
@@ -2838,12 +2841,13 @@
 !
 ! ==========================================
 !
-!   VEGCtation   < [height]  [diamtr]  [nstems]  [drag] >
+!   VEGCtation   < [iveg]  [height]  [diamtr]  [nstems]  [drag] >
 !
 ! ==========================================
 !
       IF (KEYWIS ('VEGC')) THEN
-        IVEG  = 2
+!        IVEG  = 3
+        CALL ININTG ('IVEG', IVEG, 'STA', 1)                        
         ILMAX = 0
         FRSTV%H = 0.
         FRSTV%D = 0.
diff --git a/SWAN/Src/swanpre2.F b/SWAN/Src/swanpre2.F
index 44a6730626d02f46b854eb0619c6490d7e4f3e62..ac7d17b85169a0a4d05cb86d0e16c46a75f5dd3e 100755
--- a/SWAN/Src/swanpre2.F
+++ b/SWAN/Src/swanpre2.F
@@ -1686,6 +1686,16 @@
                    JDSXC  = MCMVAR                                        40.65
                    ALOCMP = .TRUE.                                        40.97
                 ENDIF                                                     40.65
+                IF (IVTYPE.EQ.81 .AND. JDSKC.LE.1) THEN                   40.65
+                   MCMVAR = MCMVAR+1                                      40.65
+                   JDSKC  = MCMVAR                                        40.65
+                   ALOCMP = .TRUE.                                        40.97
+                ENDIF                                                     40.65
+                IF (IVTYPE.EQ.82 .AND. JDQQQ.LE.1) THEN                   40.65
+                   MCMVAR = MCMVAR+1                                      40.65
+                   JDQQQ  = MCMVAR                                        40.65
+                   ALOCMP = .TRUE.                                        40.97
+                ENDIF                                                     40.65
                 IF (IVTYPE.EQ.60 .AND. JGENR.LE.1) THEN                   40.85
                    MCMVAR = MCMVAR+1                                      40.85
                    JGENR  = MCMVAR                                        40.85
@@ -2009,6 +2019,16 @@
                JDSXC  = MCMVAR                                            40.65
                ALOCMP = .TRUE.                                            40.97
             ENDIF                                                         40.65
+            IF (IVTYPE.EQ.81 .AND. JDSKC.LE.1) THEN                       40.65
+               MCMVAR = MCMVAR+1                                          40.65
+               JDSKC  = MCMVAR                                            40.65
+               ALOCMP = .TRUE.                                            40.97
+            ENDIF                                                         40.65
+            IF (IVTYPE.EQ.82 .AND. JDQQQ.LE.1) THEN                       40.65
+               MCMVAR = MCMVAR+1                                          40.65
+               JDQQQ  = MCMVAR                                            40.65
+               ALOCMP = .TRUE.                                            40.97
+            ENDIF                                                         40.65
             IF (IVTYPE.EQ.60 .AND. JGENR.LE.1) THEN                       40.85
                MCMVAR = MCMVAR+1                                          40.85
                JGENR  = MCMVAR                                            40.85
@@ -2205,14 +2225,10 @@
         NREF = 0
 !       --- append node number to FILENM in case of                       40.30
 !           parallel computing                                            40.30
-#if defined SWAN_SINGLESPEC
-!jcw prevent additon of -001 etc to spec filename 
-#else
         IF ( PARLL ) THEN                                                 40.30
            ILPOS = INDEX ( FILENM, ' ' )-1                                40.30
            WRITE(FILENM(ILPOS+1:ILPOS+4),33) INODE                        40.30
         END IF                                                            40.30
-#endif
 !PUN        FILENM = TRIM(LOCALDIR)//DIRCH2//TRIM(FILENM)                     40.95
         ORQTMP%PSNAME = PSNAME                                            40.31
         ORQTMP%OQI(1) = NREF                                              40.31
diff --git a/SWAN/Src/swmod1.F b/SWAN/Src/swmod1.F
index d23313737eb41bbea484f37276b700d203385786..e612a57c68d78cb988023b403a8962470d0d7b0f 100755
--- a/SWAN/Src/swmod1.F
+++ b/SWAN/Src/swmod1.F
@@ -612,7 +612,7 @@
 !                     DISSRF and DISWCP added)
 !     40.61, Sep. 06: NMOVAR increased by 1 (quantity DISMUD added)
 !     40.61, Sep. 06: NMOVAR increased by 1 (quantity DISVEG added)
-!TSK  40.61, Sep. 06: NMOVAR increased by 1 (quantity DICDVG added)
+!TSK   40.61, Sep. 06: NMOVAR increased by 1 (quantity DICDVG added)
 !     40.64, Apr. 07: NMOVAR increased by 2 (quantities QP and BFI added)
 !     40.85, Aug. 08: NMOVAR increased by 10 (quantities GENE, GENW, REDI, REDQ, REDT,
 !                                             PROPA, PROPX, PROPT, PROPS and RADS added)
@@ -732,8 +732,6 @@
 !  (78) [  'HICE']
 #ifdef SWAN_MODEL
 !  (79) ['DIRBOT']
-!  (80) ['']
-!  (81) ['SWAN_SPEC'] Spectrum_stokes
 #endif
 ! OVSNAM(NMOVAR)   short name of output quantity
 !  ( 1) [    'Xp']
@@ -812,8 +810,6 @@
 !  (78) ['  hice']
 #ifdef SWAN_MODEL
 !  (79) ['Dirbot']
-!  (80) ['']
-!  (81) ['SWAN_SPEC'] Spectrum_stokes
 #endif
 ! OVLNAM(NMOVAR)   long name of output quantity
 !  ( 1) [                        'X user coordinate']
@@ -892,8 +888,6 @@
 !  (78) [                            'ice thickness']
 #ifdef SWAN_MODEL
 !  (79) [            'Average bottom wave direction']
-!  (80) ['']
-!  (81) ['SWAN_SPEC'] Spectrum_stokes
 #endif
 ! OVUNIT(NMOVAR)   unit of of output quantity
 !  ( 1) [CALCULAT] =UL
@@ -972,8 +966,6 @@
 !  (78) [     'm']
 #ifdef SWAN_MODEL
 !  (79) [CALCULAT] =UDI
-!  (80) ['']
-!  (81) ['SWAN_SPEC'] Spectrum_stokes
 #endif
 ! SNAME [        ] name of output point set
 ! UAP   [  'W/m2'] unit of dissipation
@@ -1109,8 +1101,6 @@
 !  (78)  [  1]   hice
 #ifdef SWAN_MODEL
 !  (79)  [  2]   DirBot
-!  (80) ['']
-!  (81) ['SWAN_SPEC'] Spectrum_stokes
 #endif
 ! OVEXCV(NMOVAR)   exception value for output quantity
 !  ( 1)  [ -1.E10] Xp
@@ -1187,8 +1177,6 @@
 !  (78)  [    -9.] hice
 #ifdef SWAN_MODEL
 !  (79)  [  -999.] DirBot
-!  (80) ['']
-!  (81) ['SWAN_SPEC'] Spectrum_stokes
 #endif
 ! OVLEXP(NMOVAR)   lower expected limit of output quantity
 !  ( 1)  [ -1.E10] Xp
@@ -1265,8 +1253,6 @@
 !  (78)  [     0.] hice
 #ifdef SWAN_MODEL
 !  (79)  [     0.] DirBot
-!  (80)  [     0.] Veg
-!  (81)  [     0.] Spectrum_stokes
 #endif
 ! OVLLIM(NMOVAR)   lower limit of validity of output quantity
 !  ( 1)  [ -1.E10] Xp
@@ -1343,8 +1329,6 @@
 !  (78)  [     0.] hice
 #ifdef SWAN_MODEL
 !  (79)  [     0.] DirBot
-!  (80) ['']
-!  (81) ['SWAN_SPEC'] Spectrum_stokes
 #endif
 ! OVHEXP(NMOVAR)   upper expected limit of output quantity
 !  ( 1)  [  1.E10] Xp
@@ -1421,8 +1405,6 @@
 !  (78)  [    10.] hice
 #ifdef SWAN_MODEL
 !  (79)  [   360.] DirBot
-!  (80) ['']
-!  (81) ['SWAN_SPEC'] Spectrum_stokes
 #endif
 ! OVULIM(NMOVAR)   upper limit of validity
 !  ( 1)  [  1.E10] Xp
@@ -1499,8 +1481,6 @@
 !  (78)  [   100.] hice
 #ifdef SWAN_MODEL
 !  (79)  [   360.] DirBot
-!  (80) ['']
-!  (81) ['SWAN_SPEC'] Spectrum_stokes
 #endif
 ! SINCQ  [       ] sin of ALCQ
 ! SINPQ  [       ] sin of ALPQ
@@ -2008,7 +1988,7 @@
 ! JDISS  [  2] total dissipation within array COMPDA
 ! JDPSAV [  1] saved depth (for setup) within array COMPDA
 ! JDP1   [  7] old depth within array COMPDA
-! JDP2   [  9] new depth within array COMPDA - TSK
+! Tsk JDP2   [  9] new depth within array COMPDA
 ! JDP3   [ 15] last read depth within array COMPDA
 ! JDSXB  [  1] bottom friction dissipation within array COMPDA
 !              set by command TABLE/BLOCK
@@ -2177,7 +2157,8 @@
       INTEGER             JDSXS
       INTEGER             JDSXW
       INTEGER             JDSXM
-      INTEGER             JDSXV,       JDSXC,       JDSXT
+      INTEGER             JDSXV,       JDSXC,       JDSKC,       JDQQQ  
+      INTEGER             JDSXT
       INTEGER             JGSXW,       JGENR
       INTEGER             JRSXQ,       JRSXT,       JREDS,       JRADS
       INTEGER             JTSXG,       JTSXT,       JTSXS,       JTRAN
@@ -2486,6 +2467,11 @@
 ! IVEG   [     0] indicates vegetation dissipation term:
 !                 =0; no dissipation due to vegetation
 !                 =1; dissipation due to vegetation according to Dalrymple (1984)
+!                 =2; dissipation due to vegetation according to Jacobsen et al. (2019)
+!                 =3; dissipation due to veg-COAWST according to
+!                 Dalrymple (1984)
+!                 =4; dissipation due to veg-COAWST(1984) according to
+!                 Jacobsen et al. (2019) 
 ! ITURBV [     0] indicates if turbulent viscosity is activated           40.35
 !                 =0; not activated                                       40.35
 !                 =1; activated                                           40.35
diff --git a/SWAN/Src/swmod2.F b/SWAN/Src/swmod2.F
index 0e97ddfd3cb660b238cac3d27ad9e0ff382fb6b8..8798055eb1a6a4b576c89ce5abd1c3ea1bc2ac49 100755
--- a/SWAN/Src/swmod2.F
+++ b/SWAN/Src/swmod2.F
@@ -1473,4 +1473,4 @@
 !PUN!/impi      MODULE MPI
 !PUN!/impi      INCLUDE 'mpif.h'
 !PUN!/impi      END MODULE MPI
-#endif
\ No newline at end of file
+#endif
diff --git a/SWAN/Src/waves_coupler.F b/SWAN/Src/waves_coupler.F
index a6ea9903315fa4ddf647fdfdf7e379fb0d1c62a0..80319f9652b1227bd602579f61f16d9e072f909a 100755
--- a/SWAN/Src/waves_coupler.F
+++ b/SWAN/Src/waves_coupler.F
@@ -183,17 +183,8 @@
       integer :: cid, cad, MXCGLc
       character (len=70)  :: nc_name
       character (len=20)  :: to_add
-# ifdef SPECTRUM_STOKES
-      character (len=520) :: wostring
-      character (len=140) :: owstring
-# else
       character (len=120) :: wostring
       character (len=120) :: owstring
-# endif
-# ifdef SPECTRUM_STOKES
-      character (len=5)   :: LSTCODE
-      integer             :: IZ
-# endif
 
       real :: cff
 
@@ -535,8 +526,10 @@
       write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
       cid=cid+cad
 !
-#if defined VEGETATION && defined VEG_SWAN_COUPLING \
-      && defined VEG_STREAMING
+!#if defined VEGETATION && defined VEG_SWAN_COUPLING \
+!      && defined VEG_STREAMING
+#if defined VEGETATION && defined VEG_SWAN_COUPLING 
+!      && defined VEG_STREAMING
       to_add=':DISVEG'
       cad=LEN_TRIM(to_add)
       write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
@@ -546,43 +539,6 @@
       cad=LEN_TRIM(to_add)
       write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
       cid=cid+cad
-#endif
-#ifdef SPECTRUM_STOKES
-      DO IZ=1,MSC
-        IF (IZ.LT.10) THEN
-         WRITE(LSTCODE,"(A4,I1)") ":KS0",IZ
-        ELSE
-         WRITE(LSTCODE,"(A3,I2)") ":KS",IZ 
-        END IF
-        to_add=LSTCODE
-        cad=LEN_TRIM(to_add)
-        write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
-        cid=cid+cad  
-      END DO
-
-      DO IZ=1,MSC
-        IF (IZ.LT.10) THEN
-         WRITE(LSTCODE,"(A4,I1)") ":US0",IZ
-        ELSE
-         WRITE(LSTCODE,"(A3,I2)") ":US",IZ 
-        END IF
-        to_add=LSTCODE
-        cad=LEN_TRIM(to_add)
-        write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
-        cid=cid+cad  
-      END DO
-
-      DO IZ=1,MSC
-        IF (IZ.LT.10) THEN
-         WRITE(LSTCODE,"(A4,I1)") ":VS0",IZ
-        ELSE
-         WRITE(LSTCODE,"(A3,I2)") ":VS",IZ 
-        END IF
-        to_add=LSTCODE
-        cad=LEN_TRIM(to_add)
-        write(wostring(cid:cid+cad-1),'(a)') to_add(1:cad)
-        cid=cid+cad  
-      END DO
 #endif  
 !
 !  Finalize and remove trailing spaces from the wostring
@@ -770,9 +726,6 @@
 
 # ifdef ROMS_COUPLING
       SUBROUTINE WAV2OCN_COUPLING (MIP, NVOQP, VOQR, VOQ, IRQ,          &
-# ifdef SPECTRUM_STOKES
-     &                             KSS, USS, VSS,                       &
-# endif
      &                             IVTYPE, COMPDA, Numcouple, ng, io)
 !
 !=======================================================================
@@ -824,11 +777,6 @@
 
       real :: COMPDA(MCGRD,MCMVAR)
       real :: VOQ(MIP,NVOQP)
-# ifdef SPECTRUM_STOKES
-      real :: KSS(MIP,MSC)
-      real :: USS(MIP,MSC)
-      real :: VSS(MIP,MSC)
-# endif
 !
 !  Local variable declarations.
 !
@@ -837,10 +785,7 @@
       integer :: Istr, Iend, Jstr, Jend, MXCGLc, start
       integer :: Isize, Jsize, INDXG, NPROCS, OFFSET
       integer, pointer :: points(:)
-# ifdef SPECTRUM_STOKES
-      character (len=4) :: SCODE
-      integer :: IZ
-# endif
+
       real(m8), pointer :: avdata(:)
       real(m8), pointer :: DIRE(:)
       real(m8), pointer :: DIRN(:)
@@ -992,62 +937,17 @@
         CALL AttrVect_importRAttr (AttrVect_G(ng)%wav2ocn_AV,           &
      &                             "WQP",avdata)
 !
-# if defined VEGETATION && defined VEG_SWAN_COUPLING \
-      && defined VEG_STREAMING      
+!#if defined VEGETATION && defined VEG_SWAN_COUPLING \
+!      && defined VEG_STREAMING      
+#if defined VEGETATION && defined VEG_SWAN_COUPLING 
+!      && defined VEG_STREAMING      
       ELSE IF (IVTYPE.eq.57) THEN
         CALL AttrVect_importRAttr (AttrVect_G(ng)%wav2ocn_AV,           &
      &                             "DISVEG",avdata)
       ELSE IF (IVTYPE.eq.80) THEN
         CALL AttrVect_importRAttr (AttrVect_G(ng)%wav2ocn_AV,           &
      &                             "DICDVG",avdata)
-# endif   
-# ifdef SPECTRUM_STOKES
-      ELSE IF (IVTYPE.eq.81) THEN
-        DO IZ=1,MSC
-          IF (IZ.LT.10) THEN
-            WRITE(SCODE,"(A3,I1)") "KS0",IZ
-          ELSE
-            WRITE(SCODE,"(A2,I2)") "KS", IZ
-          END IF
-
-          DO IP=1,gsmsize
-            avdata(IP)= REAL(KSS(points(IP),IZ), m8)
-          END DO
-
-          CALL AttrVect_importRAttr (AttrVect_G(ng)%wav2ocn_AV,         &
-     &                             SCODE,avdata)
-        END DO
-
-        DO IZ=1,MSC
-          IF (IZ.LT.10) THEN
-            WRITE(SCODE,"(A3,I1)") "US0",IZ
-          ELSE
-            WRITE(SCODE,"(A2,I2)") "US", IZ
-          END IF
-
-          DO IP=1,gsmsize
-            avdata(IP)= REAL(USS(points(IP),IZ), m8)
-          END DO
-
-          CALL AttrVect_importRAttr (AttrVect_G(ng)%wav2ocn_AV,         &
-     &                             SCODE,avdata)
-        END DO
-
-        DO IZ=1,MSC
-          IF (IZ.LT.10) THEN
-            WRITE(SCODE,"(A3,I1)") "VS0",IZ
-          ELSE
-            WRITE(SCODE,"(A2,I2)") "VS", IZ
-          END IF
-
-          DO IP=1,gsmsize
-            avdata(IP)= REAL(VSS(points(IP),IZ), m8)
-          END DO
-
-          CALL AttrVect_importRAttr (AttrVect_G(ng)%wav2ocn_AV,         &
-     &                             SCODE,avdata)
-        END DO
-# endif
+#endif   
       END IF
 !
       IF (IRQ.eq.Numcouple) THEN
@@ -1322,7 +1222,10 @@
       integer :: idruf, idice
       integer, dimension(MPI_STATUS_SIZE,4) :: status
 
-      real :: cff, cff2
+      real :: cff
+# if !defined BBL_MODEL
+      real :: cff2
+# endif
       real, parameter ::  Large = 1.0E+20
       real, dimension(2) :: range
       real(m8), pointer :: avdata(:)
@@ -1604,7 +1507,7 @@
         CALL SWREDUCE(range(1),1,SWREAL,SWMIN)
         CALL SWREDUCE(range(2),1,SWREAL,SWMAX)
         IF (MyRank.eq.0) THEN
-          write(SCREEN,40) 'ROMStoSWAN Min/Max VEGDENS (indv/m-2): ',   &
+          write(SCREEN,40) 'ROMStoSWAN Min/Max VEGDENS(stms/m2): ',     &
      &                      range(1),range(2)
         END IF
 !
@@ -1953,7 +1856,7 @@
 ! If not using a BBL model, then determine the non-spatially 
 ! varying friction coeff to enter into the JFRC2 array.
 !
-# if defined WAVES_OCEAN
+# if defined WAVES_OCEAN && !defined BBL_MODEL
         cff2=0.05_m8                  ! default
         IF (IBOT.EQ.1) THEN           ! Jonswap
           cff2=PBOT(3)
@@ -1962,20 +1865,20 @@
         ELSE IF (IBOT.EQ.3) THEN      ! Madsen
           cff2=PBOT(5)
         END IF
+# endif
         IP=0
         DO IY=1,MYC
           DO IX=1,MXC
             IP=IP+1
             INDX = KGRPNT(IX,IY)
             IF (INDX.GT.1) THEN
+# ifdef WAVES_OCEAN
 #  ifdef BBL_MODEL
               IF (io.eq.1) THEN
-!               COMPDA(INDX,JFRC2)=MAX(REAL(TEMPMCT(IP,idruf)), 0.0001)
-                COMPDA(INDX,JFRC2)=MAX(REAL(TEMPMCT(IP,idruf)), cff2)
+                COMPDA(INDX,JFRC2)=MAX(REAL(TEMPMCT(IP,idruf)), 0.0001)
               ELSE
                 COMPDA(INDX,JFRC2)=COMPDA(INDX,JFRC2)+                  &
-     &                             MAX(REAL(TEMPMCT(IP,idruf)), cff2)
-!    &                             MAX(REAL(TEMPMCT(IP,idruf)), 0.0001)
+     &                             MAX(REAL(TEMPMCT(IP,idruf)), 0.0001)
               END IF
 #  else
               COMPDA(INDX,JFRC2)=cff2
diff --git a/SWAN/Src/waves_coupler.ftn b/SWAN/Src/waves_coupler.ftn
index c5416b75ad9676bc31fa220dbdb3d6a9ac684cc3..4ba1493ac4ab22906ead8cd4641b23f55583e1eb 100755
--- a/SWAN/Src/waves_coupler.ftn
+++ b/SWAN/Src/waves_coupler.ftn
@@ -1385,7 +1385,7 @@
         CALL SWREDUCE(range(1),1,SWREAL,SWMIN)
         CALL SWREDUCE(range(2),1,SWREAL,SWMAX)
         IF (MyRank.eq.0) THEN
-          write(SCREEN,40) 'ROMStoSWAN Min/Max VEGDENS (indv/m-2): ',   &
+          write(SCREEN,40) 'ROMStoSWAN Min/Max VEGDENS(stms/m2): ',    &
      &                      range(1),range(2)
         END IF
 #  endif