DevelopmentActivities/CMIP6/DevelopmentsCMIP6/soil_physic: NewHydroVert.py

import numpy as np
from matplotlib.backends.backend_pdf import PdfPages
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
#
# Compute layer thickness
#
def CompThickness(x):
    nbx=len(x)
    dx=0.0
    for i in range(1,nbx):
        dx = np.append(dx,(x[i-1]+x[i])/2.0)
    # Last layer
    dx=np.append(dx,(x[nbx-1]-dx[nbx-1])*2+dx[nbx-1])
    #
    tx=[]
    for i in range(len(dx)-1):
        tx = np.append(tx,dx[i+1]-dx[i])
    #
    return dx,tx
#
# Ploting function
#
def discretizationplot(wlev, tlev, wdepth_max, tsol_max, overalltitle):
    #
    # Compute thicknesses
    #
    wy,wz=CompThickness(wlev)
    ty,tz=CompThickness(tlev)
    #
    x=np.array([0,1])

    f, (ax1, ax2) = plt.subplots(1, 2, sharey=False)

    zmax=max([np.max(wz),np.max(tz)])
    cNorm = mpl.colors.Normalize(vmin=0, vmax=zmax)
    cmap = plt.get_cmap('Dark2')    

    ztmp=np.reshape(np.repeat(wz,2), (len(wy)-1, len(x)) )
    ax1.set_yscale('symlog')
    ax1.yaxis.set_major_formatter(mpl.ticker.ScalarFormatter())
    ax1.yaxis.set_major_formatter(mpl.ticker.FormatStrFormatter('%d'))
    ax1.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, pos: str(int(round(x)))))
#    ax1.set_yticks([-0.02, -0.2, -2])
    datamap=ax1.pcolor(x,-wy, ztmp, zorder=0, norm=cNorm, cmap=cmap)
    # Add tine lines to show the middle of the layer where SM is positioned
    for i in range(len(wlev)):
        ax1.plot(x,np.array([-wlev[i],-wlev[i]]), color='k', linewidth=0.4)
#
    ax1.text(0.5, -wdepth_max, "dpu_max", fontweight='bold',ha='center',
             va='center')
    ax1.set_xlabel("Min thickness="+str(round(np.min(wz),5))+"[m]")
    ax1.set_title("nslm="+str(len(wlev)))
#
    ztmp=np.reshape(np.repeat(tz, 2), (len(ty)-1,len(x)))
    datamap = ax2.pcolor(x,-ty, ztmp, zorder=0, norm=cNorm, cmap=cmap)
    # Add tine lines to show the middle of the layer where T is positioned
    for i in range(len(tlev)):
        ax2.plot(x,np.array([-tlev[i],-tlev[i]]), color='k', linewidth=0.4)

    ax2.text(0.5, -tsol_max, "tdepth_max", fontweight='bold',ha='center',
                 va='center')
    ax2.set_title("ngrnd="+str(len(tlev)))
    ax2.set_xlabel("Max thickness="+str(round(np.max(tz),5))+"[m]")
#
    cb=plt.colorbar(datamap, cmap=cmap, shrink=0.5, norm=cNorm)
    cb.set_label("Layer Thickness [m]")
#
    plt.suptitle(overalltitle)
#
#    plt.show()
    pdf.savefig()
#
#
# As a reminder, the levels of the hydrology as in the code today 
#
dpu_max=2.0
nslm=11
zz=np.zeros(nslm, float)
dz=np.zeros(nslm+1, float)
thickness=np.zeros(nslm, float)
depth=np.zeros(nslm+1, float)

for i in range(1,nslm):
    zz[i] = dpu_max*((2**i)-1)/((2**(nslm-1))-1)
    dz[i] = zz[i]-zz[i-1]

diaglev=np.zeros(nslm, float)
for i in range(nslm-1):
    diaglev[i] = dpu_max/(2**(nslm-1) -1) * ( ( 2**(i) -1) + ( 2**(i+1) -1) )/2.0
diaglev[nslm-1]=dpu_max
print "== Old values for reference =========================================="
print "zz=", zz
print "dz=",dz
print "diagdev=",diaglev
print "======================================================================"
print "  "
print "  "
#
#
# New Methode for computing levels with the new parameters to control.
#
# All units in meters
#
# Thickness of first Layer 
FirstLayer=9.77517107e-04
try :
    FirstLayer=float(raw_input("Depth of top hydrology layer [9.77517107e-04 m] >"))
except ValueError:
    FirstLayer=9.77517107e-04
#
# Maximum depth of soil moisture
dpu_max=2.0
try :
    dpu_max=float(raw_input("Depth of hydrology [2.0m] >"))
except ValueError:
    dpu_max=2.0
#
# The Linear layer thickness starts at :
LinearFromDepth=0.75
try :
    LinearFromDepth=float(raw_input("Depth at which linear hydrology layers start [0.75 m] >"))
except ValueError:
    LinearFromDepth=0.75
#
# Maximum depth of the soil thermodynamics
tdepth_max = 8.0
try :
    tdepth_max=float(raw_input("Depth of thermodynamics [8.0m] >"))
except ValueError:
    tdepth_max = 8.0
# Depth over which the thermal profile needs to be filterd
# It could be estimate based on the wave velocity of heat in a 
# dry soil and the Courant-Levy criteria.
TFilterDepth=0.05
#
# It is coded FORTRAN style so that translations is easier
#
nltmp=500
ztmp=np.zeros(nltmp, float)
dtmp=np.zeros(nltmp+1, float)
for i in range(1,nltmp):
    if ( ztmp[i-1] < LinearFromDepth ) :
        ztmp[i] = FirstLayer*2.0*((2**i)-1)
        dtmp[i] = ztmp[i]-ztmp[i-1]
    else :
        ztmp[i] = ztmp[i-1]+dtmp[i-1]
        dtmp[i] = dtmp[i-1]
nslm=np.argmin(np.abs(ztmp-dpu_max))+1
#
# Extract now the values we need to reach dpu_max
#
zz=ztmp[0:nslm]
dz=dtmp[0:nslm+1]
dz[nslm]=0.0
diaglev=np.zeros(nslm, float)
diaglev[0]=dz[1]/2.0
for i in range(1,nslm):
    diaglev[i] = diaglev[i-1]+(dz[i]+dz[i+1])/2.0
print "== New values ========================================================"
print "Number of layers in hydrol = ", nslm
print "zz=", zz
print "dz=",dz
print "diagdev=",diaglev
#
# Extension for the thermodynamics
#
ztmp=np.zeros(nltmp, float)
for i in range(nslm):
    ztmp[i]=diaglev[i]
for i in range(nslm,nltmp):
    ztmp[i]=ztmp[i-1]+2.0*(ztmp[i-1]-ztmp[i-2])
ngrnd=np.argmin(np.abs(ztmp-tdepth_max))+1
#
tdepth=np.zeros(ngrnd, float)
tdepth=ztmp[0:ngrnd]
print "Number of layers in thermosoil = ", ngrnd
print "tdepth=",tdepth
print "======================================================================"
#
# Do the plot
#
pdf = PdfPages('VerticalLayers.pdf')
discretizationplot(diaglev, tdepth, dpu_max, tdepth_max, "mixing water and temperature")
pdf.close()