Output of viscosity coefficient from Ziriv model added

setup/cuette_setup.dat

@@ -145,7 +145,7 @@ massflux 0
 ifrad    0
 ifbubble 1
 sedim    0
-dyn_pgrad 0
+dyn_pgrad 6
 pgrad     0.
 nManning  2.5E-2
 horvisc  0.

source/model/init.f90

@@ -136,7 +136,10 @@
  allocate (S(1:M),Gen(1:M),F(1:M),TKE_turb_trans(1:M),Gen_seiches(1:M))
  allocate (KT(1:M+1),k2(1:M+1),u1(1:M+1),v1(1:M+1), &
  & E1(1:M+1), eps1(1:M+1), wArea(0:M+1), w(0:M+1))
- if (dyn_pgrad%par == 5) allocate (eps_ziriv(1:M))
+ if (dyn_pgrad%par == 5 .or. dyn_pgrad%par == 6) then 
+   allocate (eps_ziriv     (1:M))
+   allocate (viscturb_ziriv(1:M))
+ endif
  allocate (WU_(1:M+1),WV_(1:M+1),GAMT(1:M+1),GAMU(1:M+1), &
  &  GAMV(1:M+1),TF(1:M+1),KLengu(1:M+1)) 
  allocate (num(1:M+1)) 

source/model/lake.f90

@@ -2490,9 +2490,9 @@ surfrad1 = surfrad
 h2_out = h2
 
 !Diagnoctics
-if (dyn_pgrad%par == 5) then !diagnostic computation of TKE dissipation in a river scenario
+if (dyn_pgrad%par == 5 .or. dyn_pgrad%par == 6) then !diagnostic computation of turbulence in a river scenario
   call ZIRIVMODEL(M,h1,sign(1._ireals,tauxw)*sqrt(abs(tauxw)/row0), &
-  & g*tan(pi*alphax/180.),z_half,eps_ziriv)
+  & g*tan(pi*alphax/180.),z_half,eps_ziriv,viscturb_ziriv)
 endif
 
 
@@ -2787,7 +2787,7 @@ if (flag_print) then ! The output in ASCII files
     &  eps1, z_half, M, &
     &  zgrid_out, ngrid_out%par, & !ngrid_out
     &  outpath, 'epsilon', i, ndec, .false.)
-    if (dyn_pgrad%par == 5) then
+    if (dyn_pgrad%par == 5 .or. dyn_pgrad%par == 6) then
       ndec = -1
       i = i + 2
       call PROFILE_OUTPUT &
@@ -2797,6 +2797,15 @@ if (flag_print) then ! The output in ASCII files
       &  eps_ziriv, z_half, M, &
       &  zgrid_out, ngrid_out%par, & !ngrid_out
       &  outpath, 'epsziriv', i, ndec, .false.)
+      ndec = -1
+      i = i + 2
+      call PROFILE_OUTPUT &
+      & (ix, iy, dnx, dny, &
+      &  year, month, day, hour, &
+      &  time, dt_out%par, &
+      &  viscturb_ziriv, z_half, M, &
+      &  zgrid_out, ngrid_out%par, & !ngrid_out
+      &  outpath, 'viscturbziriv', i, ndec, .false.)
     endif
     ndec = -1
     i = i + 2

source/model/lake_modules.f90

@@ -568,7 +568,7 @@ real(kind=ireals), allocatable :: E2(:),eps2(:), &
 & Eeps(:),Ecent(:),Epscent(:),res_E(:), &
 & res_eps(:),dresd(:,:,:,:),dres2dE(:),dres2deps(:), &
 & veg_sw(:) 
-real(kind=ireals), allocatable :: eps_ziriv(:)
+real(kind=ireals), allocatable :: eps_ziriv(:), viscturb_ziriv(:)
 real(kind=ireals), allocatable, target :: row(:), row2(:), rowc(:)
 real(kind=ireals), allocatable :: WU_(:),WV_(:),GAMT(:),GAMU(:),GAMV(:), TF(:),KLengu(:)
 real(kind=ireals), allocatable :: PEMF    (:) , PDENS_DOWN (:), &

source/model/out_mod.f90

@@ -417,7 +417,7 @@ real(kind=ireals), allocatable :: accum_var_scalar(:,:,:)
 real(kind=ireals), external :: DZETA
 
 ! Integers
-integer(kind=iintegers), parameter :: n_var = 26
+integer(kind=iintegers), parameter :: n_var = 28
 integer(kind=iintegers), parameter :: n_var_scalar = 10
 integer(kind=iintegers), allocatable :: hour_old(:,:)
 integer(kind=iintegers), allocatable :: tsteps(:,:)
