Remove some parameters and hardwired those options: ibtrack_openbnd=1…

…; iwindoff=0; dzb_decay=0

Also removed option for negative roughness

sample_inputs/param.nml

@@ -271,14 +271,12 @@
 !           nchi=1: bottom roughness (in meters) specified in rough.gr3 (and in this case, negative
 !           or 0 depths in rough.gr3 indicate time-independent Cd, not roughness!).
 !           Cd is calculated using the log law, when dzb>=dzb_min; when dzb<dzb_min,
-!           Cd=Cdmax*exp[dzb_decay*(1-dzb/dzb_min)], where Cdmax=Cd(dzb=dzb_min),
-!           and dzb_decay (<=0) is a decay const specified below. We recommend dzb_decay=0
-!           and may remove this constant in the future.
+!           Cd=Cdmax, where Cdmax=Cd(dzb=dzb_min).
 !           If iwbl/=0, nchi must =1.
 !-----------------------------------------------------------------------
   nchi = 0 
   dzb_min = 0.5 !needed if nchi=1; min. bottom boundary layer thickness [m].
-  dzb_decay = 0. !needed if nchi=1; a decay const. [-]. should =0
+!  dzb_decay = 0. !needed if nchi=1; a decay const. [-]. should =0
   hmin_man = 1. !needed if nchi=-1: min. depth in Manning's formulation [m]
 
 !-----------------------------------------------------------------------
@@ -377,7 +375,7 @@
 ! Behavior when trajectory hits open bnd. If ibtrack_openbnd=0, slide with
 ! tangential vel; otherwise, stop and exit btrack 
 !-----------------------------------------------------------------------
-  ibtrack_openbnd = 1 
+!  ibtrack_openbnd = 1 !hardwired
 
 !-----------------------------------------------------------------------
 ! Wetting and drying. 
@@ -475,7 +473,6 @@
 ! variables are read in from wind.th. If nws=2, atmos. variables are
 ! read in from sflux_ files.
 ! If nws=4, ascii format is used for wind and atmos. pressure at each node (see source code).
-! If nws>0, 'iwindoff' can be used to scale wind speed (with windfactor.gr3).
 !
 ! Stress calculation:
 ! If nws=2, ihconsv=1 and iwind_form=0, the stress is calculated from heat exchange
@@ -489,7 +486,7 @@
   wtiminc = 150. !time step for atmos. forcing. Default: same as dt
   nrampwind = 1 !ramp-up option for atmos. forcing
   drampwind = 1. !needed if nrampwind/=0; ramp-up period in days
-  iwindoff = 0 !needed only if nws/=0; '1': needs windfactor.gr3
+!  iwindoff = 0 !needed only if nws/=0; '1': needs windfactor.gr3
   iwind_form = -1 !needed if nws/=0
   !If impose_net_flux/=0 and nws=2, read in net _surface_ heat flux as var 'dlwrf' 
   !(Downward Long Wave) in sflux_rad (solar radiation is still used separately),

src/Core/schism_glbl.F90

@@ -84,14 +84,14 @@ module schism_glbl
   integer,parameter :: mntracers=30 !max # of tracers, used only for dimensioning btrack arrays. Must >=ntracers
 
   !Parameters from param.nml
-  integer,save :: ipre,ipre2,indvel,imm,ihot,ics,iwbl,iharind,nws,impose_net_flux,iwindoff, &
+  integer,save :: ipre,ipre2,indvel,imm,ihot,ics,iwbl,iharind,nws,impose_net_flux, &
                   &ibc,nrampbc,nrampwind,nrampwafo,nramp,nramp_ss,ibdef,ihorcon,nstep_wwm,icou_elfe_wwm, &
                   &fwvor_advxy_stokes,fwvor_advz_stokes,fwvor_gradpress,fwvor_breaking, &
                   &cur_wwm,wafo_obcramp,iwind_form,irec_nu,itur,ihhat,inu_elev, &
                   &inu_uv,ibcc_mean,iflux,iout_sta,nspool_sta,nhot,nhot_write, &
                   &moitn0,mxitn0,nchi,ibtrack_test,nramp_elev,islip,ibtp,inunfl,shorewafo, &
                   &inv_atm_bnd,ieos_type,ieos_pres,iupwind_mom,inter_mom,ishapiro,isav, &
