1. Bugs fixed in surface momentum flux usage in k-eps closure 2. Cuette flow...

1. Bugs fixed in surface momentum flux usage in k-eps closure 2. Cuette flow b.c.s corrected 3. Output corrected

setup/cuette_driver.dat

@@ -94,10 +94,10 @@ month0      5
 day0        1
 hour0       0.
 #
-tinteg      35
+tinteg      15
 spinup_times  0
 spinup_period 0
-dt          10.
+dt          2.
 call_Flake  0
 cp_period   0.
 control_point 0
@@ -139,8 +139,8 @@ alphay      0.0
 a_veg       1.
 c_veg       1.e-3
 h_veg       0.
-#kor         0.
-kor         -999.
+kor         0.
+#kor         -999.
 # 1.e-4
 phi         60.
 lam         24.25
@@ -171,7 +171,7 @@ fetch       410.
 #-----------------------------------------------------------------------------------------
 #
 l10         0.
-h10         10.
+h10         2.
 #h10         3.
 #select_h10  3
 #1  1.2

setup/cuette_setup.dat

@@ -122,7 +122,7 @@ d_bot  1.E-2
 varalb   1
 PBLpar   0
 waveenh  0
-momflxpart 1
+momflxpart 0
 c_d      -999
 kwe      100.
 relwind  0
@@ -130,7 +130,7 @@ eos      1
 lindens  0
 nmeltpoint 1
 Turbpar  2
-stabfunc 2
+stabfunc 1
 kepsbc   1
 soiltype 5
 soil_depth  10.
@@ -167,15 +167,10 @@ deadvol 0.
 #                           INITIAL CONDITIONS FOR TEMPERATURE
 #----------------------------------------------------------------------------------------
 #
-T_profile 8
+T_profile 3
 0.2   3.410 0. 1.03E-3  1.76E-1 2.77E-1 0. 0. 
 1.    3.520 0. 1.18E-3  1.95E-1 2.81E-1 0. 0.
-3.    3.540 0. 1.41E-3  1.80E-1 2.77E-1 0. 0.
-5.    3.660 0. 1.72E-3  1.87E-1 2.75E-1 0. 0.
-7.    3.540 0. 2.07E-3  1.87E-1 2.79E-1 0. 0.
-9.    3.480 0. 1.91E-3  1.89E-1 2.73E-1 0. 0.
-11.   3.490 0. 7.32E-3  1.99E-1 2.66E-1 0. 0.
-12.   3.970 0. 14.88E-3 2.04E-1 2.63E-1 0. 0.
+2.    3.540 0. 1.41E-3  1.80E-1 2.77E-1 0. 0.
 T_soilprofile 2
 0. 4.
 10. 4.
@@ -202,171 +197,6 @@ tribheat       0
 N_triblev      -1
 N_tribin       -1
 iefflloc       -1
