Context Navigation

iniaqua.py @ 463

Last change on this file since 463 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

Rev	Line
[180]	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
[216]	19	filekinds = ['CF', 'LMDZ', 'WRF']
[180]	20
[224]	21	## e.g. # iniaqua.py -d 32,32,39 -p hybdrid -o WRF -t 19791201000000 -z tropo -q 2
[214]	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
[180]	82
[214]	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	#
[216]	245	xlon = np.zeros((dx+1), dtype=np.float32)
[214]	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
[215]	1397	return aire, apoln, apols, airesurg, rlatu, rlatv, cu, cv
[214]	1398
[212]	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
[214]	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
[212]	1414
	1415	return ssumv
	1416
[215]	1417	def SCOPY(ddx,ddy,ddz,sx,incx,sy,incy):
[212]	1418	""" Obsolete function to copy matrix values
	1419	from dyn3d/cray.F
	1420	"""
[215]	1421	fname = 'SCOPY'
[212]	1422
[215]	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
[212]	1438
	1439	return sy
	1440
[215]	1441	def exner_hyb (dx, dy, dz, psv, pv, aire, apoln, apols):
[212]	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
[215]	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)
[212]	1474
[215]	1475	ppn = np.zeros((dy+1, dx+1), dtype=np.float)
	1476	pps = np.zeros((dy+1, dx+1), dtype=np.float)
[214]	1477
[215]	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)
[214]	1480
[215]	1481	if dz == 1:
[212]	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
[213]	1494	#### General case:
[212]	1495
	1496	unpl2k = 1.+ 2.* kappa
	1497
[215]	1498	for j in range(dy+1):
	1499	for i in range(dx+1):
	1500	pksv[j,i] = cpp * ( psv[j,i]/preff ) ** kappa
[212]	1501
[215]	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
[214]	1510	for j in range(dy):
	1511	for i in range(dx):
[215]	1512	pksv[j,i] = xpn[i]
	1513	pksv[dy-1,i] = xps[i]
[212]	1514	#
	1515	#
	1516	# .... Calcul des coeff. alpha et beta pour la couche l = llm ..
	1517	#
[215]	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
[212]	1522
	1523	#
	1524	# ... Calcul des coeff. alpha et beta pour l = llm-1 a l = 2 ...
	1525	#
[215]	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
[213]	1532
[212]	1533	# ***********************************************************************
	1534	# ..... Calcul de pk pour la couche 1 , pres du sol ....
	1535	#
[215]	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] )
[212]	1540
	1541	#
	1542	# ..... Calcul de pk(ij,l) , pour l = 2 a l = llm ........
	1543	#
[215]	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]
[212]	1548	#
	1549	#
[215]	1550	pkfv = SCOPY ( dx+1, dy+1, dz, pkv, 1, pkfv, 1 )
[212]	1551
	1552	# We do not filter for iniaqua
	1553	# CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
	1554
[215]	1555	return pksv, pkv, pkfv, alphav, betav
[212]	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
[213]	1588	ppn = np.zeros((iim), dtype=np.float)
	1589	ppn = np.zeros((iim), dtype=np.float)
[212]	1590
[214]	1591	ip1jm = (dx+1)*dy
	1592
[212]	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:
[213]	1612	if llm == 1:
[214]	1613
[212]	1614	# Compute pks(:),pk(:),pkf(:)
	1615
[215]	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]
[212]	1619
[215]	1620	pkf = SCOPY(dx,dy,dz, pk, 1, pkf, 1 )
[212]	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
[213]	1625	return pksv, pkv, pkfv
[212]	1626
	1627	#### General case:
	1628
	1629	# -------------
	1630	# Calcul de pks
	1631	# -------------
	1632
[215]	1633	for j,i in range(ngrid):
	1634	pks[j,i] = cpp * ( ps[j,i]/preff ) ** kappa
[212]	1635
[215]	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]
[212]	1639
	1640	xpn = SSUM(iim,ppn,1) /apoln
	1641	xps = SSUM(iim,pps,1) /apols
	1642
[215]	1643	for j,i in range(iip1):
	1644	pks[j,i] = xpn
	1645	pks[j,i+ip1jm] = xps
[212]	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):
[215]	1654	for j,i in range(ngrid):
	1655	pk[j,i,l] = dum1 * (p[j,i,l] + p[j,i,l+1])**kappa
[212]	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
[215]	1661	for j,i in range(ngrid):
	1662	pk[j,i,llm] = pk[j,i,llm-1]**2 / pk[j,i,llm-2]
[212]	1663
	1664	# calcul de pkf
	1665	# -------------
[215]	1666	pkf = SCOPY( dx,dy,dz, pk, 1, pkf, 1 )
[212]	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):
[215]	1674	for j,i in range(ngrid):
	1675	beta[j,i,l] = pk[j,i,l] / pk[j,i,l-1]
[212]	1676
	1677	return pksv, pkv, pkfv
	1678
[213]	1679	def pression(dx, dy, dz, apv, bpv, psv):
[212]	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
[215]	1693	press = np.zeros((dz+1, dy+1, dx+1), dtype=np.float)
[212]	1694
	1695	for l in range(dz+1):
[214]	1696	press[l,:,:] = apv[l] + bpv[l]*psv
[212]	1697
	1698	return press
	1699
[180]	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
[212]	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)
[214]	1770	apv = np.zeros((dz+1), dtype=np.float)
	1771	bpv = np.zeros((dz+1), dtype=np.float)
[212]	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
[214]	1814	bpv[0:dz] = np.exp(1.-1./sig[0:dz]**2)
	1815	bpv[llmp1-1] = 0.
[212]	1816
[214]	1817	apv = pa * (sig - bp)
[212]	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
[214]	1829	bpv[0]=1.
	1830	bpv[1:dz] = np.exp(1.-1./sig[1:dz]**2)
	1831	bpv[llmp1-1] = 0.
[212]	1832
[214]	1833	apv[0] = 0.
	1834	apv[1:dz+1] = pa*(sig[1:dz+1]-bpv[1:dz+1])
[212]	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
[214]	1857	bpv[0] = 1.
	1858	bpv[1:dz] = np.exp(1.-1./sig[1:dz]**2)
	1859	bpv[llmp1-1] = 0.
[212]	1860
[214]	1861	apv[0] = 0.
	1862	apv[1:dz+1] = pa*(sig[1:dz+1] - bpv[1:dz+1])
[212]	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())
[214]	1877	apv[l] = np.float(values[0])
	1878	bpv[l] = np.float(values[0])
[212]	1879
	1880	sfobj.close()
[214]	1881	if apv[0] == 0. or apv[dz+1] == 0. or bpv[0] == 1. or bpv[dz+1] == 0.:
[212]	1882	print errormsg
	1883	print ' ' + fname + ': bad ap or bp values !!'
	1884	print ' k ap bp ___________'
	1885	for k in range(dz+1):
[214]	1886	print k, apv[k], bpv[k]
[212]	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):
[214]	1899	print k, bpv[k], apv[k]
[212]	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):
[214]	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 )
[212]	1909	print ' PRESNIVS(', l, ')=', pnivs[l], ' Z ~ ', \
	1910	np.log(preff/pnivs[l])*scaleheight, ' DZ ~ ', \
[214]	1911	scaleheightnp.log((apv[l]+bpv[l]preff)/ \
	1912	np.max([apv[l+1]+bpv[l+1]*preff, 1.e-10]))
[212]	1913
	1914	print ' ' + fname + ': PRESNIVS [Pa]:', pnivs
	1915
[214]	1916	return pnivs, apv, bpv
[212]	1917
	1918
	1919	def presnivs_calc_noterre(hybrid, dz):
[180]	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
[212]	1925	fname = 'presnivs_calc_noterre'
[180]	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)
[214]	1938	apv = np.zeros((dz+1), dtype=np.float)
	1939	bpv = np.zeros((dz+1), dtype=np.float)
