Numeric-analytical solution for speed in Cuette-Poiseuille flow added

setup/cuette_setup.dat

@@ -146,7 +146,7 @@ ifrad    0
 ifbubble 1
 sedim    0
 dyn_pgrad 6
-pgrad     0.
+pgrad     0.1
 nManning  2.5E-2
 horvisc  0.
 backdiff 0

source/model/init.f90

@@ -139,6 +139,7 @@
  if (dyn_pgrad%par == 5 .or. dyn_pgrad%par == 6) then 
    allocate (eps_ziriv     (1:M  ))
    allocate (viscturb_ziriv(1:M  ))
+   allocate (speed_ziriv   (1:M+1))
  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)) 

source/model/lake.f90

@@ -2490,7 +2490,7 @@ h2_out = h2
 !Diagnoctics
 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,viscturb_ziriv)
+  & g*tan(pi*alphax/180.),z_half,z_full,eps_ziriv,viscturb_ziriv,speed_ziriv)
 endif
 
 
@@ -2885,30 +2885,30 @@ if (runmode%par == 1)  then
     if (year<1000) i = year+1000
     write (*,'(3i8,2f8.2)') i, month, day, hour, time/(month_sec)
 
-    write(*,'(a)')
-    write(*,*) 'Methane generation, transport and oxidation rates'
-    write(*,*) 'Total bottom diff, Total surface diff, Total bottom ebul, Total surface ebul'
-    write(*,'(a,4(2x,e12.5))') 'Summer ', fdifftot(1), fdiff_lake_surftot(1), febulbottot(1), febultot(1)
-    write(*,'(a,4(2x,e12.5))') 'Winter ', fdifftot(2), fdiff_lake_surftot(2), febulbottot(2), febultot(2)
-    write(*,*) 'Total prod from "new" org, Total prod from "old" org, Total oxidation in water, &
-    & Total oxidation in ice (for winter)'
-    write(*,'(a,3(2x,e12.5))') 'Summer ', &
-    & rprod_total_newC_integr(1), rprod_total_oldC_integr(1), methoxidwat(1)
-    write(*,'(a,4(2x,e12.5))') 'Winter ', &
-    & rprod_total_newC_integr(2), rprod_total_oldC_integr(2), methoxidwat(2), ice_meth_oxid
-    write(*,'(a,(2x,e12.5))') 'Change of integral methane amount in water column and ice cover', &
-    & totmethc - totmeth0 + tot_ice_meth_bubbles*mf
-    write(*,'(a)')
-    write(*,*) 'Checking methane balance residuals:'   
-    write(*,'(a,1(2x,e12.5))') 'Total production in soil - total bottom flux ', &
-    & sum(rprod_total_newC_integr(1:2)) + sum(rprod_total_oldC_integr(1:2)) - &
-    & sum(febulbottot(1:2)) - sum(fdifftot(1:2))
-    write(*,'(a,1(2x,e12.5))') 'Total bottom flux - total surface flux - total oxidation - & 
-    & methane total column amount change', &
-    & sum(febulbottot(1:2)) + sum(fdifftot(1:2)) - &
-    & sum(febultot(1:2)) - sum(fdiff_lake_surftot(1:2)) - &
-    & sum(methoxidwat(1:2)) - ice_meth_oxid - &
-    & (totmethc - totmeth0 + tot_ice_meth_bubbles*mf)
+    !write(*,'(a)')
+    !write(*,*) 'Methane generation, transport and oxidation rates'
+    !write(*,*) 'Total bottom diff, Total surface diff, Total bottom ebul, Total surface ebul'
+    !write(*,'(a,4(2x,e12.5))') 'Summer ', fdifftot(1), fdiff_lake_surftot(1), febulbottot(1), febultot(1)
+    !write(*,'(a,4(2x,e12.5))') 'Winter ', fdifftot(2), fdiff_lake_surftot(2), febulbottot(2), febultot(2)
+    !write(*,*) 'Total prod from "new" org, Total prod from "old" org, Total oxidation in water, &
+    !& Total oxidation in ice (for winter)'
+    !write(*,'(a,3(2x,e12.5))') 'Summer ', &
+    !& rprod_total_newC_integr(1), rprod_total_oldC_integr(1), methoxidwat(1)
+    !write(*,'(a,4(2x,e12.5))') 'Winter ', &
+    !& rprod_total_newC_integr(2), rprod_total_oldC_integr(2), methoxidwat(2), ice_meth_oxid
+    !write(*,'(a,(2x,e12.5))') 'Change of integral methane amount in water column and ice cover', &
+    !& totmethc - totmeth0 + tot_ice_meth_bubbles*mf
+    !write(*,'(a)')
+    !write(*,*) 'Checking methane balance residuals:'   
+    !write(*,'(a,1(2x,e12.5))') 'Total production in soil - total bottom flux ', &
+    !& sum(rprod_total_newC_integr(1:2)) + sum(rprod_total_oldC_integr(1:2)) - &
+    !& sum(febulbottot(1:2)) - sum(fdifftot(1:2))
+    !write(*,'(a,1(2x,e12.5))') 'Total bottom flux - total surface flux - total oxidation - & 
+    !& methane total column amount change', &
+    !& sum(febulbottot(1:2)) + sum(fdifftot(1:2)) - &
+    !& sum(febultot(1:2)) - sum(fdiff_lake_surftot(1:2)) - &
+    !& sum(methoxidwat(1:2)) - ice_meth_oxid - &
+    !& (totmethc - totmeth0 + tot_ice_meth_bubbles*mf)
 
 
     open(12345,file=outpath(1:len_trim(outpath))//'lastscrout.dat')

source/model/lake_modules.f90

@@ -567,7 +567,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(:), viscturb_ziriv(:)
+real(kind=ireals), allocatable :: eps_ziriv(:), viscturb_ziriv(:), speed_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 = 29
+integer(kind=iintegers), parameter :: n_var = 30
 integer(kind=iintegers), parameter :: n_var_scalar = 10
 integer(kind=iintegers), allocatable :: hour_old(:,:)
 integer(kind=iintegers), allocatable :: tsteps(:,:)
@@ -508,6 +508,7 @@ var(26,ix,iy,1:M+1)  = qwater(1:M+1,1)*molm3tomcM
 var(27,ix,iy,1:M)   = k2             (1:M  )
 var(28,ix,iy,1:M)   = viscturb_ziriv (1:M  )
 var(29,ix,iy,1:M)   = eps_ziriv      (1:M  )
+var(30,ix,iy,1:M+1) = speed_ziriv    (1:M+1)
 
 var_scalar(1,ix,iy) = S_integr_positive
 var_scalar(2,ix,iy) = S_integr_negative
@@ -556,8 +557,9 @@ if (int(hour)/=hour_old(ix,iy)) then
   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'
   write (out_unit,*) '20 - dissipation rate from Ziriv model, m**2/s**3'
+  write (out_unit,*) '21 - longitudinal speed from Ziriv model, m/s'
   do i=1,M 
-    write (out_unit,'(f14.6,f12.5,e10.3,f10.3,3e10.3,4e16.7,3f11.5,3e16.7,3e10.3)')      &
+    write (out_unit,'(f14.6,f12.5,e10.3,f10.3,3e10.3,4e16.7,3f11.5,3e16.7,3e10.3,f11.5)')      &
     & -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), &
