Context Navigation

source: trunk/forcagequimarche.pro @ 30

Last change on this file since 30 was 30, checked in by pinsard, 17 years ago
add step2_diff.pro (not ok yes) and restore divfred call
File size: 11.4 KB

Rev	Line
[30]	1	;+
[28]	2	;
	3	; @file_comments
	4	; calcul d'une matrice 3D de conductivité sigma=f(T,S) en fonction de T et S
	5	; de ORCA, sur la grille T
	6	;
	7	; @categories
	8	; orca grid
	9	;
	10	; @param orcares {in}{required}{type=string}
	11	; code of ORCA grid to use for output results
	12	; must be 'ORCA2'
	13	;
	14	; @param iyear {in}{required}{type=integer}
	15	; real year
	16	;
	17	; @param ian {in}{required}{type=integer}
	18	; simulation year
	19	;
	20	; @keyword DRAKKAR_EXP {type=string}
	21	; code for Drakkar experiment
	22	; only used when orcares = ORCA025
	23	; must be G42 ++ G70
	24	;
	25	; @restrictions
	26	; - Requires SAXO tools
	27	;
	28	; @todo
	29	; validation of ORCA2 processing (now some nan in output files !)
	30	; add ORCA025
	31	; provide tools to plot output files
	32	; produce a NetCDF GDT or CF compliant
[30]	33	; add divchoice parameter to choose either div or divfred
[28]	34	;
	35	; @pre
	36	; see geomag_profile.sh
	37	; be sure to have ++ in the directory defined in ${GEOMAG_ID}/
	38	;
	39	; @post
	40	; see geomag_profile.sh
	41	; ++ in the directory defined in ${GEOMAG_OD}/
	42	;
	43	; @examples
	44	; for ORCA2
	45	; IDL> .run forcagequimarche
	46	; IDL> forcagequimarche, 'ORCA2', '1993', '01'
	47	;
	48	; @history
[30]	49	; fplod 2007-07-30T10:12:38Z aedon.locean-ipsl.upmc.fr (Darwin)
	50	; restore divfred call instead of div ("official" saxo divergence module)
[28]	51	; fplod 2007-07-27T10:29:18Z cerbere.locean-ipsl.upmc.fr (Linux)
[29]	52	; add orcares in the name of output files
	53	; fplod 2007-07-27T10:29:18Z cerbere.locean-ipsl.upmc.fr (Linux)
[28]	54	; add orcares parameter
	55	; add optional keyword drakkar_exp
	56	; add header for idldoc
	57	; add check of io directories
	58	; replace io directories user dependant to $GEOMAG_ID and $GEOMAG_OD
	59	; replace call of divfred by call of div
	60	;
	61	; @version
	62	; $Id$
	63	;
	64	;-
	65	;
	66	PRO forcagequimarche,orcares,iyear,ian,DRAKKAR_EXP = drakkar_exp
	67	;
	68	;----
	69	; check input parameters
	70	;----
	71	;
	72	; check orcares definition
	73	;
	74	CASE orcares OF
	75	'ORCA2': BEGIN
	76	msg = 'iii : valid orcares parameter = ' + orcares
	77	ras = report(msg)
	78	filename_oce = 'meshmask_bab.nc'
	79	IF keyword_set(DRAKKAR_EXP) THEN BEGIN
	80	msg = 'www : unused DRAKKAR_EXP keyword = ' + drakkar_exp
	81	ras = report(msg)
	82	END
	83	END
	84	ELSE : BEGIN
	85	msg = 'eee : invalid orcares parameter = ' + orcares
	86	ras = report(msg)
	87	msg = 'eee : orcares must be ORCA2 or ORCA025++'
	88	ras = report(msg)
	89	RETURN
	90	END
	91	ENDCASE
[2]	92
[28]	93	; ++ check iyear
	94	; ++ check ian
[2]	95
[28]	96	;++@init2
	97	; init grid
	98	initocemeshmask, orcares, DRAKKAR_EXP = drakkar_exp
[2]	99
	100	@common
