! This is part of the netCDF package. ! Copyright 2006 University Corporation for Atmospheric Research/Unidata. ! See COPYRIGHT file for conditions of use. ! This is a simple example which reads a small dummy array, from a ! netCDF data file created by the companion program simple_xy_wr.f90. ! This is intended to illustrate the use of the netCDF fortran 77 ! API. This example program is part of the netCDF tutorial, which can ! be found at: ! http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-tutorial ! Full documentation of the netCDF Fortran 90 API can be found at: ! http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-f90 ! $Id: simple_xy_rd.f90,v 1.7 2006/12/09 18:44:58 russ Exp $ program compute_xco2 use netcdf implicit none ! This is the name of the data file we will read. ! We are reading 2D data, a 6 x 12 grid. integer, parameter :: NT = 1, NX = 442, NY = 265, NZ = 47, NC = 19 !CONUS ! integer, parameter :: NT = 1, NX = 285, NY = 234, NZ = 47, NC = 19 !China ! integer, parameter :: J=101 I=217 from NCL integer, parameter :: J=101+1, I=217+1 !Lamount plus 1 for fortran index integer :: II, rec_xhu ! lat = 36.6 ; ! lon = -97.49 ; ! integer :: data_in(NT, NY, NX) real,dimension(NX, NY, NT) :: PSFC! xhu no idea why would the order totally reversed. real,dimension(NX, NY, NZ, NT) :: CO2_BIO, CO2_ANT, CO2_BCK, P, PB, pressure real,dimension(NX, NY, NZ, NT) :: CO2_OCE real,dimension(NX, NY, NZ) :: pressure_3D, CO2_3D, CO2_BCK_3D, CO2_ANT_3D, CO2_BIO_3D real,dimension(NX, NY, NZ) :: CO2_OCE_3D real,dimension(NX, NY, NZ-1) :: pressure_3D_edgeT, pressure_3D_edgeB real,dimension(NX, NY) :: XCO2, XCO2_BIO, XCO2_BCK, XCO2_ANT real,dimension(NX, NY) :: XCO2_OCE character (len = *), parameter :: Xcord_NAME = "west_east" character (len = *), parameter :: Ycord_NAME = "south_north" character (len = *), parameter :: Ccord_NAME = "DateStrLen" character (len = *), parameter :: REC_NAME = "time" character (len = *), parameter :: xco2_NAME="XCO2" character (len = *), parameter :: xco2_ant_NAME="xco2_ant" character (len = *), parameter :: xco2_bio_NAME="xco2_bio" character (len = *), parameter :: xCO2_OCE_NAME="XCO2_OCE" character (len = *), parameter :: Time_NAME="timeUTC" character (len = NC), dimension(NT) :: WRFtime integer :: x_dimid, y_dimid, rec_dimid, char_dimid ! This will be the netCDF ID for the file and data variable. integer :: ncid, varid integer :: time_varid integer :: xco2_varid, xco2_bio_varid,xco2_ant_varid integer :: xCO2_OCE_varid integer :: nc_write_id integer :: ilayer, ihour integer :: countXCO2(3), start_xco2(3), dimids(3) integer :: x, y, t integer :: reason,NWRFOUTFiles,iStation character(LEN=100), dimension(:), allocatable :: WRFOUTFileNames character(LEN=100) :: wrfout_name character(len = 100) :: output_FILE_NAME_nc character(len = 2) :: ihour_char do ihour = 17, 21 ! do ihour = 18,18 if (ihour.lt.10) then write(ihour_char,'(A1,I1)'), "0",ihour else write(ihour_char,'(I2)'), ihour end if output_FILE_NAME_nc = "simulated_XCO2_domainWide_Fortran_"//ihour_char//"UTC.nc" ! get the files call system("rm fileContentsSingleHour"//ihour_char//".txt") call system("ls ../wrfout*"//ihour_char//":00:00 > fileContentsSingleHour"//ihour_char//".txt") call system("rm "//trim(output_FILE_NAME_nc)) open(31,FILE="fileContentsSingleHour"//ihour_char//".txt",action="read") !how many II = 0 do read(31,FMT="(a)",iostat=reason) if (reason/=0) EXIT II = II+1 end do NWRFOUTFiles =II ! write(verb,'(a,I0)') "Number of WRFOUT files: " , NWRFOUTFiles if (ALLOCATED (WRFOUTFileNames)) DEALLOCATE (WRFOUTFileNames) allocate(WRFOUTFileNames(NWRFOUTFiles)) open(33,File = 'simulated_XCO2_domainWide_Fortran.txt',action="write") WRITE(33,"(A)") "Time_UTC, XCO2, XCO2_BCK, XCO2_BIO, XCO2_ANT" rewind(31) call check( nf90_create(output_FILE_NAME_nc, nf90_clobber, nc_write_id) ) call check( nf90_def_dim(nc_write_id, Xcord_NAME, NX, x_dimid) ) call check( nf90_def_dim(nc_write_id, Ycord_NAME, NY, y_dimid) ) call check( nf90_def_dim(nc_write_id, Ccord_NAME, NC, char_dimid) ) call check( nf90_def_dim(nc_write_id, REC_NAME, NF90_UNLIMITED, rec_dimid) ) dimids = (/x_dimid, y_dimid, rec_dimid/) call check( nf90_def_var(nc_write_id, time_NAME, NF90_CHAR, (/char_dimid,rec_dimid/), time_varid) ) call check( nf90_def_var(nc_write_id, xco2_NAME, NF90_REAL, dimids, xco2_varid) ) call check( nf90_put_att(nc_write_id, xco2_varid, "unit", "ppmv") ) call check( nf90_put_att(nc_write_id, xco2_varid, "description", "column-averaged CO2 concentration") ) call check( nf90_def_var(nc_write_id, xco2_bio_NAME, NF90_REAL, dimids, xco2_bio_varid) ) call check( nf90_put_att(nc_write_id, xco2_bio_varid, "unit", "ppmv") ) call check( nf90_put_att(nc_write_id, xco2_bio_varid, "description", "column-averaged CO2_bio concentration") ) call check( nf90_def_var(nc_write_id, xCO2_OCE_NAME, NF90_REAL, dimids, xCO2_OCE_varid) ) call check( nf90_put_att(nc_write_id, xCO2_OCE_varid, "unit", "ppmv") ) call check( nf90_put_att(nc_write_id, xCO2_OCE_varid, "description", "column-averaged CO2_OCE concentration") ) call check( nf90_def_var(nc_write_id, xco2_ant_NAME, NF90_REAL, dimids, xco2_ant_varid) ) call check( nf90_put_att(nc_write_id, xco2_ant_varid, "unit", "ppmv") ) call check( nf90_put_att(nc_write_id, xco2_ant_varid, "description", "column-averaged CO2_ant concentration") ) ! End define mode. call check( nf90_enddef(nc_write_id) ) start_xco2 = (/ 1, 1, 1 /) do II = 1,NWRFOUTFiles read(31,'(a)') WRFOUTFileNames(II) wrfout_name = WRFOUTFileNames(iI) ! Open the file. NF90_NOWRITE tells netCDF we want read-only access to ! the file. call check( nf90_open(wrfout_name, NF90_NOWRITE, ncid) ) ! Get the varid of the data variable, based on its name. call check( nf90_inq_varid(ncid, "CO2_BIO", varid) ) call check( nf90_get_var(ncid, varid, CO2_BIO) ) call check( nf90_inq_varid(ncid, "CO2_OCE", varid) ) call check( nf90_get_var(ncid, varid, CO2_OCE) ) call check( nf90_inq_varid(ncid, "CO2_ANT", varid) ) call check( nf90_get_var(ncid, varid, CO2_ANT) ) call check( nf90_inq_varid(ncid, "CO2_BCK", varid) ) call check( nf90_get_var(ncid, varid, CO2_BCK) ) call check( nf90_inq_varid(ncid, "P", varid) ) call check( nf90_get_var(ncid, varid, P) ) call check( nf90_inq_varid(ncid, "PB", varid) ) call check( nf90_get_var(ncid, varid, PB) ) call check( nf90_inq_varid(ncid, "PSFC", varid) ) call check( nf90_get_var(ncid, varid, PSFC) ) call check( nf90_inq_varid(ncid, "Times", varid) ) call check( nf90_get_var(ncid, varid, WRFtime) ) pressure = ( P + PB ) * 0.01 ! Close the file, freeing all resources. call check( nf90_close(ncid) ) pressure_3D = pressure(:,:,:,1) pressure_3D_edgeT = (pressure_3D(:,:,1:NZ-1) + pressure_3D(:,:,2:NZ))/2 ! get top edge pressure, collape by 1 pressure_3D_edgeB = pressure_3D_edgeT pressure_3D_edgeB(:,:,1) = (PSFC(:,:,1))/100. pressure_3D_edgeB(:,:,2:) = pressure_3D_edgeT(:,:,1:NZ-2) CO2_3D = CO2_BIO(:,:,:,1) + CO2_ANT(:,:,:,1) - CO2_BCK(:,:,:,1) + CO2_OCE(:,:,:,1) - CO2_BCK(:,:,:,1) CO2_BCK_3D = CO2_BCK(:,:,:,1) CO2_ANT_3D = CO2_ANT(:,:,:,1) - CO2_BCK(:,:,:,1) CO2_BIO_3D = CO2_BIO(:,:,:,1) - CO2_BCK(:,:,:,1) CO2_OCE_3D = CO2_OCE(:,:,:,1) - CO2_BCK(:,:,:,1) XCO2 = 0.0 XCO2_BIO = 0.0 XCO2_OCE = 0.0 XCO2_BCK = 0.0 XCO2_ANT = 0.0 do ilayer = 1, NZ-1 XCO2 = XCO2 + (pressure_3D_edgeB(:,:,ilayer)- pressure_3D_edgeT(:,:,ilayer) ) * CO2_3D(:,:,ilayer) XCO2_BIO = XCO2_BIO + (pressure_3D_edgeB(:,:,ilayer)- pressure_3D_edgeT(:,:,ilayer) ) * CO2_BIO_3D(:,:,ilayer) XCO2_OCE = XCO2_OCE + (pressure_3D_edgeB(:,:,ilayer)- pressure_3D_edgeT(:,:,ilayer) ) * CO2_OCE_3D(:,:,ilayer) XCO2_BCK = XCO2_BCK + (pressure_3D_edgeB(:,:,ilayer)- pressure_3D_edgeT(:,:,ilayer) ) * CO2_BCK_3D(:,:,ilayer) XCO2_ANT = XCO2_ANT + (pressure_3D_edgeB(:,:,ilayer)- pressure_3D_edgeT(:,:,ilayer) ) * CO2_ANT_3D(:,:,ilayer) end do ; XCO2 = XCO2 / (pressure_3D_edgeB(:,:,1)-pressure_3D_edgeT(:,:,NZ-1)) XCO2_BIO = XCO2_BIO/ (pressure_3D_edgeB(:,:,1)-pressure_3D_edgeT(:,:,NZ-1)) XCO2_OCE = XCO2_OCE/ (pressure_3D_edgeB(:,:,1)-pressure_3D_edgeT(:,:,NZ-1)) XCO2_ANT = XCO2_ANT/ (pressure_3D_edgeB(:,:,1)-pressure_3D_edgeT(:,:,NZ-1)) XCO2_BCK = XCO2_BCK/ (pressure_3D_edgeB(:,:,1)-pressure_3D_edgeT(:,:,NZ-1)) WRITE(33,"(A,2(A,F10.3),2(A,E15.3))") trim(wrfout_name(4:)), " ",XCO2(I,J), " ",XCO2_BCK(I,J), " ",XCO2_BIO(I,J), " ",XCO2_ANT(I,J) WRITE(*,"(A,2(A,F10.3),2(A,E15.3))") trim(wrfout_name(4:)), " ",XCO2(I,J), " ",XCO2_BCK(I,J), " ",XCO2_BIO(I,J), " ",XCO2_ANT(I,J) !!!!!!!!!!!!!!!!! write XCO2 into NetCDF file rec_xhu = II WRITE(*,'(A,I4,1x,6A)'), "writing record",rec_xhu,WRFtime(1), " in ",output_FILE_NAME_nc, " !!!!!!!!!!!!!" start_xco2(3) = rec_xhu countXCO2= (/NX, NY,1 /) call check( nf90_put_var(nc_write_id, XCO2_varid, XCO2, start = start_xco2, & count = countXCO2) ) call check( nf90_put_var(nc_write_id, XCO2_bio_varid, XCO2_bio, start = start_xco2, & count = countXCO2) ) call check( nf90_put_var(nc_write_id, XCO2_OCE_varid, XCO2_OCE, start = start_xco2, & count = countXCO2) ) call check( nf90_put_var(nc_write_id, XCO2_ant_varid, XCO2_ant, start = start_xco2, & count = countXCO2) ) call check( nf90_put_var(nc_write_id, time_varid, WRFtime(1), start = (/1,rec_xhu/),& count = (/NC,1/))) end do ! loop over files call check( nf90_close(nc_write_id) ) print*,"finish hour ", ihour end do ! ihour = 19, 21 contains subroutine check(status) integer, intent ( in) :: status if(status /= nf90_noerr) then print *, trim(nf90_strerror(status)) stop "Stopped" end if end subroutine check end program