@@ -579,7 +581,8 @@ if (int(hour)/=hour_old(ix,iy)) then
     & accum_var(26,ix,iy,i),                                              &
     & accum_var(27,ix,iy,i),                                              &
     & accum_var(28,ix,iy,i),                                              &
-    & accum_var(29,ix,iy,i)
+    & accum_var(29,ix,iy,i),                                              &
+    & accum_var(30,ix,iy,i)
   enddo 
   close (out_unit)
   

source/model/turb_mod.f90

@@ -1838,7 +1838,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,viscturb_ziriv)
+SUBROUTINE ZIRIVMODEL(M,h1,u_star_surf,gradp,z_half,z_full,eps_ziriv,viscturb_ziriv,speed_ziriv)
 use PHYS_CONSTANTS, only : kappa
 use NUMERIC_PARAMS, only : pi
 implicit none
@@ -1847,31 +1847,43 @@ integer(kind=iintegers), intent(in) :: M      !> Number of layers
 real(kind=ireals), intent(in ) :: h1          !> Riverflow depth, m
 real(kind=ireals), intent(in ) :: u_star_surf !> Friction velocity in water at the surface, m/s
 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(in ) :: z_half(1:M) !> A set of depths where dissipation and viscosity 
+                                              !! are to be evaluated
+real(kind=ireals), intent(in ) :: z_full(1:M+1) !> A set of depths where speed 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
+real(kind=ireals), intent(out) :: speed_ziriv(1:M+1) !> Longitudinal speed, m/s
 !Local variables
 integer(kind=iintegers) :: i !Loop index
 real(kind=ireals) :: u_star, z_Obukhov, h_max
 