-#inflowprof
-#1  100. 0. 8. 
-#2  100. 0. 8.
-#3  100. 0. 8.
-#3  100. 0. 8.
-#5  100. 0. 8.
-#6  100. 0. 8.
-#7  100. 0. 8.
-#8  100. 0. 8.
-#9  100. 0. 8.
-#10 100. 0. 8.
-#11 100. 0. 8.
-#12 100. 0. 8.
-#13 100. 0. 8.
-#14 100. 0. 8.
-#15 100. 0. 8.
-#16 100. 0. 8.
-#17 100. 0. 8.
-#18 100. 0. 8.
-#19 100. 0. 8.
-#20 100. 0. 8.
-#21 100. 0. 8.
-#22 100. 0. 8.
-#23 100. 0. 8.
-#24 100. 0. 8.
-#25 100. 0. 8.
-#26 100. 0. 8.
-#27 100. 0. 8.
-#28 100. 0. 8.
-#29 100. 0. 8.
-#30 100. 0. 8.
-#31 100. 0. 8.
-#32 100. 0. 8.
-#33 100. 0. 8.
-#34 100. 0. 8.
-#35 100. 0. 8.
-#36 100. 0. 8.
-#37 100. 0. 8.
-#38 100. 0. 8.
-#39 100. 0. 8.
-#40 100. 0. 8.
-#41 100. 0. 8.
-#42 100. 0. 8. 
-#43 100. 0. 8.
-#44 100. 0. 8.
-#45 100. 0. 8.
-#46 100. 0. 8.
-#47 100. 0. 8.
-#48 100. 0. 8.
-#49 100. 0. 8.
-#50 100. 0. 8.
-#51 100. 0. 8.
-#52 100. 0. 8.
-#53 100. 0. 8.
-#54 100. 0. 8.
-#55 100. 0. 8.
-#56 100. 0. 8.
-#57 100. 0. 8.
-#58 100. 0. 8.
-#59 100. 0. 8.
-#60 100. 0. 8.
-#61 100. 0. 8.
-#62 100. 0. 8.
-#63 100. 0. 8.
-#64 100. 0. 8.
-#65 100. 0. 8.
-#66 100. 0. 8.
-#67 100. 0. 8.
-#68 100. 0. 8.
-#69 100. 0. 8.
-#70 100. 0. 8.
-#71 100. 0. 8.
-#72 100. 0. 8.
-#73 100. 0. 8.
-#74 100. 0. 8.
-#75 100. 0. 8.
-#76 100. 0. 8.
-#77 100. 0. 8.
-#78 100. 0. 8.
-#79 100. 0. 8.
-#80 100. 0. 8.
-#81 100. 0. 8.
-#
-#outflowprof
-#1  100. 0. 
-#2  100. 0.
-#3  100. 0.
-#4  100. 0.
-#5  100. 0.
-#6  100. 0.
-#7  100. 0.
-#8  100. 0.
-#9  100. 0.
-#10 100. 0.
-#11 100. 0.
-#12 100. 0.
-#13 100. 0.
-#14 100. 0.
-#15 100. 0.
-#16 100. 0.
-#17 100. 0.
-#18 100. 0.
-#19 100. 0.
-#20 100. 0.
-#21 100. 0.
-#22 100. 0.
-#23 100. 0.
-#24 100. 0.
-#25 100. 0.
-#26 100. 0.
-#27 100. 0.
-#28 100. 0.
-#29 100. 0.
-#30 100. 0.
-#31 100. 0.
-#32 100. 0.
-#33 100. 0.
-#34 100. 0.
-#35 100. 0.
-#36 100. 0.
-#37 100. 0.
-#38 100. 0.
-#39 100. 0.
-#40 100. 0.
-#41 100. 0.
-#42 100. 0. 
-#43 100. 0.
-#44 100. 0.
-#45 100. 0.
-#46 100. 0.
-#47 100. 0.
-#48 100. 0.
-#49 100. 0.
-#50 100. 0.
-#51 100. 0.
-#52 100. 0.
-#53 100. 0.
-#54 100. 0.
-#55 100. 0.
-#56 100. 0.
-#57 100. 0.
-#58 100. 0.
-#59 100. 0.
-#60 100. 0.
-#61 100. 0.
-#62 100. 0.
-#63 100. 0.
-#64 100. 0.
-#65 100. 0.
-#66 100. 0.
-#67 100. 0.
-#68 100. 0.
-#69 100. 0.
-#70 100. 0.
-#71 100. 0.
-#72 100. 0.
-#73 100. 0.
-#74 100. 0.
-#75 100. 0.
-#76 100. 0.
-#77 100. 0.
-#78 100. 0.
-#79 100. 0.
-#80 100. 0.
-#81 100. 0.
 #
 #----------------------------------------------------------------------------------------
 #            DATA ASSIMILATION CONTROLS (NOT OPERATIONAL: PUT EVERYTHING TO 0)
@@ -422,19 +252,18 @@ zserout -999.
 rtemp 1
 -999. -999. -999.
 #
-ngrid_out 12
-0.5
+ngrid_out 11
+0.
+0.2
+0.4
+0.6
+0.8
 1.
+1.2
+1.4
+1.6
+1.8
 2.
-3.
-4.
-5.
-6.
-7.
-8.
-9.
-9.5
-10.
 #
 ngridsoil_out 11
 0.