[180]	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)
[214]	2037	bpv[l] = np.exp(1.-1./(newsig**2))
	2038	apv[l] = pa * (newsig - bp[l] )
[180]	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):
[214]	2049	apv[l] = 0.
	2050	bpv[l] = sig[l]
[180]	2051
[214]	2052	apv[llmp1-1] = 0.
[180]	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):
[214]	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])
[180]	2070
	2071	if hybrid:
	2072	aps[dz-1] = aps[dz-2]**2 / aps[dz-3]
[214]	2073	bps[dz-1] = 0.5*(bpv[dz-1] + bpv[dz])
[180]	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
[215]	2115	def massdair(dx,dy,dz,p,airesurg):
[212]	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
[215]	2132	masse = np.zeros((dz+1,dy+1,dx+1), dtype=np.float)
[212]	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
[215]	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] )
[212]	2187
[215]	2188	for j in range(1,dy):
	2189	masse[l,j,dx] = masse[l,j-1,dx]
[212]	2190
	2191	return masse
	2192
[215]	2193	def geopot(dx, dy, dz, teta, pk, pks):
[212]	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
[215]	2215	phis = np.zeros((dy+1,dx+1), dtype=np.float)
	2216	phi = np.zeros((dz+1,dy+1,dx+1), dtype=np.float)
[212]	2217
[215]	2218	print ' Lluis in ' + fname + ' shapes phis:',phis.shape,'phi:',phi.shape, \
	2219	'teta:',teta.shape,'pks:',pks.shape,'pk:',pk.shape
	2220
[212]	2221	#-----------------------------------------------------------------------
	2222	# calcul de phi au niveau 1 pres du sol .....
	2223
[215]	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] )
[212]	2227
	2228	# calcul de phi aux niveaux superieurs .......
	2229
[215]	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] )
[212]	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
[215]	2274	def ugeostr(dx,dy,dz,phis,phi,rlatu,rlatv,cu):
[212]	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'
[215]	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)
[212]	2286
[215]	2287	print ' Lluis in ' + fname + ': shapes phis:',phis.shape,'phi:',phi.shape,'u:',u.shape
	2288
	2289	for j in range(dy):
[212]	2290	if np.abs(np.sin(rlatv[j])) < 1.e-4:
	2291	zlat = 1.e-4
	2292	else:
[215]	2293	zlat=rlatv[j]
[212]	2294
	2295	fact = np.cos(zlat)
	2296	fact = fact*fact
	2297	fact = fact*fact
	2298	fact = fact*fact
[215]	2299	fact = (1.-fact)/ (2.omegnp.sin(zlat)*(rlatu[j+1]-rlatu[j]))
[212]	2300	fact = -fact/rad
[215]	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)
[212]	2305
[215]	2306	um = dump2d(dz,dy,'Vent-u geostrophique')
[212]	2307
	2308	# calcul des champ de vent:
	2309
[215]	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]
[212]	2317
[215]	2318	ucov[l,j,dx]=ucov[l,j,0]
[212]	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
[215]	2343	IFF = 0
[212]	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
[213]	2363	if J > Nvals or J < 1: quit()
[212]	2364	ran1=R[J]
	2365	R[J]=(np.float(IX1)+np.float(IX2)RM2)RM1
	2366
	2367	return ran1
	2368
