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
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; add divchoice parameter to choose either div or divfred
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
49; fplod 2007-07-30T10:12:38Z aedon.locean-ipsl.upmc.fr (Darwin)
50; restore divfred call instead of div ("official" saxo divergence module)
51; fplod 2007-07-27T10:29:18Z cerbere.locean-ipsl.upmc.fr (Linux)
52; add orcares in the name of output files
53; fplod 2007-07-27T10:29:18Z cerbere.locean-ipsl.upmc.fr (Linux)
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;
66PRO 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
92
93; ++ check iyear
94; ++ check ian
95
96;++@init2
97; init grid
98initocemeshmask, orcares, DRAKKAR_EXP = drakkar_exp
99
100@common
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;
151e_exp='ESS'
152key_portrait = 0
153; stockage des fichiers brut
154  ioDATA=geomag_id_env
155
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;
166; title
167     t_exp= e_exp
168     t_bt  = 'bar_transp'
169ioORLN2 = geomag_id_env
170;
171;facteur d'echelle vertical  for partial steps
172e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
173e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
174SIGMAsed=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)
176BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br,/cont_nofill)
177;
178; vertical integration:
179e3t3d=make_array(jpi,jpj,jpk)
180        for k=0, jpk-1 do begin              &$
181          for j=0,jpj-1 do begin              &$
182            for i=0,jpi-1 do begin             &$
183              e3t3d(i,j,k) = e3t(k)    &$
184            endfor                     &$
185          endfor                      &$
186        endfor
187jpt = 73
188;
189;vud = make_array(jpi,jpj,jpt)
190;vvd = make_array(jpi, jpj, jpt)
191divBustar = make_array(jpi, jpj, jpt)
192diver3=replicate(0.,182,149,1,73)
193sigma3=replicate(0.,182,149,1,73)
194;
195FOR jt = 0, jpt-1 DO BEGIN &$
196;
197; ouverture des fichiers et stockage en memoire partial steps
198  vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U)  &$
199;stop
200  vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V)  &$
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.00273147624197*temp + 0.00035133182334*temp*temp + 0.09139808809909*salin + 0.00241425798890*salin*temp -0.00023998958774*salin*salin
208  mask_t=where(conduct GT 1.e+19)
209  conduct(mask_t)=0.
210; Somme conduct au point T
211;
212; Calcul SIGMA
213;
214  SIGMAoc=total(conduct*e3t3d*tmask,3)
215  SIGMA=SIGMAsed+SIGMAoc
216;
217  SIGMA_u=(SIGMA+shift(SIGMA,-1,0))/2.
218  SIGMA_v=(SIGMA+shift(SIGMA,0,1))/2.
219;
220; Calcul B en points u et v
221;
222  BR_u=(BR+shift(BR,-1,0))/2.
223  BR_v=(BR+shift(BR,0,1))/2.
224;
225; Calcul integrale conduct
226;
227  conduct_u=(conduct+shift(conduct,-1,0,0))/2.
228  conduct_v=(conduct+shift(conduct,0,1,0))/2.
229;
230  u_cond_u=total( vu*conduct_u*e3u3d*umask(),3)
231  v_cond_v=total( vv*conduct_v*e3v3d*vmask(),3)
232;
233  Bu_star= BR_u*u_cond_u/SIGMA_u
234  Bv_star= BR_v*v_cond_v/SIGMA_v
235;
236; Divergence du champ
237  Diver=divfred(Bu_star,Bv_star)
238
239Diver=Diver*1e-6
240;stop
241lecontinent=where(Diver GT 1.E08)
242Diver(lecontinent)=0.
243
244;stop
245;bande de recouvrement::
246
247diver3(1:180,0:147,*,jt)=Diver(*,*)
248diver3(0,*,*,*)=diver3(180,*,*,*)
249diver3(181,*,*,*)=diver3(1,*,*,*)
250
251sigma3(1:180,0:147,*,jt)=SIGMA(*,*)
252sigma3(0,*,*,*)=sigma3(180,*,*,*)
253sigma3(181,*,*,*)=sigma3(1,*,*,*)
254
255print,  jt
256
257ENDFOR
258; on ferme le NetCDF
259;NCDF_CLOSE,id3
260;NCDF_CLOSE,id4
261
262temps=fltarr(73)
263temps(0)=0.
264for jt=0,71 do begin &$
265temps(jt+1)=temps(jt) +5*86400. &$
266endfor
267print,temps
268
269   vargrid = 'T'
270   iodir = geomag_od_env
271; Nom
272   idout = NCDF_CREATE(iodir+'DivBustar_5d_'+iyear+'_grid_T_'+orcares+'.nc',/clobber)
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
319; Creation de la profondeur
320   NCDF_VARPUT, idout, id2, gdept
321; Creation du calendrier
322
323   NCDF_VARPUT, idout, id3, temps
324
325
326; Ecriture des donnees
327
328; ecriture des glam_8
329   NCDF_VARPUT, idout, id4 , diver3
330
331   NCDF_CLOSE, idout
332;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
333;;;;meme topo pour SIGMA
334;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
335
336  idout = NCDF_CREATE(iodir+'Sigma_5d_'+iyear+'_grid_T_'+orcares+'.nc',/clobber)
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
383; Creation de la profondeur
384   NCDF_VARPUT, idout, id2, gdept
385; Creation du calendrier
386
387   NCDF_VARPUT, idout, id3, temps
388
389
390; Ecriture des donnees
391
392; ecriture des glam_8
393   NCDF_VARPUT, idout, id4 , sigma3
394
395   NCDF_CLOSE, idout
396
397
398END
Note: See TracBrowser for help on using the repository browser.