@@ -504,6 +504,10 @@ do i = 1, M+1
 enddo
 var(26,ix,iy,1:M+1)  = qwater(1:M+1,1)*molm3tomcM
 
+!Coefficients of turbulent viscosity
+var(27,ix,iy,1:M)  = k2             (1:M)
+var(28,ix,iy,1:M)  = viscturb_ziriv (1:M)
+
 var_scalar(1,ix,iy) = S_integr_positive
 var_scalar(2,ix,iy) = S_integr_negative
 var_scalar(3,ix,iy) = Gen_integr
@@ -536,7 +540,7 @@ if (int(hour)/=hour_old(ix,iy)) then
   write (out_unit,*) '3 - salinity, kg/kg' 
   write (out_unit,*) '4 - water density, kg/m**3'  
   write (out_unit,*) '5 - turbulent kinetic energy, normalized'
-  write (out_unit,*) '6 - disspation rate, normalized'
+  write (out_unit,*) '6 - dissipation rate, normalized'
   write (out_unit,*) '7 - eddy viscosity for TKE, m**2/s'
   write (out_unit,*) '8  - mass flux, m/s'
   write (out_unit,*) '9  - downdraft temperature, C' 
@@ -548,8 +552,10 @@ if (int(hour)/=hour_old(ix,iy)) then
   write (out_unit,*) '15 - oxygen, mg/l'
   write (out_unit,*) '16 - co2, mcmol/l'
   write (out_unit,*) '17 - ch4, mcmol/l'
+  write (out_unit,*) '18 - eddy viscosity from k-epsilon model, m**2/s'
+  write (out_unit,*) '19 - eddy viscosity from Ziriv model, m**2/s'
   do i=1,M 
-    write (out_unit,'(f14.6,f12.5,e10.3,f10.3,3e10.3,4e14.7,3f11.5,3e14.7)')      &
+    write (out_unit,'(f14.6,f12.5,e10.3,f10.3,3e10.3,4e14.7,3f11.5,3e14.7,2e10.3)')      &
     & -dzeta_int(i)*h1                         /accum_var(12,ix,iy,1),  &
     & (accum_var(8,ix,iy,i) - scale_output%par*maxval(accum_var(8,ix,iy,:)) ) &
     &  /accum_var(15,ix,iy,1), &
@@ -568,7 +574,9 @@ if (int(hour)/=hour_old(ix,iy)) then
     & accum_var(23,ix,iy,i),                                              &
     & accum_var(24,ix,iy,i),                                              &
     & accum_var(25,ix,iy,i),                                              &
-    & accum_var(26,ix,iy,i)
+    & accum_var(26,ix,iy,i),                                              &
+    & accum_var(27,ix,iy,i),                                              &
+    & accum_var(28,ix,iy,i)
   enddo 
   close (out_unit)
   

source/model/turb_mod.f90

@@ -1848,7 +1848,7 @@ enddo
 END SUBROUTINE ED_TEMP_HOSTETLER2
 
 !> An analytical model for neutrally stratified river flow using Obukhov coordinate
-SUBROUTINE ZIRIVMODEL(M,h1,u_star_surf,gradp,z_half,eps_ziriv)
+SUBROUTINE ZIRIVMODEL(M,h1,u_star_surf,gradp,z_half,eps_ziriv,viscturb_ziriv)
 use PHYS_CONSTANTS, only : kappa
 use NUMERIC_PARAMS, only : pi
 implicit none
@@ -1859,6 +1859,7 @@ real(kind=ireals), intent(in ) :: u_star_surf !> Friction velocity in water at t
 real(kind=ireals), intent(in ) :: gradp       !> Longitudinal pressure gradient, Pa/m
 real(kind=ireals), intent(in ) :: z_half(1:M) !> A set of depths where analytical solution is to be evaluated
 real(kind=ireals), intent(out) :: eps_ziriv(1:M) !> TKE dissipation rate, m**2/s**3
+real(kind=ireals), intent(out) :: viscturb_ziriv(1:M) !> Coefficient of turbulent viscosity, m**2/s
 !Local variables
 integer(kind=iintegers) :: i !Loop index
 real(kind=ireals) :: u_star, z_Obukhov, h_max
@@ -1870,6 +1871,7 @@ if (gradp*u_star_surf > 0.) then
     z_Obukhov = sin(pi*z_half(i)/h1)*h1/pi
     u_star = sqrt(abs(gradp)*z_half(i) + u_star_surf**2)
     eps_ziriv(i) = u_star**3/(kappa*z_Obukhov)