[28]	101	;
	102	; check for input files
	103	;
	104	; test if ${GEOMAG_ID} defined
	105	geomag_id_env=GETENV('GEOMAG_ID')
	106	CASE geomag_id_env OF
	107	'' : BEGIN
	108	msg = 'eee : ${GEOMAG_ID} is not defined'
	109	ras = report(msg)
	110	RETURN
	111	END
	112	ELSE: BEGIN
	113	msg = 'iii : ${GEOMAG_ID} is ' + geomag_id_env
	114	ras = report(msg)
	115	END
	116	ENDCASE
	117	;
	118	iodirin = isadirectory(geomag_id_env)
	119	;
	120	; existence and protection of ${GEOMAG_ID}
	121	IF (FILE_TEST(iodirin, /DIRECTORY, /EXECUTABLE, /READ) EQ 0) THEN BEGIN
	122	msg = 'eee : the directory' + iodirin + ' is not accessible.'
	123	ras = report(msg)
	124	RETURN
	125	ENDIF
	126	;
	127	; test if ${GEOMAG_OD} defined
	128	geomag_od_env=GETENV('GEOMAG_OD')
	129	CASE geomag_od_env OF
	130	'' : BEGIN
	131	msg = 'eee : ${GEOMAG_OD} is not defined'
	132	ras = report(msg)
	133	RETURN
	134	END
	135	ELSE: BEGIN
	136	msg = 'iii : ${GEOMAG_OD} is ' + geomag_od_env
	137	ras = report(msg)
	138	END
	139	ENDCASE
	140	;
	141	; check if output data will be possible
	142	iodirout = isadirectory(geomag_od_env)
	143	;
	144	; existence and protection
	145	IF (FILE_TEST(iodirout, /DIRECTORY, /WRITE) EQ 0) THEN BEGIN
	146	msg = 'eee : the directory' + iodirout + ' was not found.'
	147	ras = report(msg)
	148	RETURN
	149	ENDIF
	150	;
[2]	151	e_exp='ESS'
	152	key_portrait = 0
[28]	153	; stockage des fichiers brut
	154	ioDATA=geomag_id_env
[2]	155
[28]	156	CASE orcares OF
	157	'ORCA2': BEGIN
	158	file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc'
	159	file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc'
	160	file_T=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_T.nc'
	161	file_Sed= geomag_id_env+'cond_sed_ORCA2.nc'
	162	file_Br= geomag_id_env+'Br_ORCA2.nc'
	163	END
	164	ENDCASE
	165	;
[2]	166	; title
	167	t_exp= e_exp
	168	t_bt = 'bar_transp'
[28]	169	ioORLN2 = geomag_id_env
	170	;
[2]	171	;facteur d'echelle vertical for partial steps
	172	e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
	173	e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
	174	SIGMAsed=read_ncdf('cond_sed',0,/timestep,/nostruct,/tout,filename=file_Sed,/cont_nofill)
	175	;BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br)
	176	BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br,/cont_nofill)
[28]	177	;
[2]	178	; vertical integration:
[28]	179	e3t3d=make_array(jpi,jpj,jpk)
[2]	180	for k=0, jpk-1 do begin &$
[28]	181	for j=0,jpj-1 do begin &$
[2]	182	for i=0,jpi-1 do begin &$
	183	e3t3d(i,j,k) = e3t(k) &$
	184	endfor &$
	185	endfor &$
	186	endfor
	187	jpt = 73
[28]	188	;
	189	;vud = make_array(jpi,jpj,jpt)
[2]	190	;vvd = make_array(jpi, jpj, jpt)
	191	divBustar = make_array(jpi, jpj, jpt)
	192	diver3=replicate(0.,182,149,1,73)
	193	sigma3=replicate(0.,182,149,1,73)
[28]	194	;
	195	FOR jt = 0, jpt-1 DO BEGIN &$
	196	;
[2]	197	; ouverture des fichiers et stockage en memoire partial steps
[28]	198	vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U) &$
[2]	199	;stop
[28]	200	vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V) &$
[2]	201	;stop
	202	; lecture salinite & temperature
	203	temp= read_ncdf('votemper',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T)
	204	;stop
	205	salin=read_ncdf('vosaline',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T)
	206	;stop
	207	conduct=0.02047780622061 + 0.00273147624197temp + 0.00035133182334temptemp + 0.09139808809909salin + 0.00241425798890salintemp -0.00023998958774salinsalin
	208	mask_t=where(conduct GT 1.e+19)
	209	conduct(mask_t)=0.
	210	; Somme conduct au point T
[28]	211	;
[2]	212	; Calcul SIGMA
[28]	213	;
[2]	214	SIGMAoc=total(conducte3t3dtmask,3)
	215	SIGMA=SIGMAsed+SIGMAoc
[28]	216	;
[2]	217	SIGMA_u=(SIGMA+shift(SIGMA,-1,0))/2.
	218	SIGMA_v=(SIGMA+shift(SIGMA,0,1))/2.
[28]	219	;
[2]	220	; Calcul B en points u et v
[28]	221	;
[2]	222	BR_u=(BR+shift(BR,-1,0))/2.
	223	BR_v=(BR+shift(BR,0,1))/2.