source/model/init_var_mod.f90

@@ -15,7 +15,7 @@
   & hx1, hx2, hy1, hy2, &
   & hx1t, hx2t, hy1t, hy2t, &
   & hx1ml, hx2ml, hy1ml, hy2ml, &
-  & veg, snmelt, snowmass, cdmw2, velfrict_prev, &
+  & veg, snmelt, snowmass, cdmw, velfrict_prev, &
   & roughness, eflux0_kinem, Elatent, &
   & totalevap, totalmelt, totalprecip, totalwat, totalpen, &
   & time, dhwfsoil, dhw, dhw0, dhi, dhi0, dls0, &
@@ -89,7 +89,7 @@
   real(kind=ireals), intent(inout), allocatable :: hx1ml(:), hx2ml(:), hy1ml(:), hy2ml(:)
   real(kind=ireals), intent(out) :: veg
   real(kind=ireals), intent(out) :: snmelt, snowmass
-  real(kind=ireals), intent(out) :: cdmw2, velfrict_prev
+  real(kind=ireals), intent(out) :: cdmw, velfrict_prev
   real(kind=ireals), intent(out) :: roughness
   real(kind=ireals), intent(out) :: eflux0_kinem, Elatent
   real(kind=ireals), intent(out) :: totalevap, totalmelt, totalprecip, totalwat, totalpen
@@ -459,7 +459,7 @@
  
   snmelt = 0.e0_ireals 
   snowmass = 0.e0_ireals 
-  cdmw2 = 1.d-15
+  cdmw = 1.d-15
   velfrict_prev = 1.d-2 
   roughness = 1.d-3
   eflux0_kinem = 0.e0_ireals

source/model/lake.f90

@@ -655,7 +655,7 @@ if (init(ix,iy) == 0_iintegers) then
   & hx1, hx2, hy1, hy2, &
   & hx1t, hx2t, hy1t, hy2t, &
   & hx1ml, hx2ml, hy1ml, hy2ml, &
-  & veg, snmelt, snowmass, cdmw2, velfrict_prev, &
+  & veg, snmelt, snowmass, cdmw, velfrict_prev, &
   & roughness, eflux0_kinem, Elatent, &
   & totalevap, totalmelt, totalprecip, totalwat, totalpen, &
   & time, dhwfsoil, dhw, dhw0, dhi, dhi0, dls0, &
