Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

iniaqua.py @ 1639

Last change on this file since 1639 was 414, checked in by lfita, 10 years ago
Removing WRFV3 folder and keeping only the essential
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/branches/LMDZ_WRFmeas/tools/iniaqua.py @ 1639

Download in other formats: