/[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.3 - (hide annotations) (download)
Thu Dec 18 06:11:14 2003 UTC (21 years, 7 months ago) by dimitri
Branch: MAIN
Changes since 1.2: +8 -8 lines
bug fix to pkg/cal/cal_toseconds.F
and length-of-year change to ocean_inversion_project/*

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 dimitri 1.3 PARAMETER(nb_seconds_per_year=31556880)
39 dimitri 1.1 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 dimitri 1.3 do j=1,ny
197 dimitri 1.1 C-- convert from PSU.kg/m^2/s to m/yr
198 dimitri 1.3 H2OFLX(i,j,m)=-tmp(i,j)*nb_seconds_per_year/999.8/SAL(i,j,1,m)
199     enddo
200 dimitri 1.1 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 dimitri 1.3 do j=1,ny
322 dimitri 1.1 C-- convert from PSU.kg/m^2/s to m/yr
323 dimitri 1.3 H2OFLX(i,j,m)=-tmp(i,j)*nb_seconds_per_year/999.8/SAL(i,j,1,m)
324     enddo
325 dimitri 1.1 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 dimitri 1.3 & tmp(i,j)*nb_seconds_per_year/999.8/SAL(i,j,1,m)
482 dimitri 1.1 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

  ViewVC Help
Powered by ViewVC 1.1.22