source: lmdz_wrf/trunk/tools/iniaqua.py @ 195

## L. Fita, LMD October 2014. 
# Generation of initial conditions for an aqua-planet from the LMDZ model (r1818) 'iniacademic.F90' program
#   Author:    Frederic Hourdin      original: 15/01/93
# The forcing defined here is from Held and Suarez, 1994, Bulletin
# of the American Meteorological Society, 75, 1825.

from optparse import OptionParser
import numpy as np
from netCDF4 import Dataset as NetCDFFile
import os
import re
import nc_var_tools as ncvar
from lmdz_const import *

main = 'iniaqua.py'
errormsg='ERROR -- error -- ERROR -- error'
warnmsg='WARNING -- warning -- WARNING -- warning'

filekinds = ['CF', 'WRF']

## g.e. # iniaqua.py -d 32,32,39 -o WRF -t 19791201000000 -y true

def sig_hybrid(sig,pa,preff):
    """ Function utilisee pour calculer des valeurs de sigma modifie
     pour conserver les coordonnees verticales decrites dans
     esasig.def/z2sig.def lors du passage en coordonnees hybrides
     F. Forget 2002
       sig= sigma level
       pa=
       preff= reference pressure
     Connaissant sig (niveaux "sigma" ou on veut mettre les couches)
     L'objectif est de calculer newsig telle que
       (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig
     Cela ne se resoud pas analytiquement: 
     => on resoud par iterration bourrine 
     ----------------------------------------------
     Information  : where exp(1-1./x**2) become << x
           x      exp(1-1./x**2) /x
           1           1
           0.68       0.5
           0.5        1.E-1
           0.391      1.E-2
           0.333      1.E-3
           0.295      1.E-4
           0.269      1.E-5
           0.248      1.E-6
        => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25
    """
    fname = 'sig_hybrid'

# maximum number of iterations
    maxiter = 9999

    nsig = sig
    x1=0
    x2=1
    if sig >= 1.:
        nsig = sig
    elif sig*preff/pa >= 0.25:
        for j in range(maxiter):
            F = ((1. -pa/preff)*np.exp(1.-1./nsig**2)+(pa/preff)*nsig)/sig
#            print J,'nsig=', newsig, 'F=', F
            if F > 1:
                X2 = newsig
                nsig = (X1+nsig)*0.5
            else:
                X1 = newsig
                nsig = (X2+nsig)*0.5

#       Test : on arete lorsque on approxime sig a moins de 0.01 m pres 
#         (en pseudo altitude) :
            if np.abs(10.*np.log(F)) < 1.e-5: break
    else:
        nsig= sig*preff/pa

    return nsig

def presnivs_calc(hybrid, dz):
    """ From dyn3d/disvert_noterre.F calculation of vertical pressure levels
      hybrid= whether hydbrid coordinates have to be used
      dz= numbef of vertical layers
    """

    fname = 'presnivs_calc'

#-----------------------------------------------------------------------
#    ....  Calculs  de ap(l) et de bp(l)  ....
#    .........................................
#
#   .....  pa et preff sont lus  sur les fichiers start par dynetat0 .....
#-----------------------------------------------------------------------
#
    llmp1 = dz + 1
    pnivs = np.zeros((dz), dtype=np.float)
    s = np.zeros((dz), dtype=np.float)
    sig = np.zeros((dz+1), dtype=np.float)
    ap = np.zeros((dz+1), dtype=np.float)
    bp = np.zeros((dz+1), dtype=np.float)

# Ouverture possible de fichiers typiquement E.T.

    if os.path.isfile('easig.def'):
#-----------------------------------------------------------------------
#   cas 1 on lit les options dans esasig.def:
#   ----------------------------------------

        ofet = open('esasig.def', 'r')
#        Lecture de esasig.def :
#        Systeme peu souple, mais qui respecte en theorie
#        La conservation de l'energie (conversion Energie potentielle
#        <-> energie cinetique, d'apres la note de Frederic Hourdin...

        print '*****************************'
        print "WARNING reading 'esasig.def'"
        print '*****************************'
        for line in ofet:
            linevalues = ncvar.reduce_spaces(line)
            scaleheight = np.float(linevalues[0])
            dz0 = np.float(linevalues[1])
            dz1 = np.float(linevalues[2])
            nhaut = np.float(linevalues[3])

        ofet.close()
        dz0 = dz0/scaleheight
        dz1 = dz1/scaleheight

        sig1=(1.-dz1)/tanh(.5*(dz-1)/nhaut)

        esig=1.
        for l in range(19):
            esig=-np.log((1./sig1-1.)*np.exp(-dz0)/esig)/(dz-1.)

        csig=(1./sig1-1.)/(np.exp(esig)-1.)

        for l in range(1,dz):
            zz=csig*(np.exp(esig*(l-1.))-1.)
            sig[l] = 1./(1.+zz)*np.tanh(.5*(dz+1-l)/nhaut)

        sig[0] = 1.
        sig[dz] = 0.
        quoi = 1. + 2.* kappa
        s[dz-1]  = 1.
        s[dz-2] = quoi
        if dz > 1:
            for ll in range(1, dz-2):
                l = dz+1 - ll
                quand = sig[l+1]/sig[l]
                s[l-1] = quoi*(1.-quand)*s[l] + quand*s[l+1]
#
        snorm=(1.-.5*sig[1]+kappa*(1.-sig[1]))*s[0]+.5*sig[1]*s[1]
        for l in range(dz):
            s[l] = s[l]/ snorm
    elif os.path.isfile('z2sig.def'):
        fet = open('z2sig.def', 'r')
#-----------------------------------------------------------------------
#   cas 2 on lit les options dans z2sig.def:
#   ----------------------------------------
        print '****************************'
        print 'Reading z2sig.def'
        print '****************************'

        for line in ofet:
            linevalues = ncvar.reduce_spaces(line)
            scaleheight = np.float(linevalues[0])
            for l in range(dz):
                zsig[l] = linevalues[l+1]

        ofet.close()

        sig[0] = 1.
        for l in range(1,dz):
            sig[l] = 0.5 * ( np.exp(-zsig[l]/scaleheight) +                          \
              np.exp(-zsig[l-1]/scaleheight) )

        sig[dz+1] = 0.

#-----------------------------------------------------------------------
    else:
        print errormsg
        print '  ' + fname + ": didn't you forget something ???"
        print "    We need file  'z2sig.def' ! (OR 'esasig.def')"
        quit(-1)
#-----------------------------------------------------------------------

    nivsigs = np.arange(dz)*1.
    nivsig = np.arange(llmp1)*1.

    if hybrid:
# use hybrid coordinates
        print "*********************************"
        print "Using hybrid vertical coordinates"
        print 
#        Coordonnees hybrides avec mod
        for l in range(dz):
            newsig = sig_hybrid(sig[l],pa,preff)
            bp[l] = np.exp(1.-1./(newsig**2))
            ap[l] = pa * (newsig - bp[l] )

        bp[llmp1-1] = 0.
        ap[llmp1-1] = 0.
    else:
# use sigma coordinates
        print "********************************"
        print "Using sigma vertical coordinates"
        print
#       Pour ne pas passer en coordonnees hybrides
        for l in range(dz):
            ap[l] = 0.
            bp[l] = sig[l]

        ap[llmp1-1] = 0.

    bp[llmp1-1] = 0.

    print 'BP________ ', bp
    print 'AP________ ', ap

#   Calcul au milieu des couches  (llm = dz):
#   WARNING : le choix de placer le milieu des couches au niveau de
#   pression intermediaire est arbitraire et pourrait etre modifie.
#   Le calcul du niveau pour la derniere couche 
#   (on met la meme distance (en log pression)  entre P(llm)
#   et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
#   Specifique.  Ce choix est specifie ici ET dans exner_milieu.F

    for l in range(dz-1):
        aps[0:dz-1] = 0.5*( ap[0:dz-1]+ap[1:dz]) 
        bps[0:dz-1] = 0.5*( bp[0:dz-1]+bp[1:dz]) 

    if hybrid:
        aps[dz-1] = aps[dz-2]**2 / aps[dz-3] 
        bps[dz-1] = 0.5*(bp[dz-1] + bp[dz])
    else:
        bps[dz-1] = bps[dz-2]**2 / bps[dz-3] 
        aps[dz-1] = 0.

    print 'BPs_______ ', bps
    print 'APs_______ ', aps

    for l in range(dz):
        psnivs[l] = aps[l]+bps[l]*preff
        psalt[l] = -scaleheight*np.log(presnivs[l]/preff)

    return psnivs, psalt

def global_lonlat(dx,dy):
    """ Function to generate 2D matrices with the global longitude, latitudes
      dx, dy: dimension of the desired matrix
    >>> global_lonlat(5,5)
    array([[  36.,  108.,  180.,  252.,  324.],
          [  36.,  108.,  180.,  252.,  324.],
          [  36.,  108.,  180.,  252.,  324.],
          [  36.,  108.,  180.,  252.,  324.],
          [  36.,  108.,  180.,  252.,  324.]]), array([[-72., -72., -72., -72., -72.],
          [-36., -36., -36., -36., -36.],
          [  0.,   0.,   0.,   0.,   0.],
          [ 36.,  36.,  36.,  36.,  36.],
          [ 72.,  72.,  72.,  72.,  72.]]))
    """ 

    fname = 'global_lonlat'

    longitude = np.zeros((dy,dx), dtype=np.float)
    latitude = np.zeros((dy,dx), dtype=np.float)

    for ix in range(dx):
        longitude[:,ix] = 360.*(1./2. + ix )/(dx)

    for iy in range(dy):
        latitude[iy,:] =  180.*(1./2. + iy )/(dy) - 90.

    return longitude, latitude

####### ###### ##### #### ### ## #

filekindsnames = "'" + ncvar.numVector_String(filekinds, "', '") + "'"

parser = OptionParser()
parser.add_option("-o", "--outputkind", type='choice', dest="okind", 
  choices=filekinds, help="kind of output: " + filekindsnames, metavar="KIND")
parser.add_option("-d", "--dimensions", dest="dims", 
  help="dimensions of the initial conditions: dimx,dimy,dimz", metavar="VALUES")
parser.add_option("-t", "--time", dest="time", 
  help="time of the initial conditions ([YYYY][MM][DD][HH][MI][SS] format)", metavar="VALUE")
parser.add_option("-y", "--hybrid", dest="hybrid",
  help="whether vertical levels have to compute hydbid or not (true/false)", metavar="VAR")

(opts, args) = parser.parse_args()

#######    #######
## MAIN
    #######

#  dynamic variables
#  vcov, ucov: covariant winds
#  teta: potential temperature
#  q: advecting fields (humidity species)
#  ps: surface pressure
#  masse: air mass
#  phis: surface geopotential

#  Local:
#  p: pressure at the interface between layers (half-sigma)
#  pks: Exner function at the surface
#  pk: Exner functino at the half-sigma layers
#  pkf: Filtred Exner function at the half-sigma layers
#  phi: geopotential height
#  ddsin,zsig,tetapv,w_pv: auxiliar variables
#  tetastrat: potential temporeature in the stratosphere (K)
#  teta0,ttp,delt_y,delt_z,eps: constants for the T profile
#  k_f,k_c_a,k_c_s: calling constants
#  ok_geost: Initialisation geostrohic wind
#  ok_pv: Polar Vortex
#  phi_pv,dphi_pv,gam_pv: polar vortex constants

dimx = int(opts.dims.split(',')[0])
dimy = int(opts.dims.split(',')[1])
dimz = int(opts.dims.split(',')[2])

if opts.hybrid == 'true':
    hybrid_calc = True
else:
    hybrid_calc = False

ofile = 'iniaqua.nc'

print 'dimensions: ',dimx,', ',dimy,', ',dimz

if opts.okind == 'CF':
    varnames = ['zg', 'ta', 'ua', 'va', 'pres', 'r']
# Reference time from 1949-12-01 00:00:00 UTC
    timev = ncvar.datetimeStr_conversion(opts.time,'YmdHMS','cfTime')
    dimnames = ['time', 'z', 'y', 'x']