[216]	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'
[212]	2378
[216]	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	"""
[224]	2557	fname = 'generic_NetCDFile'
[216]	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]
[224]	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])
[216]	2580
[224]	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))
[216]	2586
	2587	# Attributes
[224]	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])
[216]	2593
	2594	# Extra attributes
[224]	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()
[216]	2625
	2626	return
	2627
[180]	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")
[213]	2637	parser.add_option("-p", "--pressure_exner", dest="pexner",
[212]	2638	help="how as to b computed Exner pressure ('hybrid': hybrid coordinates, 'middle': middle layer)",
	2639	metavar="VALUE")
[215]	2640	parser.add_option("-q", "--NWaterSpecies", dest="nqtot",
	2641	help="Number of water species", metavar="VALUE")
[180]	2642	parser.add_option("-t", "--time", dest="time",
	2643	help="time of the initial conditions ([YYYY][MM][DD][HH][MI][SS] format)", metavar="VALUE")
[212]	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")
[180]	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])
[215]	2680	nqtot = int(opts.nqtot)
[180]	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']
[216]	2688	timev = float(opts.time)
[180]	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
[212]	2733	# For extra-terrestrial planets
	2734	#presnivs, pseudoalt = presnivs_calc_noterre(opts.znivs, dimz)
[214]	2735	presnivs, ap, bp = presnivs_calc(opts.znivs, dimz)
[212]	2736	lon, lat = global_lonlat(dimy,dimx)
	2737	lonu, latu = global_lonlat(dimy,dimx+1)
	2738	lonv, latv = global_lonlat(dimy+1,dimx)
[180]	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)
[212]	2745	clat4 = np.zeros((dimy+1, dimx+1), dtype=np.float)
[180]	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):
[212]	2755	clat4[j,:]=np.cos(latu[j,0])**4
[180]	2756
	2757	# Vertical profile
[212]	2758	tetajl = np.zeros((dimz, dimy+1, dimx), dtype=np.float)
[224]	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)
[180]	2763
[212]	2764	for l in range (dimz):
[180]	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) / \
[212]	2775	((latu3.14159eps+0.0001)*2))-delt_znp.log(zsig)
[180]	2776
	2777	# Profile stratospheric isotherm (+vortex)
[212]	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])
[180]	2783
[212]	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]
[215]	2789	thetarappel = tetajl.copy()
[180]	2790
[212]	2791	# 3. Initialize fields (if necessary)
	2792	# surface pressure
[180]	2793
[212]	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."
[215]	2798	press = np.ones((dimy+1, dimx+1), dtype=np.float)*101080.
[212]	2799	else:
	2800	# use reference surface pressure
[215]	2801	press = np.ones((dimy+1, dimx+1), dtype=np.float)*preff
[212]	2802
	2803	# ground geopotential
[215]	2804	phiss = np.zeros((dimy+1, dimx+1), dtype=np.float)
[180]	2805
[215]	2806	pres = pression(dimx, dimy, dimz, ap, bp, press)
[212]	2807
[215]	2808	aire, apolnorth, apolsouth, airesurge, rlatitudeu, rlatitudev, cuwind, cvwind = \
	2809	inigeom(dimx, dimy)
[214]	2810
[212]	2811	if opts.pexner == 'hybdrid':
[215]	2812	pks, pk, pkf, alpha, beta = exner_hyb(dimx, dimy, dimz, press, pres, aire, \
	2813	apolnorth, apolsouth)
[212]	2814	else:
[216]	2815	pks, pk, pkf = exner_milieu(dimx, dimy, dimz, press, pres, beta)
[212]	2816
[215]	2817	masse = massdair(dimx,dimy,dimz,pres,airesurge)
[212]	2818
	2819	# bulk initialization of temperature
[215]	2820	theta = thetarappel.copy()
[212]	2821
	2822	# geopotential
[215]	2823	phisfc = np.zeros((dimy+1, dimx+1), dtype=np.float)
	2824	phiall = np.zeros((dimz+1, dimy+1, dimx+1), dtype=np.float)
[212]	2825
[215]	2826	phisfc, phiall = geopot(dimx,dimy,dimz,theta,pk,pks)
[212]	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:
[215]	2833	ucov = ugeostr(dimx,dimy,dimz,phisfc,phiall,rlatitudeu,rlatitudev,cuwind)
[212]	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
[215]	2847	zz = RAN1(idum,97)
[212]	2848	idum = 0
[215]	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))
[212]	2853
	2854	# maintain periodicity in longitude
[215]	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]
[212]	2858
[180]	2859	ncf = NetCDFFile(ofile, 'w')
	2860
[216]	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)
[180]	2865
[216]	2866	# File filling
	2867	for varn in ncf.variables.keys():
	2868	varobj = ncf.variables[varn]
[224]	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]
[216]	2887
[180]	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 @ 463

Download in other formats: