source: trunk/forcage.pro @ 2

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

initial import from /usr/work/fvi/OPA/geomag/

File size: 6.4 KB
Line 
1PRO forcagequimarche,iyear,ian
2@init2
3@initorca2_bab
4
5
6
7
8@common
9
10
11
12;ian='01'
13;iyear='1993'
14
15e_exp='ESS'
16rep_fred='/usr/work/sur/fvi/OPA/geomag/'
17key_portrait = 0
18; stockage des fichiers brut
19  ioDATA='/usr/work/sur/fvi/OPA/ORCA2/'
20  file_U=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_U.nc'
21  print, file_U
22  file_V=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_V.nc'
23  print, file_V
24  file_T=e_exp+'_5d_'+ian+'0101_'+ian+'1231_grid_T.nc'
25  print, file_T
26  file_Sed= rep_fred+'cond_sed_ORCA2.nc'
27  file_Br= rep_fred+'Br_ORCA2.nc'
28
29
30; title
31     t_exp= e_exp
32;     t_exp0= e_exp0
33     t_bt  = 'bar_transp'
34ioORLN2 = '/usr/work/sur/fvi/OPA/ORCA2'
35;facteur d'echelle vertical  for partial steps
36e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
37e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
38SIGMAsed=read_ncdf('cond_sed',0,/timestep,/nostruct,/tout,filename=file_Sed,/cont_nofill)
39;BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br)
40BR=read_ncdf('Br',0,/timestep,/nostruct,/tout,filename=file_Br,/cont_nofill)
41
42
43
44
45; vertical integration:
46e3t3d=make_array(jpi,jpj,jpk)   
47        for k=0, jpk-1 do begin              &$
48          for j=0,jpj-1 do begin              &$     
49            for i=0,jpi-1 do begin             &$
50              e3t3d(i,j,k) = e3t(k)    &$
51            endfor                     &$
52          endfor                      &$
53        endfor
54jpt = 73
55
56;vud = make_array(jpi,jpj,jpt)   
57;vvd = make_array(jpi, jpj, jpt)
58divBustar = make_array(jpi, jpj, jpt)
59diver2=fltarr(jpi,jpj,1,73)
60
61; ouverture des fichiers dans lesquels on va écrire
62;id3=NCDF_OPEN('/usr/work/sur/fvi/OPA/geomag/U_5d_'+iyear+'_grid_T.nc',/write)
63;id4=NCDF_OPEN('/usr/work/sur/fvi/OPA/geomag/V_5d_'+iyear+'_grid_T.nc',/write)
64;id4=NCDF_CREATE('/usr/work/sur/fvi/OPA/ORCA2/DivBustar_5d_'+iyear+'_grid_T.nc',/clobber)
65
66;id4=NCDF_OPEN('/usr/work/sur/fvi/OPA/ORCA2/DivBustar_5d_'+iyear+'_grid_T.nc',/write)
67
68FOR jt = 0, jpt-1 DO BEGIN &$   
69
70; ouverture des fichiers et stockage en memoire partial steps
71  vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U)  &$   
72;stop
73  vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V)  &$   
74;stop
75; lecture salinite & temperature
76  temp= read_ncdf('votemper',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T)
77;stop
78  salin=read_ncdf('vosaline',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T)
79;stop
80  conduct=0.02047780622061 + 0.00273147624197*temp + 0.00035133182334*temp*temp + 0.09139808809909*salin + 0.00241425798890*salin*temp -0.00023998958774*salin*salin
81  mask_t=where(conduct GT 1.e+19)
82  conduct(mask_t)=0.
83; Somme conduct au point T
84
85
86; Calcul SIGMA
87
88  SIGMAoc=total(conduct*e3t3d*tmask,3)
89  SIGMA=SIGMAsed+SIGMAoc
90
91
92
93  SIGMA_u=(SIGMA+shift(SIGMA,-1,0))/2.
94
95  SIGMA_v=(SIGMA+shift(SIGMA,0,1))/2.
96
97; Calcul B en points u et v
98
99  BR_u=(BR+shift(BR,-1,0))/2.
100
101  BR_v=(BR+shift(BR,0,1))/2.
102
103
104; Calcul integrale conduct
105
106  conduct_u=(conduct+shift(conduct,-1,0,0))/2.
107
108  conduct_v=(conduct+shift(conduct,0,1,0))/2.
109 
110  u_cond_u=total( vu*conduct_u*e3u3d*umask(),3)
111
112  v_cond_v=total( vv*conduct_v*e3v3d*vmask(),3)
113
114
115  Bu_star= BR_u*u_cond_u/SIGMA_u
116
117  Bv_star= BR_v*v_cond_v/SIGMA_v
118
119; Divergence du champ
120
121  Diver=div(Bu_star,Bv_star)
122
123Diver=Diver*1e-6
124
125;Diver=Diver*1e11
126;for jj=0,147 do begin &$
127;  for ji=0,179 do begin &$
128;   if( Diver(ji,jj) GT 1e10 ) then begin &$
129;       Diver(ji,jj) = 0. &$
130;   endif &$
131;  endfor &$
132;endfor
133
134;stop
135
136 diver2(*,*,0,jt)=Diver(*,*)
137;stop
138; Somme sur la verticale partial steps
139;                vum=total( vu*e3u3d*umask(),3 )  &$   
140;                vvm=total( vv*e3v3d*vmask(),3 )  &$   
141
142; Shift sur la grille T partial steps
143;             vut= (vum+shift(vum,1,0) )*0.5   &$   
144;             vvt= (vvm+shift(vvm,0,1) )*0.5   &$
145; Bande de recouvrement
146;             vut(0, *) = vut(jpi-2, *)
147;             vvt(*, 0) = 0.
148; stockage dans le fichier de sortie
149;NCDF_VARPUT, id3,'sossheig',vut, offset = [0, 0, jt]
150;NCDF_VARPUT, id4,'sossheig',Diver, offset = [0, 0, jt]
151
152print,  jt
153
154
155ENDFOR 
156; on ferme le NetCDF
157;NCDF_CLOSE,id3
158;NCDF_CLOSE,id4
159
160
161
162
163temps=fltarr(73)
164temps(0)=0.
165for jt=0,71 do begin &$
166temps(jt+1)=temps(jt) +5*86400. &$
167endfor
168print,temps
169
170   vargrid = 'T'
171   iodir = '/usr/work/sur/fvi/OPA/ORCA2/'
172; Nom
173   idout = NCDF_CREATE(iodir+'DivBustar_5d_'+iyear+'_grid_T.nc',/clobber)
174   print, 'Creation du fichier Netcdf'
175   NCDF_CONTROL, idout, /nofill
176; Dimension
177   xidout = NCDF_DIMDEF(idout, 'x',jpi)
178   yidout = NCDF_DIMDEF(idout, 'y',jpj)
179   didout = NCDF_DIMDEF(idout, 'deptht',1)
180   tidout = NCDF_DIMDEF(idout, 'time_counter', /unlimited)
181
182   didout1 = NCDF_DIMDEF(idout, 'deptht1',jpk)
183
184
185; Attributs globaux
186   id0  = NCDF_VARDEF(idout, 'nav_lon'     , [xidout, yidout                ], /FLOAT)
187   id1  = NCDF_VARDEF(idout, 'nav_lat'     , [xidout, yidout                ], /FLOAT)
188   id2  = NCDF_VARDEF(idout, 'deptht'      , [                didout1        ], /FLOAT)
189   id3  = NCDF_VARDEF(idout, 'time_counter', [                        tidout], /FLOAT)
190   id4  = NCDF_VARDEF(idout, 'Diver'  , [xidout, yidout, didout, tidout], /DOUBLE)
191
192; Variable 0
193   NCDF_ATTPUT, idout, id0, 'units', 'degrees_east'
194   NCDF_ATTPUT, idout, id0, 'long_name', 'Longitude'
195   NCDF_ATTPUT, idout, id0, 'nav_model', 'Default grid'
196; Variable 1
197   NCDF_ATTPUT, idout, id1, 'units', 'degrees_north'
198   NCDF_ATTPUT, idout, id1, 'long_name', 'Latitude'
199   NCDF_ATTPUT, idout, id1, 'nav_model', 'Default grid'
200; Variable 2
201   NCDF_ATTPUT, idout, id2, 'units','meters'
202   NCDF_ATTPUT, idout, id2, 'long_name','Depth'
203   NCDF_ATTPUT, idout, id2, 'nav_model','Default grid'
204; Variable3
205   NCDF_ATTPUT, idout, id3, 'units', 'seconds since 0001-01-01 00:00:00 '
206   NCDF_ATTPUT, idout, id3, 'calendar','noleap'
207   NCDF_ATTPUT, idout, id3, 'title', 'Time'
208   NCDF_ATTPUT, idout, id3, 'long_name', 'Time axis'
209   NCDF_ATTPUT, idout, id3, 'time_origin','0001-JAN-01 00:00:00'
210; Variables
211   NCDF_ATTPUT, idout, id4, 'long_name', 'Divergence'
212
213
214   NCDF_CONTROL, idout, /ENDEF
215
216; Creation de la longitude
217   NCDF_VARPUT, idout, id0, glamt
218; Creation de la latitude
219   NCDF_VARPUT, idout, id1, gphit
220; Creation de la profondeur
221   NCDF_VARPUT, idout, id2, gdept
222; Creation du calendrier
223
224   NCDF_VARPUT, idout, id3, temps
225
226   
227; Ecriture des donnees
228
229; ecriture des glam_8
230   NCDF_VARPUT, idout, id4 , diver2
231
232   NCDF_CLOSE, idout
233
234
235
236END
Note: See TracBrowser for help on using the repository browser.