elif opts.okind == 'WRF':
    varnames = ['Geopot', 'T', 'U', 'V', 'PRES', 'QVAPOR']
    timev = ncvar.datetimeStr_conversion(opts.time,'YmdHMS','WRFdatetime')
    dimnames = ['Time', 'bottom_top', 'south_north', 'west_east']

# Constants
##
llm = dimz

ok_geost = True
# Constants for Newtonian relaxation and friction
k_f = 1.
k_f = 1./(daysec*k_f)
# cooling surface
k_c_s=4.  
k_c_s=1./(daysec*k_c_s)
# cooling free atm
k_c_a=40. 
k_c_a=1./(daysec*k_c_a)
# Constants for Teta equilibrium profile
# mean Teta (S.H. 315K)
teta0=315.
# Tropopause temperature (S.H. 200K)
ttp=200.
# Deviation to N-S symmetry(~0-20K)
eps=0.
# Merid Temp. Gradient (S.H. 60K)
delt_y=60.
# Vertical Gradient (S.H. 10K)
delt_z=10.
# Polar vortex
ok_pv = False
# Latitude of edge of vortex
phi_pv=-50.
phi_pv=phi_pv*pi/180.
# Width of the edge
dphi_pv=5.
dphi_pv=dphi_pv*pi/180.
# -dT/dz vortex (in K/km)
gam_pv=4.

presnivs, pseudoalt = presnivs_calc(hybrid_calc, dimz)
lon, lat = global_lonlat(dx,dy)
lonu, latu = global_lonlat(dx+1,dy)
lonv, latv = global_lonlat(dx,dy+1)

# 2. Initialize fields towards which to relax
##

knewt_t = np.zeros((dimz), dtype=np.float)
kfrict = np.zeros((dimz), dtype=np.float)
clat4 =  np.zeros((dimy+1, dimx), dtype=np.float)

# Friction
knewt_g = k_c_a
for l in range(dimz):
    zsig=presnivs[l]/preff
    knewt_t[l]=(k_c_s-k_c_a)*np.max([0.,(zsig-0.7)/0.3])
    kfrict[l]=k_f*np.max([0.,(zsig-0.7)/0.3])

    for j in 1,range(dimy+1):
        clat4[j,:]=np.cos(rlatu[j])**4

# Vertical profile
tetajl =  np.zeros((dimy+1, dimx), dtype=np.float)
teta = np.zeros((dimy+1, dimx+1), dtype=np.float)
tetarappel = np.zeros((dimy+1, dimx+1), dtype=np.float)

for l in range (dz):
    zsig = presnivs[l]/preff
    tetastrat = ttp*zsig**(-kappa)
    tetapv = tetastrat
    if ok_pv and zsig < 0.1:
        tetapv = tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)

    ddsin = np.sin(latu)
    tetajl[l,:,:] = teta0-delt_y*ddsin+eps*ddsin-delt_z*(1.-ddsin*ddsin)*np.log(zsig)
    if planet_type == 'giant':
        tetajl[l,:,:] = teta0+(delt_y*((np.sin(latu*3.14159*eps+0.0001))**2) /       \
          ((rlatu*3.14159*eps+0.0001)**2))-delt_z*np.log(zsig)

# Profile stratospheric isotherm (+vortex)
    w_pv=(1.-np.tanh((rlatu-phi_pv)/dphi_pv))/2.
    tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
    for iy in range(dy):
        for ix in range(dx):
            tetajl[l,iy,ix]=np.max([tetajl(l,ix,iy),tetastrat])

for iz in range(dimz):
    for iy in range(dimy+1):
        tetarappel[iz,iy,0:dimx+1] = tetajl[iz,iy,:]

    tetarappel[iz,iy,0:dimx] = tetajl[iz,iy,dimx-1]


ncf = NetCDFFile(ofile, 'w')

# Dimensions creation
ncf.createDimension(dimnames[0], None)
ncf.createDimension(dimnames[1], dimz)
ncf.createDimension(dimnames[2], dimy)
ncf.createDimension(dimnames[3], dimx)

ncf.sync()
ncf.close()

print main + ": successfull writing of file '" + ofile + "' !!"