source: trunk/forcagequimarche.pro @ 29

Last change on this file since 29 was 29, checked in by pinsard, 17 years ago

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