source: trunk/forcage.pro

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

fix thanks to coding rules

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