@@ -713,7 +713,7 @@ if (init(ix,iy) /= 0_iintegers .or. &
   & oxyg(1,1), oxygsoil, &
   & DIC(1,1), DOC, POCL, POCD, phosph, &
   & snmelt, snowmass, &
-  & cdmw2, &
+  & cdmw, &
   & time, &
   & dhwfsoil, &
   & Elatent, &
@@ -933,8 +933,6 @@ endif
 ! Calculation of the whole temperature profile
 call T_SOLVER(ix,iy,dnx,dny,year,month,day,hour,phi, &
 & RadWater, RadIce, SurfTemp, fetch, dt)
-!print*, Tw1(1), Tw1(M+1)
-!print*, lamw(1), lamw(M)
 
 ! Diagnostic calculations
 do i = 1, M

source/model/lake_modules.f90

@@ -43,7 +43,7 @@ real(kind=ireals), parameter :: min_ice_thick = 0.02
 real(kind=ireals), parameter :: min_water_thick = 0.02
 real(kind=ireals), parameter :: T_phase_threshold = 1.d-5
 
-real(kind=ireals), parameter :: minwind = 1.d-2
+real(kind=ireals), parameter :: minwind = 1.d-5
 
  !   ML --- total number of layers
  !   MS --- maximal number of layers in snow
@@ -276,7 +276,6 @@ real(kind=ireals) :: surfrad
 real(kind=ireals) :: hbal
 real(kind=ireals) :: zref
 real(kind=ireals) :: botflux
-real(kind=ireals) :: cdmw2
 real(kind=ireals) :: shortwave, sabspen
 real(kind=ireals) :: precip
 real(kind=ireals) :: Sflux0

source/model/momentum_mod.f90

@@ -17,7 +17,7 @@ SUBROUTINE MOMENTUM_SOLVER(ix, iy, nx, ny, ndatamax, year, month, day, hour, &
 ! for two horizontal components od speed
 
 use ATMOS, only : &
-& cdmw2, &
+& cdmw, &
 & zref, &
 & uwind, vwind, wind, wind10, windwork, &
 & u, v, urel, vrel, &
@@ -247,7 +247,6 @@ allocate (LxLyCDL(1:M+1,1:2)); LxLyCDL(1:M+1,1:2) = missing_value
 u = u1(1)
 v = v1(1)
 
-
 if (relwind%par == 1) then
   urel=uwind-u
   vrel=vwind-v
@@ -274,12 +273,9 @@ endif
 kor2  = 0.5d0*kor*dt 
 
 wr   = sqrt(urels**2+vrels**2)
-tau_air = roa0*cdmw2*wr
+tau_air = roa0*cdmw*wr
 ufr  = sqrt(tau_air/row0)    
 
-! Sensitivity test
-! tau_air = 1.d-15
-
 ! PARAMETERIZATION OF WIND STRESS SPLIT-UP INTO TWO PARTS:
 ! Momentum exchange coefficient for wind speed at 10 m, Wuest and Lorke (2003)
 wind10 = wr*log(10./roughness0)/log(zref/roughness0)
@@ -288,7 +284,7 @@ wind10 = wr*log(10./roughness0)/log(zref/roughness0)
 !else
 !  C10 = 0.0044*wind10**(-1.15)
 !endif
-C10 = cdmw2*wr/(wind10*wind10)
+C10 = cdmw*wr/(wind10*wind10)
 ! a.PARAMETERIZATION USING SIGNIFICANT WAVE HEIGHT (NB1, pp. 8)
 ! co1=4.*kws/sqrt(C10*1000.*g) !1000 m is "average fetch" of lake Syrdakh
 ! kw=min(max(co1*wind10,kws),0.75)
@@ -1221,7 +1217,7 @@ elseif (cuette%par == 1) then
 
 endif bctop
 
-bcbot : if (cuette%par == 0) then
+bcbot : if (cuette%par == 0 .or. cuette%par == 11) then
   ! Boundary conditions for u and v at the lake bottom
   
   k2(M) = area_half(M)/area_int(M+1)*k2(M)
@@ -1254,7 +1250,7 @@ bcbot : if (cuette%par == 0) then
   
   k2(M) = area_int(M+1)/area_half(M)*k2(M)
 
-elseif (cuette%par == 1 .or. cuette%par == 11) then
+elseif (cuette%par == 1) then
 
   ! Bottom boundary
   cm_(M+1,:,:) = 0.

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 = 28
+integer(kind=iintegers), parameter :: n_var = 29
 integer(kind=iintegers), parameter :: n_var_scalar = 10
 integer(kind=iintegers), allocatable :: hour_old(:,:)
 integer(kind=iintegers), allocatable :: tsteps(:,:)
@@ -507,6 +507,7 @@ 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(29,ix,iy,1:M)  = eps_ziriv      (1:M)
 
 var_scalar(1,ix,iy) = S_integr_positive
 var_scalar(2,ix,iy) = S_integr_negative
@@ -554,8 +555,9 @@ if (int(hour)/=hour_old(ix,iy)) then
   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'
+  write (out_unit,*) '20 - dissipation rate from Ziriv model, m**2/s**3'
   do i=1,M 
-    write (out_unit,'(f14.6,f12.5,e10.3,f10.3,3e10.3,4e14.7,3f11.5,3e14.7,2e10.3)')      &
+    write (out_unit,'(f14.6,f12.5,e10.3,f10.3,3e10.3,4e16.7,3f11.5,3e16.7,3e10.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), &
@@ -576,7 +578,8 @@ if (int(hour)/=hour_old(ix,iy)) then
     & accum_var(25,ix,iy,i),                                              &
     & accum_var(26,ix,iy,i),                                              &
     & accum_var(27,ix,iy,i),                                              &
-    & accum_var(28,ix,iy,i)
+    & accum_var(28,ix,iy,i),                                              &
+    & accum_var(29,ix,iy,i)
   enddo 
   close (out_unit)
   

source/model/t_solver_mod.f90

@@ -194,6 +194,7 @@ elseif (cuette%par == 1 .or. cuette%par == 11) then
   ! B.c.s for temperature for Cuette flow
   Tsurf3 = Tw1(1)
   call T_DIFF(1_iintegers,Tsurf3,dt,snow,ice,water,deepice,nsoilcols,.false.)
+  Bal3 = HEATBALANCE(Tsurf3,surftyp(ix,iy),RadWater,RadIce,fetch,dt)
 
 endif Tbc_if
 

source/model/turb_mod.f90

@@ -734,7 +734,6 @@ allocate (rhotemp(1:M), rhosal(1:M))
     ! where the b.c. for epsilon are formulated
     
     ! dist_surf = 0._ireals !0.5d0*ddz(1)*h1
-
     if_dissbc: if (kepsbc%par <= 4) then
       !Neumann boundary conditions for TKE dissipation
       !Boundary condition at the top
@@ -995,8 +994,9 @@ allocate (rhotemp(1:M), rhosal(1:M))
   endif
 
   do i = 1, M
-    KT(i) = lam_T*max(CEt(i)*E2(i)**2/(eps2(i) + ACCk),min_diff)*cw_m_row0 !E2, eps2
+    KT(i) = max(CEt(i)*E2(i)**2/(eps2(i) + ACCk),min_diff)*cw_m_row0 !E2, eps2
   enddo
+  !print*, maxval(KT)/cw_m_row0, maxval(k2), maxval(CEt), maxval(CE)
 
   ! Diagnostic calculation
   !Turbulent transport of TKE
@@ -1106,9 +1106,6 @@ allocate (rhotemp(1:M), rhosal(1:M))
   endif
   ! End debug output
   
-  ! print*, 'S+ = ', S_integr_positive, 'S- = ', S_integr_negative
-  ! print*, 'eps_integr = ', eps_integr
-  
   ! do i = 2, M-1
   !   work(i,1) = dhw_keps/(dt*h1)*dzeta_05int(i)*(E2(i+1)-E2(i-1))/ &
   !   & (0.5d0*ddz(i-1)+ddz(i)+0.5d0*ddz(i+1))/Buoyancy0
@@ -1118,12 +1115,6 @@ allocate (rhotemp(1:M), rhosal(1:M))
 
   ! Output
   if (turb_out%par == 1) then
-    ! do i=1,M+1
-    !   write (2115,7) time/hour_sec+12.-3*24, dzeta_05int(i)*h1, &
-    !   & E2(i)*1.E9,eps2(i)*1.E9, KT(i),tau_air*1.E3,Tw1(i)
-    ! enddo
-    ! write (2116,8) time/hour_sec+12.-3*24, E2(1)*1.E9, &
-    ! & eps2(1)*1.E9, KT(1)
     if ((dE_it > 1.d-7 .or. deps_it > 1.d-9).and.(.not.semi_implicit_scheme)) then
       write (nunit_, *) nstep, nl, dE_it, deps_it
     endif
@@ -1139,12 +1130,11 @@ allocate (rhotemp(1:M), rhosal(1:M))
     E1(i) = E2(i)
     eps1(i) = eps2(i)
   enddo
-  !print*, E2(1), E2(M), eps2(1), eps2(M)
-  !print*, KT(1), KT(M)
-  !print*, tau_air, tau_gr
  
-  ! write(*,*) 'K_EPSILON finished', nstep  
-  ! STOP 
+  !Recalculation for output
+  do i = 1, M
+    k2(i) = k2(i)*CE(i)
+  enddo
      
 firstcall = .false.
 

tools/plot_contour.py

@@ -93,7 +93,7 @@ nmod_max = 12 #The last item in the list above, from the modeled data
 #Constants
 day_sec = 24.*60.*60.
 
-path = '../results/Mojai2016-2017/time_series/' # path to files
+path = '../results/cuette/time_series/' # path to files
 #pathobs = os.path.expanduser('/media/victor/main/Files/main/observdata/Iseo/')
 pathobs = os.path.expanduser('/home/victor/Files/main/observdata/Iseo/')
 #pathobs = os.path.expanduser('/media/victor/main/Files/main/observdata/Kuivajarvi/Kuiva2014sept_campaign/')

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 = 19
+ncols_hourly_file = 20
 
 year = int(input("Enter year \n"))
 nm0 = int(input("Enter month \n"))
@@ -40,7 +40,7 @@ nh0 = int(input("Enter first hour \n"))
 nhlast = int(input("Enter last hour \n"))
 #figlabel = input("Enter figure label (put 0 for none) \n")
 
-figlabel = u'г)'
+figlabel = u'' #u'г)'
 
 if nd0 == ndlast and nh0 == nhlast:
         nhstep = 0
@@ -94,9 +94,11 @@ nheader = []
 scale = []
 
 varname_ = ['nuvisc','nuviscziriv']
+#varname_ = ['eps','eps_ziriv']
 nvars = len(varname_)
 varnames_list = ['temp','o2','ch4','co2','uspeed','tke','nut', \
-                 'nuvisc','nuviscziriv','sprod','bprod','eps','tketrans']
+                 'nuvisc','nuviscziriv','sprod','bprod','eps', \
+                 'tketrans','eps_ziriv']
 
 #varname = 'Absolute value of velocity'
 #ylabel  = 'Depth, m' #$\circ C$
@@ -243,19 +245,30 @@ for k in range(nvars):
             nheader.append(7)
             scale.append(1.)
             
+#    if (nvarlist == varnames_list.index('eps')):
+#            varname.append(u'Диссипация ТКЭ (модель $k-\epsilon$)')
+#            #varname.append(u'$\epsilon$')
+#            ylabel  = u'Глубина, м' #$\circ C$
+#            xlabel  = varname[len(varname)-1] + u', м$^2$/с$^3$'
+#            basefilename.append('TKE_budget')
+#            ncol.append(5)
+#            ncols.append(6)
+#            nheader.append(7)
+#            scale.append(1.)
+
     if (nvarlist == varnames_list.index('eps')):
-            varname.append(u'Диссипация ТКЭ')
+            varname.append(u'Диссипация ТКЭ (модель $k-\epsilon$)')
             #varname.append(u'$\epsilon$')
             ylabel  = u'Глубина, м' #$\circ C$
             xlabel  = varname[len(varname)-1] + u', м$^2$/с$^3$'
-            basefilename.append('TKE_budget')
-            ncol.append(5)
-            ncols.append(6)
-            nheader.append(7)
-            scale.append(-1.)
+            basefilename.append('Profiles')
+            ncol.append(6)
+            ncols.append(ncols_hourly_file)
+            nheader.append(ncols_hourly_file)
+            scale.append(1.)
            
     if (nvarlist == varnames_list.index('tketrans')):
-            varname.append(u'Перенос ТКЭ')
+            varname.append(u'Перенос ТКЭ (модель $k-\epsilon$)')
             ylabel  = u'Глубина, м' #$\circ C$
             xlabel  = varname[len(varname)-1] + u', м$^2$/с$^3$'
             basefilename.append('TKE_budget')
@@ -264,6 +277,17 @@ for k in range(nvars):
             nheader.append(7)
             scale.append(1.)
 
+    if (nvarlist == varnames_list.index('eps_ziriv')):
+            varname.append(u'Диссипация ТКЭ (модель Ziriv)')
+            #varname.append(u'$\epsilon$')
+            ylabel  = u'Глубина, м' #$\circ C$
+            xlabel  = varname[len(varname)-1] + u', м$^2$/с$^3$'
+            basefilename.append('Profiles')
+            ncol.append(20)
+            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$' 
 
@@ -342,7 +366,7 @@ for np in range(npaths):
 				vals[j] = float(line[ncol[k]-1])
 
 			#legenditems.append(varname[k] + u' (модель)') #+ ', ' +legendtime[i]+expnames[np])
-			legenditems.append(u'модель') 
+			legenditems.append(varname[k]) 
 			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_[k])],linewidth=3,\