C program to generate netcdf output files for Gruber's C ocean inversion project C note that ECCO_MaskAreaBathy.nc is generated directly C from the model C to compile on nireas: C f77 mk_output.F write_nc_phys.F nc_util.F \ C handle_errors.F write_nc_basisfnctns.F \ C write_nc_diag_0D.F write_nc_diag_2D.F \ C -I/home/dimitri/software/netcdf/netcdf-3.5.0/include \ C -L/home/dimitri/software/netcdf/netcdf-3.5.0/lib -lnetcdf C to compile on orion: C setenv F_UFMTENDIAN big C efc -W0 -WB mk_output.F write_nc_phys.F nc_util.F \ C handle_errors.F write_nc_basisfnctns.F \ C write_nc_diag_0D.F write_nc_diag_2D.F \ C -I/u2/dmenem/software/netcdf-3.5.0/include \ C -L/u2/dmenem/software/netcdf-3.5.0/lib -lnetcdf C=========================================================== C Constants that depend on model configuration C=========================================================== C nx, ny, nz :: model domain dimensions C nb_seconds_per_year :: following your model year [s] C nb_timesteps_per_year:: following your model timestep C StationaryYears :: total number of years for C quasi-stationary integration C delz :: model thicknesses C p1 :: path for quasi-stationary integration model output C p2 :: path for time-dependent integration model output INTEGER nx , ny , nz PARAMETER(nx=360, ny=160, nz=23) INTEGER nb_seconds_per_year PARAMETER(nb_seconds_per_year=31536000) INTEGER nb_timesteps_per_year PARAMETER(nb_timesteps_per_year=8760) INTEGER StationaryYears PARAMETER(StationaryYears=3001) REAL delz(nz) DATA delz /10.,10.,15.,20.,20.,25.,35.,50.,75.,100.,150.,200., & 275.,350.,415.,450.,500.,500.,500.,500.,500.,500.,500./ CHARACTER*(38) p1 PARAMETER( p1= & '/nobackup16/menemenl/ocmip/MITgcm/exe/') C 12345678911234567892123456789312345678941234567895123456789 CHARACTER*(1) p2 PARAMETER( p2= & ' ') C & '/nobackup2/menemenl/ocmip/MITgcm/exe/') C 12345678911234567892123456789312345678941234567895123456789 C=========================================================== C Other constants and variables C=========================================================== C ndyetrac :: number of dye tracers C nrec :: number of records for time-dependent output INTEGER ndyetrac PARAMETER(ndyetrac=30) INTEGER nrec PARAMETER(nrec =56) C secs_per_month = number of seconds in each model month C PT = monthly potential temperature [C] C SAL = monthly salinity [psu] C TFLX = monthly net surface heat flux [W/m^2] C H2OFLX = monthly net surface freshwater flux [m/y] C Eul_U = monthly Eulerian Zonal Velocity [m/s] C WS_x = monthly Zonal Wind Stress [N/m^2] C Eul_V = monthly Eulerian Meridional Velocity [m/s] C WS_y = monthly Meridional Wind Stress [N/m^2] C has_eddy = logical flag, true if model has eddy induced velocities C Eddy_U = monthly Eddy Induced Zonal Velocity [m/s] C Eddy_V = monthly Eddy Induced Meridional Velocity [m/s] C dye_arr= Concentration of dye tracer [mol/cm^3] C time=time expressed as decimal years C dye_flux=dye flux for each tracer (mol/m2/s) C cum_dye_flux=cumulative dye flux for each tracer (mol/m2) C global_time=time expressed as decimal years C global_tot_dye=global total dye flux for this year for each tracer (mol) C global_cum_dye=global cumulative dye flux for each tracer (mol) C global_mean_conc= global mean dye concentration (mol/m-3) C year= year of simulation [years] C RAC = model area C hFacC = model mask REAL secs_per_month ( 12) REAL PT (nx, ny, nz, 12) REAL SAL (nx, ny, nz, 12) REAL TFLX (nx, ny, 12) REAL H2OFLX (nx, ny, 12) REAL Eul_U (nx, ny, nz, 12) REAL WS_x (nx, ny, 12) REAL Eul_V (nx, ny, nz, 12) REAL WS_y (nx, ny, 12) LOGICAL has_eddy REAL Eddy_U (nx, ny, nz, 12) REAL Eddy_V (nx, ny, nz, 12) REAL dye_arr (nx, ny, nz, ndyetrac) REAL time(nrec) REAL dye_flux (nx, ny, ndyetrac, nrec) REAL cum_dye_flux (nx, ny, ndyetrac, nrec) REAL global_time ( StationaryYears) REAL global_tot_dye (ndyetrac,StationaryYears) REAL global_cum_dye (ndyetrac,StationaryYears) REAL global_mean_conc(ndyetrac,StationaryYears) REAL RAC(nx, ny), hFacC(nx, ny, nz) INTEGER i, j, k, m, n, step, year, irec, start_step, end_step CHARACTER*(80) fn REAL tmp(nx, ny), tmp3D(nx, ny, nz) C=========================================================== print*,'generate ECCO_*_phys.nc files' C=========================================================== do m=1,12 secs_per_month(m)=nb_seconds_per_year/12 enddo C== Eddy velocity is not computed has_eddy = .FALSE. do m=1,12 do i=1,nx do j=1,ny do k=1,nz Eddy_U(i,j,k,m)=0 enddo enddo enddo do i=1,nx do j=1,ny do k=1,nz Eddy_V(i,j,k,m)=0 enddo enddo enddo enddo C=========================================================== print*,'Quasi-stationary, 1st year' C=========================================================== IF ( p1 .NE. ' ' ) THEN do m=1,12 step=m*(nb_timesteps_per_year/12) C== PT = monthly potential temperature [C] WRITE(fn,'(A,A,I10.10,A)') p1, 'Ttave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny PT (i,j,k,m)=tmp(i,j) enddo enddo enddo close(100) C== SAL = monthly salinity [psu] WRITE(fn,'(A,A,I10.10,A)') p1, 'Stave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny SAL(i,j,k,m)=tmp(i,j) enddo enddo enddo close(100) C== TFLX = monthly net surface heat flux [W/m^2] WRITE(fn,'(A,A,I10.10,A)') p1, 'tFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny TFLX(i,j,m)=-tmp(i,j) enddo enddo close(100) C== H2OFLX = monthly net surface freshwater flux [m/y] WRITE(fn,'(A,A,I10.10,A)') p1, 'sFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny C-- convert from PSU.kg/m^2/s to m/yr H2OFLX(i,j,m)=-tmp(i,j)*nb_seconds_per_year/999.8/SAL(i,j,1,m) enddo enddo close(100) C== Eul_U = monthly Eulerian Zonal Velocity [m/s] WRITE(fn,'(A,A,I10.10,A)') p1, 'uVeltave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny Eul_U(i,j,k,m)=tmp(i,j) enddo enddo enddo close(100) C== WS_x = monthly Zonal Wind Stress [N/m^2] WRITE(fn,'(A,A,I10.10,A)') p1, 'uFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny WS_x(i,j,m)=tmp(i,j) enddo enddo close(100) C== Eul_V = monthly Eulerian Meridional Velocity [m/s] WRITE(fn,'(A,A,I10.10,A)') p1, 'vVeltave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny Eul_V(i,j,k,m)=tmp(i,j) enddo enddo enddo close(100) C== WS_y = monthly Meridional Wind Stress [N/m^2] WRITE(fn,'(A,A,I10.10,A)') p1, 'vFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny WS_y(i,j,m)=tmp(i,j) enddo enddo close(100) enddo call write_nc_phys( & 'ECCO1_Stationary_1','MITgcm_checkpoint51n_post', & secs_per_month, & nx, ny, nz, & PT, SAL, TFLX, H2OFLX, & nz, & nx, ny, Eul_U, WS_x, & nx, ny, Eul_V, WS_y, & has_eddy, & nx, ny, Eddy_U, & nx, ny, Eddy_V) C=========================================================== print*,'Quasi-stationary, last year' C=========================================================== do m=1,12 step=(StationaryYears-1)*nb_timesteps_per_year+ & m*(nb_timesteps_per_year/12) C== PT = monthly potential temperature [C] WRITE(fn,'(A,A,I10.10,A)') p1, 'Ttave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny PT (i,j,k,m)=tmp(i,j) enddo enddo enddo close(100) C== SAL = monthly salinity [psu] WRITE(fn,'(A,A,I10.10,A)') p1, 'Stave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny SAL(i,j,k,m)=tmp(i,j) enddo enddo enddo close(100) C== TFLX = monthly net surface heat flux [W/m^2] WRITE(fn,'(A,A,I10.10,A)') p1, 'tFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny TFLX(i,j,m)=-tmp(i,j) enddo enddo close(100) C== H2OFLX = monthly net surface freshwater flux [m/y] WRITE(fn,'(A,A,I10.10,A)') p1, 'sFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny C-- convert from PSU.kg/m^2/s to m/yr H2OFLX(i,j,m)=-tmp(i,j)*nb_seconds_per_year/999.8/SAL(i,j,1,m) enddo enddo close(100) C== Eul_U = monthly Eulerian Zonal Velocity [m/s] WRITE(fn,'(A,A,I10.10,A)') p1, 'uVeltave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny Eul_U(i,j,k,m)=tmp(i,j) enddo enddo enddo close(100) C== WS_x = monthly Zonal Wind Stress [N/m^2] WRITE(fn,'(A,A,I10.10,A)') p1, 'uFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny WS_x(i,j,m)=tmp(i,j) enddo enddo close(100) C== Eul_V = monthly Eulerian Meridional Velocity [m/s] WRITE(fn,'(A,A,I10.10,A)') p1, 'vVeltave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny Eul_V(i,j,k,m)=tmp(i,j) enddo enddo enddo close(100) C== WS_y = monthly Meridional Wind Stress [N/m^2] WRITE(fn,'(A,A,I10.10,A)') p1, 'vFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny WS_y(i,j,m)=tmp(i,j) enddo enddo close(100) enddo call write_nc_phys( & 'ECCO1_Stationary_3001','MITgcm_checkpoint51n_post', & secs_per_month, & nx, ny, nz, & PT, SAL, TFLX, H2OFLX, & nz, & nx, ny, Eul_U, WS_x, & nx, ny, Eul_V, WS_y, & has_eddy, & nx, ny, Eddy_U, & nx, ny, Eddy_V) C IF ( p1 .NE. ' ' ) THEN ENDIF C=========================================================== print*,'Time-dependent average of last 10 years (232-241)' C=========================================================== IF ( p2 .NE. ' ' ) THEN do m=1,12 C== initialize to zero do i=1,nx do j=1,ny do k=1,nz PT (i,j,k,m)=0 SAL(i,j,k,m)=0 enddo TFLX (i,j,m)=0 H2OFLX(i,j,m)=0 enddo enddo do i=1,nx do j=1,ny do k=1,nz Eul_U(i,j,k,m)=0 enddo WS_x(i,j,m)=0 enddo enddo do i=1,nx do j=1,ny do k=1,nz Eul_V(i,j,k,m)=0 enddo WS_y(i,j,m)=0 enddo enddo n=0 do year=231,240 n=n+1 step=year*nb_timesteps_per_year+m*(nb_timesteps_per_year/12) C== PT = monthly potential temperature [C] WRITE(fn,'(A,A,I10.10,A)') p2, 'Ttave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny PT(i,j,k,m)=PT(i,j,k,m)+tmp(i,j) enddo enddo enddo close(100) C== SAL = monthly salinity [psu] WRITE(fn,'(A,A,I10.10,A)') p2, 'Stave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny SAL(i,j,k,m)=SAL(i,j,k,m)+tmp(i,j) enddo enddo enddo close(100) C== TFLX = monthly net surface heat flux [W/m^2] WRITE(fn,'(A,A,I10.10,A)') p2, 'tFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny TFLX(i,j,m)=TFLX(i,j,m)-tmp(i,j) enddo enddo close(100) C== H2OFLX = monthly net surface freshwater flux [m/y] WRITE(fn,'(A,A,I10.10,A)') p2, 'sFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny C-- convert from PSU.kg/m^2/s to m/yr H2OFLX(i,j,m)=H2OFLX(i,j,m)- & tmp(i,j)*nb_seconds_per_year/999.8/SAL(i,j,1,m) enddo enddo close(100) C== Eul_U = monthly Eulerian Zonal Velocity [m/s] WRITE(fn,'(A,A,I10.10,A)') p2, 'uVeltave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny Eul_U(i,j,k,m)=Eul_U(i,j,k,m)+tmp(i,j) enddo enddo enddo close(100) C== WS_x = monthly Zonal Wind Stress [N/m^2] WRITE(fn,'(A,A,I10.10,A)') p2, 'uFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny WS_x(i,j,m)=WS_x(i,j,m)+tmp(i,j) enddo enddo close(100) C== Eul_V = monthly Eulerian Meridional Velocity [m/s] WRITE(fn,'(A,A,I10.10,A)') p2, 'vVeltave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny Eul_V(i,j,k,m)=Eul_V(i,j,k,m)+tmp(i,j) enddo enddo enddo close(100) C== WS_y = monthly Meridional Wind Stress [N/m^2] WRITE(fn,'(A,A,I10.10,A)') p2, 'vFluxtave.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny WS_y(i,j,m)=WS_y(i,j,m)+tmp(i,j) enddo enddo close(100) enddo C== normalize do i=1,nx do j=1,ny do k=1,nz PT (i,j,k,m)=PT (i,j,k,m)/n SAL(i,j,k,m)=SAL(i,j,k,m)/n enddo TFLX (i,j,m)=TFLX (i,j,m)/n H2OFLX(i,j,m)=H2OFLX(i,j,m)/n enddo enddo do i=1,nx do j=1,ny do k=1,nz Eul_U(i,j,k,m)=Eul_U(i,j,k,m)/n enddo WS_x(i,j,m)=WS_x(i,j,m)/n enddo enddo do i=1,nx do j=1,ny do k=1,nz Eul_V(i,j,k,m)=Eul_V(i,j,k,m)/n enddo WS_y(i,j,m)=WS_y(i,j,m)/n enddo enddo enddo call write_nc_phys( & 'ECCO_Timedep','MITgcm_checkpoint51n_post', & secs_per_month, & nx, ny, nz, & PT, SAL, TFLX, H2OFLX, & nz, & nx, ny, Eul_U, WS_x, & nx, ny, Eul_V, WS_y, & has_eddy, & nx, ny, Eddy_U, & nx, ny, Eddy_V) C IF ( p2 .NE. ' ' ) THEN ENDIF C=========================================================== print*,'generate ECCO_*_BasisFNCTNS_*.nc and diag_2D files' C=========================================================== C=========================================================== print*,'Quasi-stationary annual mean basis functions' C== dye_arr= Concentration of dye tracer [mol/cm^3] C=========================================================== IF ( p1 .NE. ' ' ) THEN year = StationaryYears step = year * nb_timesteps_per_year do n=1,ndyetrac WRITE(fn,'(A,A,I2.2,A,I10.10,A)') & p1, 'PTRtave', n, '.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny dye_arr(i,j,k,n)=100*100*100*tmp(i,j) enddo enddo enddo close(100) enddo call write_nc_basisfnctns( & 'ECCO','MITgcm_checkpoint51n_post','Stationary', & nx,ny,nz,ndyetrac, & year,nb_seconds_per_year,nb_timesteps_per_year, & dye_arr) C== 2-D diagnostics C dye_flux=dye flux for each tracer (mol/m2/s) time(1)=year do n=1,ndyetrac WRITE(fn,'(A,A,I2.2,A,I10.10,A)') & p1, 'PtrFlux', n, '.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny dye_flux(i,j,n,1)=tmp(i,j) enddo enddo close(100) enddo C cum_dye_flux=cumulative dye flux for each tracer (mol/m2) do n=1,ndyetrac do i=1,nx do j=1,ny cum_dye_flux(i,j,n,1)=0. enddo enddo enddo do step=nb_timesteps_per_year, & (StationaryYears*nb_timesteps_per_year), & nb_timesteps_per_year do n=1,ndyetrac WRITE(fn,'(A,A,I2.2,A,I10.10,A)') & p1, 'PtrFlux', n, '.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny cum_dye_flux(i,j,n,1)=cum_dye_flux(i,j,n,1)+ & tmp(i,j)*nb_seconds_per_year enddo enddo close(100) enddo enddo call write_nc_diag_2D( & 'ECCO','MITgcm_checkpoint51n_post','Stationary', & nx,ny,ndyetrac, & 1, time, dye_flux,cum_dye_flux) C IF ( p1 .NE. ' ' ) THEN ENDIF C=========================================================== print*,'Time-dependent annual mean basis functions' C== 1/10-years for 1775-1965, 1/year for 1970-2005 C== dye_arr= Concentration of dye tracer [mol/cm^3] C=========================================================== IF ( p2 .NE. ' ' ) THEN n=0 do year=1775,1965,10 n=n+1 time(n)=year enddo do year=1970,2005 n=n+1 time(n)=year enddo do irec=1,nrec year=time(irec) step = (year-1764) * nb_timesteps_per_year do n=1,ndyetrac WRITE(fn,'(A,A,I2.2,A,I10.10,A)') & p2, 'PTRtave', n, '.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) do k=1,nz read(100,rec=k) tmp do i=1,nx do j=1,ny dye_arr(i,j,k,n)=100*100*100*tmp(i,j) enddo enddo enddo close(100) enddo call write_nc_basisfnctns( & 'ECCO','MITgcm_checkpoint51n_post','Timedep', & nx,ny,nz,ndyetrac, & year,nb_seconds_per_year,nb_timesteps_per_year, & dye_arr) C== 2-D diagnostics C dye_flux=dye flux for each tracer (mol/m2/s) do n=1,ndyetrac WRITE(fn,'(A,A,I2.2,A,I10.10,A)') & p2, 'PtrFlux', n, '.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny dye_flux(i,j,n,irec)=tmp(i,j) enddo enddo close(100) enddo C cum_dye_flux=cumulative dye flux for each tracer (mol/m2) do n=1,ndyetrac do i=1,nx do j=1,ny if (irec.eq.1) then cum_dye_flux(i,j,n,irec)=0. else cum_dye_flux(i,j,n,irec)=cum_dye_flux(i,j,n,irec-1) endif enddo enddo enddo start_step=nb_timesteps_per_year if (irec.gt.1) & start_step=nb_timesteps_per_year+ & (time(irec-1)-1764)*nb_timesteps_per_year end_step= (year-1764) * nb_timesteps_per_year do step=start_step,end_step,nb_timesteps_per_year do n=1,ndyetrac WRITE(fn,'(A,A,I2.2,A,I10.10,A)') & p2, 'PtrFlux', n, '.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp do i=1,nx do j=1,ny cum_dye_flux(i,j,n,irec)=cum_dye_flux(i,j,n,irec)+ & tmp(i,j)*nb_seconds_per_year enddo enddo close(100) enddo enddo enddo call write_nc_diag_2D( & 'ECCO','MITgcm_checkpoint51n_post','Timedep', & nx,ny,ndyetrac, & nrec, time, dye_flux,cum_dye_flux) C IF ( p2 .NE. ' ' ) THEN ENDIF C=========================================================== print*,'write_nc_diag_0D quasi-stationary diagnostics' C=========================================================== IF ( p1 .NE. ' ' ) THEN WRITE(fn,'(A,A)') p1, 'RAC.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) RAC close(100) WRITE(fn,'(A,A)') p1, 'hFacC.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*nz*4) read(100,rec=1) hFacC close(100) irec=0 do year=1,StationaryYears irec=irec+1 global_time(irec)=year step=year*nb_timesteps_per_year do n=1,ndyetrac C global_tot_dye=global total dye flux for this year for each tracer (mol) WRITE(fn,'(A,A,I2.2,A,I10.10,A)') & p1, 'PtrFlux', n, '.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp close(100) global_tot_dye(n,irec)=0. do i=1,nx do j=1,ny global_tot_dye(n,irec)=global_tot_dye(n,irec)+ & RAC(i,j)*hFacC(i,j,1)*tmp(i,j)*nb_seconds_per_year enddo enddo C global_cum_dye=global cumulative dye flux for each tracer (mol) global_cum_dye(n,irec)=global_tot_dye(n,irec) if (irec.gt.1) global_cum_dye(n,irec) = & global_cum_dye(n,irec) + global_cum_dye(n,irec-1) C global_mean_conc= global mean dye concentration (mol/m-3) WRITE(fn,'(A,A,I2.2,A,I10.10,A)') & p1, 'PTRtave', n, '.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*nz*4) read(100,rec=1) tmp3D close(100) global_mean_conc(n,irec)=0. do i=1,nx do j=1,ny do k=1,nz global_mean_conc(n,irec)=global_mean_conc(n,irec)+ & RAC(i,j)*hFacC(i,j,k)*tmp3D(i,j,k) enddo enddo enddo enddo enddo call write_nc_diag_0D( & 'ECCO','MITgcm_checkpoint51n_post','Stationary', & StationaryYears, global_time, ndyetrac, & global_tot_dye, global_cum_dye, global_mean_conc) C IF ( p1 .NE. ' ' ) THEN ENDIF C=========================================================== print*,'write_nc_diag_0D time-dependent diagnostics' C=========================================================== IF ( p2 .NE. ' ' ) THEN WRITE(fn,'(A,A)') p2, 'RAC.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) RAC close(100) WRITE(fn,'(A,A)') p2, 'hFacC.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*nz*4) read(100,rec=1) hFacC close(100) irec=0 do year=1765,2005 irec=irec+1 global_time(irec)=year step=(year-1764)*nb_timesteps_per_year do n=1,ndyetrac C global_tot_dye=global total dye flux for this year for each tracer (mol) WRITE(fn,'(A,A,I2.2,A,I10.10,A)') & p2, 'PtrFlux', n, '.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*4) read(100,rec=1) tmp close(100) global_tot_dye(n,irec)=0. do i=1,nx do j=1,ny global_tot_dye(n,irec)=global_tot_dye(n,irec)+ & RAC(i,j)*hFacC(i,j,1)*tmp(i,j)*nb_seconds_per_year enddo enddo C global_cum_dye=global cumulative dye flux for each tracer (mol) global_cum_dye(n,irec)=global_tot_dye(n,irec) if (irec.gt.1) global_cum_dye(n,irec) = & global_cum_dye(n,irec) + global_cum_dye(n,irec-1) C global_mean_conc= global mean dye concentration (mol/m-3) WRITE(fn,'(A,A,I2.2,A,I10.10,A)') & p2, 'PTRtave', n, '.', step, '.data' open(100,file=fn,status='old',access='direct', & recl=nx*ny*nz*4) read(100,rec=1) tmp3D close(100) global_mean_conc(n,irec)=0. do i=1,nx do j=1,ny do k=1,nz global_mean_conc(n,irec)=global_mean_conc(n,irec)+ & RAC(i,j)*hFacC(i,j,k)*tmp3D(i,j,k) enddo enddo enddo enddo enddo call write_nc_diag_0D( & 'ECCO','MITgcm_checkpoint51n_post','Timedep', & irec, global_time, ndyetrac, & global_tot_dye, global_cum_dye, global_mean_conc) C IF ( p2 .NE. ' ' ) THEN ENDIF stop end