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

Last change on this file since 186 was 180, checked in by lfita, 11 years ago

Adding python script to generate same initial conditions as 'iniaqua'
File size: 13.5 KB
Line 
1## L. Fita, LMD October 2014.
2# Generation of initial conditions for an aqua-planet from the LMDZ model (r1818) 'iniacademic.F90' program
3#   Author:    Frederic Hourdin      original: 15/01/93
4# The forcing defined here is from Held and Suarez, 1994, Bulletin
5# of the American Meteorological Society, 75, 1825.
6
7from optparse import OptionParser
8import numpy as np
9from netCDF4 import Dataset as NetCDFFile
10import os
11import re
12import nc_var_tools as ncvar
13from lmdz_const import *
14
15main = 'iniaqua.py'
16errormsg='ERROR -- error -- ERROR -- error'
17warnmsg='WARNING -- warning -- WARNING -- warning'
18
19filekinds = ['CF', 'WRF']
20
21## g.e. # iniaqua.py -d 32,32,39 -o WRF -t 19791201000000 -y true
22
23def sig_hybrid(sig,pa,preff):
24    """ Function utilisee pour calculer des valeurs de sigma modifie
25     pour conserver les coordonnees verticales decrites dans
26     esasig.def/z2sig.def lors du passage en coordonnees hybrides
27     F. Forget 2002
28       sig= sigma level
29       pa=
30       preff= reference pressure
31     Connaissant sig (niveaux "sigma" ou on veut mettre les couches)
32     L'objectif est de calculer newsig telle que
33       (1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig = sig
34     Cela ne se resoud pas analytiquement:
35     => on resoud par iterration bourrine
36     ----------------------------------------------
37     Information  : where exp(1-1./x**2) become << x
38           x      exp(1-1./x**2) /x
39           1           1
40           0.68       0.5
41           0.5        1.E-1
42           0.391      1.E-2
43           0.333      1.E-3
44           0.295      1.E-4
45           0.269      1.E-5
46           0.248      1.E-6
47        => on peut utiliser newsig = sig*preff/pa si sig*preff/pa < 0.25
48    """
49    fname = 'sig_hybrid'
50
51# maximum number of iterations
52    maxiter = 9999
53
54    nsig = sig
55    x1=0
56    x2=1
57    if sig >= 1.:
58        nsig = sig
59    elif sig*preff/pa >= 0.25:
60        for j in range(maxiter):
61            F = ((1. -pa/preff)*np.exp(1.-1./nsig**2)+(pa/preff)*nsig)/sig
62#            print J,'nsig=', newsig, 'F=', F
63            if F > 1:
64                X2 = newsig
65                nsig = (X1+nsig)*0.5
66            else:
67                X1 = newsig
68                nsig = (X2+nsig)*0.5
69
70#       Test : on arete lorsque on approxime sig a moins de 0.01 m pres
71#         (en pseudo altitude) :
72            if np.abs(10.*np.log(F)) < 1.e-5: break
73    else:
74        nsig= sig*preff/pa
75
76    return nsig
77
78def presnivs_calc(hybrid, dz):
79    """ From dyn3d/disvert_noterre.F calculation of vertical pressure levels
80      hybrid= whether hydbrid coordinates have to be used
81      dz= numbef of vertical layers
82    """
83
84    fname = 'presnivs_calc'
85
86#-----------------------------------------------------------------------
87#    ....  Calculs  de ap(l) et de bp(l)  ....
88#    .........................................
89#
90#   .....  pa et preff sont lus  sur les fichiers start par dynetat0 .....
91#-----------------------------------------------------------------------
92#
93    llmp1 = dz + 1
94    pnivs = np.zeros((dz), dtype=np.float)
95    s = np.zeros((dz), dtype=np.float)
96    sig = np.zeros((dz+1), dtype=np.float)
97    ap = np.zeros((dz+1), dtype=np.float)
98    bp = np.zeros((dz+1), dtype=np.float)
99
100# Ouverture possible de fichiers typiquement E.T.
101
102    if os.path.isfile('easig.def'):
103#-----------------------------------------------------------------------
104#   cas 1 on lit les options dans esasig.def:
105#   ----------------------------------------
106
107        ofet = open('esasig.def', 'r')
108#        Lecture de esasig.def :
109#        Systeme peu souple, mais qui respecte en theorie
110#        La conservation de l'energie (conversion Energie potentielle
111#        <-> energie cinetique, d'apres la note de Frederic Hourdin...
112
113        print '*****************************'
114        print "WARNING reading 'esasig.def'"
115        print '*****************************'
116        for line in ofet:
117            linevalues = ncvar.reduce_spaces(line)
118            scaleheight = np.float(linevalues[0])
119            dz0 = np.float(linevalues[1])
120            dz1 = np.float(linevalues[2])
121            nhaut = np.float(linevalues[3])
122
123        ofet.close()
124        dz0 = dz0/scaleheight
125        dz1 = dz1/scaleheight
126
127        sig1=(1.-dz1)/tanh(.5*(dz-1)/nhaut)
128
129        esig=1.
130        for l in range(19):
131            esig=-np.log((1./sig1-1.)*np.exp(-dz0)/esig)/(dz-1.)
132
133        csig=(1./sig1-1.)/(np.exp(esig)-1.)
134
135        for l in range(1,dz):
136            zz=csig*(np.exp(esig*(l-1.))-1.)
137            sig[l] = 1./(1.+zz)*np.tanh(.5*(dz+1-l)/nhaut)
138
139        sig[0] = 1.
140        sig[dz] = 0.
141        quoi = 1. + 2.* kappa
142        s[dz-1]  = 1.
143        s[dz-2] = quoi
144        if dz > 1:
145            for ll in range(1, dz-2):
146                l = dz+1 - ll
147                quand = sig[l+1]/sig[l]
148                s[l-1] = quoi*(1.-quand)*s[l] + quand*s[l+1]
149#
150        snorm=(1.-.5*sig[1]+kappa*(1.-sig[1]))*s[0]+.5*sig[1]*s[1]
151        for l in range(dz):
152            s[l] = s[l]/ snorm
153    elif os.path.isfile('z2sig.def'):
154        fet = open('z2sig.def', 'r')
155#-----------------------------------------------------------------------
156#   cas 2 on lit les options dans z2sig.def:
157#   ----------------------------------------
158        print '****************************'
159        print 'Reading z2sig.def'
160        print '****************************'
161
162        for line in ofet:
163            linevalues = ncvar.reduce_spaces(line)
164            scaleheight = np.float(linevalues[0])
165            for l in range(dz):
166                zsig[l] = linevalues[l+1]
167
168        ofet.close()
169
170        sig[0] = 1.
171        for l in range(1,dz):
172            sig[l] = 0.5 * ( np.exp(-zsig[l]/scaleheight) +                          \
173              np.exp(-zsig[l-1]/scaleheight) )
174
175        sig[dz+1] = 0.
176
177#-----------------------------------------------------------------------
178    else:
179        print errormsg
180        print '  ' + fname + ": didn't you forget something ???"
181        print "    We need file  'z2sig.def' ! (OR 'esasig.def')"
182        quit(-1)
183#-----------------------------------------------------------------------
184
185    nivsigs = np.arange(dz)*1.
186    nivsig = np.arange(llmp1)*1.
187
188    if hybrid:
189# use hybrid coordinates
190        print "*********************************"
191        print "Using hybrid vertical coordinates"
192        print 
193#        Coordonnees hybrides avec mod
194        for l in range(dz):
195            newsig = sig_hybrid(sig[l],pa,preff)
196            bp[l] = np.exp(1.-1./(newsig**2))
197            ap[l] = pa * (newsig - bp[l] )
198
199        bp[llmp1-1] = 0.
200        ap[llmp1-1] = 0.
201    else:
202# use sigma coordinates
203        print "********************************"
204        print "Using sigma vertical coordinates"
205        print
206#       Pour ne pas passer en coordonnees hybrides
207        for l in range(dz):
208            ap[l] = 0.
209            bp[l] = sig[l]
210
211        ap[llmp1-1] = 0.
212
213    bp[llmp1-1] = 0.
214
215    print 'BP________ ', bp
216    print 'AP________ ', ap
217
218#   Calcul au milieu des couches  (llm = dz):
219#   WARNING : le choix de placer le milieu des couches au niveau de
220#   pression intermediaire est arbitraire et pourrait etre modifie.
221#   Le calcul du niveau pour la derniere couche
222#   (on met la meme distance (en log pression)  entre P(llm)
223#   et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
224#   Specifique.  Ce choix est specifie ici ET dans exner_milieu.F
225
226    for l in range(dz-1):
227        aps[0:dz-1] = 0.5*( ap[0:dz-1]+ap[1:dz]) 
228        bps[0:dz-1] = 0.5*( bp[0:dz-1]+bp[1:dz]) 
229
230    if hybrid:
231        aps[dz-1] = aps[dz-2]**2 / aps[dz-3] 
232        bps[dz-1] = 0.5*(bp[dz-1] + bp[dz])
233    else:
234        bps[dz-1] = bps[dz-2]**2 / bps[dz-3] 
235        aps[dz-1] = 0.
236
237    print 'BPs_______ ', bps
238    print 'APs_______ ', aps
239
240    for l in range(dz):
241        psnivs[l] = aps[l]+bps[l]*preff
242        psalt[l] = -scaleheight*np.log(presnivs[l]/preff)
243
244    return psnivs, psalt
245
246def global_lonlat(dx,dy):
247    """ Function to generate 2D matrices with the global longitude, latitudes
248      dx, dy: dimension of the desired matrix
249    >>> global_lonlat(5,5)
250    array([[  36.,  108.,  180.,  252.,  324.],
251          [  36.,  108.,  180.,  252.,  324.],
252          [  36.,  108.,  180.,  252.,  324.],
253          [  36.,  108.,  180.,  252.,  324.],
254          [  36.,  108.,  180.,  252.,  324.]]), array([[-72., -72., -72., -72., -72.],
255          [-36., -36., -36., -36., -36.],
256          [  0.,   0.,   0.,   0.,   0.],
257          [ 36.,  36.,  36.,  36.,  36.],
258          [ 72.,  72.,  72.,  72.,  72.]]))
259    """ 
260
261    fname = 'global_lonlat'
262
263    longitude = np.zeros((dy,dx), dtype=np.float)
264    latitude = np.zeros((dy,dx), dtype=np.float)
265
266    for ix in range(dx):
267        longitude[:,ix] = 360.*(1./2. + ix )/(dx)
268
269    for iy in range(dy):
270        latitude[iy,:] =  180.*(1./2. + iy )/(dy) - 90.
271
272    return longitude, latitude
273
274####### ###### ##### #### ### ## #
275
276filekindsnames = "'" + ncvar.numVector_String(filekinds, "', '") + "'"
277
278parser = OptionParser()
279parser.add_option("-o", "--outputkind", type='choice', dest="okind", 
280  choices=filekinds, help="kind of output: " + filekindsnames, metavar="KIND")
281parser.add_option("-d", "--dimensions", dest="dims", 
282  help="dimensions of the initial conditions: dimx,dimy,dimz", metavar="VALUES")
283parser.add_option("-t", "--time", dest="time", 
284  help="time of the initial conditions ([YYYY][MM][DD][HH][MI][SS] format)", metavar="VALUE")
285parser.add_option("-y", "--hybrid", dest="hybrid",
286  help="whether vertical levels have to compute hydbid or not (true/false)", metavar="VAR")
287
288(opts, args) = parser.parse_args()
289
290#######    #######
291## MAIN
292    #######
293
294#  dynamic variables
295#  vcov, ucov: covariant winds
296#  teta: potential temperature
297#  q: advecting fields (humidity species)
298#  ps: surface pressure
299#  masse: air mass
300#  phis: surface geopotential
301
302#  Local:
303#  p: pressure at the interface between layers (half-sigma)
304#  pks: Exner function at the surface
305#  pk: Exner functino at the half-sigma layers
306#  pkf: Filtred Exner function at the half-sigma layers
307#  phi: geopotential height
308#  ddsin,zsig,tetapv,w_pv: auxiliar variables
309#  tetastrat: potential temporeature in the stratosphere (K)
310#  teta0,ttp,delt_y,delt_z,eps: constants for the T profile
311#  k_f,k_c_a,k_c_s: calling constants
312#  ok_geost: Initialisation geostrohic wind
313#  ok_pv: Polar Vortex
314#  phi_pv,dphi_pv,gam_pv: polar vortex constants
315
316dimx = int(opts.dims.split(',')[0])
317dimy = int(opts.dims.split(',')[1])
318dimz = int(opts.dims.split(',')[2])
319
320if opts.hybrid == 'true':
321    hybrid_calc = True
322else:
323    hybrid_calc = False
324
325ofile = 'iniaqua.nc'
326
327print 'dimensions: ',dimx,', ',dimy,', ',dimz
328
329if opts.okind == 'CF':
330    varnames = ['zg', 'ta', 'ua', 'va', 'pres', 'r']
331# Reference time from 1949-12-01 00:00:00 UTC
332    timev = ncvar.datetimeStr_conversion(opts.time,'YmdHMS','cfTime')
333    dimnames = ['time', 'z', 'y', 'x']
334elif opts.okind == 'WRF':
335    varnames = ['Geopot', 'T', 'U', 'V', 'PRES', 'QVAPOR']
336    timev = ncvar.datetimeStr_conversion(opts.time,'YmdHMS','WRFdatetime')
337    dimnames = ['Time', 'bottom_top', 'south_north', 'west_east']
338
339# Constants
340##
341llm = dimz
342
343ok_geost = True
344# Constants for Newtonian relaxation and friction
345k_f = 1.
346k_f = 1./(daysec*k_f)
347# cooling surface
348k_c_s=4. 
349k_c_s=1./(daysec*k_c_s)
350# cooling free atm
351k_c_a=40. 
352k_c_a=1./(daysec*k_c_a)
353# Constants for Teta equilibrium profile
354# mean Teta (S.H. 315K)
355teta0=315.
356# Tropopause temperature (S.H. 200K)
357ttp=200.
358# Deviation to N-S symmetry(~0-20K)
359eps=0.
360# Merid Temp. Gradient (S.H. 60K)
361delt_y=60.
362# Vertical Gradient (S.H. 10K)
363delt_z=10.
364# Polar vortex
365ok_pv = False
366# Latitude of edge of vortex
367phi_pv=-50.
368phi_pv=phi_pv*pi/180.
369# Width of the edge
370dphi_pv=5.
371dphi_pv=dphi_pv*pi/180.
372# -dT/dz vortex (in K/km)
373gam_pv=4.
374
375presnivs, pseudoalt = presnivs_calc(hybrid_calc, dimz)
376lon, lat = global_lonlat(dx,dy)
377lonu, latu = global_lonlat(dx+1,dy)
378lonv, latv = global_lonlat(dx,dy+1)
379
380# 2. Initialize fields towards which to relax
381##
382
383knewt_t = np.zeros((dimz), dtype=np.float)
384kfrict = np.zeros((dimz), dtype=np.float)
385clat4 =  np.zeros((dimy+1, dimx), dtype=np.float)
386
387# Friction
388knewt_g = k_c_a
389for l in range(dimz):
390    zsig=presnivs[l]/preff
391    knewt_t[l]=(k_c_s-k_c_a)*np.max([0.,(zsig-0.7)/0.3])
392    kfrict[l]=k_f*np.max([0.,(zsig-0.7)/0.3])
393
394    for j in 1,range(dimy+1):
395        clat4[j,:]=np.cos(rlatu[j])**4
396
397# Vertical profile
398tetajl =  np.zeros((dimy+1, dimx), dtype=np.float)
399teta = np.zeros((dimy+1, dimx+1), dtype=np.float)
400tetarappel = np.zeros((dimy+1, dimx+1), dtype=np.float)
401
402for l in range (dz):
403    zsig = presnivs[l]/preff
404    tetastrat = ttp*zsig**(-kappa)
405    tetapv = tetastrat
406    if ok_pv and zsig < 0.1:
407        tetapv = tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
408
409    ddsin = np.sin(latu)
410    tetajl[l,:,:] = teta0-delt_y*ddsin+eps*ddsin-delt_z*(1.-ddsin*ddsin)*np.log(zsig)
411    if planet_type == 'giant':
412        tetajl[l,:,:] = teta0+(delt_y*((np.sin(latu*3.14159*eps+0.0001))**2) /       \
413          ((rlatu*3.14159*eps+0.0001)**2))-delt_z*np.log(zsig)
414
415# Profile stratospheric isotherm (+vortex)
416    w_pv=(1.-np.tanh((rlatu-phi_pv)/dphi_pv))/2.
417    tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
418    for iy in range(dy):
419        for ix in range(dx):
420            tetajl[l,iy,ix]=np.max([tetajl(l,ix,iy),tetastrat])
421
422for iz in range(dimz):
423    for iy in range(dimy+1):
424        tetarappel[iz,iy,0:dimx+1] = tetajl[iz,iy,:]
425
426    tetarappel[iz,iy,0:dimx] = tetajl[iz,iy,dimx-1]
427
428
429ncf = NetCDFFile(ofile, 'w')
430
431# Dimensions creation
432ncf.createDimension(dimnames[0], None)
433ncf.createDimension(dimnames[1], dimz)
434ncf.createDimension(dimnames[2], dimy)
435ncf.createDimension(dimnames[3], dimx)
436
437ncf.sync()
438ncf.close()
439
440print main + ": successfull writing of file '" + ofile + "' !!"
Note: See TracBrowser for help on using the repository browser.