+    viscturb_ziriv(i) = kappa*u_star*z_Obukhov
   enddo
 else
   !Wind forcing opposite to the pressure gradient, 
@@ -1878,6 +1880,7 @@ else
     z_Obukhov = sin(pi*z_half(i)/h1)*h1/pi
     u_star = sqrt(abs(abs(gradp)*z_half(i) - u_star_surf**2))
     eps_ziriv(i) = u_star**3/(kappa*z_Obukhov)
+    viscturb_ziriv(i) = kappa*u_star*z_Obukhov
   enddo
 endif
 

tools/plot_hourly.py

@@ -28,7 +28,9 @@ plt.rc('text.latex',preamble='\usepackage[russian]{babel}')
 #        '../results/windforc_kor/hourly/', \
 #        '../results/windforc_dynpgrad4/hourly/'] # path to files
 #path = ['../results/Mojai2016-2017_damw5/hourly/']#,'../results/Mojai2016-2017_damtribheat=0/hourly/']
-path = ['../results/river/hourly/']
+path = ['../results/cuette/hourly/']
+
+ncols_hourly_file = 19
 
 year = int(input("Enter year \n"))
 nm0 = int(input("Enter month \n"))
@@ -91,10 +93,10 @@ ncols = []
 nheader = []
 scale = []
 
-varname_ = 'nut'
+varname_ = ['nuvisc','nuviscziriv']
+nvars = len(varname_)
 varnames_list = ['temp','o2','ch4','co2','uspeed','tke','nut', \
-                 'sprod','bprod','eps','tketrans']
-nvarlist = varnames_list.index(varname_)
+                 'nuvisc','nuviscziriv','sprod','bprod','eps','tketrans']
 
 #varname = 'Absolute value of velocity'
 #ylabel  = 'Depth, m' #$\circ C$
@@ -105,6 +107,8 @@ nvarlist = varnames_list.index(varname_)
 #nheader = 14 #Number of lines in the header
 
 logscale = False
+for k in range(nvars):
+    nvarlist = varnames_list.index(varname_[k])
     if (nvarlist == varnames_list.index('temp')):
             #varname.append(u'Temperature, model')
             #ylabel  = u'Depth, m' #$\circ C$
@@ -113,11 +117,12 @@ if (nvarlist == varnames_list.index('temp')):
             xlabel  = varname[len(varname)-1] + u',~$^{\circ}$C'
             basefilename.append('Profiles')
             ncol.append(2) #The column number to display, 2 for temperature
-        ncols.append(17) #Total number of columns
-        nheader.append(17) #Number of lines in the header
+            ncols.append(ncols_hourly_file) #Total number of columns
+            nheader.append(ncols_hourly_file) #Number of lines in the header
             scale.append(1.)
             xmin = 6.
             xmax = 26.
+    
     if (nvarlist == varnames_list.index('o2')):
             #varname.append(u'O$_2$, model')
             #ylabel  = u'Depth, m' #$\circ C$
@@ -127,11 +132,12 @@ if (nvarlist == varnames_list.index('o2')):
             xlabel  = varname[len(varname)-1] + u',~мг/л'
             basefilename.append('Profiles')
             ncol.append(15) #The column number to display, 2 for temperature
-        ncols.append(17) #Total number of columns
-        nheader.append(17) #Number of lines in the header
+            ncols.append(ncols_hourly_file) #Total number of columns
+            nheader.append(ncols_hourly_file) #Number of lines in the header
             scale.append(1.)
             xmin = 0.
             xmax = 15.
+    
     if (nvarlist == varnames_list.index('ch4')):
             #varname.append(u'CH$_4$, model')
             #ylabel  = u'Depth, m' #$\circ C$
@@ -141,12 +147,13 @@ if (nvarlist == varnames_list.index('ch4')):
             xlabel  = varname[len(varname)-1] + u',~ мкмоль/л'
             basefilename.append('Profiles')
             ncol.append(17) #The column number to display, 2 for temperature
-        ncols.append(17) #Total number of columns
-        nheader.append(17) #Number of lines in the header
+            ncols.append(ncols_hourly_file) #Total number of columns
+            nheader.append(ncols_hourly_file) #Number of lines in the header
             scale.append(1.)
             logscale = True
             xmin = 8.E-2
             xmax = 120.
+    
     if (nvarlist == varnames_list.index('co2')):
             #varname.append(u'CO$_2$, model')
             #ylabel  = u'Depth, m' #$\circ C$
@@ -156,19 +163,20 @@ if (nvarlist == varnames_list.index('co2')):
             xlabel  = varname[len(varname)-1] + u',~мкмоль/л'
             basefilename.append('Profiles')
             ncol.append(16) #The column number to display, 2 for temperature
-        ncols.append(17) #Total number of columns
-        nheader.append(17) #Number of lines in the header
+            ncols.append(ncols_hourly_file) #Total number of columns
+            nheader.append(ncols_hourly_file) #Number of lines in the header
             scale.append(1.)
             #xmin = -8.E-2
             #xmax = 8.E-2
+    
     if (nvarlist == varnames_list.index('uspeed')):
             varname.append(u'Компонента скорости по оси $x$')
             ylabel  = u'Глубина, м' #$\circ C$
             xlabel  = varname[len(varname)-1] + u',~м/с'
             basefilename.append('Profiles')
             ncol.append(12) #The column number to display, 2 for temperature
-        ncols.append(17) #Total number of columns
-        nheader.append(17) #Number of lines in the header
+            ncols.append(ncols_hourly_file) #Total number of columns
+            nheader.append(ncols_hourly_file) #Number of lines in the header
             scale.append(1.)
             #xmin = -8.E-2
             #xmax = 8.E-2
@@ -179,8 +187,8 @@ if (nvarlist == varnames_list.index('tke')):
             xlabel  = varname[len(varname)-1] + u',~м$^2$/с$^2$'
             basefilename.append('Profiles')
             ncol.append(5) #The column number to display, 2 for temperature
-        ncols.append(17) #Total number of columns
-        nheader.append(17) #Number of lines in the header
+            ncols.append(ncols_hourly_file) #Total number of columns
+            nheader.append(ncols_hourly_file) #Number of lines in the header
             scale.append(1.)
     
     if (nvarlist == varnames_list.index('nut')):
@@ -189,8 +197,28 @@ if (nvarlist == varnames_list.index('nut')):
             xlabel  = varname[len(varname)-1] + u',~м$^2$/с'
             basefilename.append('Profiles')
             ncol.append(7) #The column number to display, 2 for temperature
-        ncols.append(17) #Total number of columns
-        nheader.append(17) #Number of lines in the header
+            ncols.append(ncols_hourly_file) #Total number of columns
+            nheader.append(ncols_hourly_file) #Number of lines in the header
+            scale.append(1.)
+    
+    if (nvarlist == varnames_list.index('nuvisc')):
+            varname.append(u'Коэффициент вязкости по модели $k-\epsilon$')
+            ylabel  = u'Глубина, м' #$\circ C$
+            xlabel  = varname[len(varname)-1] + u',~м$^2$/с'
+            basefilename.append('Profiles')
+            ncol.append(18) #The column number to display, 2 for temperature
+            ncols.append(ncols_hourly_file) #Total number of columns
+            nheader.append(ncols_hourly_file) #Number of lines in the header
+            scale.append(1.)
+    
+    if (nvarlist == varnames_list.index('nuviscziriv')):
+            varname.append(u'Коэффициент вязкости по модели Zi')
+            ylabel  = u'Глубина, м' #$\circ C$
+            xlabel  = varname[len(varname)-1] + u',~м$^2$/с'
+            basefilename.append('Profiles')
+            ncol.append(19) #The column number to display, 2 for temperature
+            ncols.append(ncols_hourly_file) #Total number of columns
+            nheader.append(ncols_hourly_file) #Number of lines in the header
             scale.append(1.)
     
     if (nvarlist == varnames_list.index('sprod')):
@@ -239,12 +267,11 @@ if (nvarlist == varnames_list.index('tketrans')):
 #xlabel = u'Температура, $^{\circ}C$' #u'Слагаемые уравнения баланса ТКЭ, м$^2$/с$^3$'
 #xlabel = u'Methane, mcmol/l$' 
 
-nvars = len(varname)
-
 #fileout = 'uspeed'
 #fileout = 'TKE_terms'
 #fileout = 'temp'
-fileout = 'ch4'
+#fileout = 'ch4'
+fileout = 'viscosities'
 
 #Files names suffixes (timestamps)
 nfiles = 0 #number of files
@@ -318,7 +345,7 @@ for np in range(npaths):
 			legenditems.append(u'модель') 
 			for m in range(len(vals)) : vals[m] = vals[m]*scale[k]
 			plt.plot(vals,depths,ls=linestyles[np], \
-                        color=linecolors[varnames_list.index(varname_)],linewidth=3,\
+                        color=linecolors[varnames_list.index(varname_[k])],linewidth=3,\
                         markersize=markersize)
 	#print(depths)
 	#print(vals)