update EPA MBM model setup

Author: wzhengui

src/Utility/Pre-Processing/EPA_MBM/param/param.nml

@@ -157,7 +157,7 @@
 ! For applications that include eddying regime, refer to the manual.
 !-----------------------------------------------------------------------
   ihorcon = 2
-  hvis_coef0 = 0.008 !const. diffusion # if ihorcon/=0; <=0.025 for ihorcon=2, <=0.125 for ihorcon=1
+  hvis_coef0 = 0.015 !const. diffusion # if ihorcon/=0; <=0.025 for ihorcon=2, <=0.125 for ihorcon=1
 !  cdh = 0.01 !needed only if ihorcon/=0; land friction coefficient - not active yet
 
 !-----------------------------------------------------------------------
@@ -606,9 +606,9 @@
 ! stemp_stc=0 would turn off this effect.
 !-----------------------------------------------------------------------
   istemp =1
-  #stemp_stc  = 5.0   !heat transfer coefficient W.m-2.K-1
-  #stemp_dz(1)= 2.0   !equivalent sediment buffer depth (m) for heat into sediment
-  #stemp_dz(2)= 0.1   !equivalent sediment buffer depth (m) for heat out of sediment
+  !stemp_stc  = 5.0   !heat transfer coefficient W.m-2.K-1
+  !stemp_dz(1)= 2.0   !equivalent sediment buffer depth (m) for heat into sediment
+  !stemp_dz(2)= 0.1   !equivalent sediment buffer depth (m) for heat out of sediment
 
 !-----------------------------------------------------------------------
 ! Turbulence closure.
@@ -800,7 +800,7 @@
 ! 0.0 <= vclose_surf_frac < 1.0
 ! 1: fully from surface (i.e. no correction as before); 0: fully from bottom 
 !-----------------------------------------------------------------------
-  vclose_surf_frac=1.0
+!  vclose_surf_frac=1.0
 
 !-----------------------------------------------------------------------
 ! Option to enforce strict mass conservation for each tracer model (only works with itr_met=3,4)

src/Utility/Pre-Processing/EPA_MBM/param/param.npz

@@ -44,6 +44,7 @@
 p.flag['diffmin.gr3']          = 0 #hydro
 p.flag['diffmax.gr3']          = 0 #hydro
 p.flag['shapiro.gr3']          = 0 #hydro
+p.flag['soil_*.gr3']           = 0 #hydro
 p.flag['tvd.prop']             = 0 #hydro
 p.flag['windrot_geo2proj.gr3'] = 0 #hydro
 p.flag['veg_*.gr3']            = 0 #hydro: veg module (SAV drag effect)
@@ -156,7 +157,7 @@ def pattr(x,a): p.flag[x]=a
    if p.flag['ICM_nu.nc']==1: pm['inu_tr(7)']=2; pm['vnf1']=0.01
    if p.flag['SED']==1 and p.flag['WWM']==1: #for offline transport mode
       for i in [13,21,27,29,30]:  pm['iof_hydro({})'.format(i)]=1
-   s=num2date(p.StartT); pm['start_year']=s.year; pm['start_month']=s.month; pm['start_day']=s.day; pm['rnday']=int(p.EndT-p.StartT)
+   s=num2date(p.StartT); pm['start_year']=s.year; pm['start_month']=s.month; pm['start_day']=s.day; pm['rnday']=int(p.EndT-p.StartT-1)
 if p.flag['ICM'] in [10,20]: #offline ICM mode
    pm['nspool']=int(6*3600/p.dt_offline); pm['ihfskip']=int(720*3600/p.dt_offline)
    pm['nhot_write']=pm['ihfskip']; pm['dt']=p.dt_offline; pm['wtiminc']=pm['dt']
@@ -327,8 +328,7 @@ def pattr(x,a): p.flag[x]=a
 fname='diffmin.gr3'
 if p.flag[fname]==1:
    print('writing '+fname)
-   vi=ones(gd.np)*1e-6; vi[read(p.region+'lowerbay.reg').inside(gd.xy)]=8e-5
-   gd.save(fname,value=vi)
+   gd.save(fname,value=1e-6)
 
 fname='shapiro.gr3'
 if p.flag[fname]==1:
@@ -337,6 +337,12 @@ def pattr(x,a): p.flag[x]=a
    slope=gd.compute_gradient(fmt=2)[2]; vi=shapiro_max*tanh(2*slope/threshold_slope); vi[slope<0.005]=0
    gd.save('shapiro.gr3',value=vi)
 
+fname='soil_*.gr3'
+if p.flag[fname]==1:
+   print('writing '+fname)
+   gd.save('soil_conductivity.gr3',value=5)
+   vi=ones(gd.np); vi[gd.z>10]=0.1; gd.save('soil_thick.gr3',value=vi)
+
 fname='veg_*.gr3'
 if p.flag[fname]==1:
    print('writing '+fname)