[28]	224	;
[2]	225	; Calcul integrale conduct
[28]	226	;
[2]	227	conduct_u=(conduct+shift(conduct,-1,0,0))/2.
	228	conduct_v=(conduct+shift(conduct,0,1,0))/2.
[28]	229	;
[2]	230	u_cond_u=total( vuconduct_ue3u3d*umask(),3)
	231	v_cond_v=total( vvconduct_ve3v3d*vmask(),3)
[28]	232	;
[2]	233	Bu_star= BR_u*u_cond_u/SIGMA_u
	234	Bv_star= BR_v*v_cond_v/SIGMA_v
[28]	235	;
[2]	236	; Divergence du champ
[30]	237	Diver=divfred(Bu_star,Bv_star)
[2]	238
	239	Diver=Diver*1e-6
	240	;stop
	241	lecontinent=where(Diver GT 1.E08)
	242	Diver(lecontinent)=0.
	243
	244	;stop
	245	;bande de recouvrement::
	246
	247	diver3(1:180,0:147,,jt)=Diver(,*)
	248	diver3(0,,,)=diver3(180,,,)
	249	diver3(181,,,)=diver3(1,,,)
	250
	251	sigma3(1:180,0:147,,jt)=SIGMA(,*)
	252	sigma3(0,,,)=sigma3(180,,,)
	253	sigma3(181,,,)=sigma3(1,,,)
	254
	255	print, jt
	256
[28]	257	ENDFOR
[2]	258	; on ferme le NetCDF
	259	;NCDF_CLOSE,id3
	260	;NCDF_CLOSE,id4
	261
	262	temps=fltarr(73)
	263	temps(0)=0.
	264	for jt=0,71 do begin &$
	265	temps(jt+1)=temps(jt) +5*86400. &$
	266	endfor
	267	print,temps
	268
	269	vargrid = 'T'
[28]	270	iodir = geomag_od_env
[2]	271	; Nom
[29]	272	idout = NCDF_CREATE(iodir+'DivBustar_5d_'+iyear+'_grid_T_'+orcares+'.nc',/clobber)
[2]	273	print, 'Creation du fichier Netcdf'
	274	NCDF_CONTROL, idout, /nofill
	275	; Dimension
	276	xidout = NCDF_DIMDEF(idout, 'x',jpiglo)
	277	yidout = NCDF_DIMDEF(idout, 'y',jpjglo)
	278	didout = NCDF_DIMDEF(idout, 'deptht',1)
	279	tidout = NCDF_DIMDEF(idout, 'time_counter', /unlimited)
	280
	281	didout1 = NCDF_DIMDEF(idout, 'deptht1',jpk)
	282
	283
	284	; Attributs globaux
	285	id0 = NCDF_VARDEF(idout, 'nav_lon' , [xidout, yidout ], /FLOAT)
	286	id1 = NCDF_VARDEF(idout, 'nav_lat' , [xidout, yidout ], /FLOAT)
	287	id2 = NCDF_VARDEF(idout, 'deptht' , [ didout1 ], /FLOAT)
	288	id3 = NCDF_VARDEF(idout, 'time_counter', [ tidout], /FLOAT)
	289	id4 = NCDF_VARDEF(idout, 'Diver' , [xidout, yidout, didout, tidout], /DOUBLE)
	290
	291	; Variable 0
	292	NCDF_ATTPUT, idout, id0, 'units', 'degrees_east'
	293	NCDF_ATTPUT, idout, id0, 'long_name', 'Longitude'
	294	NCDF_ATTPUT, idout, id0, 'nav_model', 'Default grid'
	295	; Variable 1
	296	NCDF_ATTPUT, idout, id1, 'units', 'degrees_north'
	297	NCDF_ATTPUT, idout, id1, 'long_name', 'Latitude'
	298	NCDF_ATTPUT, idout, id1, 'nav_model', 'Default grid'
	299	; Variable 2
	300	NCDF_ATTPUT, idout, id2, 'units','meters'
	301	NCDF_ATTPUT, idout, id2, 'long_name','Depth'
	302	NCDF_ATTPUT, idout, id2, 'nav_model','Default grid'
	303	; Variable3
	304	NCDF_ATTPUT, idout, id3, 'units', 'seconds since 0001-01-01 00:00:00 '
	305	NCDF_ATTPUT, idout, id3, 'calendar','noleap'
	306	NCDF_ATTPUT, idout, id3, 'title', 'Time'
	307	NCDF_ATTPUT, idout, id3, 'long_name', 'Time axis'
	308	NCDF_ATTPUT, idout, id3, 'time_origin','0001-JAN-01 00:00:00'
	309	; Variables
	310	NCDF_ATTPUT, idout, id4, 'long_name', 'Divergence'
	311
	312
	313	NCDF_CONTROL, idout, /ENDEF
	314
	315	; Creation de la longitude
	316	NCDF_VARPUT, idout, id0, glamt
	317	; Creation de la latitude
	318	NCDF_VARPUT, idout, id1, gphit