-                  &nstep_ice,niter_shap,iunder_deep,ibtrack_openbnd,flag_fib,ielm_transport,max_subcyc, &
+                  &nstep_ice,niter_shap,iunder_deep,flag_fib,ielm_transport,max_subcyc, &
                   &itransport_only,meth_sink,iloadtide,nc_out
   integer,save :: ntrs(natrm),nnu_pts(natrm),mnu_pts
   integer,save,dimension(:),allocatable :: iof_hydro,iof_wwm,iof_gen,iof_age,iof_sed,iof_eco, &
@@ -102,7 +102,7 @@ module schism_glbl
                       &vdmin_pp1,tdmin_pp1,vdmax_pp2,vdmin_pp2,tdmin_pp2, &
                       &h1_pp,h2_pp,dtb_min,dtb_max,thetai,theta2,rtol0, &
                       &vnh1,vnh2,vnf1,vnf2,rnday,btrack_nudge,hmin_man, &
-                      &prmsl_ref,hmin_radstress,dzb_decay,eos_a,eos_b,eps1_tvd_imp,eps2_tvd_imp, &
+                      &prmsl_ref,hmin_radstress,eos_a,eos_b,eps1_tvd_imp,eps2_tvd_imp, &
                       &xlsc0,rearth_pole,rearth_eq,hvis_coef0,disch_coef(10),hw_depth,hw_ratio, &
                       &slr_rate,rho0,shw,gen_wsett,turbinj,h1_bcc,h2_bcc,vclose_surf_frac, &
                       &hmin_airsea_ex
@@ -442,7 +442,7 @@ module schism_glbl
                                   !WARNING: airt[12] are in C not K. The
                                   !original air T in sflux_air*.nc is in K but
                                   !get_wind() converts it to C
-                                  &tau(:,:),tau_bot_node(:,:),windfactor(:),pr1(:),airt1(:), &
+                                  &tau(:,:),tau_bot_node(:,:),pr1(:),airt1(:), &
                                   &shum1(:),pr2(:),airt2(:),shum2(:),pr(:), &
                                   &sflux(:),srad(:),tauxz(:),tauyz(:),fluxsu(:), &
                                   &fluxlu(:),hradu(:),hradd(:),cori(:), & !chi(:)

src/Hydro/bktrk_subs.F90

@@ -1441,7 +1441,8 @@ subroutine quicksearch(idx,itr,l_ns,ipsgb,gcor0,frame0,time,x0,y0,z0,nnel, &
 
         vtan=-su2(jlev,isd)*sny(isd)+sv2(jlev,isd)*snx(isd)
         !If open bnd is hit, optionally stop with nfl=1
-        if(ibtrack_openbnd/=0.and.isbs(isd)>0) vtan=0._rkind
+        !if(ibtrack_openbnd/=0.and.isbs(isd)>0) vtan=0._rkind
+        if(isbs(isd)>0) vtan=0._rkind
         xvel=-vtan*sny(isd)
         yvel=vtan*snx(isd)
 

src/Hydro/schism_init.F90

@@ -170,9 +170,9 @@ subroutine schism_init(iorder,indir,iths,ntime)
       namelist /OPT/ gen_wsett,flag_fib,ics,rearth_pole,rearth_eq,indvel, &
      &imm,ibdef,ihot,ihydraulics,izonal5,slam0,sfea0,iupwind_mom,ihorcon, &
      &hvis_coef0,ishapiro,shapiro0,niter_shap,ihdif,thetai,nrampbc,drampbc, &
-     &nramp,dramp,nadv,dtb_min,dtb_max,h0,nchi,dzb_min,dzb_decay, &
+     &nramp,dramp,nadv,dtb_min,dtb_max,h0,nchi,dzb_min, &
      &hmin_man,ncor,rlatitude,coricoef,nws,impose_net_flux,wtiminc,iwind_form,nrampwind, &
-     &drampwind,iwindoff,ihconsv,isconsv,itur,dfv0,dfh0,h1_pp,h2_pp,vdmax_pp1, &
+     &drampwind,ihconsv,isconsv,itur,dfv0,dfh0,h1_pp,h2_pp,vdmax_pp1, &
      &vdmax_pp2,vdmin_pp1,vdmin_pp2,tdmin_pp1,tdmin_pp2,mid,stab,xlsc0, &
      &ibcc_mean,flag_ic,start_year,start_month,start_day,start_hour,utc_start, &
      &itr_met,h_tvd,eps1_tvd_imp,eps2_tvd_imp,ip_weno, &
@@ -185,7 +185,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
      &fwvor_advxy_stokes,fwvor_advz_stokes,fwvor_gradpress,fwvor_breaking,wafo_obcramp, &
      &iwbl,cur_wwm,if_source,nramp_ss,dramp_ss,ieos_type,ieos_pres,eos_a,eos_b,slr_rate, &
      &rho0,shw,isav,nstep_ice,iunder_deep,h1_bcc,h2_bcc,hw_depth,hw_ratio, &
-     &ibtrack_openbnd,level_age,vclose_surf_frac,iadjust_mass_consv0,ipre2, &
+     &level_age,vclose_surf_frac,iadjust_mass_consv0,ipre2, &
      &ielm_transport,max_subcyc,i_hmin_airsea_ex,hmin_airsea_ex,itransport_only,meth_sink, &
      &iloadtide
 
@@ -279,8 +279,6 @@ subroutine schism_init(iorder,indir,iths,ntime)
       if(ibtp/=0.and.ibtp/=1) call parallel_abort('Unknown ibtp')
       if(ibc==0) then
         if(myrank==0) write(16,*)'You are using baroclinic model'
-!        call get_param('param.in','nrampbc',1,nrampbc,tmp,stringvalue)
-!        if(nrampbc/=0) call get_param('param.in','drampbc',2,itmp,drampbc,stringvalue)
       else !ibc=1
         if(ibtp==0) then
           if(myrank==0) write(16,*)'Barotropic model without ST calculation'
@@ -437,10 +435,10 @@ subroutine schism_init(iorder,indir,iths,ntime)
       gen_wsett=real(0.d0,rkind); flag_fib=1; ics=1; rearth_pole=6378206.4_rkind; rearth_eq=6378206.4_rkind; 
       imm=0; ibdef=10; ihot=0; ihydraulics=0; izonal5=0; slam0=-124._rkind; sfea0=45._rkind; 
       ihdif=0; thetai=0.6_rkind; nrampbc=0; drampbc=1._rkind;  
-      nramp=1; dramp=1._rkind; nadv=1; dtb_min=10._rkind; dtb_max=30._rkind; h0=0.01_rkind; nchi=0; dzb_min=0.5_rkind; dzb_decay=0._rkind;  
+      nramp=1; dramp=1._rkind; nadv=1; dtb_min=10._rkind; dtb_max=30._rkind; h0=0.01_rkind; nchi=0; dzb_min=0.5_rkind 
       hmin_man=1._rkind; ncor=0; rlatitude=46._rkind; coricoef=0._rkind; 
       nws=0; impose_net_flux=0; wtiminc=dt; iwind_form=-1; nrampwind=1; 
-      drampwind=1; iwindoff=0; ihconsv=0; isconsv=0; i_hmin_airsea_ex=2; itur=0; dfv0=0.01_rkind; dfh0=real(1.d-4,rkind); 
+      drampwind=1; ihconsv=0; isconsv=0; i_hmin_airsea_ex=2; itur=0; dfv0=0.01_rkind; dfh0=real(1.d-4,rkind); 
       h1_pp=20._rkind; h2_pp=50._rkind; vdmax_pp1=0.01_rkind; vdmax_pp2=0.01_rkind
       vdmin_pp1=real(1.d-5,rkind); vdmin_pp2=vdmin_pp1; tdmin_pp1=vdmin_pp1; tdmin_pp2=vdmin_pp1
       mid='KL'; stab='KC'; xlsc0=0.1_rkind;  
@@ -459,7 +457,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
       fwvor_advxy_stokes=1; fwvor_advz_stokes=1; fwvor_gradpress=1; fwvor_breaking=1; wafo_obcramp=0;
       iwbl=0; cur_wwm=0; if_source=0; nramp_ss=1; dramp_ss=2._rkind; ieos_type=0; ieos_pres=0; eos_a=-0.1_rkind; eos_b=1001._rkind;
       slr_rate=120._rkind; rho0=1000._rkind; shw=4184._rkind; isav=0; nstep_ice=1; h1_bcc=50._rkind; h2_bcc=100._rkind
-      hw_depth=1.d6; hw_ratio=0.5d0; iunder_deep=0; ibtrack_openbnd=1; level_age=-999;
+      hw_depth=1.d6; hw_ratio=0.5d0; iunder_deep=0; level_age=-999;
       !vclose_surf_frac \in [0,1]: correction factor for vertical vel & flux. 1: no correction
       vclose_surf_frac=1.0
       iadjust_mass_consv0=0 !Enforce mass conservation for a tracer 
@@ -663,9 +661,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
 
       if(nchi==1) then
 !       dzb_min: min. bottom boundary layer thickness [m]
-!        call get_param('param.in','dzb_min',2,itmp,dzb_min,stringvalue)
-!        call get_param('param.in','dzb_decay',2,itmp,dzb_decay,stringvalue)
-        if(dzb_min<=0._rkind.or.dzb_decay>0._rkind) call parallel_abort('INIT: dzb_min<=0 or dzb_decay>0')
+        if(dzb_min<=0._rkind) call parallel_abort('INIT: dzb_min<=0') 
       endif
       if(nchi==-1) then
 !       Min depth used in Manning formulation
@@ -723,7 +719,6 @@ subroutine schism_init(iorder,indir,iths,ntime)
 !      if(nws>0) then
 !        call get_param('param.in','nrampwind',1,nrampwind,tmp,stringvalue)
 !        call get_param('param.in','drampwind',2,itmp,drampwind,stringvalue)
-!        call get_param('param.in','iwindoff',1,iwindoff,tmp,stringvalue)
 !      endif !nws
 
 !     Heat and salt conservation flags
@@ -1712,7 +1707,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
 !     All other arrays
 !      allocate(sdbt(2+ntracers,nvrt,nsa), & !webt(nvrt,nea), bubt(2,nea), & 
        allocate(windx1(npa),windy1(npa),windx2(npa),windy2(npa),windx(npa),windy(npa), &
-         &  tau(2,npa),tau_bot_node(3,npa),iadv(npa),windfactor(npa),pr1(npa),airt1(npa),shum1(npa), &
+         &  tau(2,npa),tau_bot_node(3,npa),iadv(npa),pr1(npa),airt1(npa),shum1(npa), &
          &  pr2(npa),airt2(npa),shum2(npa),pr(npa),sflux(npa),srad(npa),tauxz(npa),tauyz(npa), &
          &  fluxsu(npa),fluxlu(npa),hradu(npa),hradd(npa),cori(nsa),Cd(nsa), &
          &  Cdp(npa),rmanning(npa),rough_p(npa),dfv(nvrt,npa),elev_nudge(npa),uv_nudge(npa), &
@@ -3202,6 +3197,10 @@ subroutine schism_init(iorder,indir,iths,ntime)
      &call parallel_abort('Check rough.gr3')
         do i=1,np_global
           read(32,*)j,xtmp,ytmp,tmp
+          if(tmp<0.d0) then
+            write(errmsg,*)'INIT: negative rough at node ',i,tmp
+            call parallel_abort(errmsg)
+          endif
           if(ipgl(i)%rank==myrank) rough_p(ipgl(i)%id)=tmp
         enddo !i
         close(32)
@@ -3273,25 +3272,25 @@ subroutine schism_init(iorder,indir,iths,ntime)
 #endif
       endif
 
-      windfactor=1 !intialize for default
-      if(nws>0) then
-        if(iwindoff/=0) then
-          open(32,file=in_dir(1:len_in_dir)//'windfactor.gr3',status='old')
-          read(32,*)
-          read(32,*) itmp1,itmp2
-          if(itmp1/=ne_global.or.itmp2/=np_global) &
-     &call parallel_abort('Check windfactor.gr3')
-          do i=1,np_global
-            read(32,*)j,xtmp,ytmp,tmp
-            if(tmp<0.d0) then
-              write(errmsg,*)'Wind scaling factor must be positive:',i,tmp
-              call parallel_abort(errmsg)
-            endif
-            if(ipgl(i)%rank==myrank) windfactor(ipgl(i)%id)=tmp
-          enddo !i
-          close(32)
-        endif
-      endif !nws>0
+!      windfactor=1 !intialize for default
+!      if(nws>0) then
+!        if(iwindoff/=0) then
+!          open(32,file=in_dir(1:len_in_dir)//'windfactor.gr3',status='old')
+!          read(32,*)
+!          read(32,*) itmp1,itmp2
+!          if(itmp1/=ne_global.or.itmp2/=np_global) &
+!     &call parallel_abort('Check windfactor.gr3')
+!          do i=1,np_global
+!            read(32,*)j,xtmp,ytmp,tmp
+!            if(tmp<0.d0) then
+!              write(errmsg,*)'Wind scaling factor must be positive:',i,tmp
+!              call parallel_abort(errmsg)
+!            endif
+!            if(ipgl(i)%rank==myrank) windfactor(ipgl(i)%id)=tmp
+!          enddo !i
+!          close(32)
+!        endif
+!      endif !nws>0
 
 !     Alloc. the large array for nws=4-6 option (may consider changing to unformatted binary read)
 !      if(nws==4) then

src/Hydro/schism_step.F90

@@ -183,7 +183,7 @@ subroutine schism_step(it)
                      &tot_salt_gb,dav_mag,tvol,tmass,tpe,tkne,enerf,ener_ob, &
                      &av_dep,vel_m1,vel_m2,xtmp,ytmp,ftmp,tvol12,fluxbnd, &
                      &fluxchan,fluxchan1,fluxchan2,tot_s,flux_s,ah,ubm,ramp_ss,Cdmax, &
-                     &wmag_e,wmag_factor, & !wmag_e and wmag_facotr addedy by wangzg
+                     &wmag_e, & !wmag_e and wmag_facotr addedy by wangzg
                      &bthick_ori,big_ubstar,big_vbstar,zsurf,tot_bedmass,w1,w2,slr_elev, &
                      &i34inv,av_cff1,av_cff2,av_cff3,av_cff2_chi,av_cff3_chi, &
                      &sav_cfk,sav_cfpsi,sav_h_sd,sav_alpha_sd,sav_nv_sd,sav_c,beta_bar, &
@@ -585,21 +585,21 @@ subroutine schism_step(it)
       endif !nws>=2
 
 !...  Re-scale wind
-      if(nws>0) then; if(iwindoff/=0) then
-        do i=1,npa
-          windx(i)=windx(i)*windfactor(i)
-          windy(i)=windy(i)*windfactor(i)
-        enddo !i
-      endif; endif
+!      if(nws>0) then; if(iwindoff/=0) then
+!        do i=1,npa
+!          windx(i)=windx(i)*windfactor(i)
+!          windy(i)=windy(i)*windfactor(i)
+!        enddo !i
+!      endif; endif
 
 #ifdef USE_ICM
 ! calculating WMS used for reareation,added by wangzg
       if(irea==1) then
-!$OMP parallel do default(shared) private(i,wmag_e,wmag_factor)
+!$OMP parallel do default(shared) private(i,wmag_e)
         do i=1,nea
           if(idry_e(i)==1) cycle
           wmag_e=sum(sqrt(windx(elnode(1:i34(i),i))**2.d0+windy(elnode(1:i34(i),i))**2.d0))/real(i34(i),rkind)
-          wmag_factor=sum(windfactor(elnode(1:i34(i),i)))/real(i34(i),rkind)
+!          wmag_factor=sum(windfactor(elnode(1:i34(i),i)))/real(i34(i),rkind)
           WMS(i)=wmag_e !*wmag_factor !no windfactor for DO reareation
         enddo !i
 !$OMP end parallel do
@@ -759,8 +759,8 @@ subroutine schism_step(it)
             tau(1,i)=0.d0
             tau(2,i)=0.d0
           else !rescale as well
-            tau(1,i)=-tauxz(i)/rho0*rampwind*windfactor(i)**2.d0 !sign and scale difference between stresses tauxz and tau
-            tau(2,i)=-tauyz(i)/rho0*rampwind*windfactor(i)**2.d0
+            tau(1,i)=-tauxz(i)/rho0*rampwind !*windfactor(i)**2.d0 !sign and scale difference between stresses tauxz and tau
+            tau(2,i)=-tauyz(i)/rho0*rampwind !*windfactor(i)**2.d0
           endif
         else !if(nws==1.or.nws>=4.or.nws>=2.and.ihconsv==0.or.iwind_form==-1) then
           wmag=sqrt(windx(i)**2.d0+windy(i)**2.d0)
@@ -1897,8 +1897,12 @@ subroutine schism_step(it)
 
 !         Wet node
           htot=dp(i)+eta2(i)
-          if(rough_p(i)<=0.d0) then !time-independent Cd
-            Cdp(i)=abs(rough_p(i))
+          if(rough_p(i)<0.d0) then 
+            !Cdp(i)=abs(rough_p(i))
+            write(errmsg,*)'STEP: rough_p<0 at node ',iplg(i),rough_p(i)
+            call parallel_abort(errmsg)
+          else if(rough_p(i)==0.d0) then 
+            Cdp(i)=0.d0
           else !roughness >0 
             bthick_ori=znl(kbp(i)+1,i)-znl(kbp(i),i)  !thickness of bottom bnd layer
             bthick=max(dzb_min,bthick_ori)
@@ -1914,9 +1918,9 @@ subroutine schism_step(it)
             else
               Cdp(i)=1.d0/(2.5d0*log(bthick/rough_p(i)))**2.d0
 
-              if(dzb_decay/=0.d0.and.bthick_ori<bthick) then !dzb_decay=0 leads to no decay
-                Cdp(i)=Cdp(i)*exp(dzb_decay*(1.d0-bthick_ori/bthick))
-              endif
+!              if(dzb_decay/=0.d0.and.bthick_ori<bthick) then !dzb_decay=0 leads to no decay
+!                Cdp(i)=Cdp(i)*exp(dzb_decay*(1.d0-bthick_ori/bthick))
+!              endif
               !WBL
 #ifdef USE_WWM
               ! Quantities used in both formulations for the WBL
@@ -7169,10 +7173,9 @@ subroutine schism_step(it)
         if(nwild(i)==1) then !marsh elem
           if(smax>0.5d0) then !drowned
             imarsh(i)=0
-!Error: OMP race
-            Cdp(elnode(1:i34(i),i))=0.001d0
-            Cd(elside(1:i34(i),i))=0.001d0
-            rough_p(elnode(1:i34(i),i))=1.d-4
+!            Cdp(elnode(1:i34(i),i))=0.001d0
+!            Cd(elside(1:i34(i),i))=0.001d0
+!            rough_p(elnode(1:i34(i),i))=1.d-4
           endif !smax
         else !non-marsh elem @last step
           if(smax<=0.d0.and.smin>=-1.d0) then !create marsh
@@ -7192,27 +7195,62 @@ subroutine schism_step(it)
       enddo !i=1,ne
 !$OMP end do
 
-      !Set Cd etc for marsh
+      !Set Cd etc for marsh and also drowned marsh
 !$OMP workshare
       sav_di=0.d0; sav_h=0.d0; sav_nv=0.d0; sav_alpha=0.d0
 !$OMP end workshare
 !$OMP do 
-      do i=1,ne
-        if(imarsh(i)==1) then
-          if(isav==0) then
-            Cdp(elnode(1:i34(i),i))=0.05d0
-            Cd(elside(1:i34(i),i))=0.05d0
-            rough_p(elnode(1:i34(i),i))=1.d-2
-          else
-            sav_di(elnode(1:i34(i),i))=sav_di0
-            sav_h(elnode(1:i34(i),i))=sav_h0
-            sav_nv(elnode(1:i34(i),i))=sav_nv0
-            sav_alpha(elnode(1:i34(i),i))=sav_di0*sav_nv0*sav_cd/2.d0
-          endif !isav
-        endif
+      do i=1,np
+        do j=1,nne(i)
+          ie=indel(j,i)
+          if(imarsh(ie)==1) then
+            if(isav==0) then
+              Cdp(i)=0.05d0
+              rough_p(i)=1.d-2
+            else !isav/=0
+              sav_di(i)=sav_di0
+              sav_h(i)=sav_h0
+              sav_nv(i)=sav_nv0
+              sav_alpha(i)=sav_di0*sav_nv0*sav_cd0/2.d0
+            endif !isav 
+          endif !imarsh
+
+          !drowned marsh
+          if(nwild(ie)==1.and.imarsh(ie)==0) then
+            Cdp(i)=0.001d0
+            rough_p(i)=1.d-4
+          endif
+        enddo !j
       enddo !i
 !$OMP end do
 
+!$OMP do
+      do i=1,ns
+        do j=1,2
+          ie=isdel(j,i)
+          if(isav==0.and.imarsh(ie)==1) Cd(i)=0.05d0
+          if(nwild(ie)==1.and.imarsh(ie)==0) Cd(i)=0.001d0
+        enddo !j
+      enddo !i
+!$OMP end do
+
+
+!      do i=1,ne
+!        if(imarsh(i)==1) then
+!          if(isav==0) then
+!            Cdp(elnode(1:i34(i),i))=0.05d0
+!            Cd(elside(1:i34(i),i))=0.05d0
+!            rough_p(elnode(1:i34(i),i))=1.d-2
+!          else
+!            sav_di(elnode(1:i34(i),i))=sav_di0
+!            sav_h(elnode(1:i34(i),i))=sav_h0
+!            sav_nv(elnode(1:i34(i),i))=sav_nv0
+!            sav_alpha(elnode(1:i34(i),i))=sav_di0*sav_nv0*sav_cd0/2.d0
+!          endif !isav
+!        endif
+!      enddo !i
+!$OMP end do
+
 !$OMP master
       call exchange_p2d(Cdp)
       call exchange_p2d(rough_p)