@@ -604,10 +610,10 @@ def pattr(x,a): p.flag[x]=a
    if p.flag['ICM'] in [2,20]: vnames,svars,rats=vnames[:-4],svars[:-4],rats[:-4]
 
    #interp
-   fpt=(M.CB74_time>=p.StartT)*(M.CB74_time<p.EndT+2); mti=M.CB74_time[fpt]; cz=M.CB74_depth
+   i1=abs(M.CB74_time-p.StartT).argmin(); i1=abs(M.CB74_time-p.EndT-2).argmin(); mti=M.CB74_time[i1:i2]; cz=M.CB74_depth
    bind=pindex(read('ICM_nudge.gr3').z!=0); lz=abs(compute_zcor(vd.sigma[bind],gd.dp[bind])).T.clip(cz.min(),cz.max())
    nobn=bind.size; ntr=len(svars); nt=len(mti); nvrt=vd.nvrt
-   mys=array([M.attr('CB74_'+svar)[:,fpt]*rat for svar, rat in zip(svars,rats)]).transpose([1,0,2]) #interp in time
+   mys=array([M.attr('CB74_'+svar)[:,i1:i2]*rat for svar, rat in zip(svars,rats)]).transpose([1,0,2]) #interp in time
 
    #interp vertically, and write ICM_nu.nc 
    from netCDF4 import Dataset; dt=30; irec=0
@@ -642,11 +648,11 @@ def pattr(x,a): p.flag[x]=a
    #read data
    if p.flag['ICM'] in [2,20]: svars,rats=svars[:-4],srats[:-4]
    S=read(p.atmdep,1); S.NO3=S.wetno3+S.dryno3; S.NH4=S.wetnh3+S.drynh3; S.PO4=S.wetpo4; S.ORGN=S.wetorn; S.ORGP=S.wetorp
-   fps=S.sind[near_pts(gd.xy,S.xy[:])]; fpt=(S.time>=p.StartT)*(S.time<=(p.EndT+2))
-   nt=sum(fpt); ntr=len(svars); trs=zeros([ntr,gd.np,nt],'float32')
+   fps=S.sind[near_pts(gd.xy,S.xy[:])]; i1=abs(S.time-p.StartT).argmin(); i2=abs(S.time-p.EndT-2).argmin()
+   nt=(i2-i1); ntr=len(svars); trs=zeros([ntr,gd.np,nt],'float32')
    for m,[svar,srat] in enumerate(zip(svars,srats)):
        if srat==0: continue
-       trs[m]=S.attr(svar)[fps][:,fpt]*srat
+       trs[m]=S.attr(svar)[fps][:,i1:i2]*srat
 
    #write netcdf
    C=zdata(); C.dimname=['node','time','ntr','one']; C.dims=[gd.np,nt,ntr,1]
@@ -810,7 +816,7 @@ def pattr(x,a): p.flag[x]=a
    #-----------------------------------
    #watershed loading
    #-----------------------------------
-   fpn=sort(pindex((C.time>=p.StartT)*(C.time<(p.EndT+2)))); mtime=C.time[fpn]
+   fpn=sort(pindex((C.time>=(p.StartT-1e-4))*(C.time<(p.EndT+2)))); mtime=C.time[fpn]
    nt,ntr,nsource=len(mtime),len(svars),len(sinde); flow=zeros([nt,nsource]); trs=zeros([nt,ntr,nsource])
    sflow=C.flow[C.river=='susquehanna'].sum(axis=0)[fpn] #susquehanna river flow
 
@@ -873,7 +879,7 @@ def pattr(x,a): p.flag[x]=a
    #--------------------------------------------------------------------------
    evars0=['clay1','clay2','silt', 'sand', 'RPOC', 'SRPOC', 'RPON', 'SRPON', 'RPOP', 'SRPOP', 'PIP']
    evars= ['clay' ,'clay' ,'silt', 'sand', 'RPOC', 'G3OC',  'RPON', 'G3ON',  'RPOP', 'G3OP',  'PIP']
-   fps=sort(pindex((C.sho_time>=p.StartT)*(C.sho_time<(p.EndT+2)))); fpt=near_pts(C.sho_time[fps],mtime)
+   fps=sort(pindex((C.sho_time>=(p.StartT-1e-4))*(C.sho_time<(p.EndT+2)))); fpt=near_pts(C.sho_time[fps],mtime)
    for i, id in enumerate(near_pts(C.sho_sxy,gd.exy[sinde])):
        for evar0,evar in zip(evars0,evars):
            rat=0.5 if evar=='clay' else 1