source: trunk/forcageforMOED.pro

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

fix thanks to coding rules
File size: 6.8 KB
Line 
1PRO forcageforMOED
2    @init2
3    @initorca2_bab
4
5    @common
6
7    ian='01'
8    iyear='1993'
9
10    e_exp='ESS'
11    rep_fred='/usr/work/sur/fvi/OPA/geomag/'
12    key_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
28    t_bt  = 'bar_transp'
29    ioORLN2 = '/usr/work/sur/fvi/OPA/ORCA2'
30    ;facteur d'echelle vertical for partial steps
31    e3v3d=read_ncdf('e3v_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
32    e3u3d=read_ncdf('e3u_ps',0,/timestep,iodir=ioORLN2,/nostruct,/tout,filename='meshmask_bab.nc')
33    SIGMAsed=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,/cont_nofill)
35
36    ; vertical integration:
37    e3t3d=make_array(jpi,jpj,jpk)
38    for k=0, jpk-1 do begin              &$
39        for j=0,jpj-1 do begin              &$
40            for i=0,jpi-1 do begin             &$
41                e3t3d(i,j,k) = e3t(k)    &$
42            endfor                     &$
43        endfor                      &$
44    endfor
45    jpt = 73
46
47    ;vud = make_array(jpi,jpj,jpt)
48    ;vvd = make_array(jpi, jpj, jpt)
49    divBustar = make_array(jpi, jpj, jpt)
50    diver3=replicate(0.,182,149,1,73)
51    sigma3=replicate(0.,182,149,1,73)
52    Jpu=replicate(0.,182,149,1,73)
53    Jpv=replicate(0.,182,149,1,73)
54
55
56    FOR jt = 0, jpt-1 DO BEGIN &$
57
58        ; ouverture des fichiers et stockage en memoire partial steps
59          vu=read_ncdf('vozocrtx',jt,jt, /timestep, iodir=ioDATA,/nostruct,/TOUT,filename=file_U)  &$
60        ;stop
61          vv=read_ncdf('vomecrty',jt,jt, /timestep,iodir=ioDATA,/nostruct,/TOUT,filename=file_V)  &$
62        ;stop
63        ; lecture salinite & temperature
64          temp= read_ncdf('votemper',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T)
65        ;stop
66          salin=read_ncdf('vosaline',jt,jt,/timestep,iodir=ioDATA,/nostruct,/tout,filename=file_T)
67        ;stop
68          conduct=0.02047780622061 + 0.00273147624197*temp + 0.00035133182334*temp*temp + 0.09139808809909*salin + 0.00241425798890*salin*temp -0.00023998958774*salin*salin
69          mask_t=where(conduct GT 1.e+19)
70          conduct(mask_t)=0.
71        ; Somme conduct au point T
72
73
74        ; Calcul SIGMA
75
76          SIGMAoc=total(conduct*e3t3d*tmask,3)
77          SIGMA=SIGMAsed+SIGMAoc
78
79          SIGMA_u=(SIGMA+shift(SIGMA,-1,0))/2.
80
81          SIGMA_v=(SIGMA+shift(SIGMA,0,1))/2.
82
83        ; Calcul B en points u et v
84
85          BR_u=(BR+shift(BR,-1,0))/2.
86
87          BR_v=(BR+shift(BR,0,1))/2.
88
89
90        ; Calcul integrale conduct
91
92          conduct_u=(conduct+shift(conduct,-1,0,0))/2.
93
94          conduct_v=(conduct+shift(conduct,0,1,0))/2.
95
96          u_cond_u=total( vu*conduct_u*e3u3d*umask(),3)
97
98          v_cond_v=total( vv*conduct_v*e3v3d*vmask(),3)
99
100
101        ;calcul de Jp= u_star x Fz k
102        Jpv_star= -1*BR_u*u_cond_u
103
104        Jpu_star= BR_v*v_cond_v
105
106        ; Divergence du champ
107
108        Diver=divfred(Jpu_star,Jpv_star)
109
110        Diver=Diver*1e-6
111        ;stop
112        lecontinent=where(Diver GT 1.E08)
113        Diver(lecontinent)=0.
114
115
116
117        ;stop
118        ;bande de recouvrement:
119
120        diver3(1:180,0:147,*,jt)=Diver(*,*)
121        diver3(0,*,*,*)=diver3(180,*,*,*)
122        diver3(181,*,*,*)=diver3(1,*,*,*)
123
124        sigma3(1:180,0:147,*,jt)=SIGMA(*,*)
125        sigma3(0,*,*,*)=sigma3(180,*,*,*)
126        sigma3(181,*,*,*)=sigma3(1,*,*,*)
127
128        Jpu(1:180,0:147,*,jt)=Jpu_star(*,*)
129        Jpu(0,*,*,*)=Jpu(180,*,*,*)
130        Jpu(181,*,*,*)=Jpu(1,*,*,*)
131
132        Jpv(1:180,0:147,*,jt)=Jpv_star(*,*)
133        Jpv(0,*,*,*)=Jpv(180,*,*,*)
134        Jpv(181,*,*,*)=Jpv(1,*,*,*)
135
136        print,  jt
137
138    ENDFOR
139
140    temps=fltarr(73)
141    temps(0)=0.
142    for jt=0,71 do begin &$
143        temps(jt+1)=temps(jt) +5*86400. &$
144    endfor
145    print,temps
146
147    vargrid = 'T'
148    iodir = '/usr/work/sur/fvi/OPA/ORCA2/'
149    ; Nom
150    idout = NCDF_CREATE(iodir+'ForcageMOED_5d_'+iyear+'_grid_T.nc',/clobber)
151    print, 'Creation du fichier Netcdf'
152    NCDF_CONTROL, idout, /nofill
153    ; Dimension
154    xidout = NCDF_DIMDEF(idout, 'x',jpiglo)
155    yidout = NCDF_DIMDEF(idout, 'y',jpjglo)
156    didout = NCDF_DIMDEF(idout, 'deptht',1)
157    tidout = NCDF_DIMDEF(idout, 'time_counter', /unlimited)
158
159    didout1 = NCDF_DIMDEF(idout, 'deptht1',jpk)
160
161
162; Attributs globaux
163    id0  = NCDF_VARDEF(idout, 'nav_lon'     , [xidout, yidout                ], /FLOAT)
164    id1  = NCDF_VARDEF(idout, 'nav_lat'     , [xidout, yidout                ], /FLOAT)
165    id2  = NCDF_VARDEF(idout, 'deptht'      , [                didout1        ], /FLOAT)
166    id3  = NCDF_VARDEF(idout, 'time_counter', [                        tidout], /FLOAT)
167    id4  = NCDF_VARDEF(idout, 'DivJp'  , [xidout, yidout, didout, tidout], /DOUBLE)
168    id5  = NCDF_VARDEF(idout, 'Sigma'  , [xidout, yidout, didout, tidout], /DOUBLE)
169    id6  = NCDF_VARDEF(idout, 'Jpu'  , [xidout, yidout, didout, tidout], /DOUBLE)
170    id7  = NCDF_VARDEF(idout, 'Jpv'  , [xidout, yidout, didout, tidout], /DOUBLE)
171
172; Variable 0
173    NCDF_ATTPUT, idout, id0, 'units', 'degrees_east'
174    NCDF_ATTPUT, idout, id0, 'long_name', 'Longitude'
175    NCDF_ATTPUT, idout, id0, 'nav_model', 'Default grid'
176; Variable 1
177    NCDF_ATTPUT, idout, id1, 'units', 'degrees_north'
178    NCDF_ATTPUT, idout, id1, 'long_name', 'Latitude'
179    NCDF_ATTPUT, idout, id1, 'nav_model', 'Default grid'
180; Variable 2
181    NCDF_ATTPUT, idout, id2, 'units','meters'
182    NCDF_ATTPUT, idout, id2, 'long_name','Depth'
183    NCDF_ATTPUT, idout, id2, 'nav_model','Default grid'
184; Variable3
185    NCDF_ATTPUT, idout, id3, 'units', 'seconds since 0001-01-01 00:00:00 '
186    NCDF_ATTPUT, idout, id3, 'calendar','noleap'
187    NCDF_ATTPUT, idout, id3, 'title', 'Time'
188    NCDF_ATTPUT, idout, id3, 'long_name', 'Time axis'
189    NCDF_ATTPUT, idout, id3, 'time_origin','0001-JAN-01 00:00:00'
190; Variables
191    NCDF_ATTPUT, idout, id4, 'long_name', 'Divergence Jp'
192; Variables
193    NCDF_ATTPUT, idout, id5, 'long_name', 'Total Conductance'
194    NCDF_ATTPUT, idout, id6, 'long_name', 'Jpu=u_star x B (x comp)'
195    NCDF_ATTPUT, idout, id7, 'long_name', 'Jpv=u_star x B (y comp)'
196
197
198    NCDF_CONTROL, idout, /ENDEF
199
200; Creation de la longitude
201    NCDF_VARPUT, idout, id0, glamt
202; Creation de la latitude
203    NCDF_VARPUT, idout, id1, gphit
204; Creation de la profondeur
205    NCDF_VARPUT, idout, id2, gdept
206; Creation du calendrier
207
208    NCDF_VARPUT, idout, id3, temps
209
210
211    ; Ecriture des donnees
212
213    ; ecriture des glam_8
214    NCDF_VARPUT, idout, id4 , diver3
215    NCDF_VARPUT, idout, id5 , sigma3
216    NCDF_VARPUT, idout, id6 , Jpu
217    NCDF_VARPUT, idout, id7 , Jpv
218
219    NCDF_CLOSE, idout
220
221END
Note: See TracBrowser for help on using the repository browser.