1 |
dimitri |
1.1 |
C program to generate netcdf output files for Gruber's |
2 |
|
|
C ocean inversion project |
3 |
|
|
|
4 |
|
|
C note that ECCO_MaskAreaBathy.nc is generated directly |
5 |
|
|
C from the model |
6 |
|
|
|
7 |
|
|
C to compile on nireas: |
8 |
|
|
C f77 mk_output.F write_nc_phys.F nc_util.F \ |
9 |
|
|
C handle_errors.F write_nc_basisfnctns.F \ |
10 |
|
|
C write_nc_diag_0D.F write_nc_diag_2D.F \ |
11 |
|
|
C -I/home/dimitri/software/netcdf/netcdf-3.5.0/include \ |
12 |
|
|
C -L/home/dimitri/software/netcdf/netcdf-3.5.0/lib -lnetcdf |
13 |
|
|
|
14 |
|
|
C to compile on orion: |
15 |
|
|
C setenv F_UFMTENDIAN big |
16 |
|
|
C efc -W0 -WB mk_output.F write_nc_phys.F nc_util.F \ |
17 |
|
|
C handle_errors.F write_nc_basisfnctns.F \ |
18 |
|
|
C write_nc_diag_0D.F write_nc_diag_2D.F \ |
19 |
|
|
C -I/u2/dmenem/software/netcdf-3.5.0/include \ |
20 |
|
|
C -L/u2/dmenem/software/netcdf-3.5.0/lib -lnetcdf |
21 |
|
|
|
22 |
|
|
C=========================================================== |
23 |
|
|
C Constants that depend on model configuration |
24 |
|
|
C=========================================================== |
25 |
|
|
|
26 |
|
|
C nx, ny, nz :: model domain dimensions |
27 |
|
|
C nb_seconds_per_year :: following your model year [s] |
28 |
|
|
C nb_timesteps_per_year:: following your model timestep |
29 |
|
|
C StationaryYears :: total number of years for |
30 |
|
|
C quasi-stationary integration |
31 |
|
|
C delz :: model thicknesses |
32 |
|
|
C p1 :: path for quasi-stationary integration model output |
33 |
|
|
C p2 :: path for time-dependent integration model output |
34 |
|
|
|
35 |
|
|
INTEGER nx , ny , nz |
36 |
|
|
PARAMETER(nx=90, ny=40, nz=15) |
37 |
|
|
INTEGER nb_seconds_per_year |
38 |
|
|
PARAMETER(nb_seconds_per_year=31557600) |
39 |
|
|
INTEGER nb_timesteps_per_year |
40 |
|
|
PARAMETER(nb_timesteps_per_year=180) |
41 |
|
|
INTEGER StationaryYears |
42 |
|
|
PARAMETER(StationaryYears=3001) |
43 |
|
|
REAL delz(nz) |
44 |
|
|
DATA delz /50, 70, 100, 140, 190, 240, 290, 340, 390, |
45 |
|
|
& 440, 490, 540, 590, 640, 690/ |
46 |
|
|
CHARACTER*(52) p1 |
47 |
|
|
PARAMETER( p1= |
48 |
|
|
& '/tmp1/dmenem/checkpoint51n_branch/exe4x4steadystate/') |
49 |
|
|
CHARACTER*(54) p2 |
50 |
|
|
PARAMETER( p2= |
51 |
|
|
& '/tmp1/dmenem/checkpoint51n_branch/exe4x4timedependent/') |
52 |
|
|
|
53 |
|
|
C=========================================================== |
54 |
|
|
C Other constants and variables |
55 |
|
|
C=========================================================== |
56 |
|
|
|
57 |
|
|
C ndyetrac :: number of dye tracers |
58 |
|
|
C nrec :: number of records fro time-dependent output |
59 |
|
|
|
60 |
|
|
INTEGER ndyetrac |
61 |
|
|
PARAMETER(ndyetrac=30) |
62 |
|
|
INTEGER nrec |
63 |
|
|
PARAMETER(nrec =56) |
64 |
|
|
|
65 |
|
|
C secs_per_month = number of seconds in each model month |
66 |
|
|
C PT = monthly potential temperature [C] |
67 |
|
|
C SAL = monthly salinity [psu] |
68 |
|
|
C TFLX = monthly net surface heat flux [W/m^2] |
69 |
|
|
C H2OFLX = monthly net surface freshwater flux [m/y] |
70 |
|
|
C Eul_U = monthly Eulerian Zonal Velocity [m/s] |
71 |
|
|
C WS_x = monthly Zonal Wind Stress [N/m^2] |
72 |
|
|
C Eul_V = monthly Eulerian Meridional Velocity [m/s] |
73 |
|
|
C WS_y = monthly Meridional Wind Stress [N/m^2] |
74 |
|
|
C has_eddy = logical flag, true if model has eddy induced velocities |
75 |
|
|
C Eddy_U = monthly Eddy Induced Zonal Velocity [m/s] |
76 |
|
|
C Eddy_V = monthly Eddy Induced Meridional Velocity [m/s] |
77 |
|
|
C dye_arr= Concentration of dye tracer [mol/cm^3] |
78 |
|
|
C time=time expressed as decimal years |
79 |
|
|
C dye_flux=dye flux for each tracer (mol/m2/s) |
80 |
|
|
C cum_dye_flux=cumulative dye flux for each tracer (mol/m2) |
81 |
|
|
C global_time=time expressed as decimal years |
82 |
|
|
C global_tot_dye=global total dye flux for this year for each tracer (mol) |
83 |
|
|
C global_cum_dye=global cumulative dye flux for each tracer (mol) |
84 |
|
|
C global_mean_conc= global mean dye concentration (mol/m-3) |
85 |
|
|
C year= year of simulation [years] |
86 |
|
|
C RAC = model area |
87 |
|
|
C hFacC = model mask |
88 |
|
|
|
89 |
|
|
REAL secs_per_month ( 12) |
90 |
|
|
REAL PT (nx, ny, nz, 12) |
91 |
|
|
REAL SAL (nx, ny, nz, 12) |
92 |
|
|
REAL TFLX (nx, ny, 12) |
93 |
|
|
REAL H2OFLX (nx, ny, 12) |
94 |
|
|
REAL Eul_U (nx, ny, nz, 12) |
95 |
|
|
REAL WS_x (nx, ny, 12) |
96 |
|
|
REAL Eul_V (nx, ny, nz, 12) |
97 |
|
|
REAL WS_y (nx, ny, 12) |
98 |
|
|
LOGICAL has_eddy |
99 |
|
|
REAL Eddy_U (nx, ny, nz, 12) |
100 |
|
|
REAL Eddy_V (nx, ny, nz, 12) |
101 |
|
|
REAL dye_arr (nx, ny, nz, ndyetrac) |
102 |
|
|
REAL time(nrec) |
103 |
|
|
REAL dye_flux (nx, ny, ndyetrac, nrec) |
104 |
|
|
REAL cum_dye_flux (nx, ny, ndyetrac, nrec) |
105 |
|
|
REAL global_time ( StationaryYears) |
106 |
|
|
REAL global_tot_dye (ndyetrac,StationaryYears) |
107 |
|
|
REAL global_cum_dye (ndyetrac,StationaryYears) |
108 |
|
|
REAL global_mean_conc(ndyetrac,StationaryYears) |
109 |
|
|
INTEGER year |
110 |
|
|
REAL RAC(nx, ny), hFacC(nx, ny, nz) |
111 |
|
|
|
112 |
|
|
INTEGER i, j, k, m, n, step, year, irec, start_step, end_step |
113 |
|
|
CHARACTER*(80) fn |
114 |
|
|
REAL tmp(nx, ny), tmp3D(nx, ny, nz) |
115 |
|
|
|
116 |
|
|
C=========================================================== |
117 |
|
|
print*,'generate ECCO_*_phys.nc files' |
118 |
|
|
C=========================================================== |
119 |
|
|
|
120 |
|
|
do m=1,12 |
121 |
|
|
secs_per_month(m)=nb_seconds_per_year/12 |
122 |
|
|
enddo |
123 |
|
|
|
124 |
|
|
C== Eddy velocity is not computed |
125 |
|
|
has_eddy = .FALSE. |
126 |
|
|
do m=1,12 |
127 |
|
|
do i=1,nx |
128 |
|
|
do j=1,ny |
129 |
|
|
do k=1,nz |
130 |
|
|
Eddy_U(i,j,k,m)=0 |
131 |
|
|
enddo |
132 |
|
|
enddo |
133 |
|
|
enddo |
134 |
|
|
do i=1,nx |
135 |
|
|
do j=1,ny |
136 |
|
|
do k=1,nz |
137 |
|
|
Eddy_V(i,j,k,m)=0 |
138 |
|
|
enddo |
139 |
|
|
enddo |
140 |
|
|
enddo |
141 |
|
|
enddo |
142 |
|
|
|
143 |
|
|
C=========================================================== |
144 |
|
|
print*,'Quasi-stationary, 1st year' |
145 |
|
|
C=========================================================== |
146 |
|
|
|
147 |
|
|
do m=1,12 |
148 |
|
|
step=m*15 |
149 |
|
|
|
150 |
|
|
C== PT = monthly potential temperature [C] |
151 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'Ttave.', step, '.data' |
152 |
|
|
open(100,file=fn,status='old',access='direct', |
153 |
|
|
& recl=nx*ny*4) |
154 |
|
|
do k=1,nz |
155 |
|
|
read(100,rec=k) tmp |
156 |
|
|
do i=1,nx |
157 |
|
|
do j=1,ny |
158 |
|
|
PT (i,j,k,m)=tmp(i,j) |
159 |
|
|
enddo |
160 |
|
|
enddo |
161 |
|
|
enddo |
162 |
|
|
close(100) |
163 |
|
|
|
164 |
|
|
C== SAL = monthly salinity [psu] |
165 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'Stave.', step, '.data' |
166 |
|
|
open(100,file=fn,status='old',access='direct', |
167 |
|
|
& recl=nx*ny*4) |
168 |
|
|
do k=1,nz |
169 |
|
|
read(100,rec=k) tmp |
170 |
|
|
do i=1,nx |
171 |
|
|
do j=1,ny |
172 |
|
|
SAL(i,j,k,m)=tmp(i,j) |
173 |
|
|
enddo |
174 |
|
|
enddo |
175 |
|
|
enddo |
176 |
|
|
close(100) |
177 |
|
|
|
178 |
|
|
C== TFLX = monthly net surface heat flux [W/m^2] |
179 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'tFluxtave.', step, '.data' |
180 |
|
|
open(100,file=fn,status='old',access='direct', |
181 |
|
|
& recl=nx*ny*4) |
182 |
|
|
read(100,rec=1) tmp |
183 |
|
|
do i=1,nx |
184 |
|
|
do j=1,ny |
185 |
|
|
TFLX(i,j,m)=-tmp(i,j) |
186 |
|
|
enddo |
187 |
|
|
enddo |
188 |
|
|
close(100) |
189 |
|
|
|
190 |
|
|
C== H2OFLX = monthly net surface freshwater flux [m/y] |
191 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'sFluxtave.', step, '.data' |
192 |
|
|
open(100,file=fn,status='old',access='direct', |
193 |
|
|
& recl=nx*ny*4) |
194 |
|
|
read(100,rec=1) tmp |
195 |
|
|
do i=1,nx |
196 |
|
|
do j=1,ny |
197 |
|
|
C-- convert from PSU.kg/m^2/s to m/yr |
198 |
|
|
H2OFLX(i,j,m)=-tmp(i,j)*365.25*24*3600/999.8/SAL(i,j,1,m) |
199 |
|
|
enddo |
200 |
|
|
enddo |
201 |
|
|
close(100) |
202 |
|
|
|
203 |
|
|
C== Eul_U = monthly Eulerian Zonal Velocity [m/s] |
204 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'uVeltave.', step, '.data' |
205 |
|
|
open(100,file=fn,status='old',access='direct', |
206 |
|
|
& recl=nx*ny*4) |
207 |
|
|
do k=1,nz |
208 |
|
|
read(100,rec=k) tmp |
209 |
|
|
do i=1,nx |
210 |
|
|
do j=1,ny |
211 |
|
|
Eul_U(i,j,k,m)=tmp(i,j) |
212 |
|
|
enddo |
213 |
|
|
enddo |
214 |
|
|
enddo |
215 |
|
|
close(100) |
216 |
|
|
|
217 |
|
|
C== WS_x = monthly Zonal Wind Stress [N/m^2] |
218 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'uFluxtave.', step, '.data' |
219 |
|
|
open(100,file=fn,status='old',access='direct', |
220 |
|
|
& recl=nx*ny*4) |
221 |
|
|
read(100,rec=1) tmp |
222 |
|
|
do i=1,nx |
223 |
|
|
do j=1,ny |
224 |
|
|
WS_x(i,j,m)=tmp(i,j) |
225 |
|
|
enddo |
226 |
|
|
enddo |
227 |
|
|
close(100) |
228 |
|
|
|
229 |
|
|
C== Eul_V = monthly Eulerian Meridional Velocity [m/s] |
230 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'vVeltave.', step, '.data' |
231 |
|
|
open(100,file=fn,status='old',access='direct', |
232 |
|
|
& recl=nx*ny*4) |
233 |
|
|
do k=1,nz |
234 |
|
|
read(100,rec=k) tmp |
235 |
|
|
do i=1,nx |
236 |
|
|
do j=1,ny |
237 |
|
|
Eul_V(i,j,k,m)=tmp(i,j) |
238 |
|
|
enddo |
239 |
|
|
enddo |
240 |
|
|
enddo |
241 |
|
|
close(100) |
242 |
|
|
|
243 |
|
|
C== WS_y = monthly Meridional Wind Stress [N/m^2] |
244 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'vFluxtave.', step, '.data' |
245 |
|
|
open(100,file=fn,status='old',access='direct', |
246 |
|
|
& recl=nx*ny*4) |
247 |
|
|
read(100,rec=1) tmp |
248 |
|
|
do i=1,nx |
249 |
|
|
do j=1,ny |
250 |
|
|
WS_y(i,j,m)=tmp(i,j) |
251 |
|
|
enddo |
252 |
|
|
enddo |
253 |
|
|
close(100) |
254 |
|
|
enddo |
255 |
|
|
|
256 |
|
|
call write_nc_phys( |
257 |
|
|
& 'ECCO_Stationary_1','MIT GCM Release 1', |
258 |
|
|
& secs_per_month, |
259 |
|
|
& nx, ny, nz, |
260 |
|
|
& PT, SAL, TFLX, H2OFLX, |
261 |
|
|
& nz, |
262 |
|
|
& nx, ny, Eul_U, WS_x, |
263 |
|
|
& nx, ny, Eul_V, WS_y, |
264 |
|
|
& has_eddy, |
265 |
|
|
& nx, ny, Eddy_U, |
266 |
|
|
& nx, ny, Eddy_V) |
267 |
|
|
|
268 |
|
|
C=========================================================== |
269 |
|
|
print*,'Quasi-stationary, last year' |
270 |
|
|
C=========================================================== |
271 |
|
|
|
272 |
|
|
do m=1,12 |
273 |
|
|
step=540000+m*15 |
274 |
|
|
|
275 |
|
|
C== PT = monthly potential temperature [C] |
276 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'Ttave.', step, '.data' |
277 |
|
|
open(100,file=fn,status='old',access='direct', |
278 |
|
|
& recl=nx*ny*4) |
279 |
|
|
do k=1,nz |
280 |
|
|
read(100,rec=k) tmp |
281 |
|
|
do i=1,nx |
282 |
|
|
do j=1,ny |
283 |
|
|
PT (i,j,k,m)=tmp(i,j) |
284 |
|
|
enddo |
285 |
|
|
enddo |
286 |
|
|
enddo |
287 |
|
|
close(100) |
288 |
|
|
|
289 |
|
|
C== SAL = monthly salinity [psu] |
290 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'Stave.', step, '.data' |
291 |
|
|
open(100,file=fn,status='old',access='direct', |
292 |
|
|
& recl=nx*ny*4) |
293 |
|
|
do k=1,nz |
294 |
|
|
read(100,rec=k) tmp |
295 |
|
|
do i=1,nx |
296 |
|
|
do j=1,ny |
297 |
|
|
SAL(i,j,k,m)=tmp(i,j) |
298 |
|
|
enddo |
299 |
|
|
enddo |
300 |
|
|
enddo |
301 |
|
|
close(100) |
302 |
|
|
|
303 |
|
|
C== TFLX = monthly net surface heat flux [W/m^2] |
304 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'tFluxtave.', step, '.data' |
305 |
|
|
open(100,file=fn,status='old',access='direct', |
306 |
|
|
& recl=nx*ny*4) |
307 |
|
|
read(100,rec=1) tmp |
308 |
|
|
do i=1,nx |
309 |
|
|
do j=1,ny |
310 |
|
|
TFLX(i,j,m)=-tmp(i,j) |
311 |
|
|
enddo |
312 |
|
|
enddo |
313 |
|
|
close(100) |
314 |
|
|
|
315 |
|
|
C== H2OFLX = monthly net surface freshwater flux [m/y] |
316 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'sFluxtave.', step, '.data' |
317 |
|
|
open(100,file=fn,status='old',access='direct', |
318 |
|
|
& recl=nx*ny*4) |
319 |
|
|
read(100,rec=1) tmp |
320 |
|
|
do i=1,nx |
321 |
|
|
do j=1,ny |
322 |
|
|
C-- convert from PSU.kg/m^2/s to m/yr |
323 |
|
|
H2OFLX(i,j,m)=-tmp(i,j)*365.25*24*3600/999.8/SAL(i,j,1,m) |
324 |
|
|
enddo |
325 |
|
|
enddo |
326 |
|
|
close(100) |
327 |
|
|
|
328 |
|
|
C== Eul_U = monthly Eulerian Zonal Velocity [m/s] |
329 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'uVeltave.', step, '.data' |
330 |
|
|
open(100,file=fn,status='old',access='direct', |
331 |
|
|
& recl=nx*ny*4) |
332 |
|
|
do k=1,nz |
333 |
|
|
read(100,rec=k) tmp |
334 |
|
|
do i=1,nx |
335 |
|
|
do j=1,ny |
336 |
|
|
Eul_U(i,j,k,m)=tmp(i,j) |
337 |
|
|
enddo |
338 |
|
|
enddo |
339 |
|
|
enddo |
340 |
|
|
close(100) |
341 |
|
|
|
342 |
|
|
C== WS_x = monthly Zonal Wind Stress [N/m^2] |
343 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'uFluxtave.', step, '.data' |
344 |
|
|
open(100,file=fn,status='old',access='direct', |
345 |
|
|
& recl=nx*ny*4) |
346 |
|
|
read(100,rec=1) tmp |
347 |
|
|
do i=1,nx |
348 |
|
|
do j=1,ny |
349 |
|
|
WS_x(i,j,m)=tmp(i,j) |
350 |
|
|
enddo |
351 |
|
|
enddo |
352 |
|
|
close(100) |
353 |
|
|
|
354 |
|
|
C== Eul_V = monthly Eulerian Meridional Velocity [m/s] |
355 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'vVeltave.', step, '.data' |
356 |
|
|
open(100,file=fn,status='old',access='direct', |
357 |
|
|
& recl=nx*ny*4) |
358 |
|
|
do k=1,nz |
359 |
|
|
read(100,rec=k) tmp |
360 |
|
|
do i=1,nx |
361 |
|
|
do j=1,ny |
362 |
|
|
Eul_V(i,j,k,m)=tmp(i,j) |
363 |
|
|
enddo |
364 |
|
|
enddo |
365 |
|
|
enddo |
366 |
|
|
close(100) |
367 |
|
|
|
368 |
|
|
C== WS_y = monthly Meridional Wind Stress [N/m^2] |
369 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p1, 'vFluxtave.', step, '.data' |
370 |
|
|
open(100,file=fn,status='old',access='direct', |
371 |
|
|
& recl=nx*ny*4) |
372 |
|
|
read(100,rec=1) tmp |
373 |
|
|
do i=1,nx |
374 |
|
|
do j=1,ny |
375 |
|
|
WS_y(i,j,m)=tmp(i,j) |
376 |
|
|
enddo |
377 |
|
|
enddo |
378 |
|
|
close(100) |
379 |
|
|
enddo |
380 |
|
|
|
381 |
|
|
call write_nc_phys( |
382 |
|
|
& 'ECCO_Stationary_3001','MIT GCM Release 1', |
383 |
|
|
& secs_per_month, |
384 |
|
|
& nx, ny, nz, |
385 |
|
|
& PT, SAL, TFLX, H2OFLX, |
386 |
|
|
& nz, |
387 |
|
|
& nx, ny, Eul_U, WS_x, |
388 |
|
|
& nx, ny, Eul_V, WS_y, |
389 |
|
|
& has_eddy, |
390 |
|
|
& nx, ny, Eddy_U, |
391 |
|
|
& nx, ny, Eddy_V) |
392 |
|
|
|
393 |
|
|
C=========================================================== |
394 |
|
|
print*,'Time-dependent average of last 10 years (232-241)' |
395 |
|
|
C=========================================================== |
396 |
|
|
|
397 |
|
|
do m=1,12 |
398 |
|
|
|
399 |
|
|
C== initialize to zero |
400 |
|
|
do i=1,nx |
401 |
|
|
do j=1,ny |
402 |
|
|
do k=1,nz |
403 |
|
|
PT (i,j,k,m)=0 |
404 |
|
|
SAL(i,j,k,m)=0 |
405 |
|
|
enddo |
406 |
|
|
TFLX (i,j,m)=0 |
407 |
|
|
H2OFLX(i,j,m)=0 |
408 |
|
|
enddo |
409 |
|
|
enddo |
410 |
|
|
do i=1,nx |
411 |
|
|
do j=1,ny |
412 |
|
|
do k=1,nz |
413 |
|
|
Eul_U(i,j,k,m)=0 |
414 |
|
|
enddo |
415 |
|
|
WS_x(i,j,m)=0 |
416 |
|
|
enddo |
417 |
|
|
enddo |
418 |
|
|
do i=1,nx |
419 |
|
|
do j=1,ny |
420 |
|
|
do k=1,nz |
421 |
|
|
Eul_V(i,j,k,m)=0 |
422 |
|
|
enddo |
423 |
|
|
WS_y(i,j,m)=0 |
424 |
|
|
enddo |
425 |
|
|
enddo |
426 |
|
|
|
427 |
|
|
n=0 |
428 |
|
|
do year=231,240 |
429 |
|
|
n=n+1 |
430 |
|
|
step=year*nb_timesteps_per_year+m*15 |
431 |
|
|
|
432 |
|
|
C== PT = monthly potential temperature [C] |
433 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p2, 'Ttave.', step, '.data' |
434 |
|
|
open(100,file=fn,status='old',access='direct', |
435 |
|
|
& recl=nx*ny*4) |
436 |
|
|
do k=1,nz |
437 |
|
|
read(100,rec=k) tmp |
438 |
|
|
do i=1,nx |
439 |
|
|
do j=1,ny |
440 |
|
|
PT(i,j,k,m)=PT(i,j,k,m)+tmp(i,j) |
441 |
|
|
enddo |
442 |
|
|
enddo |
443 |
|
|
enddo |
444 |
|
|
close(100) |
445 |
|
|
|
446 |
|
|
C== SAL = monthly salinity [psu] |
447 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p2, 'Stave.', step, '.data' |
448 |
|
|
open(100,file=fn,status='old',access='direct', |
449 |
|
|
& recl=nx*ny*4) |
450 |
|
|
do k=1,nz |
451 |
|
|
read(100,rec=k) tmp |
452 |
|
|
do i=1,nx |
453 |
|
|
do j=1,ny |
454 |
|
|
SAL(i,j,k,m)=SAL(i,j,k,m)+tmp(i,j) |
455 |
|
|
enddo |
456 |
|
|
enddo |
457 |
|
|
enddo |
458 |
|
|
close(100) |
459 |
|
|
|
460 |
|
|
C== TFLX = monthly net surface heat flux [W/m^2] |
461 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p2, 'tFluxtave.', step, '.data' |
462 |
|
|
open(100,file=fn,status='old',access='direct', |
463 |
|
|
& recl=nx*ny*4) |
464 |
|
|
read(100,rec=1) tmp |
465 |
|
|
do i=1,nx |
466 |
|
|
do j=1,ny |
467 |
|
|
TFLX(i,j,m)=TFLX(i,j,m)-tmp(i,j) |
468 |
|
|
enddo |
469 |
|
|
enddo |
470 |
|
|
close(100) |
471 |
|
|
|
472 |
|
|
C== H2OFLX = monthly net surface freshwater flux [m/y] |
473 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p2, 'sFluxtave.', step, '.data' |
474 |
|
|
open(100,file=fn,status='old',access='direct', |
475 |
|
|
& recl=nx*ny*4) |
476 |
|
|
read(100,rec=1) tmp |
477 |
|
|
do i=1,nx |
478 |
|
|
do j=1,ny |
479 |
|
|
C-- convert from PSU.kg/m^2/s to m/yr |
480 |
|
|
H2OFLX(i,j,m)=H2OFLX(i,j,m)- |
481 |
|
|
& tmp(i,j)*365.25*24*3600/999.8/SAL(i,j,1,m) |
482 |
|
|
enddo |
483 |
|
|
enddo |
484 |
|
|
close(100) |
485 |
|
|
|
486 |
|
|
C== Eul_U = monthly Eulerian Zonal Velocity [m/s] |
487 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p2, 'uVeltave.', step, '.data' |
488 |
|
|
open(100,file=fn,status='old',access='direct', |
489 |
|
|
& recl=nx*ny*4) |
490 |
|
|
do k=1,nz |
491 |
|
|
read(100,rec=k) tmp |
492 |
|
|
do i=1,nx |
493 |
|
|
do j=1,ny |
494 |
|
|
Eul_U(i,j,k,m)=Eul_U(i,j,k,m)+tmp(i,j) |
495 |
|
|
enddo |
496 |
|
|
enddo |
497 |
|
|
enddo |
498 |
|
|
close(100) |
499 |
|
|
|
500 |
|
|
C== WS_x = monthly Zonal Wind Stress [N/m^2] |
501 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p2, 'uFluxtave.', step, '.data' |
502 |
|
|
open(100,file=fn,status='old',access='direct', |
503 |
|
|
& recl=nx*ny*4) |
504 |
|
|
read(100,rec=1) tmp |
505 |
|
|
do i=1,nx |
506 |
|
|
do j=1,ny |
507 |
|
|
WS_x(i,j,m)=WS_x(i,j,m)+tmp(i,j) |
508 |
|
|
enddo |
509 |
|
|
enddo |
510 |
|
|
close(100) |
511 |
|
|
|
512 |
|
|
C== Eul_V = monthly Eulerian Meridional Velocity [m/s] |
513 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p2, 'vVeltave.', step, '.data' |
514 |
|
|
open(100,file=fn,status='old',access='direct', |
515 |
|
|
& recl=nx*ny*4) |
516 |
|
|
do k=1,nz |
517 |
|
|
read(100,rec=k) tmp |
518 |
|
|
do i=1,nx |
519 |
|
|
do j=1,ny |
520 |
|
|
Eul_V(i,j,k,m)=Eul_V(i,j,k,m)+tmp(i,j) |
521 |
|
|
enddo |
522 |
|
|
enddo |
523 |
|
|
enddo |
524 |
|
|
close(100) |
525 |
|
|
|
526 |
|
|
C== WS_y = monthly Meridional Wind Stress [N/m^2] |
527 |
|
|
WRITE(fn,'(A,A,I10.10,A)') p2, 'vFluxtave.', step, '.data' |
528 |
|
|
open(100,file=fn,status='old',access='direct', |
529 |
|
|
& recl=nx*ny*4) |
530 |
|
|
read(100,rec=1) tmp |
531 |
|
|
do i=1,nx |
532 |
|
|
do j=1,ny |
533 |
|
|
WS_y(i,j,m)=WS_y(i,j,m)+tmp(i,j) |
534 |
|
|
enddo |
535 |
|
|
enddo |
536 |
|
|
close(100) |
537 |
|
|
enddo |
538 |
|
|
|
539 |
|
|
C== normalize |
540 |
|
|
do i=1,nx |
541 |
|
|
do j=1,ny |
542 |
|
|
do k=1,nz |
543 |
|
|
PT (i,j,k,m)=PT (i,j,k,m)/n |
544 |
|
|
SAL(i,j,k,m)=SAL(i,j,k,m)/n |
545 |
|
|
enddo |
546 |
|
|
TFLX (i,j,m)=TFLX (i,j,m)/n |
547 |
|
|
H2OFLX(i,j,m)=H2OFLX(i,j,m)/n |
548 |
|
|
enddo |
549 |
|
|
enddo |
550 |
|
|
do i=1,nx |
551 |
|
|
do j=1,ny |
552 |
|
|
do k=1,nz |
553 |
|
|
Eul_U(i,j,k,m)=Eul_U(i,j,k,m)/n |
554 |
|
|
enddo |
555 |
|
|
WS_x(i,j,m)=WS_x(i,j,m)/n |
556 |
|
|
enddo |
557 |
|
|
enddo |
558 |
|
|
do i=1,nx |
559 |
|
|
do j=1,ny |
560 |
|
|
do k=1,nz |
561 |
|
|
Eul_V(i,j,k,m)=Eul_V(i,j,k,m)/n |
562 |
|
|
enddo |
563 |
|
|
WS_y(i,j,m)=WS_y(i,j,m)/n |
564 |
|
|
enddo |
565 |
|
|
enddo |
566 |
|
|
enddo |
567 |
|
|
|
568 |
|
|
call write_nc_phys( |
569 |
|
|
& 'ECCO_Timedep','MIT GCM Release 1', |
570 |
|
|
& secs_per_month, |
571 |
|
|
& nx, ny, nz, |
572 |
|
|
& PT, SAL, TFLX, H2OFLX, |
573 |
|
|
& nz, |
574 |
|
|
& nx, ny, Eul_U, WS_x, |
575 |
|
|
& nx, ny, Eul_V, WS_y, |
576 |
|
|
& has_eddy, |
577 |
|
|
& nx, ny, Eddy_U, |
578 |
|
|
& nx, ny, Eddy_V) |
579 |
|
|
|
580 |
|
|
C=========================================================== |
581 |
|
|
print*,'generate ECCO_*_BasisFNCTNS_*.nc and diag_2D files' |
582 |
|
|
C=========================================================== |
583 |
|
|
|
584 |
|
|
C=========================================================== |
585 |
|
|
print*,'Quasi-stationary annual mean basis functions' |
586 |
|
|
C== dye_arr= Concentration of dye tracer [mol/cm^3] |
587 |
|
|
C=========================================================== |
588 |
|
|
|
589 |
|
|
year = StationaryYears |
590 |
|
|
step = year * nb_timesteps_per_year |
591 |
|
|
do n=1,ndyetrac |
592 |
|
|
WRITE(fn,'(A,A,I2.2,A,I10.10,A)') |
593 |
|
|
& p1, 'PTRtave', n, '.', step, '.data' |
594 |
|
|
open(100,file=fn,status='old',access='direct', |
595 |
|
|
& recl=nx*ny*4) |
596 |
|
|
do k=1,nz |
597 |
|
|
read(100,rec=k) tmp |
598 |
|
|
do i=1,nx |
599 |
|
|
do j=1,ny |
600 |
|
|
dye_arr(i,j,k,n)=100*100*100*tmp(i,j) |
601 |
|
|
enddo |
602 |
|
|
enddo |
603 |
|
|
enddo |
604 |
|
|
close(100) |
605 |
|
|
enddo |
606 |
|
|
call write_nc_basisfnctns( |
607 |
|
|
& 'ECCO','MIT GCM Release 1','Stationary', |
608 |
|
|
& nx,ny,nz,ndyetrac, |
609 |
|
|
& year,nb_seconds_per_year,nb_timesteps_per_year, |
610 |
|
|
& dye_arr) |
611 |
|
|
|
612 |
|
|
C== 2-D diagnostics |
613 |
|
|
|
614 |
|
|
C dye_flux=dye flux for each tracer (mol/m2/s) |
615 |
|
|
time(1)=year |
616 |
|
|
do n=1,ndyetrac |
617 |
|
|
WRITE(fn,'(A,A,I2.2,A,I10.10,A)') |
618 |
|
|
& p1, 'PtrFlux', n, '.', step, '.data' |
619 |
|
|
open(100,file=fn,status='old',access='direct', |
620 |
|
|
& recl=nx*ny*4) |
621 |
|
|
read(100,rec=1) tmp |
622 |
|
|
do i=1,nx |
623 |
|
|
do j=1,ny |
624 |
|
|
dye_flux(i,j,n,1)=tmp(i,j) |
625 |
|
|
enddo |
626 |
|
|
enddo |
627 |
|
|
close(100) |
628 |
|
|
enddo |
629 |
|
|
|
630 |
|
|
C cum_dye_flux=cumulative dye flux for each tracer (mol/m2) |
631 |
|
|
do n=1,ndyetrac |
632 |
|
|
do i=1,nx |
633 |
|
|
do j=1,ny |
634 |
|
|
cum_dye_flux(i,j,n,1)=0. |
635 |
|
|
enddo |
636 |
|
|
enddo |
637 |
|
|
enddo |
638 |
|
|
do step=nb_timesteps_per_year,540180,nb_timesteps_per_year |
639 |
|
|
do n=1,ndyetrac |
640 |
|
|
WRITE(fn,'(A,A,I2.2,A,I10.10,A)') |
641 |
|
|
& p1, 'PtrFlux', n, '.', step, '.data' |
642 |
|
|
open(100,file=fn,status='old',access='direct', |
643 |
|
|
& recl=nx*ny*4) |
644 |
|
|
read(100,rec=1) tmp |
645 |
|
|
do i=1,nx |
646 |
|
|
do j=1,ny |
647 |
|
|
cum_dye_flux(i,j,n,1)=cum_dye_flux(i,j,n,1)+ |
648 |
|
|
& tmp(i,j)*nb_seconds_per_year |
649 |
|
|
enddo |
650 |
|
|
enddo |
651 |
|
|
close(100) |
652 |
|
|
enddo |
653 |
|
|
enddo |
654 |
|
|
|
655 |
|
|
call write_nc_diag_2D( |
656 |
|
|
& 'ECCO','MIT GCM Release 1','Stationary', |
657 |
|
|
& nx,ny,ndyetrac, |
658 |
|
|
& 1, time, dye_flux,cum_dye_flux) |
659 |
|
|
|
660 |
|
|
C=========================================================== |
661 |
|
|
print*,'Time-dependent annual mean basis functions' |
662 |
|
|
C== 1/10-years for 1775-1965, 1/year for 1970-2005 |
663 |
|
|
C== dye_arr= Concentration of dye tracer [mol/cm^3] |
664 |
|
|
C=========================================================== |
665 |
|
|
|
666 |
|
|
n=0 |
667 |
|
|
do year=1775,1965,10 |
668 |
|
|
n=n+1 |
669 |
|
|
time(n)=year |
670 |
|
|
enddo |
671 |
|
|
do year=1970,2005 |
672 |
|
|
n=n+1 |
673 |
|
|
time(n)=year |
674 |
|
|
enddo |
675 |
|
|
|
676 |
|
|
do irec=1,nrec |
677 |
|
|
year=time(irec) |
678 |
|
|
step = (year-1764) * nb_timesteps_per_year |
679 |
|
|
do n=1,ndyetrac |
680 |
|
|
WRITE(fn,'(A,A,I2.2,A,I10.10,A)') |
681 |
|
|
& p2, 'PTRtave', n, '.', step, '.data' |
682 |
|
|
open(100,file=fn,status='old',access='direct', |
683 |
|
|
& recl=nx*ny*4) |
684 |
|
|
do k=1,nz |
685 |
|
|
read(100,rec=k) tmp |
686 |
|
|
do i=1,nx |
687 |
|
|
do j=1,ny |
688 |
|
|
dye_arr(i,j,k,n)=100*100*100*tmp(i,j) |
689 |
|
|
enddo |
690 |
|
|
enddo |
691 |
|
|
enddo |
692 |
|
|
close(100) |
693 |
|
|
enddo |
694 |
|
|
call write_nc_basisfnctns( |
695 |
|
|
& 'ECCO','MIT GCM Release 1','Timedep', |
696 |
|
|
& nx,ny,nz,ndyetrac, |
697 |
|
|
& year,nb_seconds_per_year,nb_timesteps_per_year, |
698 |
|
|
& dye_arr) |
699 |
|
|
|
700 |
|
|
C== 2-D diagnostics |
701 |
|
|
|
702 |
|
|
C dye_flux=dye flux for each tracer (mol/m2/s) |
703 |
|
|
do n=1,ndyetrac |
704 |
|
|
WRITE(fn,'(A,A,I2.2,A,I10.10,A)') |
705 |
|
|
& p2, 'PtrFlux', n, '.', step, '.data' |
706 |
|
|
open(100,file=fn,status='old',access='direct', |
707 |
|
|
& recl=nx*ny*4) |
708 |
|
|
read(100,rec=1) tmp |
709 |
|
|
do i=1,nx |
710 |
|
|
do j=1,ny |
711 |
dimitri |
1.2 |
dye_flux(i,j,n,irec)=tmp(i,j) |
712 |
dimitri |
1.1 |
enddo |
713 |
|
|
enddo |
714 |
|
|
close(100) |
715 |
|
|
enddo |
716 |
|
|
|
717 |
|
|
C cum_dye_flux=cumulative dye flux for each tracer (mol/m2) |
718 |
|
|
do n=1,ndyetrac |
719 |
|
|
do i=1,nx |
720 |
|
|
do j=1,ny |
721 |
|
|
if (irec.eq.1) then |
722 |
|
|
cum_dye_flux(i,j,n,irec)=0. |
723 |
|
|
else |
724 |
|
|
cum_dye_flux(i,j,n,irec)=cum_dye_flux(i,j,n,irec-1) |
725 |
|
|
endif |
726 |
|
|
enddo |
727 |
|
|
enddo |
728 |
|
|
enddo |
729 |
|
|
start_step=nb_timesteps_per_year |
730 |
|
|
if (irec.gt.1) |
731 |
|
|
& start_step=nb_timesteps_per_year+ |
732 |
|
|
& (time(irec-1)-1764)*nb_timesteps_per_year |
733 |
dimitri |
1.2 |
end_step= (year-1764) * nb_timesteps_per_year |
734 |
dimitri |
1.1 |
do step=start_step,end_step,nb_timesteps_per_year |
735 |
|
|
do n=1,ndyetrac |
736 |
|
|
WRITE(fn,'(A,A,I2.2,A,I10.10,A)') |
737 |
|
|
& p2, 'PtrFlux', n, '.', step, '.data' |
738 |
|
|
open(100,file=fn,status='old',access='direct', |
739 |
|
|
& recl=nx*ny*4) |
740 |
|
|
read(100,rec=1) tmp |
741 |
|
|
do i=1,nx |
742 |
|
|
do j=1,ny |
743 |
dimitri |
1.2 |
cum_dye_flux(i,j,n,irec)=cum_dye_flux(i,j,n,irec)+ |
744 |
dimitri |
1.1 |
& tmp(i,j)*nb_seconds_per_year |
745 |
|
|
enddo |
746 |
|
|
enddo |
747 |
|
|
close(100) |
748 |
|
|
enddo |
749 |
|
|
enddo |
750 |
|
|
enddo |
751 |
|
|
|
752 |
|
|
call write_nc_diag_2D( |
753 |
|
|
& 'ECCO','MIT GCM Release 1','Timedep', |
754 |
|
|
& nx,ny,ndyetrac, |
755 |
|
|
& nrec, time, dye_flux,cum_dye_flux) |
756 |
|
|
|
757 |
|
|
C=========================================================== |
758 |
|
|
print*,'write_nc_diag_0D quasi-stationary diagnostics' |
759 |
|
|
C=========================================================== |
760 |
|
|
|
761 |
|
|
WRITE(fn,'(A,A)') p1, 'RAC.data' |
762 |
|
|
open(100,file=fn,status='old',access='direct', |
763 |
|
|
& recl=nx*ny*4) |
764 |
|
|
read(100,rec=1) RAC |
765 |
|
|
close(100) |
766 |
|
|
WRITE(fn,'(A,A)') p1, 'hFacC.data' |
767 |
|
|
open(100,file=fn,status='old',access='direct', |
768 |
|
|
& recl=nx*ny*nz*4) |
769 |
|
|
read(100,rec=1) hFacC |
770 |
|
|
close(100) |
771 |
|
|
|
772 |
|
|
irec=0 |
773 |
dimitri |
1.2 |
do year=1,StationaryYears |
774 |
dimitri |
1.1 |
irec=irec+1 |
775 |
|
|
global_time(irec)=year |
776 |
|
|
step=year*nb_timesteps_per_year |
777 |
|
|
do n=1,ndyetrac |
778 |
|
|
|
779 |
|
|
C global_tot_dye=global total dye flux for this year for each tracer (mol) |
780 |
|
|
WRITE(fn,'(A,A,I2.2,A,I10.10,A)') |
781 |
|
|
& p1, 'PtrFlux', n, '.', step, '.data' |
782 |
|
|
open(100,file=fn,status='old',access='direct', |
783 |
|
|
& recl=nx*ny*4) |
784 |
|
|
read(100,rec=1) tmp |
785 |
|
|
close(100) |
786 |
|
|
global_tot_dye(n,irec)=0. |
787 |
|
|
do i=1,nx |
788 |
|
|
do j=1,ny |
789 |
|
|
global_tot_dye(n,irec)=global_tot_dye(n,irec)+ |
790 |
|
|
& RAC(i,j)*hFacC(i,j,1)*tmp(i,j)*nb_seconds_per_year |
791 |
|
|
enddo |
792 |
|
|
enddo |
793 |
|
|
|
794 |
|
|
C global_cum_dye=global cumulative dye flux for each tracer (mol) |
795 |
|
|
global_cum_dye(n,irec)=global_tot_dye(n,irec) |
796 |
|
|
if (irec.gt.1) global_cum_dye(n,irec) = |
797 |
|
|
& global_cum_dye(n,irec) + global_cum_dye(n,irec-1) |
798 |
|
|
|
799 |
|
|
C global_mean_conc= global mean dye concentration (mol/m-3) |
800 |
|
|
WRITE(fn,'(A,A,I2.2,A,I10.10,A)') |
801 |
|
|
& p1, 'PTRtave', n, '.', step, '.data' |
802 |
|
|
open(100,file=fn,status='old',access='direct', |
803 |
|
|
& recl=nx*ny*nz*4) |
804 |
|
|
read(100,rec=1) tmp3D |
805 |
|
|
close(100) |
806 |
|
|
global_mean_conc(n,irec)=0. |
807 |
|
|
do i=1,nx |
808 |
|
|
do j=1,ny |
809 |
|
|
do k=1,nz |
810 |
|
|
global_mean_conc(n,irec)=global_mean_conc(n,irec)+ |
811 |
|
|
& RAC(i,j)*hFacC(i,j,k)*tmp3D(i,j,k) |
812 |
|
|
enddo |
813 |
|
|
enddo |
814 |
|
|
enddo |
815 |
|
|
enddo |
816 |
|
|
enddo |
817 |
|
|
|
818 |
|
|
call write_nc_diag_0D( |
819 |
|
|
& 'ECCO','MIT GCM Release 1','Stationary', |
820 |
|
|
& StationaryYears, global_time, ndyetrac, |
821 |
|
|
& global_tot_dye, global_cum_dye, global_mean_conc) |
822 |
|
|
|
823 |
|
|
C=========================================================== |
824 |
|
|
print*,'write_nc_diag_0D time-dependent diagnostics' |
825 |
|
|
C=========================================================== |
826 |
|
|
|
827 |
|
|
irec=0 |
828 |
|
|
do year=1765,2005 |
829 |
|
|
irec=irec+1 |
830 |
|
|
global_time(irec)=year |
831 |
|
|
step=(year-1764)*nb_timesteps_per_year |
832 |
|
|
do n=1,ndyetrac |
833 |
|
|
|
834 |
|
|
C global_tot_dye=global total dye flux for this year for each tracer (mol) |
835 |
|
|
WRITE(fn,'(A,A,I2.2,A,I10.10,A)') |
836 |
|
|
& p2, 'PtrFlux', n, '.', step, '.data' |
837 |
|
|
open(100,file=fn,status='old',access='direct', |
838 |
|
|
& recl=nx*ny*4) |
839 |
|
|
read(100,rec=1) tmp |
840 |
|
|
close(100) |
841 |
|
|
global_tot_dye(n,irec)=0. |
842 |
|
|
do i=1,nx |
843 |
|
|
do j=1,ny |
844 |
|
|
global_tot_dye(n,irec)=global_tot_dye(n,irec)+ |
845 |
|
|
& RAC(i,j)*hFacC(i,j,1)*tmp(i,j)*nb_seconds_per_year |
846 |
|
|
enddo |
847 |
|
|
enddo |
848 |
|
|
|
849 |
|
|
C global_cum_dye=global cumulative dye flux for each tracer (mol) |
850 |
|
|
global_cum_dye(n,irec)=global_tot_dye(n,irec) |
851 |
|
|
if (irec.gt.1) global_cum_dye(n,irec) = |
852 |
|
|
& global_cum_dye(n,irec) + global_cum_dye(n,irec-1) |
853 |
|
|
|
854 |
|
|
C global_mean_conc= global mean dye concentration (mol/m-3) |
855 |
|
|
WRITE(fn,'(A,A,I2.2,A,I10.10,A)') |
856 |
|
|
& p2, 'PTRtave', n, '.', step, '.data' |
857 |
|
|
open(100,file=fn,status='old',access='direct', |
858 |
|
|
& recl=nx*ny*nz*4) |
859 |
|
|
read(100,rec=1) tmp3D |
860 |
|
|
close(100) |
861 |
|
|
global_mean_conc(n,irec)=0. |
862 |
|
|
do i=1,nx |
863 |
|
|
do j=1,ny |
864 |
|
|
do k=1,nz |
865 |
|
|
global_mean_conc(n,irec)=global_mean_conc(n,irec)+ |
866 |
|
|
& RAC(i,j)*hFacC(i,j,k)*tmp3D(i,j,k) |
867 |
|
|
enddo |
868 |
|
|
enddo |
869 |
|
|
enddo |
870 |
|
|
enddo |
871 |
|
|
enddo |
872 |
|
|
|
873 |
|
|
call write_nc_diag_0D( |
874 |
|
|
& 'ECCO','MIT GCM Release 1','Timedep', |
875 |
|
|
& irec, global_time, ndyetrac, |
876 |
|
|
& global_tot_dye, global_cum_dye, global_mean_conc) |
877 |
|
|
|
878 |
|
|
stop |
879 |
|
|
end |