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