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

  ViewVC Help
Powered by ViewVC 1.1.22