/[MITgcm]/MITgcm_contrib/ocean_inversion_project/write_netCDF/mk_output.F
ViewVC logotype

Annotation of /MITgcm_contrib/ocean_inversion_project/write_netCDF/mk_output.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.4 - (hide annotations) (download)
Mon Nov 1 20:03:35 2004 UTC (20 years, 8 months ago) by dimitri
Branch: MAIN
Changes since 1.3: +92 -59 lines 
Tfreezing=-3 in model/src/freeze.F model/src/ini_theta.F
Modified ocean_inversion_project/README
         ocean_inversion_project/input_ecco1x1/data.1765-2005.[1-2]
         ocean_inversion_project/write_netCDF/mk_output.F
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 dimitri 1.4 PARAMETER(nx=360, ny=160, nz=23)
37 dimitri 1.1 INTEGER nb_seconds_per_year
38 dimitri 1.4 PARAMETER(nb_seconds_per_year=31536000)
39 dimitri 1.1 INTEGER nb_timesteps_per_year
40 dimitri 1.4 PARAMETER(nb_timesteps_per_year=8760)
41 dimitri 1.1 INTEGER StationaryYears
42     PARAMETER(StationaryYears=3001)
43     REAL delz(nz)
44 dimitri 1.4 DATA delz /10.,10.,15.,20.,20.,25.,35.,50.,75.,100.,150.,200.,
45     & 275.,350.,415.,450.,500.,500.,500.,500.,500.,500.,500./
46 dimitri 1.1 CHARACTER*(52) p1
47     PARAMETER( p1=
48 dimitri 1.4 & ' ')
49 dimitri 1.1 CHARACTER*(54) p2
50     PARAMETER( p2=
51 dimitri 1.4 & '/nobackup19/menemenl/ocmip/MITgcm/exe_freeze/')
52 dimitri 1.1
53     C===========================================================
54     C Other constants and variables
55     C===========================================================
56    
57     C ndyetrac :: number of dye tracers
58 dimitri 1.4 C nrec :: number of records for time-dependent output
59 dimitri 1.1
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 dimitri 1.4 IF ( p1 .NE. ' ' ) THEN
148    
149     do m=1,12
150     step=m*(nb_timesteps_per_year/12)
151 dimitri 1.1
152     C== PT = monthly potential temperature [C]
153     WRITE(fn,'(A,A,I10.10,A)') p1, 'Ttave.', step, '.data'
154     open(100,file=fn,status='old',access='direct',
155     & recl=nx*ny*4)
156     do k=1,nz
157     read(100,rec=k) tmp
158     do i=1,nx
159     do j=1,ny
160     PT (i,j,k,m)=tmp(i,j)
161     enddo
162     enddo
163     enddo
164     close(100)
165    
166     C== SAL = monthly salinity [psu]
167     WRITE(fn,'(A,A,I10.10,A)') p1, 'Stave.', step, '.data'
168     open(100,file=fn,status='old',access='direct',
169     & recl=nx*ny*4)
170     do k=1,nz
171     read(100,rec=k) tmp
172     do i=1,nx
173     do j=1,ny
174     SAL(i,j,k,m)=tmp(i,j)
175     enddo
176     enddo
177     enddo
178     close(100)
179    
180     C== TFLX = monthly net surface heat flux [W/m^2]
181     WRITE(fn,'(A,A,I10.10,A)') p1, 'tFluxtave.', step, '.data'
182     open(100,file=fn,status='old',access='direct',
183     & recl=nx*ny*4)
184     read(100,rec=1) tmp
185     do i=1,nx
186     do j=1,ny
187     TFLX(i,j,m)=-tmp(i,j)
188     enddo
189     enddo
190     close(100)
191    
192     C== H2OFLX = monthly net surface freshwater flux [m/y]
193     WRITE(fn,'(A,A,I10.10,A)') p1, 'sFluxtave.', step, '.data'
194     open(100,file=fn,status='old',access='direct',
195     & recl=nx*ny*4)
196     read(100,rec=1) tmp
197     do i=1,nx
198 dimitri 1.3 do j=1,ny
199 dimitri 1.1 C-- convert from PSU.kg/m^2/s to m/yr
200 dimitri 1.3 H2OFLX(i,j,m)=-tmp(i,j)*nb_seconds_per_year/999.8/SAL(i,j,1,m)
201     enddo
202 dimitri 1.1 enddo
203     close(100)
204    
205     C== Eul_U = monthly Eulerian Zonal Velocity [m/s]
206     WRITE(fn,'(A,A,I10.10,A)') p1, 'uVeltave.', step, '.data'
207     open(100,file=fn,status='old',access='direct',
208     & recl=nx*ny*4)
209     do k=1,nz
210     read(100,rec=k) tmp
211     do i=1,nx
212     do j=1,ny
213     Eul_U(i,j,k,m)=tmp(i,j)
214     enddo
215     enddo
216     enddo
217     close(100)
218    
219     C== WS_x = monthly Zonal Wind Stress [N/m^2]
220     WRITE(fn,'(A,A,I10.10,A)') p1, 'uFluxtave.', step, '.data'
221     open(100,file=fn,status='old',access='direct',
222     & recl=nx*ny*4)
223     read(100,rec=1) tmp
224     do i=1,nx
225     do j=1,ny
226     WS_x(i,j,m)=tmp(i,j)
227     enddo
228     enddo
229     close(100)
230    
231     C== Eul_V = monthly Eulerian Meridional Velocity [m/s]
232     WRITE(fn,'(A,A,I10.10,A)') p1, 'vVeltave.', step, '.data'
233     open(100,file=fn,status='old',access='direct',
234     & recl=nx*ny*4)
235     do k=1,nz
236     read(100,rec=k) tmp
237     do i=1,nx
238     do j=1,ny
239     Eul_V(i,j,k,m)=tmp(i,j)
240     enddo
241     enddo
242     enddo
243     close(100)
244    
245     C== WS_y = monthly Meridional Wind Stress [N/m^2]
246     WRITE(fn,'(A,A,I10.10,A)') p1, 'vFluxtave.', step, '.data'
247     open(100,file=fn,status='old',access='direct',
248     & recl=nx*ny*4)
249     read(100,rec=1) tmp
250     do i=1,nx
251     do j=1,ny
252     WS_y(i,j,m)=tmp(i,j)
253     enddo
254     enddo
255     close(100)
256 dimitri 1.4 enddo
257 dimitri 1.1
258 dimitri 1.4 call write_nc_phys(
259     & 'ECCO1_Stationary_1','MITgcm_checkpoint51n_post',
260 dimitri 1.1 & secs_per_month,
261     & nx, ny, nz,
262     & PT, SAL, TFLX, H2OFLX,
263     & nz,
264     & nx, ny, Eul_U, WS_x,
265     & nx, ny, Eul_V, WS_y,
266     & has_eddy,
267     & nx, ny, Eddy_U,
268     & nx, ny, Eddy_V)
269    
270     C===========================================================
271     print*,'Quasi-stationary, last year'
272     C===========================================================
273    
274 dimitri 1.4 do m=1,12
275     step=(StationaryYears-1)*nb_timesteps_per_year+
276     & m*(nb_timesteps_per_year/12)
277 dimitri 1.1
278     C== PT = monthly potential temperature [C]
279     WRITE(fn,'(A,A,I10.10,A)') p1, 'Ttave.', step, '.data'
280     open(100,file=fn,status='old',access='direct',
281     & recl=nx*ny*4)
282     do k=1,nz
283     read(100,rec=k) tmp
284     do i=1,nx
285     do j=1,ny
286     PT (i,j,k,m)=tmp(i,j)
287     enddo
288     enddo
289     enddo
290     close(100)
291    
292     C== SAL = monthly salinity [psu]
293     WRITE(fn,'(A,A,I10.10,A)') p1, 'Stave.', step, '.data'
294     open(100,file=fn,status='old',access='direct',
295     & recl=nx*ny*4)
296     do k=1,nz
297     read(100,rec=k) tmp
298     do i=1,nx
299     do j=1,ny
300     SAL(i,j,k,m)=tmp(i,j)
301     enddo
302     enddo
303     enddo
304     close(100)
305    
306     C== TFLX = monthly net surface heat flux [W/m^2]
307     WRITE(fn,'(A,A,I10.10,A)') p1, 'tFluxtave.', step, '.data'
308     open(100,file=fn,status='old',access='direct',
309     & recl=nx*ny*4)
310     read(100,rec=1) tmp
311     do i=1,nx
312     do j=1,ny
313     TFLX(i,j,m)=-tmp(i,j)
314     enddo
315     enddo
316     close(100)
317    
318     C== H2OFLX = monthly net surface freshwater flux [m/y]
319     WRITE(fn,'(A,A,I10.10,A)') p1, 'sFluxtave.', step, '.data'
320     open(100,file=fn,status='old',access='direct',
321     & recl=nx*ny*4)
322     read(100,rec=1) tmp
323     do i=1,nx
324 dimitri 1.3 do j=1,ny
325 dimitri 1.1 C-- convert from PSU.kg/m^2/s to m/yr
326 dimitri 1.3 H2OFLX(i,j,m)=-tmp(i,j)*nb_seconds_per_year/999.8/SAL(i,j,1,m)
327     enddo
328 dimitri 1.1 enddo
329     close(100)
330    
331     C== Eul_U = monthly Eulerian Zonal Velocity [m/s]
332     WRITE(fn,'(A,A,I10.10,A)') p1, 'uVeltave.', step, '.data'
333     open(100,file=fn,status='old',access='direct',
334     & recl=nx*ny*4)
335     do k=1,nz
336     read(100,rec=k) tmp
337     do i=1,nx
338     do j=1,ny
339     Eul_U(i,j,k,m)=tmp(i,j)
340     enddo
341     enddo
342     enddo
343     close(100)
344    
345     C== WS_x = monthly Zonal Wind Stress [N/m^2]
346     WRITE(fn,'(A,A,I10.10,A)') p1, 'uFluxtave.', step, '.data'
347     open(100,file=fn,status='old',access='direct',
348     & recl=nx*ny*4)
349     read(100,rec=1) tmp
350     do i=1,nx
351     do j=1,ny
352     WS_x(i,j,m)=tmp(i,j)
353     enddo
354     enddo
355     close(100)
356    
357     C== Eul_V = monthly Eulerian Meridional Velocity [m/s]
358     WRITE(fn,'(A,A,I10.10,A)') p1, 'vVeltave.', step, '.data'
359     open(100,file=fn,status='old',access='direct',
360     & recl=nx*ny*4)
361     do k=1,nz
362     read(100,rec=k) tmp
363     do i=1,nx
364     do j=1,ny
365     Eul_V(i,j,k,m)=tmp(i,j)
366     enddo
367     enddo
368     enddo
369     close(100)
370    
371     C== WS_y = monthly Meridional Wind Stress [N/m^2]
372     WRITE(fn,'(A,A,I10.10,A)') p1, 'vFluxtave.', step, '.data'
373     open(100,file=fn,status='old',access='direct',
374     & recl=nx*ny*4)
375     read(100,rec=1) tmp
376     do i=1,nx
377     do j=1,ny
378     WS_y(i,j,m)=tmp(i,j)
379     enddo
380     enddo
381     close(100)
382 dimitri 1.4 enddo
383 dimitri 1.1
384 dimitri 1.4 call write_nc_phys(
385     & 'ECCO1_Stationary_3001','MITgcm_checkpoint51n_post',
386 dimitri 1.1 & secs_per_month,
387     & nx, ny, nz,
388     & PT, SAL, TFLX, H2OFLX,
389     & nz,
390     & nx, ny, Eul_U, WS_x,
391     & nx, ny, Eul_V, WS_y,
392     & has_eddy,
393     & nx, ny, Eddy_U,
394     & nx, ny, Eddy_V)
395    
396 dimitri 1.4 C IF ( p1 .NE. ' ' ) THEN
397     ENDIF
398    
399 dimitri 1.1 C===========================================================
400     print*,'Time-dependent average of last 10 years (232-241)'
401     C===========================================================
402    
403 dimitri 1.4 IF ( p2 .NE. ' ' ) THEN
404    
405     do m=1,12
406 dimitri 1.1
407     C== initialize to zero
408     do i=1,nx
409     do j=1,ny
410     do k=1,nz
411     PT (i,j,k,m)=0
412     SAL(i,j,k,m)=0
413     enddo
414     TFLX (i,j,m)=0
415     H2OFLX(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_U(i,j,k,m)=0
422     enddo
423     WS_x(i,j,m)=0
424     enddo
425     enddo
426     do i=1,nx
427     do j=1,ny
428     do k=1,nz
429     Eul_V(i,j,k,m)=0
430     enddo
431     WS_y(i,j,m)=0
432     enddo
433     enddo
434    
435     n=0
436     do year=231,240
437     n=n+1
438 dimitri 1.4 step=year*nb_timesteps_per_year+m*(nb_timesteps_per_year/12)
439 dimitri 1.1
440     C== PT = monthly potential temperature [C]
441     WRITE(fn,'(A,A,I10.10,A)') p2, 'Ttave.', step, '.data'
442     open(100,file=fn,status='old',access='direct',
443     & recl=nx*ny*4)
444     do k=1,nz
445     read(100,rec=k) tmp
446     do i=1,nx
447     do j=1,ny
448     PT(i,j,k,m)=PT(i,j,k,m)+tmp(i,j)
449     enddo
450     enddo
451     enddo
452     close(100)
453    
454     C== SAL = monthly salinity [psu]
455     WRITE(fn,'(A,A,I10.10,A)') p2, 'Stave.', step, '.data'
456     open(100,file=fn,status='old',access='direct',
457     & recl=nx*ny*4)
458     do k=1,nz
459     read(100,rec=k) tmp
460     do i=1,nx
461     do j=1,ny
462     SAL(i,j,k,m)=SAL(i,j,k,m)+tmp(i,j)
463     enddo
464     enddo
465     enddo
466     close(100)
467    
468     C== TFLX = monthly net surface heat flux [W/m^2]
469     WRITE(fn,'(A,A,I10.10,A)') p2, 'tFluxtave.', step, '.data'
470     open(100,file=fn,status='old',access='direct',
471     & recl=nx*ny*4)
472     read(100,rec=1) tmp
473     do i=1,nx
474     do j=1,ny
475     TFLX(i,j,m)=TFLX(i,j,m)-tmp(i,j)
476     enddo
477     enddo
478     close(100)
479    
480     C== H2OFLX = monthly net surface freshwater flux [m/y]
481     WRITE(fn,'(A,A,I10.10,A)') p2, 'sFluxtave.', step, '.data'
482     open(100,file=fn,status='old',access='direct',
483     & recl=nx*ny*4)
484     read(100,rec=1) tmp
485     do i=1,nx
486     do j=1,ny
487     C-- convert from PSU.kg/m^2/s to m/yr
488     H2OFLX(i,j,m)=H2OFLX(i,j,m)-
489 dimitri 1.3 & tmp(i,j)*nb_seconds_per_year/999.8/SAL(i,j,1,m)
490 dimitri 1.1 enddo
491     enddo
492     close(100)
493    
494     C== Eul_U = monthly Eulerian Zonal Velocity [m/s]
495     WRITE(fn,'(A,A,I10.10,A)') p2, 'uVeltave.', step, '.data'
496     open(100,file=fn,status='old',access='direct',
497     & recl=nx*ny*4)
498     do k=1,nz
499     read(100,rec=k) tmp
500     do i=1,nx
501     do j=1,ny
502     Eul_U(i,j,k,m)=Eul_U(i,j,k,m)+tmp(i,j)
503     enddo
504     enddo
505     enddo
506     close(100)
507    
508     C== WS_x = monthly Zonal Wind Stress [N/m^2]
509     WRITE(fn,'(A,A,I10.10,A)') p2, 'uFluxtave.', step, '.data'
510     open(100,file=fn,status='old',access='direct',
511     & recl=nx*ny*4)
512     read(100,rec=1) tmp
513     do i=1,nx
514     do j=1,ny
515     WS_x(i,j,m)=WS_x(i,j,m)+tmp(i,j)
516     enddo
517     enddo
518     close(100)
519    
520     C== Eul_V = monthly Eulerian Meridional Velocity [m/s]
521     WRITE(fn,'(A,A,I10.10,A)') p2, 'vVeltave.', step, '.data'
522     open(100,file=fn,status='old',access='direct',
523     & recl=nx*ny*4)
524     do k=1,nz
525     read(100,rec=k) tmp
526     do i=1,nx
527     do j=1,ny
528     Eul_V(i,j,k,m)=Eul_V(i,j,k,m)+tmp(i,j)
529     enddo
530     enddo
531     enddo
532     close(100)
533    
534     C== WS_y = monthly Meridional Wind Stress [N/m^2]
535     WRITE(fn,'(A,A,I10.10,A)') p2, 'vFluxtave.', step, '.data'
536     open(100,file=fn,status='old',access='direct',
537     & recl=nx*ny*4)
538     read(100,rec=1) tmp
539     do i=1,nx
540     do j=1,ny
541     WS_y(i,j,m)=WS_y(i,j,m)+tmp(i,j)
542     enddo
543     enddo
544     close(100)
545     enddo
546    
547     C== normalize
548     do i=1,nx
549     do j=1,ny
550     do k=1,nz
551     PT (i,j,k,m)=PT (i,j,k,m)/n
552     SAL(i,j,k,m)=SAL(i,j,k,m)/n
553     enddo
554     TFLX (i,j,m)=TFLX (i,j,m)/n
555     H2OFLX(i,j,m)=H2OFLX(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_U(i,j,k,m)=Eul_U(i,j,k,m)/n
562     enddo
563     WS_x(i,j,m)=WS_x(i,j,m)/n
564     enddo
565     enddo
566     do i=1,nx
567     do j=1,ny
568     do k=1,nz
569     Eul_V(i,j,k,m)=Eul_V(i,j,k,m)/n
570     enddo
571     WS_y(i,j,m)=WS_y(i,j,m)/n
572     enddo
573     enddo
574 dimitri 1.4 enddo
575 dimitri 1.1
576 dimitri 1.4 call write_nc_phys(
577 dimitri 1.1 & 'ECCO_Timedep','MIT GCM Release 1',
578     & secs_per_month,
579     & nx, ny, nz,
580     & PT, SAL, TFLX, H2OFLX,
581     & nz,
582     & nx, ny, Eul_U, WS_x,
583     & nx, ny, Eul_V, WS_y,
584     & has_eddy,
585     & nx, ny, Eddy_U,
586     & nx, ny, Eddy_V)
587    
588 dimitri 1.4 C IF ( p2 .NE. ' ' ) THEN
589     ENDIF
590    
591 dimitri 1.1 C===========================================================
592     print*,'generate ECCO_*_BasisFNCTNS_*.nc and diag_2D files'
593     C===========================================================
594    
595     C===========================================================
596     print*,'Quasi-stationary annual mean basis functions'
597     C== dye_arr= Concentration of dye tracer [mol/cm^3]
598     C===========================================================
599    
600 dimitri 1.4 IF ( p1 .NE. ' ' ) THEN
601    
602     year = StationaryYears
603     step = year * nb_timesteps_per_year
604     do n=1,ndyetrac
605 dimitri 1.1 WRITE(fn,'(A,A,I2.2,A,I10.10,A)')
606     & p1, 'PTRtave', n, '.', step, '.data'
607     open(100,file=fn,status='old',access='direct',
608     & recl=nx*ny*4)
609     do k=1,nz
610     read(100,rec=k) tmp
611     do i=1,nx
612     do j=1,ny
613     dye_arr(i,j,k,n)=100*100*100*tmp(i,j)
614     enddo
615     enddo
616     enddo
617     close(100)
618 dimitri 1.4 enddo
619     call write_nc_basisfnctns(
620 dimitri 1.1 & 'ECCO','MIT GCM Release 1','Stationary',
621     & nx,ny,nz,ndyetrac,
622     & year,nb_seconds_per_year,nb_timesteps_per_year,
623     & dye_arr)
624    
625     C== 2-D diagnostics
626    
627     C dye_flux=dye flux for each tracer (mol/m2/s)
628 dimitri 1.4 time(1)=year
629     do n=1,ndyetrac
630 dimitri 1.1 WRITE(fn,'(A,A,I2.2,A,I10.10,A)')
631     & p1, 'PtrFlux', n, '.', step, '.data'
632     open(100,file=fn,status='old',access='direct',
633     & recl=nx*ny*4)
634     read(100,rec=1) tmp
635     do i=1,nx
636     do j=1,ny
637     dye_flux(i,j,n,1)=tmp(i,j)
638     enddo
639     enddo
640     close(100)
641 dimitri 1.4 enddo
642 dimitri 1.1
643     C cum_dye_flux=cumulative dye flux for each tracer (mol/m2)
644 dimitri 1.4 do n=1,ndyetrac
645 dimitri 1.1 do i=1,nx
646     do j=1,ny
647     cum_dye_flux(i,j,n,1)=0.
648     enddo
649     enddo
650 dimitri 1.4 enddo
651     do step=nb_timesteps_per_year,
652     & (StationaryYears*nb_timesteps_per_year),
653     & nb_timesteps_per_year
654 dimitri 1.1 do n=1,ndyetrac
655     WRITE(fn,'(A,A,I2.2,A,I10.10,A)')
656     & p1, 'PtrFlux', n, '.', step, '.data'
657     open(100,file=fn,status='old',access='direct',
658     & recl=nx*ny*4)
659     read(100,rec=1) tmp
660     do i=1,nx
661     do j=1,ny
662     cum_dye_flux(i,j,n,1)=cum_dye_flux(i,j,n,1)+
663     & tmp(i,j)*nb_seconds_per_year
664     enddo
665     enddo
666     close(100)
667     enddo
668 dimitri 1.4 enddo
669 dimitri 1.1
670 dimitri 1.4 call write_nc_diag_2D(
671 dimitri 1.1 & 'ECCO','MIT GCM Release 1','Stationary',
672     & nx,ny,ndyetrac,
673     & 1, time, dye_flux,cum_dye_flux)
674    
675 dimitri 1.4 C IF ( p1 .NE. ' ' ) THEN
676     ENDIF
677    
678 dimitri 1.1 C===========================================================
679     print*,'Time-dependent annual mean basis functions'
680     C== 1/10-years for 1775-1965, 1/year for 1970-2005
681     C== dye_arr= Concentration of dye tracer [mol/cm^3]
682     C===========================================================
683    
684 dimitri 1.4 IF ( p2 .NE. ' ' ) THEN
685    
686     n=0
687     do year=1775,1965,10
688 dimitri 1.1 n=n+1
689     time(n)=year
690 dimitri 1.4 enddo
691     do year=1970,2005
692 dimitri 1.1 n=n+1
693     time(n)=year
694 dimitri 1.4 enddo
695 dimitri 1.1
696 dimitri 1.4 do irec=1,nrec
697 dimitri 1.1 year=time(irec)
698     step = (year-1764) * nb_timesteps_per_year
699     do n=1,ndyetrac
700     WRITE(fn,'(A,A,I2.2,A,I10.10,A)')
701     & p2, 'PTRtave', n, '.', step, '.data'
702     open(100,file=fn,status='old',access='direct',
703     & recl=nx*ny*4)
704     do k=1,nz
705     read(100,rec=k) tmp
706     do i=1,nx
707     do j=1,ny
708     dye_arr(i,j,k,n)=100*100*100*tmp(i,j)
709     enddo
710     enddo
711     enddo
712     close(100)
713     enddo
714     call write_nc_basisfnctns(
715     & 'ECCO','MIT GCM Release 1','Timedep',
716     & nx,ny,nz,ndyetrac,
717     & year,nb_seconds_per_year,nb_timesteps_per_year,
718     & dye_arr)
719    
720     C== 2-D diagnostics
721    
722     C dye_flux=dye flux for each tracer (mol/m2/s)
723     do n=1,ndyetrac
724     WRITE(fn,'(A,A,I2.2,A,I10.10,A)')
725     & p2, 'PtrFlux', n, '.', step, '.data'
726     open(100,file=fn,status='old',access='direct',
727     & recl=nx*ny*4)
728     read(100,rec=1) tmp
729     do i=1,nx
730     do j=1,ny
731 dimitri 1.2 dye_flux(i,j,n,irec)=tmp(i,j)
732 dimitri 1.1 enddo
733     enddo
734     close(100)
735     enddo
736    
737     C cum_dye_flux=cumulative dye flux for each tracer (mol/m2)
738     do n=1,ndyetrac
739     do i=1,nx
740     do j=1,ny
741     if (irec.eq.1) then
742     cum_dye_flux(i,j,n,irec)=0.
743     else
744     cum_dye_flux(i,j,n,irec)=cum_dye_flux(i,j,n,irec-1)
745     endif
746     enddo
747     enddo
748     enddo
749     start_step=nb_timesteps_per_year
750     if (irec.gt.1)
751     & start_step=nb_timesteps_per_year+
752     & (time(irec-1)-1764)*nb_timesteps_per_year
753 dimitri 1.2 end_step= (year-1764) * nb_timesteps_per_year
754 dimitri 1.1 do step=start_step,end_step,nb_timesteps_per_year
755     do n=1,ndyetrac
756     WRITE(fn,'(A,A,I2.2,A,I10.10,A)')
757     & p2, 'PtrFlux', n, '.', step, '.data'
758     open(100,file=fn,status='old',access='direct',
759     & recl=nx*ny*4)
760     read(100,rec=1) tmp
761     do i=1,nx
762     do j=1,ny
763 dimitri 1.2 cum_dye_flux(i,j,n,irec)=cum_dye_flux(i,j,n,irec)+
764 dimitri 1.1 & tmp(i,j)*nb_seconds_per_year
765     enddo
766     enddo
767     close(100)
768     enddo
769     enddo
770 dimitri 1.4 enddo
771 dimitri 1.1
772 dimitri 1.4 call write_nc_diag_2D(
773 dimitri 1.1 & 'ECCO','MIT GCM Release 1','Timedep',
774     & nx,ny,ndyetrac,
775     & nrec, time, dye_flux,cum_dye_flux)
776    
777 dimitri 1.4 C IF ( p2 .NE. ' ' ) THEN
778     ENDIF
779    
780 dimitri 1.1 C===========================================================
781     print*,'write_nc_diag_0D quasi-stationary diagnostics'
782     C===========================================================
783    
784 dimitri 1.4 IF ( p1 .NE. ' ' ) THEN
785    
786     WRITE(fn,'(A,A)') p1, 'RAC.data'
787     open(100,file=fn,status='old',access='direct',
788 dimitri 1.1 & recl=nx*ny*4)
789 dimitri 1.4 read(100,rec=1) RAC
790     close(100)
791     WRITE(fn,'(A,A)') p1, 'hFacC.data'
792     open(100,file=fn,status='old',access='direct',
793 dimitri 1.1 & recl=nx*ny*nz*4)
794 dimitri 1.4 read(100,rec=1) hFacC
795     close(100)
796 dimitri 1.1
797 dimitri 1.4 irec=0
798     do year=1,StationaryYears
799 dimitri 1.1 irec=irec+1
800     global_time(irec)=year
801     step=year*nb_timesteps_per_year
802     do n=1,ndyetrac
803    
804     C global_tot_dye=global total dye flux for this year for each tracer (mol)
805     WRITE(fn,'(A,A,I2.2,A,I10.10,A)')
806     & p1, 'PtrFlux', n, '.', step, '.data'
807     open(100,file=fn,status='old',access='direct',
808     & recl=nx*ny*4)
809     read(100,rec=1) tmp
810     close(100)
811     global_tot_dye(n,irec)=0.
812     do i=1,nx
813     do j=1,ny
814     global_tot_dye(n,irec)=global_tot_dye(n,irec)+
815     & RAC(i,j)*hFacC(i,j,1)*tmp(i,j)*nb_seconds_per_year
816     enddo
817     enddo
818    
819     C global_cum_dye=global cumulative dye flux for each tracer (mol)
820     global_cum_dye(n,irec)=global_tot_dye(n,irec)
821     if (irec.gt.1) global_cum_dye(n,irec) =
822     & global_cum_dye(n,irec) + global_cum_dye(n,irec-1)
823    
824     C global_mean_conc= global mean dye concentration (mol/m-3)
825     WRITE(fn,'(A,A,I2.2,A,I10.10,A)')
826     & p1, 'PTRtave', n, '.', step, '.data'
827     open(100,file=fn,status='old',access='direct',
828     & recl=nx*ny*nz*4)
829     read(100,rec=1) tmp3D
830     close(100)
831     global_mean_conc(n,irec)=0.
832     do i=1,nx
833     do j=1,ny
834     do k=1,nz
835     global_mean_conc(n,irec)=global_mean_conc(n,irec)+
836     & RAC(i,j)*hFacC(i,j,k)*tmp3D(i,j,k)
837     enddo
838     enddo
839     enddo
840     enddo
841 dimitri 1.4 enddo
842 dimitri 1.1
843 dimitri 1.4 call write_nc_diag_0D(
844 dimitri 1.1 & 'ECCO','MIT GCM Release 1','Stationary',
845     & StationaryYears, global_time, ndyetrac,
846     & global_tot_dye, global_cum_dye, global_mean_conc)
847    
848 dimitri 1.4 C IF ( p1 .NE. ' ' ) THEN
849     ENDIF
850    
851 dimitri 1.1 C===========================================================
852     print*,'write_nc_diag_0D time-dependent diagnostics'
853     C===========================================================
854    
855 dimitri 1.4 IF ( p2 .NE. ' ' ) THEN
856    
857     irec=0
858     do year=1765,2005
859 dimitri 1.1 irec=irec+1
860     global_time(irec)=year
861     step=(year-1764)*nb_timesteps_per_year
862     do n=1,ndyetrac
863    
864     C global_tot_dye=global total dye flux for this year for each tracer (mol)
865     WRITE(fn,'(A,A,I2.2,A,I10.10,A)')
866     & p2, 'PtrFlux', n, '.', step, '.data'
867     open(100,file=fn,status='old',access='direct',
868     & recl=nx*ny*4)
869     read(100,rec=1) tmp
870     close(100)
871     global_tot_dye(n,irec)=0.
872     do i=1,nx
873     do j=1,ny
874     global_tot_dye(n,irec)=global_tot_dye(n,irec)+
875     & RAC(i,j)*hFacC(i,j,1)*tmp(i,j)*nb_seconds_per_year
876     enddo
877     enddo
878    
879     C global_cum_dye=global cumulative dye flux for each tracer (mol)
880     global_cum_dye(n,irec)=global_tot_dye(n,irec)
881     if (irec.gt.1) global_cum_dye(n,irec) =
882     & global_cum_dye(n,irec) + global_cum_dye(n,irec-1)
883    
884     C global_mean_conc= global mean dye concentration (mol/m-3)
885     WRITE(fn,'(A,A,I2.2,A,I10.10,A)')
886     & p2, 'PTRtave', n, '.', step, '.data'
887     open(100,file=fn,status='old',access='direct',
888     & recl=nx*ny*nz*4)
889     read(100,rec=1) tmp3D
890     close(100)
891     global_mean_conc(n,irec)=0.
892     do i=1,nx
893     do j=1,ny
894     do k=1,nz
895     global_mean_conc(n,irec)=global_mean_conc(n,irec)+
896     & RAC(i,j)*hFacC(i,j,k)*tmp3D(i,j,k)
897     enddo
898     enddo
899     enddo
900     enddo
901 dimitri 1.4 enddo
902 dimitri 1.1
903 dimitri 1.4 call write_nc_diag_0D(
904 dimitri 1.1 & 'ECCO','MIT GCM Release 1','Timedep',
905     & irec, global_time, ndyetrac,
906     & global_tot_dye, global_cum_dye, global_mean_conc)
907    
908 dimitri 1.4 C IF ( p2 .NE. ' ' ) THEN
909     ENDIF
910    
911 dimitri 1.1 stop
912     end
  ViewVC Help
Powered by ViewVC 1.1.22