source: trunk/forcagequimarche.pro @ 36

Last change on this file since 36 was 36, checked in by pinsard, 16 years ago

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