-if (gradp*u_star_surf > 0.) then 
+if (gradp*u_star_surf >= 0.) then 
   !Wind forcing along the pressure gradient, 
   !monotonic decrease of flow speed with depth
-  do i = 1, M
+  speed_ziriv(M+1) = 0.
+  do i = M, 1, -1
     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
+    speed_ziriv(i) = speed_ziriv(i+1) + &
+    & u_star/(kappa*z_Obukhov)*(z_full(i+1) - z_full(i))
   enddo
+  speed_ziriv(:) = speed_ziriv(:)*sign(1._ireals,u_star_surf)
 else
   !Wind forcing opposite to the pressure gradient, 
-  !a maxumum of flow speed between top and bottom boundary
-  do i = 1, M
+  !a maxumum of flow speed locates between top and bottom boundary
+  h_max = u_star_surf**2/abs(gradp)
+  do i = M, 1, -1
     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
+    speed_ziriv(i) = speed_ziriv(i+1) + &
+    & u_star/(kappa*z_Obukhov)*(z_full(i+1) - z_full(i)) * &
+    & sign(1._ireals,z_half(i)-h_max)
   enddo
+  speed_ziriv(:) = speed_ziriv(:)*sign(1._ireals, - u_star_surf)
 endif
 
 END SUBROUTINE ZIRIVMODEL

tools/plot_hourly.py

@@ -30,7 +30,7 @@ plt.rc('text.latex',preamble='\usepackage[russian]{babel}')
 #path = ['../results/Mojai2016-2017_damw5/hourly/']#,'../results/Mojai2016-2017_damtribheat=0/hourly/']
 path = ['../results/cuette/hourly/']
 
-ncols_hourly_file = 20
+ncols_hourly_file = 21
 
 year = int(input("Enter year \n"))
 nm0 = int(input("Enter month \n"))
@@ -93,12 +93,13 @@ ncols = []
 nheader = []
 scale = []
 
-varname_ = ['nuvisc','nuviscziriv']
+#varname_ = ['nuvisc','nuviscziriv']
 #varname_ = ['eps','eps_ziriv']
+varname_ = ['uspeed','speed_ziriv']
 nvars = len(varname_)
 varnames_list = ['temp','o2','ch4','co2','uspeed','tke','nut', \
                  'nuvisc','nuviscziriv','sprod','bprod','eps', \
-                 'tketrans','eps_ziriv']
+                 'tketrans','eps_ziriv','speed_ziriv']
 
 #varname = 'Absolute value of velocity'
 #ylabel  = 'Depth, m' #$\circ C$
@@ -288,6 +289,17 @@ for k in range(nvars):
             nheader.append(ncols_hourly_file)
             scale.append(1.)
 
+    if (nvarlist == varnames_list.index('speed_ziriv')):
+            varname.append(u'Продольная скорость (модель Ziriv)')
+            #varname.append(u'$\epsilon$')
+            ylabel  = u'Глубина, м' #$\circ C$
+            xlabel  = varname[len(varname)-1] + u', м/с'
+            basefilename.append('Profiles')
+            ncol.append(21)
+            ncols.append(ncols_hourly_file)
+            nheader.append(ncols_hourly_file)
+            scale.append(1.)
+
 #xlabel = u'Температура, $^{\circ}C$' #u'Слагаемые уравнения баланса ТКЭ, м$^2$/с$^3$'
 #xlabel = u'Methane, mcmol/l$'