[28]	319	; Creation de la profondeur
	320	NCDF_VARPUT, idout, id2, gdept
[2]	321	; Creation du calendrier
	322
	323	NCDF_VARPUT, idout, id3, temps
	324
[28]	325
[2]	326	; Ecriture des donnees
	327
[28]	328	; ecriture des glam_8
[2]	329	NCDF_VARPUT, idout, id4 , diver3
	330
	331	NCDF_CLOSE, idout
	332	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	333	;;;;meme topo pour SIGMA
	334	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	335
[29]	336	idout = NCDF_CREATE(iodir+'Sigma_5d_'+iyear+'_grid_T_'+orcares+'.nc',/clobber)
[2]	337	print, 'Creation du fichier Netcdf'
	338	NCDF_CONTROL, idout, /nofill
	339	; Dimension
	340	xidout = NCDF_DIMDEF(idout, 'x',jpiglo)
	341	yidout = NCDF_DIMDEF(idout, 'y',jpjglo)
	342	didout = NCDF_DIMDEF(idout, 'deptht',1)
	343	tidout = NCDF_DIMDEF(idout, 'time_counter', /unlimited)
	344
	345	didout1 = NCDF_DIMDEF(idout, 'deptht1',jpk)
	346
	347
	348	; Attributs globaux
	349	id0 = NCDF_VARDEF(idout, 'nav_lon' , [xidout, yidout ], /FLOAT)
	350	id1 = NCDF_VARDEF(idout, 'nav_lat' , [xidout, yidout ], /FLOAT)
	351	id2 = NCDF_VARDEF(idout, 'deptht' , [ didout1 ], /FLOAT)
	352	id3 = NCDF_VARDEF(idout, 'time_counter', [ tidout], /FLOAT)
	353	id4 = NCDF_VARDEF(idout, 'Sigma' , [xidout, yidout, didout, tidout], /DOUBLE)
	354
	355	; Variable 0
	356	NCDF_ATTPUT, idout, id0, 'units', 'degrees_east'
	357	NCDF_ATTPUT, idout, id0, 'long_name', 'Longitude'
	358	NCDF_ATTPUT, idout, id0, 'nav_model', 'Default grid'
	359	; Variable 1
	360	NCDF_ATTPUT, idout, id1, 'units', 'degrees_north'
	361	NCDF_ATTPUT, idout, id1, 'long_name', 'Latitude'
	362	NCDF_ATTPUT, idout, id1, 'nav_model', 'Default grid'
	363	; Variable 2
	364	NCDF_ATTPUT, idout, id2, 'units','meters'
	365	NCDF_ATTPUT, idout, id2, 'long_name','Depth'
	366	NCDF_ATTPUT, idout, id2, 'nav_model','Default grid'
	367	; Variable3
	368	NCDF_ATTPUT, idout, id3, 'units', 'seconds since 0001-01-01 00:00:00 '
	369	NCDF_ATTPUT, idout, id3, 'calendar','noleap'
	370	NCDF_ATTPUT, idout, id3, 'title', 'Time'
	371	NCDF_ATTPUT, idout, id3, 'long_name', 'Time axis'
	372	NCDF_ATTPUT, idout, id3, 'time_origin','0001-JAN-01 00:00:00'
	373	; Variables
	374	NCDF_ATTPUT, idout, id4, 'long_name', 'Total Conductance'
	375
	376
	377	NCDF_CONTROL, idout, /ENDEF
	378
	379	; Creation de la longitude
	380	NCDF_VARPUT, idout, id0, glamt
	381	; Creation de la latitude
	382	NCDF_VARPUT, idout, id1, gphit
[28]	383	; Creation de la profondeur
	384	NCDF_VARPUT, idout, id2, gdept
[2]	385	; Creation du calendrier
	386
	387	NCDF_VARPUT, idout, id3, temps
	388
[28]	389
[2]	390	; Ecriture des donnees
	391
[28]	392	; ecriture des glam_8
[2]	393	NCDF_VARPUT, idout, id4 , sigma3
	394
	395	NCDF_CLOSE, idout
	396
	397
	398	END

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

Download in other formats:

Original Format