Context Navigation

iniaqua.py @ 664

Last change on this file since 664 was 224, checked in by lfita, 10 years ago
Working version up to the filling. Now problems with dimensions of the variables
File size: 94.4 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
7	from optparse import OptionParser
8	import numpy as np
9	from netCDF4 import Dataset as NetCDFFile
10	import os
11	import re
12	import nc_var_tools as ncvar
13	from lmdz_const import *
14
15	main = 'iniaqua.py'
16	errormsg='ERROR -- error -- ERROR -- error'
17	warnmsg='WARNING -- warning -- WARNING -- warning'
18
19	filekinds = ['CF', 'LMDZ', 'WRF']
20
21	## e.g. # iniaqua.py -d 32,32,39 -p hybdrid -o WRF -t 19791201000000 -z tropo -q 2
22	def fxy(dx, dy):
23	"""!
24	! $Id: fxy.F 1403 2010-07-01 09:02:53Z fairhead $
25	!
26	c Auteur : P. Le Van
27	c
28	c Calcul des longitudes et des latitudes pour une fonction f(x,y)
29	c a tangente sinusoidale et eventuellement avec zoom .
30	c
31	c
32	"""
33	fname ='fxy'
34
35	#c ...... calcul des latitudes et de y' .....
36	#c
37	vrlatu = np.zeros((dy+1), dtype=np.float)
38	vyprimu = np.zeros((dy+1), dtype=np.float)
39
40	vrlatu = ffy(np.arange(dy+1)*1. + 1.)
41	vyprimu = ffy(np.arange(dy+1)*1. + 1.)
42
43	vrlatv = np.zeros((dy), dtype=np.float)
44	vrlatu1 = np.zeros((dy), dtype=np.float)
45	vrlatu2 = np.zeros((dy), dtype=np.float)
46	yprimv = np.zeros((dy), dtype=np.float)
47	vyprimu1 = np.zeros((dy), dtype=np.float)
48	vyprimu2 = np.zeros((dy), dtype=np.float)
49
50	vrlatv = ffy(np.arange(dy)*1. + 1. + 0.5)
51	vrlatu1 = ffy(np.arange(dy)*1. + 1. + 0.25)
52	vrlatu2 = ffy(np.arange(dy)*1. + 1. + 0.75)
53
54	vyprimv = fyprim(np.arange(dy)*1. + 1. + 0.5)
55	vyprimu1 = fyprim(np.arange(dy)*1. + 1. + 0.25)
56	vyprimu2 = fyprim(np.arange(dy)*1. + 1. + 0.75)
57
58	#c
59	#c ..... calcul des longitudes et de x' .....
60	#c
61	vrlonv = np.zeros((dx+1), dtype=np.float)
62	vrlonu = np.zeros((dx+1), dtype=np.float)
63	vrlonm025 = np.zeros((dx+1), dtype=np.float)
64	vrlonp025 = np.zeros((dx+1), dtype=np.float)
65	vxprimv = np.zeros((dx+1), dtype=np.float)
66	vxprimu = np.zeros((dx+1), dtype=np.float)
67	vxprimm025 = np.zeros((dx+1), dtype=np.float)
68	vxprimo025 = np.zeros((dx+1), dtype=np.float)
69
70	vrlonv = ffy(np.arange(dx+1)*1. + 1.)
71	vrlonu = ffy(np.arange(dx+1)*1. + 1. + 0.5)
72	vrlonm025 = ffy(np.arange(dx+1)*1. + 1. - 0.25)
73	vrlonp025 = ffy(np.arange(dx+1)*1. + 1. + 0.25)
74
75	vxprimv = fxprim(np.arange(dx+1)*1. + 1.)
76	vxprimu = fxprim(np.arange(dx+1)*1. + 1. + 0.5)
77	vxprimm025 = fxprim(np.arange(dx+1)*1. + 1. - 0.25)
78	vxprimp025 = fxprim(np.arange(dx+1)*1. + 1. + 0.25)
79
80	return vrlatu, vyprimu, vrlatv, vyprimv, vrlatu1, vyprimu1, vrlatu2, vyprimu2, \
81	vrlonu, vxprimu, vrlonv, vxprimv, vrlonm025, vxprimm025, vrlonp025, vxprimp025
82
83	# From grid/fxy_sin.h
84	# ................................................................
85	# ................ Fonctions in line ...........................
86	# ................................................................
87	#
88
89	def fy(rj,dy):
90	"""
91	"""
92	val = np.arcsin(1.+2.*((1.-rj)/np.float(dy)))
93
94	return val
95
96	def fyprim(rj,dy):
97	"""
98	"""
99	val = 1./np.sqrt((rj-1.)*(dy+1.-rj))
100
101	return val
102
103	def fx(ri,dx):
104	"""
105	"""
106	val = 2.np.pi/np.float(dx) ( ri - 0.5*np.float(dx) - 1. )
107
108	return val
109
110	def fxprim(ri,dx):
111	"""
112	"""
113	val = 2.*np.pi/np.float(dx)
114
115	return val
116
117	#
118	#
119	# La valeur de pi est passee par le common/const/ou /const2/ .
120	# Sinon, il faut la calculer avant d'appeler ces fonctions .
121	#
122	# ----------------------------------------------------------------
123	# Fonctions a changer eventuellement, selon x(x) et y(y) choisis .
124	# -----------------------------------------------------------------
125	#
126	# ..... ici, on a l'application particuliere suivante ........
127	#
128	# **************************************
129	# ** x = 2. * pi/iim * X **
130	# ** y = pi/jjm * Y **
131	# **************************************
132	#
133	# ..................................................................
134	# ..................................................................
135	#
136	#
137	#
138	#-----------------------------------------------------------------------
139
140	def fxysinus (ddy,ddx):
141	"""c
142	c Calcul des longitudes et des latitudes pour une fonction f(x,y)
143	c avec y = Asin( j ) .
144	c
145	c Auteur : P. Le Van
146	c
147	c
148	from: dyn3d/fxysinus.F
149	"""
150	fname = 'fxysinus'
151
152	# ...... calcul des latitudes et de y' .....
153	#
154	rlatu = np.zeros((ddy+1), dtype=np.float)
155	yprimu = np.zeros((ddy+1), dtype=np.float)
156
157	for j in range(ddy+1):
158	rlatu[j] = fy( np.float(j+1), ddy)
159	yprimu[j] = fyprim( np.float(j+1), ddy)
160
161	rlatv = np.zeros((ddy), dtype=np.float)
162	rlatu1 = np.zeros((ddy), dtype=np.float)
163	rlatu2 = np.zeros((ddy), dtype=np.float)
164	yprimv = np.zeros((ddy), dtype=np.float)
165	yprimu1 = np.zeros((ddy), dtype=np.float)
166	yprimu2 = np.zeros((ddy), dtype=np.float)
167
168	for j in range(ddy):
169	rlatv[j] = fy( np.float(j) + 0.5, ddy)
170	rlatu1[j] = fy( np.float(j) + 0.25, ddy)
171	rlatu2[j] = fy( np.float(j) + 0.75, ddy)
172
173	yprimv[j] = fyprim( np.float(j) + 0.5, ddy)
174	yprimu1[j] = fyprim( np.float(j) + 0.25, ddy)
175	yprimu2[j] = fyprim( np.float(j) + 0.75, ddy)
176
177	#
178	# ..... calcul des longitudes et de x' .....
179	#
180	rlonv = np.zeros((ddx+1), dtype=np.float)
181	rlonu = np.zeros((ddx+1), dtype=np.float)
182	rlonm025 = np.zeros((ddx+1), dtype=np.float)
183	rlonp025 = np.zeros((ddx+1), dtype=np.float)
184	xprimv = np.zeros((ddx+1), dtype=np.float)
185	xprimu = np.zeros((ddx+1), dtype=np.float)
186	xprimm025 = np.zeros((ddx+1), dtype=np.float)
187	xprimp025 = np.zeros((ddx+1), dtype=np.float)
188
189	for i in range(ddx + 1):
190	rlonv[i] = fx( np.float(i), ddx )
191	rlonu[i] = fx( np.float(i)+ 0.5, ddx )
192	rlonm025[i] = fx( np.float(i)- 0.25, ddx )
193	rlonp025[i] = fx( np.float(i)+ 0.25, ddx )
194
195	xprimv[i] = fxprim(np.float(i), ddx )
196	xprimu[i] = fxprim(np.float(i)+ 0.5, ddx )
197	xprimm025[i] = fxprim(np.float(i)- 0.25, ddx )
198	xprimp025[i] = fxprim(np.float(i)+ 0.25, ddx )
199
200	return rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu, \
201	xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025
202
203	def fxhyp (dx, dy, xzoomdeg, grossism):
204	"""
205	c Auteur : P. Le Van
206
207
208	c Calcule les longitudes et derivees dans la grille du GCM pour une
209	c fonction f(x) a tangente hyperbolique .
210	c
211	c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
212	c dzoom etant la distance totale de la zone du zoom
213	c tau la raideur de la transition de l'interieur a l'exterieur du zoom
214	c
215	c On doit avoir grossism x dzoom < pi ( radians ) , en longitude.
216	c ********************************************************************
217	"""
218	fname = 'fxhyp'
219
220	nmax = 30000
221	nmax2 = 2*nmax
222
223	scal180 = True
224
225	# scal180 = .TRUE. si on veut avoir le premier point scalaire pour
226	# une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres.
227	# sinon scal180 = .FALSE.
228
229	# ...... arguments d'entree .......
230	#
231	# REAL xzoomdeg,dzooma,tau,grossism
232
233	# ...... arguments de sortie ......
234	rlonm025 = np.zeros((dx+1), dtype=np.float)
235	xprimm025 = np.zeros((dx+1), dtype=np.float)
236	rlonv = np.zeros((dx+1), dtype=np.float)
237	xprimv = np.zeros((dx+1), dtype=np.float)
238	rlonu = np.zeros((dx+1), dtype=np.float)
239	xprimu = np.zeros((dx+1), dtype=np.float)
240	rlonp025 = np.zeros((dx+1), dtype=np.float)
241	xprimp025 = np.zeros((dx+1), dtype=np.float)
242
243	# .... variables locales ....
244	#
245	xlon = np.zeros((dx+1), dtype=np.float32)
246	xprimm = np.zeros((dx+1), dtype=np.float32)
247	xtild = np.zeros((nmax2), dtype=np.float32)
248	fhyp = np.zeros((nmax2), dtype=np.float32)
249	Xprimt = np.zeros((nmax2), dtype=np.float32)
250	Xf = np.zeros((nmax2), dtype=np.float32)
251	xxpr = np.zeros((nmax2), dtype=np.float32)
252	xvrai = np.zeros((dx+1), dtype=np.float32)
253	xxprim = np.zeros((dx+1), dtype=np.float32)
254
255	pi = np.pi
256	depi = 2. * np.pi
257	epsilon = 1.e-3
258	xzoom = xzoomdeg * pi/180.
259
260	if dx == 1:
261	rlonm025[1] = -pi/2.
262	rlonv[1] = 0.
263	rlonu[1] = pi
264	rlonp025[1] = pi/2.
265	rlonm025[2] = rlonm025[1] + depi
266	rlonv[2] = rlonv[1] + depi
267	rlonu[2] = rlonu[1] + depi
268	rlonp025[2] = rlonp025[1] + depi
269
270	xprimm025 = 1.
271	xprimv = 1.
272	xprimu = 1.
273	xprimp025 = 1.
274	champmin =depi
275	champmax = depi
276	return dzooma, tau, rlonm025, xprimm025, rlonv, xprimv, rlonu, xprimu, \
277	rlonp025, xprimp025, champmin, champmax
278
279	decalx = .75
280	if grossism == 1. and scal180:
281	decalx = 1.
282
283	print 'FXHYP scal180,decalx', scal180,decalx
284
285	if dzooma > 1.:
286	dzoom = dzooma * depi
287	elif dzooma < 25.:
288	print erormsg
289	print ' ' + fname +": Le param. dzoomx pour 'fxhyp' est trop petit dzooma:",\
290	dzooma,'!L augmenter et relancer !'
291	quit(-1)
292	else:
293	dzoom = dzooma * pi/180.
294
295
296	print ' xzoom( rad.),grossism,tau,dzoom (radians)'
297	print xzoom,grossism,tau,dzoom
298
299	xtild = - pi + np.range(namx2) * depi /nmax2
300
301	for i in range(nmax,nmax2):
302	fa = tau* ( dzoom/2. - xtild[i] )
303	fb = xtild[i] * ( pi - xtild[i] )
304
305	if 200.* fb < - fa:
306	fhyp[i] = -1.
307	elif 200. * fb < fa:
308	fhyp[i] = 1.
309	else:
310	if np.abs(fa) < 1.e-13 and np.abs(fb) < 1.e-13:
311	if 200.*fb + fa < 1.e-10:
312	fhyp[i] = -1.
313	elif 200.*fb - fa < 1.e-10:
314	fhyp[i] = 1.
315	else:
316	fhyp[i] = np.tanh(fa/fb)
317	if xtild[i] == 0.: fhyp[i] = 1.
318	if xtild[i] == pi: fhyp[i] = -1.
319
320	## .... Calcul de beta ....
321
322	ffdx = 0.
323
324	for i in range(nmax+1,nmax2):
325	xmoy = 0.5 * ( xtild[i-1] + xtild[i] )
326	fa = tau*( dzoom/2. - xmoy )
327	fb = xmoy*( pi - xmoy )
328
329	if 200.* fb < -fa:
330	fxm = -1.
331	elif 200. * fb < fa:
332	fxm = 1.
333	else:
334	if np.abs(fa) < 1.e-13 and np.abs(fb) < 1.e-13:
335	if 200.*fb + fa < 1.e-10:
336	fxm = -1.
337	elif 200.*fb - fa < 1.e-10:
338	fxm = 1.
339	else:
340	fxm = np.tanh( fa/fb )
341
342	if xmoy == 0.: fxm = 1.
343	if xmoy == pi: fxm = -1.
344
345	ffdx = ffdx + fxm * ( xtild[i] - xtild[i-1] )
346
347	beta = ( grossism * ffdx - pi ) / ( ffdx - pi )
348
349	if 2.*beta - grossism < 0.:
350	print warnmsg
351	print ' ' + fname + ': ** Attention ! La valeur beta calculee dans la ' + \
352	'routine fxhyp est mauvaise ! '
353	print ' Modifier les valeurs de grossismx ,tau ou dzoomx et relancer ! ***'
354	quit(-1)
355
356	#
357	# ..... calcul de Xprimt .....
358	#
359
360	for i in range(nmax, nmax2):
361	Xprimt[i] = beta + ( grossism - beta ) * fhyp[i]
362
363	for i in range(nmax+1, nmax2):
364	Xprimt[nmax2-i] = Xprimt[i]
365
366	# ..... Calcul de Xf ........
367
368	Xf[0] = - pi
369
370	for i in range(nmax+1, nmax2):
371	xmoy = 0.5 * ( xtild[i-1] + xtild[i] )
372	fa = tau* ( dzoom/2. - xmoy )
373	fb = xmoy * ( pi - xmoy )
374
375	if 200.* fb < -fa:
376	fxm = -1.
377	elif 200. * fb < fa:
378	fxm = 1.
379	else:
380	fxm = np.tanh( fa/fb )
381
382	if xmoy == 0.: fxm = 1.
383	if xmoy == pi: fxm = -1.
384	xxpr[i] = beta + ( grossism - beta ) * fxm
385
386	for i in range(nmax+1, nmax2):
387	xxpr[nmax2-i+1] = xxpr[i]
388
389	for i in range(nmax2):
390	Xf[i] = Xf[i-1] + xxpr[i] * ( xtild[i] - xtild[i-1] )
391
392	# *****************************************************************
393	#
394
395	# ..... xuv = 0. si calcul aux pts scalaires ........
396	# ..... xuv = 0.5 si calcul aux pts U ........
397
398	for ik in range(4):
399	if ik == 0:
400	xuv = -0.25
401	elif ik == 1:
402	xuv = 0.
403	elif ik == 2:
404	xuv = 0.50
405	elif ik == 4:
406	xuv = 0.25
407
408	xo1 = 0.
409
410	ii1=1
411	ii2=dx
412	if ik == 0 and grossism == 1.:
413	ii1 = 2
414	ii2 = dx+1
415
416	for i in range(ii1, ii2):
417	xlon2 = - pi + (np.float(i) + xuv - decalx) * depi / np.float32(dx)
418
419	Xfi = xlon2
420
421	ended = False
422	for it in range(nmax2,0,-1):
423	if Xfi > Xf[it]:
424	ended = True
425	break
426	if not ended: it = 0
427
428	# ...... Calcul de Xf(xi) ......
429	#
430	xi = xtild[it]
431
432	if it == nmax2:
433	it = nmax2 -1
434	Xf[it+1] = pi
435
436	# .....................................................................
437	#
438	# Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un
439	# polynome de degre 3 qui passe par les points (Xf(it),xtild(it) )
440	# et (Xf(it+1),xtild(it+1) )
441
442	a0, a1, a2, a3 = coefpoly( Xf[it], Xf[it+1], Xprimt[it], Xprimt[it+1], \
443	xtild[it], xtild[it+1])
444
445	Xf1 = Xf[it]
446	Xprimin = a1 + 2.* a2 * xi + 3.a3 xi *xi
447
448	for iteri in range(300):
449	xi = xi - ( Xf1 - Xfi )/ Xprimin
450
451	ended = False
452	if np.abs(xi-xo1) < epsilon:
453	ended = True
454	break
455	xo1 = xi
456	xi2 = xi * xi
457	Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi
458	Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2
459
460	if not ended:
461	print errmsg
462	print ' ' + fname + ': Pas de solution ***** ',i,xlon2,iteri
463	quit(-1)
464
465	xxprim[i] = depi/ ( np.float32(dx) * Xprimin )
466	xvrai[i] = xi + xzoom
467
468	if ik == 0 and grossism == 1:
469	xvrai[0] = xvrai[dx+1] - depi
470	xxprim[0] = xxprim[dx+1]
471
472	xlon[0:dx+1] = xvrai[0:dx+1]
473	xprimm[0:dx+1] = xxprim[0:dx+1]
474
475	for i in range(dx-1):
476	if xvrai[i+1] < xvrai[i]:
477	print errormsg
478	print ' ' + fname + ': PBS. avec rlonu(',i+1,') plus petit que rlonu(', \
479	i,')'
480	quit(-1)
481
482	#
483	# ... Reorganisation des longitudes pour les avoir entre - pi et pi ..
484	# ........................................................................
485
486	champmin = 1.e12
487	champmax = -1.e12
488	for i in range(dx):
489	champmin = np.min( [champmin,xvrai[i]] )
490	champmax = np.max( [champmax,xvrai[i]] )
491
492	# if champmin >= -pi-0.10 and champmax <= pi+0.10:
493	# GO TO 1600
494
495	#
496	# HERE -- here
497	#
498	# ELSE
499	# WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
500	# , ' et pi '
501	#c
502	# IF( xzoom.LE.0.) THEN
503	# IF( ik.EQ. 1 ) THEN
504	# DO i = 1, iim
505	# IF( xvrai(i).GE. - pi ) GO TO 80
506	# ENDDO
507	# WRITE(6,*) ' PBS. 1 ! Xvrai plus petit que - pi ! '
508	# STOP 8
509	# 80 CONTINUE
510	# is2 = i
511	# ENDIF
512	#
513	# IF( is2.NE. 1 ) THEN
514	# DO ii = is2 , iim
515	# xlon (ii-is2+1) = xvrai(ii)
516	# xprimm(ii-is2+1) = xxprim(ii)
517	# ENDDO
518	# DO ii = 1 , is2 -1
519	# xlon (ii+iim-is2+1) = xvrai(ii) + depi
520	# xprimm(ii+iim-is2+1) = xxprim(ii)
521	# ENDDO
522	# ENDIF
523	# ELSE
524	# IF( ik.EQ.1 ) THEN
525	# DO i = iim,1,-1
526	# IF( xvrai(i).LE. pi ) GO TO 90
527	# ENDDO
528	# WRITE(6,*) ' PBS. 2 ! Xvrai plus grand que pi ! '
529	# STOP 9
530	# 90 CONTINUE
531	# is2 = i
532	# ENDIF
533	# idif = iim -is2
534	# DO ii = 1, is2
535	# xlon (ii+idif) = xvrai(ii)
536	# xprimm(ii+idif) = xxprim(ii)
537	# ENDDO
538	# DO ii = 1, idif
539	# xlon (ii) = xvrai (ii+is2) - depi
540	# xprimm(ii) = xxprim(ii+is2)
541	# ENDDO
542	# ENDIF
543	# ENDIF
544	#c
545	#c ......... Fin de la reorganisation ............................
546
547	# 1600 CONTINUE
548
549
550	# xlon ( iip1) = xlon(1) + depi
551	# xprimm( iip1 ) = xprimm (1 )
552
553	# DO i = 1, iim+1
554	# xvrai(i) = xlon(i)*180./pi
555	# ENDDO
556
557	# IF( ik.EQ.1 ) THEN
558	#c WRITE(6,*) ' XLON aux pts. V-0.25 apres ( en deg. ) '
559	#c WRITE(6,18)
560	#c WRITE(6,68) xvrai
561	#c WRITE(6,*) ' XPRIM k ',ik
562	#c WRITE(6,566) xprimm
563
564	# DO i = 1,iim +1
565	# rlonm025(i) = xlon( i )
566	# xprimm025(i) = xprimm(i)
567	# ENDDO
568	# ELSE IF( ik.EQ.2 ) THEN
569	#c WRITE(6,18)
570	#c WRITE(6,*) ' XLON aux pts. V apres ( en deg. ) '
571	#c WRITE(6,68) xvrai
572	#c WRITE(6,*) ' XPRIM k ',ik
573	#c WRITE(6,566) xprimm
574
575	# DO i = 1,iim + 1
576	# rlonv(i) = xlon( i )
577	# xprimv(i) = xprimm(i)
578	# ENDDO
579
580	# ELSE IF( ik.EQ.3) THEN
581	#c WRITE(6,18)
582	#c WRITE(6,*) ' XLON aux pts. U apres ( en deg. ) '
583	#c WRITE(6,68) xvrai
584	#c WRITE(6,*) ' XPRIM ik ',ik
585	#c WRITE(6,566) xprimm
586
587	# DO i = 1,iim + 1
588	# rlonu(i) = xlon( i )
589	# xprimu(i) = xprimm(i)
590	# ENDDO
591
592	# ELSE IF( ik.EQ.4 ) THEN
593	#c WRITE(6,18)
594	#c WRITE(6,*) ' XLON aux pts. V+0.25 apres ( en deg. ) '
595	#c WRITE(6,68) xvrai
596	#c WRITE(6,*) ' XPRIM ik ',ik
597	#c WRITE(6,566) xprimm
598
599	# DO i = 1,iim + 1
600	# rlonp025(i) = xlon( i )
601	# xprimp025(i) = xprimm(i)
602	# ENDDO
603
604	# ENDIF
605
606	#5000 CONTINUE
607	#c
608	# WRITE(6,18)
609	#c
610	#c ........... fin de la boucle do 5000 ............
611
612	# DO i = 1, iim
613	# xlon(i) = rlonv(i+1) - rlonv(i)
614	# ENDDO
615	# champmin = 1.e12
616	# champmax = -1.e12
617	# DO i = 1, iim
618	# champmin = MIN( champmin, xlon(i) )
619	# champmax = MAX( champmax, xlon(i) )
620	# ENDDO
621	# champmin = champmin * 180./pi
622	# champmax = champmax * 180./pi
623
624	#18 FORMAT(/)
625	#24 FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
626	#68 FORMAT(1x,7f9.2)
627	#566 FORMAT(1x,7f9.4)
628
629	return dzooma, tau, rlonm025, xprimm025, rlonv, xprimv, rlonu, xprimu, \
630	rlonp025, xprimp025, champmin, champmax
631
632
633	def fxyhyper (ddy, ddx, yzoom, grossy, dzoomy, tauy, xzoom, grossx, dzoomx, taux, \
634	rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2,
635	rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025):
636	"""c
637	c Auteur : P. Le Van .
638	c
639	c d'apres formulations de R. Sadourny .
640	c
641	c
642	c Ce spg calcule les latitudes( routine fyhyp ) et longitudes( fxhyp )
643	c par des fonctions a tangente hyperbolique .
644	c
645	c Il y a 3 parametres ,en plus des coordonnees du centre du zoom (xzoom
646	c et yzoom ) :
647	c
648	c a) le grossissement du zoom : grossy ( en y ) et grossx ( en x )
649	c b) l' extension du zoom : dzoomy ( en y ) et dzoomx ( en x )
650	c c) la raideur de la transition du zoom : taux et tauy
651	c
652	c N.B : Il vaut mieux avoir : grossx * dzoomx < pi ( radians )
653	c ******
654	c et grossy * dzoomy < pi/2 ( radians )
655	c
656	"""
657	fname = 'fxyhyper'
658
659	# CALL fyhyp ( yzoom, grossy, dzoomy,tauy ,
660	# , rlatu, yprimu,rlatv,yprimv,rlatu2,yprimu2,rlatu1,yprimu1 ,
661	# , dymin,dymax )
662
663	# CALL fxhyp(xzoom,grossx,dzoomx,taux,rlonm025,xprimm025,rlonv,
664	# , xprimv,rlonu,xprimu,rlonp025,xprimp025 , dxmin,dxmax )
665
666
667	for i in range(ddx+1):
668	if rlonp025[i] < rlonv[i]:
669	print errormsg
670	print ' ' + fname + ' Attention ! rlonp025 < rlonv',i
671	quit(-1)
672
673	if rlonv[i] < rlonm025[i]:
674	print errormsg
675	print ' ' + fname + ' Attention ! rlonm025 > rlonv',i
676	quit(-1)
677
678	if rlonp025[i] > rlonu[i]:
679	print errormsg
680	print ' ' + fname + ' Attention ! rlonp025 > rlonu',i
681	quit(-1)
682
683	print ' * TEST DE COHERENCE OK POUR FX ** '
684
685
686	for j in range(ddy):
687	if rlatu1[j] <= rlatu2[j]:
688	print errormsg
689	print ' ' + fname + 'Attention ! rlatu1 < rlatu2',rlatu1[j],rlatu2[j],j
690	quit(-1)
691
692	if rlatu2[j] <= rlatu[j+1]:
693	print errormsg
694	print ' ' + fname + 'Attention ! rlatu2 < rlatup1',rlatu2[j],rlatu[j+1],j
695	quit(-1)
696
697	if rlatu[j] <= rlatu1[j]:
698	print errormsg
699	print ' ' + fname + 'Attention ! rlatu < rlatu1',rlatu[j],rlatu1[j],j
700	quit(-1)
701
702	if rlatv[j] <= rlatu2[j]:
703	print errormsg
704	print ' ' + fname + 'Attention ! rlatv < rlatu2 ',rlatv[j],rlatu2[j],j
705	quit(-1)
706
707	if rlatv[j] >= rlatu1[j]:
708	print errormsg
709	print ' ' + fname + 'Attention ! rlatv > rlatu1 ',rlatv[j],rlatu1[j],j
710	quit(-1)
711
712	if rlatv[j] >= rlatu[j]:
713	print errormsg
714	print ' ' + fname + 'Attention ! rlatv > rlatu ',rlatv[j],rlatu[j],j
715	quit(-1)
716
717	print ' * TEST DE COHERENCE OK POUR FY ** '
718
719	print ' Latitudes '
720	print ' *********** '
721
722	print 'Au centre du zoom, la longueur de la maille est d environ', dymin , \
723	'degres alors que la maille en dehors de la zone du zoom est d environ', \
724	dymax, 'degres'
725	print ' Si cette derniere est trop lache , modifiez les parametres grossism , ' +\
726	'tau , dzoom pour Y et repasser ! '
727
728	print ' Longitudes '
729	print ' ************ '
730	print 'Au centre du zoom, la longueur de la maille est d environ', dxmin , \
731	'degres alors que la maille en dehors de la zone du zoom est d environ', \
732	dxmax, 'degres'
733	print ' Si cette derniere est trop lache , modifiez les parametres grossism , ' +\
734	'tau , dzoom pour Y et repasser ! '
735
736	return
737
738	def inigeom(dy,dx):
739	"""c
740	c Auteur : P. Le Van
741	c
742	c ............ Version du 01/04/2001 ........................
743	c
744	c Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4 aux memes en-
745	c endroits que les aires aireij1,..aireij4 .
746
747	c Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol.
748	c
749	c
750	"""
751	fname = 'inigeom'
752
753	cvu = np.zeros((dy+1,dx+1), dtype=np.float)
754	cuv = np.zeros((dy, dx+1), dtype=np.float)
755
756	cuij1 = np.zeros((dy+1,dx+1), dtype=np.float)
757	cuij2 = np.zeros((dy+1,dx+1), dtype=np.float)
758	cuij3 = np.zeros((dy+1,dx+1), dtype=np.float)
759	cuij4 = np.zeros((dy+1,dx+1), dtype=np.float)
760	cvij1 = np.zeros((dy+1,dx+1), dtype=np.float)
761	cvij2 = np.zeros((dy+1,dx+1), dtype=np.float)
762	cvij3 = np.zeros((dy+1,dx+1), dtype=np.float)
763	cvij4 = np.zeros((dy+1,dx+1), dtype=np.float)
764	aireij1 = np.zeros((dy+1,dx+1), dtype=np.float)
765	aireij2 = np.zeros((dy+1,dx+1), dtype=np.float)
766	aireij3 = np.zeros((dy+1,dx+1), dtype=np.float)
767	aireij4 = np.zeros((dy+1,dx+1), dtype=np.float)
768
769	aire = np.zeros((dy+1,dx+1), dtype=np.float)
770	aireu = np.zeros((dy+1,dx+1), dtype=np.float)
771	airev = np.zeros((dy+1,dx+1), dtype=np.float)
772	unsaire = np.zeros((dy+1,dx+1), dtype=np.float)
773	unsair_gam1 = np.zeros((dy+1,dx+1), dtype=np.float)
774	unsair_gam2 = np.zeros((dy+1,dx+1), dtype=np.float)
775	airesurg = np.zeros((dy+1,dx+1), dtype=np.float)
776	unsairez = np.zeros((dy+1,dx+1), dtype=np.float)
777	unsairz_gam = np.zeros((dy+1,dx+1), dtype=np.float)
778	fext = np.zeros((dy+1,dx+1), dtype=np.float)
779
780	alpha1 = np.zeros((dy+1,dx+1), dtype=np.float)
781	alpha2 = np.zeros((dy+1,dx+1), dtype=np.float)
782	alpha3 = np.zeros((dy+1,dx+1), dtype=np.float)
783	alpha4 = np.zeros((dy+1,dx+1), dtype=np.float)
784	alpha1p2 = np.zeros((dy+1,dx+1), dtype=np.float)
785	alpha1p4 = np.zeros((dy+1,dx+1), dtype=np.float)
786	alpha2p3 = np.zeros((dy+1,dx+1), dtype=np.float)
787	alpha3p4 = np.zeros((dy+1,dx+1), dtype=np.float)
788
789	rlonvv = np.zeros((dx+1), dtype=np.float)
790	rlatuu = np.zeros((dy+1), dtype=np.float)
791	rlatu1 = np.zeros((dy), dtype=np.float)
792	yprimu1 = np.zeros((dy), dtype=np.float)
793	rlatu2 = np.zeros((dy), dtype=np.float)
794	yprimu2 = np.zeros((dy), dtype=np.float)
795	yprimv = np.zeros((dy), dtype=np.float)
796	yprimu = np.zeros((dy+1), dtype=np.float)
797
798	rlonm025 = np.zeros((dx+1), dtype=np.float)
799	xprimm025 = np.zeros((dx+1), dtype=np.float)
800	rlonp025 = np.zeros((dx+1), dtype=np.float)
801	xprimp025 = np.zeros((dx+1), dtype=np.float)
802
803	cu = np.zeros((dy+1,dx+1), dtype=np.float)
804	cv = np.zeros((dy+1,dx+1), dtype=np.float)
805	unscu2 = np.zeros((dy+1,dx+1), dtype=np.float)
806	unscv2 = np.zeros((dy+1,dx+1), dtype=np.float)
807	cuvsurcv = np.zeros((dy+1,dx+1), dtype=np.float)
808	cuvsurcuv = np.zeros((dy+1,dx+1), dtype=np.float)
809	cvsurcv = np.zeros((dy+1,dx+1), dtype=np.float)
810	cvsurcuv = np.zeros((dy+1,dx+1), dtype=np.float)
811	cuvscvgam1 = np.zeros((dy+1,dx+1), dtype=np.float)
812	cuvscvgam2 = np.zeros((dy+1,dx+1), dtype=np.float)
813	cvscuvgam = np.zeros((dy+1,dx+1), dtype=np.float)
814	cvusurcu = np.zeros((dy+1,dx+1), dtype=np.float)
815	cusurcvu = np.zeros((dy+1,dx+1), dtype=np.float)
816	cvuscugam1 = np.zeros((dy+1,dx+1), dtype=np.float)
817	cvuscugam2 = np.zeros((dy+1,dx+1), dtype=np.float)
818	cuscvugam = np.zeros((dy+1,dx+1), dtype=np.float)
819	airvscu2 = np.zeros((dy+1,dx+1), dtype=np.float)
820	aivscu2gam = np.zeros((dy+1,dx+1), dtype=np.float)
821	airuscv2 = np.zeros((dy+1,dx+1), dtype=np.float)
822	aiuscv2gam = np.zeros((dy+1,dx+1), dtype=np.float)
823
824	constang = np.zeros((dy+1,dx+1), dtype=np.float)
825	#
826	#
827	# ------------------------------------------------------------------
828	# - -
829	# - calcul des coeff. ( cu, cv , 1./cu2, 1./cv2 ) -
830	# - -
831	# ------------------------------------------------------------------
832	#
833	# les coef. ( cu, cv ) permettent de passer des vitesses naturelles
834	# aux vitesses covariantes et contravariantes , ou vice-versa ...
835	#
836	#
837	# on a : u (covariant) = cu * u (naturel) , u(contrav)= u(nat)/cu
838	# v (covariant) = cv * v (naturel) , v(contrav)= v(nat)/cv
839	#
840	# on en tire : u(covariant) = cu * cu * u(contravariant)
841	# v(covariant) = cv * cv * v(contravariant)
842	#
843	#
844	# on a l'application ( x(X) , y(Y) ) avec - im/2 +1 < X < im/2
845	# = =
846	# et - jm/2 < Y < jm/2
847	# = =
848	#
849	# ...................................................
850	# ...................................................
851	# . x est la longitude du point en radians .
852	# . y est la latitude du point en radians .
853	# . .
854	# . on a : cu(i,j) = rad * COS(y) * dx/dX .
855	# . cv( j ) = rad * dy/dY .
856	# . aire(i,j) = cu(i,j) * cv(j) .
857	# . .
858	# . y, dx/dX, dy/dY calcules aux points concernes .
859	# . .
860	# ...................................................
861	# ...................................................
862	#
863	#
864	#
865	# ,
866	# cv , bien que dependant de j uniquement,sera ici indice aussi en i
867	# pour un adressage plus facile en ij .
868	#
869	#
870	#
871	# ************ aux points u et v , ***************
872	# xprimu et xprimv sont respectivement les valeurs de dx/dX
873	# yprimu et yprimv . . . . . . . . . . . dy/dY
874	# rlatu et rlatv . . . . . . . . . . .la latitude
875	# cvu et cv . . . . . . . . . . . cv
876	#
877	# ************ aux points u, v, scalaires, et z **************
878	# cu, cuv, cuscal, cuz sont respectiv. les valeurs de cu
879	#
880	#
881	#
882	# Exemple de distribution de variables sur la grille dans le
883	# domaine de travail ( X,Y ) .
884	# ................................................................
885	# DX=DY= 1
886	#
887	#
888	# + represente un point scalaire ( p.exp la pression )
889	# > represente la composante zonale du vent
890	# V represente la composante meridienne du vent
891	# o represente la vorticite
892	#
893	# ---- , car aux poles , les comp.zonales covariantes sont nulles
894	#
895	#
896	#
897	# i ->
898	#
899	# 1 2 3 4 5 6 7 8
900	# j
901	# v 1 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
902	#
903	# V o V o V o V o V o V o V o V o
904	#
905	# 2 + > + > + > + > + > + > + > + >
906	#
907	# V o V o V o V o V o V o V o V o
908	#
909	# 3 + > + > + > + > + > + > + > + >
910	#
911	# V o V o V o V o V o V o V o V o
912	#
913	# 4 + > + > + > + > + > + > + > + >
914	#
915	# V o V o V o V o V o V o V o V o
916	#
917	# 5 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + --
918	#
919	#
920	# Ci-dessus, on voit que le nombre de pts.en longitude est egal
921	# a IM = 8
922	# De meme , le nombre d'intervalles entre les 2 poles est egal
923	# a JM = 4
924	#
925	# Les points scalaires ( + ) correspondent donc a des valeurs
926	# entieres de i ( 1 a IM ) et de j ( 1 a JM +1 ) .
927	#
928	# Les vents U ( > ) correspondent a des valeurs semi-
929	# entieres de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1)
930	#
931	# Les vents V ( V ) correspondent a des valeurs entieres
932	# de i ( 1 a IM ) et semi-entieres de j ( 1 +0.5 a JM +0.5)
933	#
934	#
935	#
936
937	if nitergdiv != 2:
938	gamdi_gdiv = coefdis/ ( np.float(nitergdiv) -2. )
939	else:
940	gamdi_gdiv = 0.
941
942	if nitergrot != 2:
943	gamdi_grot = coefdis/ ( np.float(nitergrot) -2. )
944	else:
945	gamdi_grot = 0.
946
947	if niterh != 2:
948	gamdi_h = coefdis/ ( np.float(niterh) -2. )
949	else:
950	gamdi_h = 0.
951
952	print 'gamdi_gd:',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,nitergdiv,nitergrot,niterh
953
954	# ----------------------------------------------------------------
955	#
956	if not fxyhypb:
957	if ysinus:
958	print ' * Inigeom , Y = Sinus ( Latitude ) * '
959
960	# .... utilisation de f(x,y ) avec y = sinus de la latitude .....
961
962	rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu, \
963	xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025 = \
964	fxysinus(dx, dy)
965	else:
966	print '* Inigeom , Y = Latitude , der. sinusoid . *'
967
968	# .... utilisation de f(x,y) a tangente sinusoidale , y etant la latit. ...
969
970	pxo = clon *np.pi /180.
971	pyo = 2.* clat* np.pi /180.
972	#
973	# .... determination de transx ( pour le zoom ) par Newton-Raphson ...
974	#
975	itmax = 10
976	eps = .1e-7
977
978	xo1 = 0.
979	for iter in range(itmax):
980	x1 = xo1
981	f = x1+ alphax*np.sin(x1-pxo)
982	df = 1.+ alphax*np.cos(x1-pxo)
983	x1 = x1 - f/df
984	xdm = np.abs( x1- xo1 )
985	if xdm > eps: xo1 = x1
986
987	transx = xo1
988
989	itmay = 10
990	eps = .1e-7
991
992	yo1 = 0.
993	for iter in range(itmay):
994	y1 = yo1
995	f = y1 + alphay*np.sin(y1-pyo)
996	df = 1. + alphay*np.cos(y1-pyo)
997	y1 = y1 -f/df
998	ydm = np.abs(y1-yo1)
999	if ydm > eps: yo1 = y1
1000
1001	transy = yo1
1002
1003	rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, yprimu2, rlonu, \
1004	xprimu, rlonv, xprimv, rlonm025, xprimm025, rlonp025, xprimp025 = \
1005	fxy(dx,dy)
1006	else:
1007	#
1008	# .... Utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol.
1009	# .....................................................................
1010
1011	print '* Inigeom , Y = Latitude , der.tg. hyperbolique *'
1012
1013	#!!!!! Lluis
1014	#!! HERE, how to do all this without zoom!?
1015	#!
1016
1017	# CALL fxyhyper( clat, grossismy, dzoomy, tauy ,
1018	# , clon, grossismx, dzoomx, taux ,
1019	# , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2 ,
1020	# , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 )
1021
1022	# -------------------------------------------------------------------
1023
1024	rlatu[0] = np.arcsin(1.)
1025	rlatu[dy] = -rlatu[0]
1026
1027	# .... calcul aux poles ....
1028
1029	yprimu[0] = 0.
1030	yprimu[dy] = 0.
1031
1032	un4rad2 = 0.25 * rad * rad
1033
1034	#
1035	# --------------------------------------------------------------------
1036	# --------------------------------------------------------------------
1037	# - -
1038	# - calcul des aires ( aire,aireu,airev, 1./aire, 1./airez ) -
1039	# - et de fext , force de coriolis extensive . -
1040	# - -
1041	# --------------------------------------------------------------------
1042	# --------------------------------------------------------------------
1043	#
1044	#
1045	#
1046	# A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont
1047	# affectees 4 aires entourant P , calculees respectivement aux points
1048	# ( i + 1/4, j - 1/4 ) : aireij1 (i,j)
1049	# ( i + 1/4, j + 1/4 ) : aireij2 (i,j)
1050	# ( i - 1/4, j + 1/4 ) : aireij3 (i,j)
1051	# ( i - 1/4, j - 1/4 ) : aireij4 (i,j)
1052	#
1053	# ,
1054	# Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X,Y).
1055	# Chaque aire centree en 1 point scalaire P(i,j) est egale a la somme
1056	# des 4 aires aireij1,aireij2,aireij3,aireij4 qui sont affectees au
1057	# point (i,j) .
1058	# On definit en outre les coefficients alpha comme etant egaux a
1059	# (aireij / aire), c.a.d par exp. alpha1(i,j)=aireij1(i,j)/aire(i,j)
1060	#
1061	# De meme, toute aire centree en 1 point U est egale a la somme des
1062	# 4 aires aireij1,aireij2,aireij3,aireij4 entourant le point U .
1063	# Idem pour airev, airez .
1064	#
1065	# On a ,pour chaque maille : dX = dY = 1
1066	#
1067	#
1068	# . V
1069	#
1070	# aireij4 . . aireij1
1071	#
1072	# U . . P . U
1073	#
1074	# aireij3 . . aireij2
1075	#
1076	# . V
1077	#
1078	#
1079	#
1080	#
1081	#
1082	# ....................................................................
1083	#
1084	# Calcul des 4 aires elementaires aireij1,aireij2,aireij3,aireij4
1085	# qui entourent chaque aire(i,j) , ainsi que les 4 elongations elemen
1086	# taires cuij et les 4 elongat. cvij qui sont calculees aux memes
1087	# endroits que les aireij .
1088	#
1089	# ....................................................................
1090	#
1091	# ....... do 35 : boucle sur les jjm + 1 latitudes .....
1092	#
1093	#
1094	for j in range(dy+1):
1095	if j == 1:
1096	yprm = yprimu1[j]
1097	rlatm = rlatu1[j]
1098
1099	coslatm = np.cos( rlatm )
1100	radclatm = 0.5* rad * coslatm
1101
1102	for i in range(dx):
1103	xprp = xprimp025[i]
1104	xprm = xprimm025[i]
1105	aireij2[0,i] = un4rad2 * coslatm * xprp * yprm
1106	aireij3[0,i] = un4rad2 * coslatm * xprm * yprm
1107	cuij2[0,i] = radclatm * xprp
1108	cuij3[0,i] = radclatm * xprm
1109	cvij2[0,i] = 0.5* rad * yprm
1110	cvij3[0,i] = cvij2[0,i]
1111
1112	for i in range(dx):
1113	aireij1[0,i] = 0.
1114	aireij4[0,i] = 0.
1115	cuij1[0,i] = 0.
1116	cuij4[0,i] = 0.
1117	cvij1[0,i] = 0.
1118	cvij4[0,i] = 0.
1119
1120	elif j == dy:
1121	yprp = yprimu2[j-1]
1122	rlatp = rlatu2[j-1]
1123
1124	coslatp = np.cos( rlatp )
1125	radclatp = 0.5* rad * coslatp
1126
1127	for i in range(dx):
1128	xprp = xprimp025[i]
1129	xprm = xprimm025[i]
1130	aireij1[dy,i] = un4rad2 * coslatp * xprp * yprp
1131	aireij4[dy,i] = un4rad2 * coslatp * xprm * yprp
1132	cuij1[dy,i] = radclatp * xprp
1133	cuij4[dy,i] = radclatp * xprm
1134	cvij1[dy,i] = 0.5 * rad* yprp
1135	cvij4[dy,i] = cvij1[dy,i]
1136
1137	for i in range(dx):
1138	aireij2[dy,i] = 0.
1139	aireij3[dy,i] = 0.
1140	cvij2[dy,i] = 0.
1141	cvij3[dy,i] = 0.
1142	cuij2[dy,i] = 0.
1143	cuij3[dy,i] = 0.
1144
1145	else:
1146
1147	rlatp = rlatu2[j-1]
1148	yprp = yprimu2[j-1]
1149	rlatm = rlatu1[j]
1150	yprm = yprimu1[j]
1151
1152	coslatm = np.cos( rlatm )
1153	coslatp = np.cos( rlatp )
1154	radclatp = 0.5* rad * coslatp
1155	radclatm = 0.5* rad * coslatm
1156
1157	for i in range(dx):
1158	xprp = xprimp025[i]
1159	xprm = xprimm025[i]
1160
1161	ai14 = un4rad2 * coslatp * yprp
1162	ai23 = un4rad2 * coslatm * yprm
1163	aireij1[j,i] = ai14 * xprp
1164	aireij2[j,i] = ai23 * xprp
1165	aireij3[j,i] = ai23 * xprm
1166	aireij4[j,i] = ai14 * xprm
1167	cuij1[j,i] = radclatp * xprp
1168	cuij2[j,i] = radclatm * xprp
1169	cuij3[j,i] = radclatm * xprm
1170	cuij4[j,i] = radclatp * xprm
1171	cvij1[j,i] = 0.5* rad * yprp
1172	cvij2[j,i] = 0.5* rad * yprm
1173	cvij3[j,i] = cvij2[j,i]
1174	cvij4[j,i] = cvij1[j,i]
1175
1176	#
1177	# ........ periodicite ............
1178	#
1179	cvij1[j,dx] = cvij1[j,0]
1180	cvij2[j,dx] = cvij2[j,0]
1181	cvij3[j,dx] = cvij3[j,0]
1182	cvij4[j,dx] = cvij4[j,0]
1183	cuij1[j,dx] = cuij1[j,0]
1184	cuij2[j,dx] = cuij2[j,0]
1185	cuij3[j,dx] = cuij3[j,0]
1186	cuij4[j,dx] = cuij4[j,0]
1187	aireij1[j,dx] = aireij1[j,0]
1188	aireij2[j,dx] = aireij2[j,0]
1189	aireij3[j,dx] = aireij3[j,0]
1190	aireij4[j,dx] = aireij4[j,0]
1191
1192	#
1193	# ..............................................................
1194	#
1195	for j in range(dy+1):
1196	for i in range(dx):
1197	aire[j,i] = aireij1[j,i] + aireij2[j,i] + aireij3[j,i] + aireij4[j,i]
1198	alpha1[j,i] = aireij1[j,i] / aire[j,i]
1199	alpha2[j,i] = aireij2[j,i] / aire[j,i]
1200	alpha3[j,i] = aireij3[j,i] / aire[j,i]
1201	alpha4[j,i] = aireij4[j,i] / aire[j,i]
1202	alpha1p2[j,i] = alpha1 [j,i] + alpha2 [j,i]
1203	alpha1p4[j,i] = alpha1 [j,i] + alpha4 [j,i]
1204	alpha2p3[j,i] = alpha2 [j,i] + alpha3 [j,i]
1205	alpha3p4[j,i] = alpha3 [j,i] + alpha4 [j,i]
1206
1207	aire[j,dx] = aire[j,0]
1208	alpha1[j,dx] = alpha1[j,0]
1209	alpha2[j,dx] = alpha2[j,0]
1210	alpha3[j,dx] = alpha3[j,0]
1211	alpha4[j,dx] = alpha4[j,0]
1212	alpha1p2[j,dx] = alpha1p2[j,0]
1213	alpha1p4[j,dx] = alpha1p4[j,0]
1214	alpha2p3[j,dx] = alpha2p3[j,0]
1215	alpha3p4[j,dx] = alpha3p4[j,0]
1216
1217	for j in range(dy+1):
1218	for i in range(dx):
1219	aireu[j,i] = aireij1[j,i] + aireij2[j,i] + aireij4[j,i+1] + aireij3[j,i+1]
1220	unsaire[j,i] = 1./ aire[j,i]
1221	unsair_gam1[j,i] = unsaire[j,i]** ( - gamdi_gdiv )
1222	unsair_gam2[j,i] = unsaire[j,i]** ( - gamdi_h )
1223	airesurg[j,i] = aire[j,i]/ g
1224
1225	aireu[j,dx] = aireu[j,0]
1226	unsaire[j,dx] = unsaire[j,0]
1227	unsair_gam1[j,dx] = unsair_gam1[j,0]
1228	unsair_gam2[j,dx] = unsair_gam2[j,0]
1229	airesurg[j,dx] = airesurg[j,0]
1230
1231
1232	for j in range(dy):
1233	for i in range(dx):
1234	airev[j,i] = aireij2[j,i]+ aireij3[j,i]+ aireij1[j,i] + aireij4[j+1,i]
1235
1236	for i in range(dx):
1237	airez = aireij2[j,i]+aireij1[j+1,i]+aireij3[j,i+1] + aireij4[j+1,i+1]
1238	unsairez[j,i] = 1./ airez
1239	unsairz_gam[j,i]= unsairez[j,i]** ( - gamdi_grot )
1240	fext[j,i] = airez * np.sin(rlatv[j])* 2.* omeg
1241
1242	airev[j,dx] = airev[j,0]
1243	unsairez[j,dx] = unsairez[j,0]
1244	fext[j,dx] = fext[j,0]
1245	unsairz_gam[j,dx] = unsairz_gam[j,0]
1246
1247
1248	#
1249	# ..... Calcul des elongations cu,cv, cvu .........
1250	#
1251	for j in range(dy):
1252	for i in range(dx):
1253	cv[j,i] = 0.5 *( cvij2[j,i]+cvij3[j,i]+cvij1[j+1,i]+cvij4[j+1,i])
1254	cvu[j,i]= 0.5 *( cvij1[j,i]+cvij4[j,i]+cvij2[j,i]+cvij3[j,i] )
1255	cuv[j,i]= 0.5 *( cuij2[j,i]+cuij3[j,i]+cuij1[j+1,i]+cuij4[j+1,i])
1256	unscv2[j,i] = 1./ ( cv[j,i]*cv[j,i] )
1257
1258	for i in range(dx):
1259	cuvsurcv [j,i] = airev[j,i] * unscv2[j,i]
1260	cvsurcuv [j,i] = 1./cuvsurcv[j,i]
1261	cuvscvgam1[j,i] = cuvsurcv [j,i] ** ( - gamdi_gdiv )
1262	cuvscvgam2[j,i] = cuvsurcv [j,i] ** ( - gamdi_h )
1263	cvscuvgam[j,i] = cvsurcuv [j,i] ** ( - gamdi_grot )
1264
1265	cv[j,dx] = cv[j,0]
1266	cvu[j,dx] = cvu[j,0]
1267	unscv2[j,dx] = unscv2[j,0]
1268	cuv[j,dx] = cuv[j,0]
1269	cuvsurcv[j,dx] = cuvsurcv[j,0]
1270	cvsurcuv[j,dx] = cvsurcuv[j,0]
1271	cuvscvgam1[j,dx] = cuvscvgam1[j,0]
1272	cuvscvgam2[j,dx] = cuvscvgam2[j,0]
1273	cvscuvgam[j,dx] = cvscuvgam[j,0]
1274
1275
1276	for j in range(1,dy):
1277	for i in range(dx):
1278	cu[j,i] = 0.5*(cuij1[j,i]+cuij4[j,i+1]+cuij2[j,i]+cuij3[j,i+1])
1279	unscu2[j,i] = 1./ ( cu[j,i] * cu[j,i] )
1280	cvusurcu[j,i] = aireu[j,i] * unscu2[j,i]
1281	cusurcvu[j,i] = 1./ cvusurcu[j,i]
1282	cvuscugam1[j,i] = cvusurcu[j,i] ** ( - gamdi_gdiv )
1283	cvuscugam2[j,i] = cvusurcu[j,i] ** ( - gamdi_h )
1284	cuscvugam[j,i] = cusurcvu[j,i] ** ( - gamdi_grot )
1285
1286	cu[j,dx] = cu[j,0]
1287	unscu2[j,dx] = unscu2[j,0]
1288	cvusurcu[j,dx] = cvusurcu[j,0]
1289	cusurcvu[j,dx] = cusurcvu[j,0]
1290	cvuscugam1[j,dx] = cvuscugam1[j,0]
1291	cvuscugam2[j,dx] = cvuscugam2[j,0]
1292	cuscvugam[j,dx] = cuscvugam[j,0]
1293
1294	#
1295	# .... calcul aux poles ....
1296	#
1297	for i in range(dx+1):
1298	cu[0, i] = 0.
1299	unscu2[0, i] = 0.
1300	cvu[0, i] = 0.
1301
1302	cu[dy,i] = 0.
1303	unscu2[dy,i] = 0.
1304	cvu[dy,i] = 0.
1305
1306	#
1307	# ..............................................................
1308	#
1309	for j in range(dy):
1310	for i in range(dx):
1311	airvscu2[j,i] = airev[j,i]/ ( cuv[j,i] * cuv[j,i] )
1312	aivscu2gam[j,i] = airvscu2[j,i]** ( - gamdi_grot )
1313
1314	airvscu2[j,dx] = airvscu2[j,0]
1315	aivscu2gam[j,dx] = aivscu2gam[j,0]
1316
1317	for j in range(dy):
1318	for i in range(dx):
1319	airuscv2[j,i] = aireu[j,i]/ ( cvu[j,i] * cvu[j,i] )
1320	aiuscv2gam[j,i] = airuscv2[j,i]** ( - gamdi_grot )
1321
1322	airuscv2[j,dx] = airuscv2[j,0]
1323	aiuscv2gam[j,dx] = aiuscv2gam[j,0]
1324
1325	#
1326	# calcul des aires aux poles :
1327	# -----------------------------
1328	#
1329	apoln = SSUM(dx,aire[0,0],1)
1330	apols = SSUM(dx,aire[dy,0],1)
1331	unsapolnga1 = 1./ ( apoln ** ( - gamdi_gdiv ) )
1332	unsapolsga1 = 1./ ( apols ** ( - gamdi_gdiv ) )
1333	unsapolnga2 = 1./ ( apoln ** ( - gamdi_h ) )
1334	unsapolsga2 = 1./ ( apols ** ( - gamdi_h ) )
1335
1336	#
1337	#-----------------------------------------------------------------------
1338	# gtitre='Coriolis version ancienne'
1339	# gfichier='fext1'
1340	# CALL writestd(fext,iip1*jjm)
1341	#
1342	# changement F. Hourdin calcul conservatif pour fext
1343	# constang contient le produit a * cos ( latitude ) * omega
1344	#
1345	for i in range(dx):
1346	constang[0,i] = 0.
1347
1348	for j in range(dy-1):
1349	for i in range(dx):
1350	constang[j+1,i] = radomegcu[j+1,i]*np.cos(rlatu[j+1])
1351
1352	for i in range(dx):
1353	constang[dy,i] = 0.
1354
1355	#
1356	# periodicite en longitude
1357	#
1358	for j in range(dy):
1359	fext[j,dx] = fext[j,0]
1360
1361	for j in range(dy+1):
1362	constang[j,dx] = constang[j,0]
1363
1364	# fin du changement
1365
1366	#
1367	#-----------------------------------------------------------------------
1368	#
1369	print ' * Coordonnees de la grille * '
1370	print ' LONGITUDES aux pts. V ( degres ) '
1371
1372	for i in range(dx+1):
1373	rlonvv[i] = rlonv[i]*180./np.pi
1374
1375	print rlonvv
1376
1377	print ' LATITUDES aux pts. V ( degres ) '
1378
1379	for i in range(dy):
1380	rlatuu[i] = rlatv[i]*180./np.pi
1381
1382	print rlatuu[dy]
1383
1384	for i in range(dx+1):
1385	rlonvv[i]=rlonu[i]*180./np.pi
1386
1387	print ' LONGITUDES aux pts. U ( degres ) '
1388	print rlonvv
1389
1390
1391	print ' LATITUDES aux pts. U ( degres ) '
1392	for i in range(dy+1):
1393	rlatuu[i]=rlatu[i]*180./np.pi
1394
1395	print rlatuu[0:dy+2]
1396
1397	return aire, apoln, apols, airesurg, rlatu, rlatv, cu, cv
1398
1399	def SSUM(n,sx,incx):
1400	""" Obsolete version of sum for non Fortan 90 code
1401	from dyn3d/cray.F
1402	"""
1403
1404	ssumv = 0.
1405
1406	ix = 0
1407
1408	if len(sx.shape) != 0:
1409	for i in range(n):
1410	ssumv=ssumv+sx[ix]
1411	ix=ix+incx
1412	else:
1413	ssumv=ssumv+sx
1414
1415	return ssumv
1416
1417	def SCOPY(ddx,ddy,ddz,sx,incx,sy,incy):
1418	""" Obsolete function to copy matrix values
1419	from dyn3d/cray.F
1420	"""
1421	fname = 'SCOPY'
1422
1423	if len(sx.shape) == 2:
1424	for j in range(ddy-1):
1425	for i in range(ddx-1):
1426	sy[iy,ix] = sx[iy,ix]
1427	ix = ix+incx
1428	iy = iy+incy
1429	elif len(sx.shape) == 3:
1430	iy = 0
1431	for j in range(ddy):
1432	ix = 0
1433	for i in range(ddx):
1434	for l in range(ddz):
1435	sy[l,iy,ix] = sx[l,iy,ix]
1436	ix = ix+incx
1437	iy = iy+incy
1438
1439	return sy
1440
1441	def exner_hyb (dx, dy, dz, psv, pv, aire, apoln, apols):
1442	"""c
1443	c Auteurs : P.Le Van , Fr. Hourdin .
1444	c ..........
1445	c
1446	c .... ngrid, ps,p sont des argum.d'entree au sous-prog ...
1447	c .... alpha,beta, pks,pk,pkf sont des argum.de sortie au sous-prog ...
1448	c
1449	c ************************************************************************
1450	c Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des
1451	c couches . Pk(l) sera calcule aux milieux des couches l ,entre les
1452	c pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
1453	c ************************************************************************
1454	c .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont
1455	c la pression et la fonction d'Exner au sol .
1456	c
1457	c -------- z
1458	c A partir des relations ( 1 ) pdz(pk) = kappa pk*dz(p) et
1459	c ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
1460	c ( voir note de Fr.Hourdin ) ,
1461	c
1462	c on determine successivement , du haut vers le bas des couches, les
1463	c coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
1464	c puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,
1465	c pk(ij,l) donne par la relation (2), pour l = 2 a l = llm .
1466	c
1467	"""
1468
1469	fname = 'exner_hyb'
1470
1471	pksv = np.zeros((dy+1, dx+1), dtype=np.float)
1472	pkv = np.zeros((dz, dy+1, dx+1), dtype=np.float)
1473	pkfv = np.zeros((dz, dy+1, dx+1), dtype=np.float)
1474
1475	ppn = np.zeros((dy+1, dx+1), dtype=np.float)
1476	pps = np.zeros((dy+1, dx+1), dtype=np.float)
1477
1478	alphav = np.zeros((dz+1, dy+1, dx+1), dtype=np.float)
1479	betav = np.zeros((dz+1, dy+1, dx+1), dtype=np.float)
1480
1481	if dz == 1:
1482	# Compute pks(:),pk(:),pkf(:)
1483	pks = (cpp/preff)*ps
1484	pk[0,:,:] = 0.5*pks
1485	# CALL SCOPY ( ngrid * llm, pk, 1, pkf, 1 ) is the same as the next line?
1486	pkf = pk
1487
1488	# No filtering... not necessary on aquaplanet
1489	# CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
1490
1491	# our work is done, exit routine
1492	return pksv, pkv, pkfv
1493
1494	#### General case:
1495
1496	unpl2k = 1.+ 2.* kappa
1497
1498	for j in range(dy+1):
1499	for i in range(dx+1):
1500	pksv[j,i] = cpp * ( psv[j,i]/preff ) ** kappa
1501
1502	for j in range(dy+1):
1503	for i in range(dx+1):
1504	ppn[j,i] = aire[j,i] * pksv[j,i]
1505	pps[j,i] = aire[dy,i] * pksv[dy,i]
1506
1507	xpn = SSUM(dx,ppn,1) /apoln
1508	xps = SSUM(dx,pps,1) /apols
1509
1510	for j in range(dy):
1511	for i in range(dx):
1512	pksv[j,i] = xpn[i]
1513	pksv[dy-1,i] = xps[i]
1514	#
1515	#
1516	# .... Calcul des coeff. alpha et beta pour la couche l = llm ..
1517	#
1518	for j in range(dy+1):
1519	for i in range(dx+1):
1520	alphav[dz-1,j,i] = 0.
1521	betav[dz-1,j,i] = 1./ unpl2k
1522
1523	#
1524	# ... Calcul des coeff. alpha et beta pour l = llm-1 a l = 2 ...
1525	#
1526	for l in range (dz-1,1,-1):
1527	for j in range(dy+1):
1528	for i in range(dx+1):
1529	dellta = pv[l,j,i]* unpl2k + pv[l+1,j,i]* ( betav[l+1,j,i]-unpl2k )
1530	alphav[l,j,i] = -pv[l+1,j,i] / dellta * alphav[l+1,j,i]
1531	betav[l,j,i] = pv[l,j,i] / dellta
1532
1533	# ***********************************************************************
1534	# ..... Calcul de pk pour la couche 1 , pres du sol ....
1535	#
1536	for j in range(dy+1):
1537	for i in range(dx+1):
1538	pkv[0,j,i] = ( pv[0,j,i]pksv[j,i] - 0.5alphav[1,j,i]*pv[1,j,i] ) \
1539	( pv[0,j,i] (1.+kappa) + 0.5( betav[1,j,i]-unpl2k ) pv[1,j,i] )
1540
1541	#
1542	# ..... Calcul de pk(ij,l) , pour l = 2 a l = llm ........
1543	#
1544	for l in range(dz):
1545	for j in range(dy+1):
1546	for i in range(dx+1):
1547	pkv[l,j,i] = alphav[l,j,i] + betav[l,j,i] * pkv[l-1,j,i]
1548	#
1549	#
1550	pkfv = SCOPY ( dx+1, dy+1, dz, pkv, 1, pkfv, 1 )
1551
1552	# We do not filter for iniaqua
1553	# CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
1554
1555	return pksv, pkv, pkfv, alphav, betav
1556
1557	def exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf ):
1558	"""c
1559	c Auteurs : F. Forget , Y. Wanherdrick
1560	c P.Le Van , Fr. Hourdin .
1561	c ..........
1562	c
1563	c .... ngrid, ps,p sont des argum.d'entree au sous-prog ...
1564	c .... beta, pks,pk,pkf sont des argum.de sortie au sous-prog ...
1565	c
1566	c ************************************************************************
1567	c Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
1568	c couches . Pk(l) sera calcule aux milieux des couches l ,entre les
1569	c pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
1570	c ************************************************************************
1571	c .. N.B : Au sommet de l'atmosphere, p(llm+1) = 0. , et ps et pks sont
1572	c la pression et la fonction d'Exner au sol .
1573	c
1574	c WARNING : CECI est une version speciale de exner_hyb originale
1575	c Utilise dans la version martienne pour pouvoir
1576	c tourner avec des coordonnees verticales complexe
1577	c => Il ne verifie PAS la condition la proportionalite en
1578	c energie totale/ interne / potentielle (F.Forget 2001)
1579	c ( voir note de Fr.Hourdin ) ,
1580	c
1581	"""
1582	fname = 'exner_milieu'
1583
1584	pks = np.zeros((dy, dx), dtype=np.float)
1585	pk = np.zeros((dz, dy, dx), dtype=np.float)
1586	pkf = np.zeros((dz, dy, dx), dtype=np.float)
1587
1588	ppn = np.zeros((iim), dtype=np.float)
1589	ppn = np.zeros((iim), dtype=np.float)
1590
1591	ip1jm = (dx+1)*dy
1592
1593	firstcall = True
1594	modname = 'exner_milieu'
1595
1596	# Sanity check
1597	if firstcall:
1598	# sanity checks for Shallow Water case (1 vertical layer)
1599	if llm == 1:
1600	if kappa != 1:
1601	print errormsg
1602	print ' ' + fname+ ': kappa!=1 , but running in Shallow Water mode!!'
1603	quit(-1)
1604	if cpp != r:
1605	print errormsg
1606	print ' ' + fname+ ': cpp!=r , but running in Shallow Water mode!!'
1607	quit(-1)
1608
1609	firstcall = False
1610
1611	#### Specific behaviour for Shallow Water (1 vertical layer) case:
1612	if llm == 1:
1613
1614	# Compute pks(:),pk(:),pkf(:)
1615
1616	for j,i in range(ngrid):
1617	pks[j,i] = (cpp/preff) * ps[j,i]
1618	pk[j,i,1] = .5*pks[j,i]
1619
1620	pkf = SCOPY(dx,dy,dz, pk, 1, pkf, 1 )
1621	# We do not filter for iniaqua
1622	# CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
1623
1624	# our work is done, exit routine
1625	return pksv, pkv, pkfv
1626
1627	#### General case:
1628
1629	# -------------
1630	# Calcul de pks
1631	# -------------
1632
1633	for j,i in range(ngrid):
1634	pks[j,i] = cpp * ( ps[j,i]/preff ) ** kappa
1635
1636	for j,i in range(iim):
1637	ppn[j,i] = aire[j,i] * pks[j,i]
1638	pps[j,i] = aire[j,i+ip1jm] * pks[j,i+ip1jm]
1639
1640	xpn = SSUM(iim,ppn,1) /apoln
1641	xps = SSUM(iim,pps,1) /apols
1642
1643	for j,i in range(iip1):
1644	pks[j,i] = xpn
1645	pks[j,i+ip1jm] = xps
1646
1647	#
1648	#
1649	# .... Calcul de pk pour la couche l
1650	# --------------------------------------------
1651	#
1652	dum1 = cpp * (2preff)*(-kappa)
1653	for l in range(llm-1):
1654	for j,i in range(ngrid):
1655	pk[j,i,l] = dum1 * (p[j,i,l] + p[j,i,l+1])**kappa
1656
1657	# .... Calcul de pk pour la couche l = llm ..
1658	# (on met la meme distance (en log pression) entre Pk(llm)
1659	# et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
1660
1661	for j,i in range(ngrid):
1662	pk[j,i,llm] = pk[j,i,llm-1]**2 / pk[j,i,llm-2]
1663
1664	# calcul de pkf
1665	# -------------
1666	pkf = SCOPY( dx,dy,dz, pk, 1, pkf, 1 )
1667
1668	# We do not filter iniaqua
1669	# CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
1670
1671	# EST-CE UTILE ?? : calcul de beta
1672	# --------------------------------
1673	for l in range(1, llm):
1674	for j,i in range(ngrid):
1675	beta[j,i,l] = pk[j,i,l] / pk[j,i,l-1]
1676
1677	return pksv, pkv, pkfv
1678
1679	def pression(dx, dy, dz, apv, bpv, psv):
1680	"""c
1681
1682	c Auteurs : P. Le Van , Fr.Hourdin .
1683
1684	c ************************************************************************
1685	c Calcule la pression p(l) aux differents niveaux l = 1 ( niveau du
1686	c sol) a l = llm +1 ,ces niveaux correspondant aux interfaces des (llm)
1687	c couches , avec p(ij,llm +1) = 0. et p(ij,1) = ps(ij) .
1688	c ************************************************************************
1689	c
1690	"""
1691	fname = 'pression'
1692
1693	press = np.zeros((dz+1, dy+1, dx+1), dtype=np.float)
1694
1695	for l in range(dz+1):
1696	press[l,:,:] = apv[l] + bpv[l]*psv
1697
1698	return press
1699
1700	def sig_hybrid(sig,pa,preff):
1701	""" Function utilisee pour calculer des valeurs de sigma modifie
1702	pour conserver les coordonnees verticales decrites dans
1703	esasig.def/z2sig.def lors du passage en coordonnees hybrides
1704	F. Forget 2002
1705	sig= sigma level
1706	pa=
1707	preff= reference pressure
1708	Connaissant sig (niveaux "sigma" ou on veut mettre les couches)
1709	L'objectif est de calculer newsig telle que
1710	(1 -pa/preff)exp(1-1./newsig2)+(pa/preff)newsig = sig
1711	Cela ne se resoud pas analytiquement:
1712	=> on resoud par iterration bourrine
1713	----------------------------------------------
1714	Information : where exp(1-1./x**2) become << x
1715	x exp(1-1./x**2) /x
1716	1 1
1717	0.68 0.5
1718	0.5 1.E-1
1719	0.391 1.E-2
1720	0.333 1.E-3
1721	0.295 1.E-4
1722	0.269 1.E-5
1723	0.248 1.E-6
1724	=> on peut utiliser newsig = sigpreff/pa si sigpreff/pa < 0.25
1725	"""
1726	fname = 'sig_hybrid'
1727
1728	# maximum number of iterations
1729	maxiter = 9999
1730
1731	nsig = sig
1732	x1=0
1733	x2=1
1734	if sig >= 1.:
1735	nsig = sig
1736	elif sig*preff/pa >= 0.25:
1737	for j in range(maxiter):
1738	F = ((1. -pa/preff)np.exp(1.-1./nsig2)+(pa/preff)nsig)/sig
1739	# print J,'nsig=', newsig, 'F=', F
1740	if F > 1:
1741	X2 = newsig
1742	nsig = (X1+nsig)*0.5
1743	else:
1744	X1 = newsig
1745	nsig = (X2+nsig)*0.5
1746
1747	# Test : on arete lorsque on approxime sig a moins de 0.01 m pres
1748	# (en pseudo altitude) :
1749	if np.abs(10.*np.log(F)) < 1.e-5: break
1750	else:
1751	nsig= sig*preff/pa
1752
1753	return nsig
1754
1755	def presnivs_calc(vert_sampling, dz):
1756	""" From dyn3d/disvert.F calculation of vertical pressure levels
1757	vert_sampling= which kind of vertical sampling is desired: "param", "tropo",
1758	"strato" and "read"
1759	dz= numbef of vertical layers
1760	"""
1761
1762	fname = 'presnivs_calc'
1763
1764	llmp1 = dz + 1
1765	pnivs = np.zeros((dz), dtype=np.float)
1766	s = np.zeros((dz), dtype=np.float)
1767	sig = np.zeros((dz+1), dtype=np.float)
1768	dsig = np.zeros((dz), dtype=np.float)
1769	dpres = np.zeros((dz), dtype=np.float)
1770	apv = np.zeros((dz+1), dtype=np.float)
1771	bpv = np.zeros((dz+1), dtype=np.float)
1772
1773	# default scaleheight is 8km for earth
1774	scaleheight = 8.
1775
1776	# vert_sampling = merge("strato", "tropo ", ok_strato) ! default value
1777
1778	if vert_sampling == 'param':
1779	# On lit les options dans sigma.def:
1780	if not os.path.isfile('easig.def'):
1781	print errormsg
1782	print ' ' + fname + ": parameters file 'easig.def' does not exist!!"
1783	quit(-1)
1784
1785	sfobj = open('sigma.def', 'r')
1786	scaleheight = np.float( ncvar.reduce_spaces(fobj.readline()))
1787	deltaz = np.float( ncvar.reduce_spaces(fobj.readline()))
1788	beta = np.float( ncvar.reduce_spaces(fobj.readline()))
1789	k0 = np.float( ncvar.reduce_spaces(fobj.readline()))
1790	k1 = np.float( ncvar.reduce_spaces(fobj.readline()))
1791	sfobj.close()
1792
1793	alpha = deltaz/(dz*scaleheight)
1794	print ':scaleheight, alpha, k0, k1, beta', scaleheight, alpha, k0, k1, beta
1795
1796	alpha=deltaz/np.tanh(1./k0)*2.
1797	zkm1=0.
1798	sig[0]=1.
1799	for l in range(dz):
1800	sig[l+1]=(np.cosh(l/k0))*(-alphak0/scaleheight) \
1801	exp(-alpha/scaleheightnp.tanh((llm-k1)/k0) \
1802	beta*(l-(llm-k1))/np.log(beta))
1803
1804	zk=-scaleheight*np.log(sig[l+1])
1805
1806	dzk1=alpha*np.tanh(l/k0)
1807	dzk2=alphanp.tanh((llm-k1)/k0)beta**(l-(llm-k1))/np.log(beta)
1808
1809	print l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
1810	zkm1=zk
1811
1812	sig[dz-1]=0.
1813
1814	bpv[0:dz] = np.exp(1.-1./sig[0:dz]**2)
1815	bpv[llmp1-1] = 0.
1816
1817	apv = pa * (sig - bp)
1818
1819	elif vert_sampling == 'tropo':
1820	for l in range(dz):
1821	x = 2.np.arcsin(1.)(l-0.5)/(dz+1.)
1822	dsig[l] = 1.0+7.0np.sin(x)*2
1823
1824	dsig = dsig/np.sum(dsig)
1825	sig[dz] = 0.
1826	for l in range(dz-1,0,-1):
1827	sig[l] = sig[l+1] + dsig[l]
1828
1829	bpv[0]=1.
1830	bpv[1:dz] = np.exp(1.-1./sig[1:dz]**2)
1831	bpv[llmp1-1] = 0.
1832
1833	apv[0] = 0.
1834	apv[1:dz+1] = pa*(sig[1:dz+1]-bpv[1:dz+1])
1835
1836	elif vert_sampling == 'strato':
1837	if dz == 39:
1838	dsigmin = 0.3
1839	elif dz == 50:
1840	dsigmin = 1.
1841	else:
1842	print ' ATTENTION discretisation z a ajuster'
1843	dsigmin = 1.
1844
1845	print 'Discretisation verticale DSIGMIN=',dsigmin
1846
1847	for l in range(dz):
1848	x = 2.np.arcsin(1.)(l - 0.5)/(dz+1)
1849	dsig[l] =(dsigmin+7.np.sin(x)*2) \
1850	(0.5(1.-np.tanh(1.(x-np.arcsin(1.))/np.arcsin(1.))))*2
1851
1852	dsig = dsig/np.sum(dsig)
1853	sig[dz] = 0.
1854	for l in range(dz-1,0,-1):
1855	sig[l] = sig[l+1] + dsig[l]
1856
1857	bpv[0] = 1.
1858	bpv[1:dz] = np.exp(1.-1./sig[1:dz]**2)
1859	bpv[llmp1-1] = 0.
1860
1861	apv[0] = 0.
1862	apv[1:dz+1] = pa*(sig[1:dz+1] - bpv[1:dz+1])
1863
1864	elif vert_sampling == 'read':
1865	# Read "ap" and "bp". First line is skipped (title line). "ap"
1866	# should be in Pa. First couple of values should correspond to
1867	# the surface, that is : "bp" should be in descending order.
1868	if not os.path.isfile('hybrid.txt'):
1869	print errormsg
1870	print ' ' + fname + ": parameters file 'hybrid.txt' does not exist!!"
1871	quit(-1)
1872	sfobj = ('hybrid.txt', 'r')
1873	# skip title line
1874	title = sfobj.readline()
1875	for l in range(dz+1):
1876	values = ncvar.readuce_space(sfobj.readline())
1877	apv[l] = np.float(values[0])
1878	bpv[l] = np.float(values[0])
1879
1880	sfobj.close()
1881	if apv[0] == 0. or apv[dz+1] == 0. or bpv[0] == 1. or bpv[dz+1] == 0.:
1882	print errormsg
1883	print ' ' + fname + ': bad ap or bp values !!'
1884	print ' k ap bp ___________'
1885	for k in range(dz+1):
1886	print k, apv[k], bpv[k]
1887
1888	else:
1889	print errormsg
1890	print ' ' + fname + ': wrong value for vert_sampling:', vert_sampling
1891	quit(-1)
1892
1893
1894	nivsigs = np.arange(dz)*1.+1.
1895	nivsig = np.arange(llmp1)*1.+1.
1896
1897	print ' ' + fname + ': k BP AP ________'
1898	for k in range(dz+1):
1899	print k, bpv[k], apv[k]
1900
1901	print 'Niveaux de pressions approximatifs aux centres des'
1902	print 'couches calcules pour une pression de surface =', preff
1903	print 'et altitudes equivalentes pour une hauteur d echelle de '
1904	print scaleheight,' km'
1905
1906	for l in range(dz):
1907	dpres[l] = bpv[l] - bpv[l+1]
1908	pnivs[l] = 0.5( apv[l]+bpv[l]preff + apv[l+1]+bpv[l+1]*preff )
1909	print ' PRESNIVS(', l, ')=', pnivs[l], ' Z ~ ', \
1910	np.log(preff/pnivs[l])*scaleheight, ' DZ ~ ', \
1911	scaleheightnp.log((apv[l]+bpv[l]preff)/ \
1912	np.max([apv[l+1]+bpv[l+1]*preff, 1.e-10]))
1913
1914	print ' ' + fname + ': PRESNIVS [Pa]:', pnivs
1915
1916	return pnivs, apv, bpv
1917
1918
1919	def presnivs_calc_noterre(hybrid, dz):
1920	""" From dyn3d/disvert_noterre.F calculation of vertical pressure levels
1921	hybrid= whether hydbrid coordinates have to be used
1922	dz= numbef of vertical layers
1923	"""
1924
1925	fname = 'presnivs_calc_noterre'
1926
1927	#-----------------------------------------------------------------------
1928	# .... Calculs de ap(l) et de bp(l) ....
1929	# .........................................
1930	#
1931	# ..... pa et preff sont lus sur les fichiers start par dynetat0 .....
1932	#-----------------------------------------------------------------------
1933	#
1934	llmp1 = dz + 1
1935	pnivs = np.zeros((dz), dtype=np.float)
1936	s = np.zeros((dz), dtype=np.float)
1937	sig = np.zeros((dz+1), dtype=np.float)
1938	apv = np.zeros((dz+1), dtype=np.float)
1939	bpv = np.zeros((dz+1), dtype=np.float)
1940
1941	# Ouverture possible de fichiers typiquement E.T.
1942
1943	if os.path.isfile('easig.def'):
1944	#-----------------------------------------------------------------------
1945	# cas 1 on lit les options dans esasig.def:
1946	# ----------------------------------------
1947
1948	ofet = open('esasig.def', 'r')
1949	# Lecture de esasig.def :
1950	# Systeme peu souple, mais qui respecte en theorie
1951	# La conservation de l'energie (conversion Energie potentielle
1952	# <-> energie cinetique, d'apres la note de Frederic Hourdin...
1953
1954	print '*****************************'
1955	print "WARNING reading 'esasig.def'"
1956	print '*****************************'
1957	for line in ofet:
1958	linevalues = ncvar.reduce_spaces(line)
1959	scaleheight = np.float(linevalues[0])
1960	dz0 = np.float(linevalues[1])
1961	dz1 = np.float(linevalues[2])
1962	nhaut = np.float(linevalues[3])
1963
1964	ofet.close()
1965	dz0 = dz0/scaleheight
1966	dz1 = dz1/scaleheight
1967
1968	sig1=(1.-dz1)/tanh(.5*(dz-1)/nhaut)
1969
1970	esig=1.
1971	for l in range(19):
1972	esig=-np.log((1./sig1-1.)*np.exp(-dz0)/esig)/(dz-1.)
1973
1974	csig=(1./sig1-1.)/(np.exp(esig)-1.)
1975
1976	for l in range(1,dz):
1977	zz=csig(np.exp(esig(l-1.))-1.)
1978	sig[l] = 1./(1.+zz)np.tanh(.5(dz+1-l)/nhaut)
1979
1980	sig[0] = 1.
1981	sig[dz] = 0.
1982	quoi = 1. + 2.* kappa
1983	s[dz-1] = 1.
1984	s[dz-2] = quoi
1985	if dz > 1:
1986	for ll in range(1, dz-2):
1987	l = dz+1 - ll
1988	quand = sig[l+1]/sig[l]
1989	s[l-1] = quoi(1.-quand)s[l] + quand*s[l+1]
1990	#
1991	snorm=(1.-.5sig[1]+kappa(1.-sig[1]))s[0]+.5sig[1]*s[1]
1992	for l in range(dz):
1993	s[l] = s[l]/ snorm
1994	elif os.path.isfile('z2sig.def'):
1995	fet = open('z2sig.def', 'r')
1996	#-----------------------------------------------------------------------
1997	# cas 2 on lit les options dans z2sig.def:
1998	# ----------------------------------------
1999	print '****************************'
2000	print 'Reading z2sig.def'
2001	print '****************************'
2002
2003	for line in ofet:
2004	linevalues = ncvar.reduce_spaces(line)
2005	scaleheight = np.float(linevalues[0])
2006	for l in range(dz):
2007	zsig[l] = linevalues[l+1]
2008
2009	ofet.close()
2010
2011	sig[0] = 1.
2012	for l in range(1,dz):
2013	sig[l] = 0.5 * ( np.exp(-zsig[l]/scaleheight) + \
2014	np.exp(-zsig[l-1]/scaleheight) )
2015
2016	sig[dz+1] = 0.
2017
2018	#-----------------------------------------------------------------------
2019	else:
2020	print errormsg
2021	print ' ' + fname + ": didn't you forget something ???"
2022	print " We need file 'z2sig.def' ! (OR 'esasig.def')"
2023	quit(-1)
2024	#-----------------------------------------------------------------------
2025
2026	nivsigs = np.arange(dz)*1.
2027	nivsig = np.arange(llmp1)*1.
2028
2029	if hybrid:
2030	# use hybrid coordinates
2031	print "*********************************"
2032	print "Using hybrid vertical coordinates"
2033	print
2034	# Coordonnees hybrides avec mod
2035	for l in range(dz):
2036	newsig = sig_hybrid(sig[l],pa,preff)
2037	bpv[l] = np.exp(1.-1./(newsig**2))
2038	apv[l] = pa * (newsig - bp[l] )
2039
2040	bp[llmp1-1] = 0.
2041	ap[llmp1-1] = 0.
2042	else:
2043	# use sigma coordinates
2044	print "********************************"
2045	print "Using sigma vertical coordinates"
2046	print
2047	# Pour ne pas passer en coordonnees hybrides
2048	for l in range(dz):
2049	apv[l] = 0.
2050	bpv[l] = sig[l]
2051
2052	apv[llmp1-1] = 0.
2053
2054	bp[llmp1-1] = 0.
2055
2056	print 'BP________ ', bp
2057	print 'AP________ ', ap
2058
2059	# Calcul au milieu des couches (llm = dz):
2060	# WARNING : le choix de placer le milieu des couches au niveau de
2061	# pression intermediaire est arbitraire et pourrait etre modifie.
2062	# Le calcul du niveau pour la derniere couche
2063	# (on met la meme distance (en log pression) entre P(llm)
2064	# et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
2065	# Specifique. Ce choix est specifie ici ET dans exner_milieu.F
2066
2067	for l in range(dz-1):
2068	aps[0:dz-1] = 0.5*( apv[0:dz-1]+apv[1:dz])
2069	bps[0:dz-1] = 0.5*( bpv[0:dz-1]+bpv[1:dz])
2070
2071	if hybrid:
2072	aps[dz-1] = aps[dz-2]**2 / aps[dz-3]
2073	bps[dz-1] = 0.5*(bpv[dz-1] + bpv[dz])
2074	else:
2075	bps[dz-1] = bps[dz-2]**2 / bps[dz-3]
2076	aps[dz-1] = 0.
2077
2078	print 'BPs_______ ', bps
2079	print 'APs_______ ', aps
2080
2081	for l in range(dz):
2082	psnivs[l] = aps[l]+bps[l]*preff
2083	psalt[l] = -scaleheight*np.log(presnivs[l]/preff)
2084
2085	return psnivs, psalt
2086
2087	def global_lonlat(dx,dy):
2088	""" Function to generate 2D matrices with the global longitude, latitudes
2089	dx, dy: dimension of the desired matrix
2090	>>> global_lonlat(5,5)
2091	array([[ 36., 108., 180., 252., 324.],
2092	[ 36., 108., 180., 252., 324.],
2093	[ 36., 108., 180., 252., 324.],
2094	[ 36., 108., 180., 252., 324.],
2095	[ 36., 108., 180., 252., 324.]]), array([[-72., -72., -72., -72., -72.],
2096	[-36., -36., -36., -36., -36.],
2097	[ 0., 0., 0., 0., 0.],
2098	[ 36., 36., 36., 36., 36.],
2099	[ 72., 72., 72., 72., 72.]]))
2100	"""
2101
2102	fname = 'global_lonlat'
2103
2104	longitude = np.zeros((dy,dx), dtype=np.float)
2105	latitude = np.zeros((dy,dx), dtype=np.float)
2106
2107	for ix in range(dx):
2108	longitude[:,ix] = 360.*(1./2. + ix )/(dx)
2109
2110	for iy in range(dy):
2111	latitude[iy,:] = 180.*(1./2. + iy )/(dy) - 90.
2112
2113	return longitude, latitude
2114
2115	def massdair(dx,dy,dz,p,airesurg):
2116	"""c
2117	c *********************************************************************
2118	c .... Calcule la masse d'air dans chaque maille ....
2119	c *********************************************************************
2120	c
2121	c Auteurs : P. Le Van , Fr. Hourdin .
2122	c ..........
2123	c
2124	c .. p est un argum. d'entree pour le s-pg ...
2125	c .. masse est un argum.de sortie pour le s-pg ...
2126	c
2127	c .... p est defini aux interfaces des llm couches .....
2128	c
2129	"""
2130	fname = 'massdair'
2131
2132	masse = np.zeros((dz+1,dy+1,dx+1), dtype=np.float)
2133	#
2134	#
2135	# Methode pour calculer massebx et masseby .
2136	# ----------------------------------------
2137	#
2138	# A chaque point scalaire P (i,j) est affecte 4 coefficients d'aires
2139	# alpha1(i,j) calcule au point ( i+1/4,j-1/4 )
2140	# alpha2(i,j) calcule au point ( i+1/4,j+1/4 )
2141	# alpha3(i,j) calcule au point ( i-1/4,j+1/4 )
2142	# alpha4(i,j) calcule au point ( i-1/4,j-1/4 )
2143	#
2144	# Avec alpha1(i,j) = aire(i+1/4,j-1/4)/ aire(i,j)
2145	#
2146	# N.B . Pour plus de details, voir s-pg ... iniconst ...
2147	#
2148	#
2149	#
2150	# alpha4 . . alpha1 . alpha4
2151	# (i,j) (i,j) (i+1,j)
2152	#
2153	# P . U . . P
2154	# (i,j) (i,j) (i+1,j)
2155	#
2156	# alpha3 . . alpha2 .alpha3
2157	# (i,j) (i,j) (i+1,j)
2158	#
2159	# V . Z . . V
2160	# (i,j)
2161	#
2162	# alpha4 . . alpha1 .alpha4
2163	# (i,j+1) (i,j+1) (i+1,j+1)
2164	#
2165	# P . U . . P
2166	# (i,j+1) (i+1,j+1)
2167	#
2168	#
2169	#
2170	# On a :
2171	#
2172	# massebx(i,j) = masse(i ,j) * ( alpha1(i ,j) + alpha2(i,j)) +
2173	# masse(i+1,j) * ( alpha3(i+1,j) + alpha4(i+1,j) )
2174	# localise au point ... U (i,j) ...
2175	#
2176	# masseby(i,j) = masse(i,j ) * ( alpha2(i,j ) + alpha3(i,j ) +
2177	# masse(i,j+1) * ( alpha1(i,j+1) + alpha4(i,j+1)
2178	# localise au point ... V (i,j) ...
2179	#
2180	#
2181	#=======================================================================
2182
2183	for l in range (dz-1):
2184	for j in range(dy+1):
2185	for i in range(dx+1):
2186	masse[l,j,i] = airesurg[j,i] * ( p[l,j,i] - p[l+2,j,i] )
2187
2188	for j in range(1,dy):
2189	masse[l,j,dx] = masse[l,j-1,dx]
2190
2191	return masse
2192
2193	def geopot(dx, dy, dz, teta, pk, pks):
2194	"""c=======================================================================
2195	c
2196	c Auteur: P. Le Van
2197	c -------
2198	c
2199	c Objet:
2200	c ------
2201	c
2202	c *******************************************************************
2203	c .... calcul du geopotentiel aux milieux des couches .....
2204	c *******************************************************************
2205	c
2206	c .... l'integration se fait de bas en haut ....
2207	c
2208	c .. ngrid,teta,pk,pks,phis sont des argum. d'entree pour le s-pg ..
2209	c phi est un argum. de sortie pour le s-pg .
2210	c
2211	c=======================================================================
2212	"""
2213	fname = 'geopot'
2214
2215	phis = np.zeros((dy+1,dx+1), dtype=np.float)
2216	phi = np.zeros((dz+1,dy+1,dx+1), dtype=np.float)
2217
2218	print ' Lluis in ' + fname + ' shapes phis:',phis.shape,'phi:',phi.shape, \
2219	'teta:',teta.shape,'pks:',pks.shape,'pk:',pk.shape
2220
2221	#-----------------------------------------------------------------------
2222	# calcul de phi au niveau 1 pres du sol .....
2223
2224	for j in range(dy):
2225	for i in range(dx):
2226	phi[0,j,i] = phis[j,i] + teta[0,j,i] * ( pks[j,i] - pk[0,j,i] )
2227
2228	# calcul de phi aux niveaux superieurs .......
2229
2230	for l in range(1,dz):
2231	for j in range(dy):
2232	for i in range(dx):
2233	phi[l,j,i] = phi[l-1,j,i] + 0.5 * ( teta[l,j,i]+teta[l-1,j,i] ) * \
2234	( pk[l-1,j,i]-pk[l,j,i] )
2235
2236	return phis, phi
2237
2238	def dump2d(im,jm,nom_z):
2239	""" Function to create a dump 2d variable
2240	from dyn3d/dump2d.F
2241	"""
2242	fname = 'dump2d'
2243
2244	z = np.zeros((im,jm), dtype=np.float)
2245
2246	zmin = z[0,0]
2247	zllm = z[0,0]
2248	imin = 0
2249	illm = 0
2250	jmin = 0
2251	jllm = 0
2252
2253	for j in range(jm):
2254	for i in range(im):
2255	if z[i,j] > zllm:
2256	illm=i
2257	jllm=j
2258	zllm=z[i,j]
2259
2260	if z[i,j] < zmin:
2261	imin=i
2262	jmin=j
2263	zmin=z[i,j]
2264
2265	print 'MIN:',zmin
2266	print 'MAX:',zllm
2267
2268	if zllm > zmin:
2269	for j in range(jm):
2270	print int(10.*(z[:,j]-zmin)/(zllm-zmin))
2271
2272	return z
2273
2274	def ugeostr(dx,dy,dz,phis,phi,rlatu,rlatv,cu):
2275	"""! Calcul du vent covariant geostrophique a partir du champ de
2276	! geopotentiel.
2277	! We actually compute: (1 - cos^8 \phi) u_g
2278	! to have a wind going smoothly to 0 at the equator.
2279	! We assume that the surface pressure is uniform so that model
2280	! levels are pressure levels.
2281	"""
2282	fname = 'ugeostr'
2283	ucov = np.zeros((dz,dy+1,dx+1), dtype=np.float)
2284	um = np.zeros((dz,dy), dtype=np.float)
2285	u = np.zeros((dz,dy,dx+1), dtype=np.float)
2286
2287	print ' Lluis in ' + fname + ': shapes phis:',phis.shape,'phi:',phi.shape,'u:',u.shape
2288
2289	for j in range(dy):
2290	if np.abs(np.sin(rlatv[j])) < 1.e-4:
2291	zlat = 1.e-4
2292	else:
2293	zlat=rlatv[j]
2294
2295	fact = np.cos(zlat)
2296	fact = fact*fact
2297	fact = fact*fact
2298	fact = fact*fact
2299	fact = (1.-fact)/ (2.omegnp.sin(zlat)*(rlatu[j+1]-rlatu[j]))
2300	fact = -fact/rad
2301	for l in range(dz):
2302	for i in range(dx):
2303	u[l,j,i] = fact*(phi[l,j+1,i]-phi[l,j,i])
2304	um[l,j]=um[l,j]+u[l,j,i]/np.float(dx)
2305
2306	um = dump2d(dz,dy,'Vent-u geostrophique')
2307
2308	# calcul des champ de vent:
2309
2310	for l in range(dz):
2311	for i in range(dx+1):
2312	ucov[l,0,i]=0.
2313	ucov[l,dy,i]=0.
2314	for j in range(1,dy):
2315	for i in range(dx):
2316	ucov[l,j,i] = 0.5(u[l,j,i]+u[l,j-1,i])cu[j,i]
2317
2318	ucov[l,j,dx]=ucov[l,j,0]
2319
2320	return ucov
2321
2322	def RAN1(IDUM, Nvals):
2323	""" Function to generate Nvals random numbers
2324	from dyn3d/ran1.F
2325	IDUM= Random Seed
2326	Nvals= number of values
2327	"""
2328	fname = 'RAN1'
2329
2330	R = np.zeros((Nvals), dtype=np.float)
2331
2332	M1 = 259200
2333	IA1 = 7141
2334	IC1 = 54773
2335	RM1 = 3.8580247E-6
2336	M2 = 134456
2337	IA2 = 8121
2338	IC2 = 28411
2339	RM2 = 7.4373773E-6
2340	M3 = 243000
2341	IA3 = 4561
2342	IC3 = 51349
2343	IFF = 0
2344
2345	if IDUM < 0 or IFF == 0:
2346	IFF = 1
2347	IX1 = np.mod(IC1-IDUM,M1)
2348	IX1 = np.mod(IA1*IX1+IC1,M1)
2349	IX2 = np.mod(IX1,M2)
2350	IX1 = np.mod(IA1*IX1+IC1,M1)
2351	IX3 = np.mod(IX1,M3)
2352	for J in range(Nvals):
2353	IX1 = np.mod(IA1*IX1+IC1,M1)
2354	IX2 = np.mod(IA2*IX2+IC2,M2)
2355	R[J] = (np.float(IX1)+np.float(IX2)RM2)RM1
2356
2357	IDUM=1
2358
2359	IX1 = np.mod(IA1*IX1+IC1,M1)
2360	IX2 = np.mod(IA2*IX2+IC2,M2)
2361	IX3 = np.mod(IA3*IX3+IC3,M3)
2362	J = 1+(Nvals*IX3)/M3
2363	if J > Nvals or J < 1: quit()
2364	ran1=R[J]
2365	R[J]=(np.float(IX1)+np.float(IX2)RM2)RM1
2366
2367	return ran1
2368
2369	def name_variables(filekind):
2370	""" Function to provide name of the variables and their atributes as function of
2371	the output type of file
2372	filekind= kind of file
2373	'CF': CF-convention
2374	'LMDZ': LMDZ style
2375	'WRF': WRF style
2376	"""
2377	fname = 'name_variables'
2378
2379	dimnames = {}
2380	varnames = {}
2381
2382	# Standard dimensions
2383	dimn = ['x','y','z','t']
2384
2385	# Standard variables
2386	varn = ['lon', 'lat', 'lev', 'time', 'temp', 'tsfc', 'u10m', 'v10m', 'u', 'v', \
2387	'zsfc', 'geopot', 'psfc', 'pres', 'H2Ov', 'H2Ol']
2388
2389	# Standard variables' attribute names
2390	stdattrn = ['standard_name', 'long_name', 'units']
2391
2392	# Extra variables' attribute names
2393	# kextrattrn = []
2394
2395	# Dictionary with the values for each standard variable
2396	# kvarn['stdn'] = ['name', 'dims(list)', 'std_name', 'long_name', 'units', 'FillValue']
2397	# None: No value
2398
2399	kdimn = {}
2400	kvarn = {}
2401	kattrn = {}
2402	if filekind == 'CF':
2403	kdimn['x'] = 'x'
2404	kdimn['y'] = 'y'
2405	kdimn['z'] = 'z'
2406	kdimn['t'] = 'time'
2407
2408	kvarn['lon'] = ['lon',['x'],'longitude','Longitude','degrees_east',None]
2409	kvarn['lat'] = ['lat',['y'],'latitude','Latitude','degrees_north',None]
2410	kvarn['lev'] = ['lev',['z'],'levels','Levels','-',None]
2411	kvarn['time'] = ['time',['t'],'time','Time', \
2412	'monutes since 1949-12-01 00:00:00', None]
2413	kvarn['temp'] = ['ta',['t','z','y','x'],'air_temperature','Air temperature', \
2414	'K',None]
2415	kvarn['tsfc'] = ['tas',['t','y','x'],'air_temperature','Air temperature','K',None]
2416	kvarn['u10m'] = ['uas',['t','y','x'],'eastward_wind','eastward wind','ms-1',None]
2417	kvarn['v10m'] = ['vas',['t','y','x'],'northward_wind','northward wind','ms-1',None]
2418	kvarn['u'] = ['ua',['t','z','y','x'],'eastward_wind','eastward wind','ms-1',None]
2419	kvarn['v'] = ['va',['t','z','y','x'],'northward_wind','northward wind','ms-1',None]
2420	kvarn['zsfc'] = ['zgs',['t','y','x'],'surface_geopotential_height', \
2421	'surface geopotential height','m2s-2',None]
2422	kvarn['geopot'] = ['zg',['t','z','y','x'],'geopotential_height','geopotential height', \
2423	'm2s-2',None]
2424	kvarn['psfc'] = ['ps',['t','z','y','x'],'surface_air_pressure','surface pressure','Pa',None]
2425	kvarn['pres'] = ['pres',['t','z','y','x'],'air_pressure','pressure','Pa',None]
2426	kvarn['H2Ov'] = ['r',['t','z','y','x'],'water_mixing_ratio','water mixing ratio','kgkg-1', \
2427	None]
2428	kvarn['H2Ol'] = ['c',['t','z','y','x'],'condensed_water_mixing_ratio', \
2429	'condensed water mixing ratio','kgkg-1',None]
2430
2431	kattrn['standard_name'] = 'standard_name'
2432	kattrn['long_name'] = 'long_name'
2433	kattrn['units'] = 'units'
2434
2435	kextrattrn = ['']
2436
2437	elif filekind == 'LMDZ':
2438	kdimn['x'] = 'x'
2439	kdimn['y'] = 'y'
2440	kdimn['z'] = 'presnivs'
2441	kdimn['t'] = 'time_counter'
2442
2443	kvarn['lon'] = ['lon',['x'],'longitude','Longitude','degrees_east', None]
2444	kvarn['lat'] = ['lat',['y'],'latitude','Latitude','degrees_north', None]
2445	kvarn['lev'] = ['presnivs',['z'],'model_level_number','Vertical levels','Pa', \
2446	None]
2447	kvarn['time'] = ['time_counter',['t'],'time','Time', \
2448	'seconds since 1980-01-01 00:24:00', None]
2449	kvarn['temp'] = ['temp',['t','z','y','x'],'Air temperature','Air temperature','K', \
2450	9.96921e+36]
2451	kvarn['tsfc'] = ['t2m',['t','y','x'],'Temperature 2m','Temperature 2m','K',9.96921e+36]
2452	kvarn['u10m'] = ['u10m',['t','y','x'],'Vent zonal 10m','Vent zonal 10m','m/s',9.96921e+36]
2453	kvarn['v10m'] = ['v10m',['t','y','x'],'Vent meridien 10m','Vent meridien 10m','m/s', \
2454	9.96921e+36]
2455	kvarn['u'] = ['vitu',['t','z','y','x'],'Zonal wind','Zonal wind','m/s',9.96921e+36]
2456	kvarn['v'] = ['vitv',['t','z','y','x'],'Meridional wind','Meridional wind','m/s',9.96921e+36]
2457	kvarn['zsfc'] = ['phis',['t','y','x'],'Surface geop.height','Surface geop.height', \
2458	'm2/s2',9.96921e+36]
2459	kvarn['geopot'] = ['geop',['t','z','y','x'],'Geopotential height','Geopotential height', \
2460	'm2/s2',9.96921e+36]
2461	kvarn['psfc'] = ['psol',['t','y','x'],'Surface Pressure','Surface Pressure','Pa', \
2462	9.96921e+36]
2463	kvarn['pres'] = ['pres',['t','z','y','x'],'Air pressure','Air pressure','Pa',9.96921e+36]
2464	kvarn['H2Ov'] = ['ovap',['t','z','y','x'],'Specific humidity','Specific humidity','kg/kg', \
2465	9.96921e+36]
2466	kvarn['H2Ol'] = ['oliq',['t','z','y','x'],'Condensed water','Condensed water','kg/kg', \
2467	9.96921e+36]
2468
2469	kattrn['standard_name'] = 'standard_name'
2470	kattrn['long_name'] = 'long_name'
2471	kattrn['units'] = 'units'
2472
2473	kextrattrn = ['coordinates']
2474
2475	elif filekind == 'WRF':
2476	kdimn['x'] = 'west_east'
2477	kdimn['y'] = 'south_north'
2478	kdimn['z'] = 'bottom_top'
2479	kdimn['t'] = 'Time'
2480
2481	kvarn['lon'] = ['XLONG',['t','y','x'], None,'LONGITUDE, WEST IS NEGATIVE','degree_east', \
2482	None]
2483	kvarn['lat'] = ['XLAT',['t','y','x'], None,'LATITUDE, SOUTH IS NEGATIVE','degree_north', \
2484	None]
2485	kvarn['lev'] = ['ZNU',['t','z'],None,'eta values on half (mass) levels','',None]
2486	kvarn['time'] = ['Times',['t','DateStrLen'],None,None,None,None]
2487	kvarn['temp'] = ['T',['t','z','y','x'],None,'perturbation potential temperature (theta-t0)','K'\
2488	,None]
2489	kvarn['tsfc'] = ['T2',['t','y','x'],None,'TEMP at 2 M','K',None]
2490	kvarn['u10m'] = ['U10',['t','y','x'],None,'U at 10 M','m s-1',None]
2491	kvarn['v10m'] = ['V10',['t','y','x'],None,'V at 10 M','m s-1',None]
2492	kvarn['u'] = ['U',['t','z','y','x'],None,'x-wind component','m s-1',None]
2493	kvarn['v'] = ['V',['t','z','y','x'],None,'y-wind component','m s-1',None]
2494	kvarn['zsfc'] = [None,['t','y','x'],None,None,None,None]
2495	kvarn['geopot'] = ['PHB',['t','z','y','x'],None,'perturbation geopotential','m2 s-2',None]
2496	kvarn['psfc'] = [None,['t','y','x'],None,None,None,None]
2497	kvarn['pres'] = ['P',['t','z','y','x'],None,'perturbation pressure','Pa',None]
2498	kvarn['H2Ov'] = ['QVAPOR',['t','z','y','x'],None,'Water vapor mixing ratio','kg kg-1',None]
2499	kvarn['H2Ol'] = ['QCLOUD',['t','z','y','x'],None,'Cloud water mixing ratio','kg kg-1',None]
2500
2501	kattrn['standard_name'] = None
2502	kattrn['long_name'] = 'description'
2503	kattrn['units'] = 'units'
2504
2505	kextrattrn = ['FieldType','MemoryOrder','stagger']
2506
2507	# elif filekind == 'OTHER':
2508	# kdimn['x'] = ''
2509	# kdimn['y'] = ''
2510	# kdimn['z'] = ''
2511	# kdimn['t'] = ''
2512
2513	# kvarn['lon'] = ['',['t','z','y','x'],'','','',]
2514	# kvarn['lat'] = ['',['t','z','y','x'],'','','',]
2515	# kvarn['lev'] = ['',['t','z','y','x'],'','','',]
2516	# kvarn['time'] = ['',['t','z','y','x'],'','','',]
2517	# kvarn['temp'] = ['',['t','z','y','x'],'','','',]
2518	# kvarn['tsfc'] = ['',['t','y','x'],'','','',]
2519	# kvarn['u10m'] = ['',['t','y','x'],'','','',]
2520	# kvarn['v10m'] = ['',['t','y','x'],'','','',]
2521	# kvarn['u'] = ['',['t','z','y','x'],'','','',]
2522	# kvarn['v'] = ['',['t','z','y','x'],'','','',]
2523	# kvarn['zsfc'] = ['',['t','y','x'],'','','',]
2524	# kvarn['geopot'] = ['',['t','z','y','x'],'','','',]
2525	# kvarn['psfc'] = ['',['t','y','x'],'','','',]
2526	# kvarn['pres'] = ['',['t','z','y','x'],'','','',]
2527	# kvarn['H2Ov'] = ['',['t','z','y','x'],'','','',]
2528	# kvarn['H2Ol'] = ['',['t','z','y','x'],'','','',]
2529
2530	# kattrn['standard_name'] = ''
2531	# kattrn['long_name'] = ''
2532	# kattrn['units'] = ''
2533
2534	# kextrattrn = ['']
2535	else:
2536	print errormsg
2537	print ' ' + fname + ": filekind '" + filekind + "' not ready !!"
2538	quit(-1)
2539
2540	return kdimn, kvarn, kattrn, kextrattrn
2541
2542	def generic_NetCDFfile(ncobj, dims, kfile, kdimns, kvarns, stdattrns, extrattrns):
2543	""" Function to fill a generic NetCDF file
2544	ncobj= NetCDF file object to which the variables have to be created
2545	kfile= kind of file
2546	'CF': CF-convention
2547	'LMDZ': LMDZ style
2548	'WRF': WRF style
2549	dims= [dimx, dimy, dimz] dimensions of the file
2550	kdimns= dictionary for the specific names for the standard dimensions (x,y,z,t)
2551	kvarns= dictionary for the specific values for the standard variables
2552	kvarn['stdn'] = ['name', 'dims(list)', 'std_name', 'long_name', 'units', 'FillValue']
2553	None: No value
2554	stdattrns= dictionary for the specific values for the variables' standard attributes
2555	extrattrns= list for the specific values for the variables' extra attributes
2556	"""
2557	fname = 'generic_NetCDFile'
2558
2559	# Creation of dimensions
2560	##
2561	newdim = ncobj.createDimension(kdimns['x'], dims[0])
2562	newdim = ncobj.createDimension(kdimns['y'], dims[1])
2563	newdim = ncobj.createDimension(kdimns['z'], dims[2])
2564	newdim = ncobj.createDimension(kdimns['t'], None)
2565
2566	if kfile == 'WRF':
2567	newdim = ncobj.createDimension('DateStrLen', 19)
2568
2569	# Creation of variables
2570	##
2571	for varn in kvarns.keys():
2572	varvals = kvarns[varn]
2573	if varvals[0] is not None:
2574	dimsvar = []
2575	for dimn in varvals[1]:
2576	if dimn == 'DateStrLen':
2577	dimsvar.append(dimn)
2578	else:
2579	dimsvar.append(kdimns[dimn])
2580
2581	if varvals[5] is not None:
2582	nerwvar = ncobj.createVariable(varvals[0], 'f4', tuple(dimsvar), \
2583	fill_value=np.float(varvals[5]))
2584	else:
2585	newvar = ncobj.createVariable(varvals[0], 'f4', tuple(dimsvar))
2586
2587	# Attributes
2588	attrns = stdattrns.keys()
2589	for iattr in range(len(attrns)):
2590	attrn = attrns[iattr]
2591	if stdattrns[attrn] is not None and varvals[4+iattr] is not None:
2592	newvar.setncattr(stdattrns[attrn], varvals[4+iattr])
2593
2594	# Extra attributes
2595	for iattr in range(len(extrattrns)):
2596	attrn = extrattrns[iattr]
2597	if kfile == 'LMDZ':
2598	if attrn == 'coordinates':
2599	attrv = ''
2600	for din in varvals[1]:
2601	attrv = attrv + kdimns[din] + ' '
2602	elif kfile == 'WRF':
2603	if attrn == 'FieldType':
2604	attrv = '104'
2605	elif attrn == 'MemoryOrder':
2606	attrv = ''
2607	Ndims = len(varvals[1])
2608	for idim in range(Ndims-1,0,-1):
2609	if varvals[1][idim] != 't':
2610	attrv = attrv + varvals[1][idim].upper()
2611	elif attrn == 'stagger':
2612	staggeredvars = {}
2613	staggeredvars['U'] = 'X'
2614	staggeredvars['V'] = 'Y'
2615	staggeredvars['PH'] = 'Z'
2616
2617	if ncvar.searchInlist(staggeredvars.keys(),varvals[0]):
2618	attrv = staggeredvars[varvals[0]]
2619	else:
2620	attrv = ''
2621
2622	newvar.setncattr(extrattrns[iattr], attrv)
2623
2624	ncobj.sync()
2625
2626	return
2627
2628	####### ###### ##### #### ### ## #
2629
2630	filekindsnames = "'" + ncvar.numVector_String(filekinds, "', '") + "'"
2631
2632	parser = OptionParser()
2633	parser.add_option("-o", "--outputkind", type='choice', dest="okind",
2634	choices=filekinds, help="kind of output: " + filekindsnames, metavar="KIND")
2635	parser.add_option("-d", "--dimensions", dest="dims",
2636	help="dimensions of the initial conditions: dimx,dimy,dimz", metavar="VALUES")
2637	parser.add_option("-p", "--pressure_exner", dest="pexner",
2638	help="how as to b computed Exner pressure ('hybrid': hybrid coordinates, 'middle': middle layer)",
2639	metavar="VALUE")
2640	parser.add_option("-q", "--NWaterSpecies", dest="nqtot",
2641	help="Number of water species", metavar="VALUE")
2642	parser.add_option("-t", "--time", dest="time",
2643	help="time of the initial conditions ([YYYY][MM][DD][HH][MI][SS] format)", metavar="VALUE")
2644	parser.add_option("-z", "--z_levels", type='choice', dest="znivs",
2645	choices=['param', 'tropo', 'strato', 'read'],
2646	help="which kind of vertical levels have to be computed ('param', 'tropo', 'strato', 'read')",
2647	metavar="VAR")
2648
2649	(opts, args) = parser.parse_args()
2650
2651	####### #######
2652	## MAIN
2653	#######
2654
2655	# dynamic variables
2656	# vcov, ucov: covariant winds
2657	# teta: potential temperature
2658	# q: advecting fields (humidity species)
2659	# ps: surface pressure
2660	# masse: air mass
2661	# phis: surface geopotential
2662
2663	# Local:
2664	# p: pressure at the interface between layers (half-sigma)
2665	# pks: Exner function at the surface
2666	# pk: Exner functino at the half-sigma layers
2667	# pkf: Filtred Exner function at the half-sigma layers
2668	# phi: geopotential height
2669	# ddsin,zsig,tetapv,w_pv: auxiliar variables
2670	# tetastrat: potential temporeature in the stratosphere (K)
2671	# teta0,ttp,delt_y,delt_z,eps: constants for the T profile
2672	# k_f,k_c_a,k_c_s: calling constants
2673	# ok_geost: Initialisation geostrohic wind
2674	# ok_pv: Polar Vortex
2675	# phi_pv,dphi_pv,gam_pv: polar vortex constants
2676
2677	dimx = int(opts.dims.split(',')[0])
2678	dimy = int(opts.dims.split(',')[1])
2679	dimz = int(opts.dims.split(',')[2])
2680	nqtot = int(opts.nqtot)
2681
2682	ofile = 'iniaqua.nc'
2683
2684	print 'dimensions: ',dimx,', ',dimy,', ',dimz
2685
2686	if opts.okind == 'CF':
2687	varnames = ['zg', 'ta', 'ua', 'va', 'pres', 'r']
2688	timev = float(opts.time)
2689	# Reference time from 1949-12-01 00:00:00 UTC
2690	timev = ncvar.datetimeStr_conversion(opts.time,'YmdHMS','cfTime')
2691	dimnames = ['time', 'z', 'y', 'x']
2692	elif opts.okind == 'WRF':
2693	varnames = ['Geopot', 'T', 'U', 'V', 'PRES', 'QVAPOR']
2694	timev = ncvar.datetimeStr_conversion(opts.time,'YmdHMS','WRFdatetime')
2695	dimnames = ['Time', 'bottom_top', 'south_north', 'west_east']
2696
2697	# Constants
2698	##
2699	llm = dimz
2700
2701	ok_geost = True
2702	# Constants for Newtonian relaxation and friction
2703	k_f = 1.
2704	k_f = 1./(daysec*k_f)
2705	# cooling surface
2706	k_c_s=4.
2707	k_c_s=1./(daysec*k_c_s)
2708	# cooling free atm
2709	k_c_a=40.
2710	k_c_a=1./(daysec*k_c_a)
2711	# Constants for Teta equilibrium profile
2712	# mean Teta (S.H. 315K)
2713	teta0=315.
2714	# Tropopause temperature (S.H. 200K)
2715	ttp=200.
2716	# Deviation to N-S symmetry(~0-20K)
2717	eps=0.
2718	# Merid Temp. Gradient (S.H. 60K)
2719	delt_y=60.
2720	# Vertical Gradient (S.H. 10K)
2721	delt_z=10.
2722	# Polar vortex
2723	ok_pv = False
2724	# Latitude of edge of vortex
2725	phi_pv=-50.
2726	phi_pv=phi_pv*pi/180.
2727	# Width of the edge
2728	dphi_pv=5.
2729	dphi_pv=dphi_pv*pi/180.
2730	# -dT/dz vortex (in K/km)
2731	gam_pv=4.
2732
2733	# For extra-terrestrial planets
2734	#presnivs, pseudoalt = presnivs_calc_noterre(opts.znivs, dimz)
2735	presnivs, ap, bp = presnivs_calc(opts.znivs, dimz)
2736	lon, lat = global_lonlat(dimy,dimx)
2737	lonu, latu = global_lonlat(dimy,dimx+1)
2738	lonv, latv = global_lonlat(dimy+1,dimx)
2739
2740	# 2. Initialize fields towards which to relax
2741	##
2742
2743	knewt_t = np.zeros((dimz), dtype=np.float)
2744	kfrict = np.zeros((dimz), dtype=np.float)
2745	clat4 = np.zeros((dimy+1, dimx+1), dtype=np.float)
2746
2747	# Friction
2748	knewt_g = k_c_a
2749	for l in range(dimz):
2750	zsig=presnivs[l]/preff
2751	knewt_t[l]=(k_c_s-k_c_a)*np.max([0.,(zsig-0.7)/0.3])
2752	kfrict[l]=k_f*np.max([0.,(zsig-0.7)/0.3])
2753
2754	for j in 1,range(dimy+1):
2755	clat4[j,:]=np.cos(latu[j,0])**4
2756
2757	# Vertical profile
2758	tetajl = np.zeros((dimz, dimy+1, dimx), dtype=np.float)
2759	#theta = np.zeros((dimy+1, dimx+1), dtype=np.float)
2760	#thetarappel = np.zeros((dimz, dimy+1, dimx+1), dtype=np.float)
2761	theta = np.zeros((dimy, dimx), dtype=np.float)
2762	thetarappel = np.zeros((dimz, dimy, dimx), dtype=np.float)
2763
2764	for l in range (dimz):
2765	zsig = presnivs[l]/preff
2766	tetastrat = ttpzsig*(-kappa)
2767	tetapv = tetastrat
2768	if ok_pv and zsig < 0.1:
2769	tetapv = tetastrat(zsig10.)*(kappacpp*gam_pv/1000./g)
2770
2771	ddsin = np.sin(latu)
2772	tetajl[l,:,:] = teta0-delt_yddsin+epsddsin-delt_z(1.-ddsinddsin)*np.log(zsig)
2773	if planet_type == 'giant':
2774	tetajl[l,:,:] = teta0+(delt_y((np.sin(latu3.14159eps+0.0001))*2) / \
2775	((latu3.14159eps+0.0001)*2))-delt_znp.log(zsig)
2776
2777	# Profile stratospheric isotherm (+vortex)
2778	for iy in range(dimy):
2779	for ix in range(dimx):
2780	w_pv=(1.-np.tanh((latu[iy,ix]-phi_pv)/dphi_pv))/2.
2781	tetastrat=tetastrat(1.-w_pv)+tetapvw_pv
2782	tetajl[l,iy,ix]=np.max([tetajl[l,iy,ix],tetastrat])
2783
2784	#for iz in range(dimz):
2785	# for iy in range(dimy+1):
2786	# tetarappel[iz,iy,0:dimx] = tetajl[iz,iy,:]
2787	#
2788	# tetarappel[iz,iy,0:dimx] = tetajl[iz,iy,dimx-1]
2789	thetarappel = tetajl.copy()
2790
2791	# 3. Initialize fields (if necessary)
2792	# surface pressure
2793
2794	if iflag_phys > 2:
2795	# specific value for CMIP5 aqua/terra planets
2796	# "Specify the initial dry mass to be equivalent to
2797	# a global mean surface pressure (101325 minus 245) Pa."
2798	press = np.ones((dimy+1, dimx+1), dtype=np.float)*101080.
2799	else:
2800	# use reference surface pressure
2801	press = np.ones((dimy+1, dimx+1), dtype=np.float)*preff
2802
2803	# ground geopotential
2804	phiss = np.zeros((dimy+1, dimx+1), dtype=np.float)
2805
2806	pres = pression(dimx, dimy, dimz, ap, bp, press)
2807
2808	aire, apolnorth, apolsouth, airesurge, rlatitudeu, rlatitudev, cuwind, cvwind = \
2809	inigeom(dimx, dimy)
2810
2811	if opts.pexner == 'hybdrid':
2812	pks, pk, pkf, alpha, beta = exner_hyb(dimx, dimy, dimz, press, pres, aire, \
2813	apolnorth, apolsouth)
2814	else:
2815	pks, pk, pkf = exner_milieu(dimx, dimy, dimz, press, pres, beta)
2816
2817	masse = massdair(dimx,dimy,dimz,pres,airesurge)
2818
2819	# bulk initialization of temperature
2820	theta = thetarappel.copy()
2821
2822	# geopotential
2823	phisfc = np.zeros((dimy+1, dimx+1), dtype=np.float)
2824	phiall = np.zeros((dimz+1, dimy+1, dimx+1), dtype=np.float)
2825
2826	phisfc, phiall = geopot(dimx,dimy,dimz,theta,pk,pks)
2827
2828	# winds
2829	ucov = np.zeros((dimz, dimy, dimx), dtype=np.float)
2830	vcov = np.zeros((dimz, dimy, dimx), dtype=np.float)
2831
2832	if ok_geost:
2833	ucov = ugeostr(dimx,dimy,dimz,phisfc,phiall,rlatitudeu,rlatitudev,cuwind)
2834
2835	# bulk initialization of tracers
2836	q = np.zeros((dimz, dimy, dimx, nqtot), dtype=np.float)
2837
2838	if planet_type == 'earth':
2839	# Earth: first two tracers will be water
2840	for i in range(nqtot):
2841	if i == 1: q[:,:,i] = 1.e-10
2842	if i == 2: q[:,:,i] = 1.e-15
2843	if i > 2: q[:,:,i] = 0.
2844
2845	# add random perturbation to temperature
2846	idum = -1
2847	zz = RAN1(idum,97)
2848	idum = 0
2849	for l in range(dimz):
2850	for j in range(dimy):
2851	for i in range(dimx):
2852	theta[l,j,i] = theta[l,j,i](1.+0.005RAN1(idum,97))
2853
2854	# maintain periodicity in longitude
2855	for l in range(dimz):
2856	for j in range(1,dimy):
2857	theta[l,j,dimx-1]=theta[l,j-1,dimx-1]
2858
2859	ncf = NetCDFFile(ofile, 'w')
2860
2861	# File structure creation
2862	filedimns, filevarns, filevarattr, filevarxtrattr = name_variables(opts.okind)
2863	generic_NetCDFfile(ncf,[dimx,dimy,dimz], opts.okind, filedimns, filevarns, \
2864	filevarattr, filevarxtrattr)
2865
2866	# File filling
2867	for varn in ncf.variables.keys():
2868	varobj = ncf.variables[varn]
2869	if varn == filevarns['lon'][0]: varobj[:] = lon
2870	elif varn == filevarns['lat'][0]: varobj[:] = lat
2871	elif varn == filevarns['lev'][0]: varobj[:] = presnivs
2872	elif varn == filevarns['time'][0]: varobj[:] = timev
2873	elif varn == filevarns['temp'][0]:
2874	print 'Lluis shapes varobj:',varobj.shape,'theta:',theta.shape
2875	varobj[0,:,:,:] = theta[:]
2876	elif varn == filevarns['tsfc'][0]: varobj[:] = theta[0,:,:]
2877	elif varn == filevarns['u10m'][0]: varobj[:] = ucov[0,:,:]
2878	elif varn == filevarns['v10m'][0]: varobj[:] = vcov[0,:,:]
2879	elif varn == filevarns['u'][0]: varobj[0,:,:,:] = ucov[:]
2880	elif varn == filevarns['v'][0]: varobj[0,:,:,:] = vcov[:]
2881	elif varn == filevarns['zsfc'][0]: varobj[0,:,:] = phisfc[:]
2882	elif varn == filevarns['geopot'][0]: varobj[0,:,:] = phisall[:]
2883	elif varn == filevarns['psfc'][0]: varobj[0,:,:] = press[:]
2884	elif varn == filevarns['pres'][0]: varobj[0,:,:,:] = pres[:]
2885	elif varn == filevarns['H2Ov'][0]: varobj[0,:,:,:] = q[:,:,:0]
2886	elif varn == filevarns['H2Ol'][0]: varobj[0,:,:,:] = q[:,:,:1]
2887
2888	ncf.sync()
2889	ncf.close()
2890
2891	print main + ": successfull writing of file '" + ofile + "' !!"

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: