Index: trunk/tools/iniaqua.py
===================================================================
--- trunk/tools/iniaqua.py	(revision 180)
+++ trunk/tools/iniaqua.py	(revision 180)
@@ -0,0 +1,440 @@
+## 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 + "' !!"
Index: trunk/tools/lmdz_const.py
===================================================================
--- trunk/tools/lmdz_const.py	(revision 180)
+++ trunk/tools/lmdz_const.py	(revision 180)
@@ -0,0 +1,145 @@
+#
+# $Id: comconst.h 1795 2013-07-18 08:20:28Z emillour $
+#
+#----------------------------------------------------------------------
+import numpy as np
+
+# dtvr: dynamical time step (in s)
+# daysec: length (in s) of a standard day
+# pi: something like 3.14159....
+# dtphys: (s) time step for the physics
+# dtdiss: (s) time step for the dissipation
+# rad: (m) radius of the planet
+# r: Reduced Gas constant r=R/mu
+#   with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
+# cpp: Specific heat Cp (J.kg-1.K-1)
+# kappa: kappa=R/Cp 
+# cotot
+# unsim: 1./iim
+# g: (m/s2) gravity
+# omeg: (rad/s) rotation rate of the planet
+# dissip_factz,dissip_deltaz,dissip_zref
+## top_bound sponge:
+# iflag_top_bound: sponge type
+# mode_top_bound: sponge mode
+# tau_top_bound: inverse of sponge characteristic time scale (Hz)
+# daylen: length of solar day, in 'standard' day length
+# year_day: Number of standard days in a year
+# molmass: (g/mol) molar mass of the atmosphere
+# ihf: (W/m2) Intrinsic heat flux (for giant planets)
+
+daysec = 86400.
+preff = 1013250.
+pa = 500.
+
+#-----------------------------------------------------------------------
+# From suphel
+
+pi = 2.*np.arcsin(1.)
+clum = 299792458.
+hpla = 6.6260755e-34
+kbol = 1.380658e-23
+navo = 6.0221367e+23
+#
+#     ----------------------------------------------------------------
+#
+#*       2.    define astronomical constants.
+#              ------------------------------
+#
+day = 86400.
+ea = 149597870000.
+epsm = 0.409093
+#
+siyea = 365.25*day*2.*pi/6.283076
+siday = day/(1.+day/siyea)
+omega = 2.*pi/siday
+#
+#
+#     ------------------------------------------------------------------
+#
+#*       3.    define geoide.
+#              --------------
+#
+g = 9.80665
+a = 6371229.
+r1sa = np.float(np.float64(1.)/np.float64(a))
+#
+#     -----------------------------------------------------------------
+#
+#*       4.    define radiation constants.
+#              ---------------------------
+#
+# z.x.li      rsigma=2. * rpi**5 * rkbol**4 /(15.* rclum**2 * rhpla**3)
+sigma = 2.*pi**5 * (kbol/hpla)**3 * kbol/clum/clum/15.
+#
+#     -----------------------------------------------------------------
+#
+#*       5.    define thermodynamic constants, gas phase.
+#              ------------------------------------------
+#
+r = navo*kbol
+md = 28.9644
+mo3 = 47.9942
+mv = 18.0153
+rd = 1000.*r/md
+rv = 1000.*r/mv
+cpd = 3.5*rd
+cvd = cpd-rd
+cpv = 4.*rv
+cvv = cpv-rv
+kappa = rd/cpd
+etv = rv/rd-1.
+#
+#     ----------------------------------------------------------------
+#
+#*       6.    define thermodynamic constants, liquid phase.
+#              ---------------------------------------------
+#
+cw = cpv
+#     ----------------------------------------------------------------
+#
+#*       7.    define thermodynamic constants, solid phase.
+#              --------------------------------------------
+#
+cs = cpv
+#
+#     ----------------------------------------------------------------
+#
+#*       8.    define thermodynamic constants, transition of phase.
+#              ----------------------------------------------------
+#
+tt = 273.16
+lvtt = 2.5008e+6
+lstt = 2.8345e+6
+lmlt = lstt-lvtt
+atm = 100000.
+#
+#     ----------------------------------------------------------------
+#
+#*       9.    saturated vapour pressure.
+#              --------------------------
+#
+estt = 611.14
+gamw = (cw-cpv)/rv
+betw = lvtt/rv+gamw*tt
+alpw = np.log(estt)+betw/tt+gamw*np.log(tt)
+gams = (cs-cpv)/rv
+bets = lstt/rv+gams*tt
+alps = np.log(estt)+bets/tt+gams*np.log(tt)
+gamd = gams-gamw
+betd = bets-betw
+alpd = alps-alpw
+#
+#     ------------------------------------------------------------------
+#
+# calculer les constantes pour les fonctions thermodynamiques
+#
+vtmp2 = cpv/cpd-1.
+hoh2o = atm/100.
+r2es = estt*rd/rv
+r3les = 17.269
+r3ies = 21.875
+r4les = 35.86
+r4ies = 7.66
+r5les = r3les*(tt-r4les)
+r5ies = r3ies*(tt-r4ies)
