/[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.8 - (hide annotations) (download)
Wed May 4 23:31:39 2005 UTC (20 years, 2 months ago) by dimitri
Branch: MAIN
CVS Tags: HEAD
Changes since 1.7: +5 -5 lines
modified mk_output.F for steady state solution
also removed a spurious integer definition

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

  ViewVC Help
Powered by ViewVC 1.1.22