source: trunk/forcagequimarche.pro @ 48

Last change on this file since 48 was 48, checked in by pinsard, 10 years ago

fix thanks to coding rules
File size: 12.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; 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
100    initocemeshmask, 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  ;
153  e_exp='ESS'
154  key_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'
163  jpt = 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'
172  ioORLN2 = geomag_id_env
173  ;
174  ;facteur d'echelle vertical for partial steps
175  e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
176  e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
177  SIGMAsed=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)
179  BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br,/cont_nofill)
180  ;
181  ; vertical integration:
182  e3t3d=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)
193  divBustar = make_array(jpi, jpj, jpt)
194  diver3=replicate(0.,182,149,1,jpt)
195  sigma3=replicate(0.,182,149,1,jpt)
196  ;
197  FOR 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 
246  Diver=Diver*1e-6
247  ;stop
248  lecontinent=where(Diver GT 1.E08)
249  Diver(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.