! 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 simple_xy_rd use netcdf implicit none ! This is the name of the data file we will read. ! character (len = *), parameter :: FILE_NAME = "wrfout_d01_2016-01-01_00:00:00" ! We are reading 2D data, a 6 x 12 grid. integer, parameter :: NT = 1, NX = 442, NY = 265, NZ = 47, NChar = 19 integer, parameter :: Nsite = 4 ! integer, parameter :: J=101 I=217 from NCL ! Lamont,dryden,caltech,park_falls integer,dimension(Nsite), parameter :: Jsite=(/101,104,97, 189/) ! NCl index, need to integer,dimension(Nsite), parameter :: Isite=(/217,64, 60, 264/) ! plus 1 for fortran index integer :: J, I ! plus 1 for fortran index integer :: II, JJ integer :: 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) :: Tpotential, QVAPOR real,dimension(NX, NY, NZ+1, NT) :: PHB, PH real,dimension(NX, NY, NT) :: HGT real,dimension(NZ) :: pressure_profile, CO2_profile, CO2_BCK_profile, CO2_ANT_profile, CO2_BIO_profile real,dimension(NZ) :: non_stag_profile real,dimension(NZ+1) :: stag_profile real,dimension(NZ-1) :: pressure_profile_edgeT, pressure_profile_edgeB real,dimension(Nsite) :: XCO2, XCO2_BIO, XCO2_BCK, XCO2_ANT character (len = *), parameter :: LVL_NAME = "level" character (len = *), parameter :: LVL_stager_NAME = "level_stagger" character (len = *), parameter :: Xcord_NAME = "west_east" character (len = *), parameter :: Ycord_NAME = "south_north" character (len = *), parameter :: Sitecord_NAME = "n_sites" character (len = *), parameter :: Ccord_NAME = "DateStrLen" character (len = *), parameter :: REC_NAME = "time" character (len = *), parameter :: xco2_NAME="XCO2" character (len = *), parameter :: XCO2_BCK_NAME="XCO2_BCK" character (len = *), parameter :: XCO2_ANT_NAME="XCO2_ANT" character (len = *), parameter :: XCO2_BIO_NAME="XCO2_BIO" character (len = *), parameter :: PHB_prof_NAME="PHB_prof" character (len = *), parameter :: CO2_prof_NAME="CO2_prof" character (len = *), parameter :: T_potential_prof_NAME="T_potential_prof" character (len = *), parameter :: QVAPOR_prof_NAME="QVAPOR_prof" character (len = *), parameter :: height_prof_NAME="height_prof" character (len = *), parameter :: pressure_prof_NAME="pressure_prof" character (len = *), parameter :: PH_prof_NAME="PH_prof" character (len = *), parameter :: Time_NAME="timeUTC" character (len = NChar), dimension(NT) :: WRFtime ! This will be the netCDF ID for the file and data variable. integer :: ncid, varid integer :: time_varid, xco2_varid,XCO2_BIO_varid,XCO2_ANT_varid,XCO2_BCK_varid integer :: CO2_prof_varid, height_prof_varid, pressure_prof_varid integer :: QVAPOR_prof_varid, T_potential_prof_varid integer :: PH_prof_varid, PHB_prof_varid integer :: nc_write_id integer :: countXCO2(2), start_xco2(2), dimids(2) integer :: count2d(2),count2d_s(2), start2d(2) integer :: dimids_prof(2), dimids_prof_s(2) integer :: ilayer ! Loop indexes, and error handling. integer :: x, y, t integer :: reason,NWRFOUTFiles,iStation character(LEN=100), dimension(:), allocatable :: WRFOUTFileNames character(LEN=100) :: FILE_NAME character(len = *), parameter :: output_FILE_NAME_nc="simulated_XCO2_4sites_Fortran.nc" integer :: x_dimid, y_dimid, rec_dimid, char_dimid, site_dimid integer :: lvl_dimid, lvl_stager_dimid ! get the files ! call system('ls ../wrfout*:00:00 > fileContents.txt') call system('ls ../wrfout*01_00:00:00 > fileContents.txt') ! for test open(31,FILE='fileContents.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 allocate(WRFOUTFileNames(NWRFOUTFiles)) open(33,File = 'simulated_XCO2_4sites_Fortran.txt',action="write") WRITE(33,"(A15,4(A45))") "Time_UTC", "Lamont","dryden","caltech","park_falls" WRITE(33,"(A15,4(2A10,2A15))") "Time_UTC", (("XCO2", "XCO2_BCK", "XCO2_BIO", "XCO2_ANT"),JJ = 1, Nsite) call check( nf90_create(output_FILE_NAME_nc, nf90_clobber, nc_write_id) ) call check( nf90_def_dim(nc_write_id, LVL_NAME, NZ, lvl_dimid) ) call check( nf90_def_dim(nc_write_id, LVL_stager_NAME, NZ+1, lvl_stager_dimid) ) call check( nf90_def_dim(nc_write_id, Ccord_NAME, NChar, char_dimid) ) call check( nf90_def_dim(nc_write_id, Sitecord_NAME, Nsite, site_dimid) ) call check( nf90_def_dim(nc_write_id, REC_NAME, NF90_UNLIMITED, rec_dimid) ) call check( nf90_def_var(nc_write_id, time_NAME, NF90_CHAR,(/char_dimid,rec_dimid/), time_varid) ) dimids = (/site_dimid, rec_dimid/) 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 at 4 sites") ) call check( nf90_put_att(nc_write_id, xco2_varid, "sites", "Lamont,dryden,caltech,park_falls") ) call check( nf90_put_att(nc_write_id, xco2_varid, "Isite", "217,64, 60, 264") ) call check( nf90_put_att(nc_write_id, xco2_varid, "Jsite", "101,104,97, 189") ) 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_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_def_var(nc_write_id, XCO2_BCK_NAME, NF90_REAL, dimids, XCO2_BCK_varid) ) call check( nf90_put_att(nc_write_id, XCO2_BCK_varid, "unit", "ppmv") ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!1!!! define profiles !!!!!!!!!!!!!!!!!!!!!!!!!!!! dimids_prof_s = (/lvl_stager_dimid, rec_dimid/) dimids_prof = (/lvl_dimid, rec_dimid/) call check( nf90_def_var(nc_write_id, PHB_prof_NAME, NF90_REAL, dimids_prof_s, PHB_prof_varid) ) call check( nf90_put_att(nc_write_id, PHB_prof_varid, "description", "base-state geopotential at lamont") ) call check( nf90_put_att(nc_write_id, PHB_prof_varid, "unit", "m2 s-2") ) call check( nf90_def_var(nc_write_id, PH_prof_NAME, NF90_REAL, dimids_prof_s, PH_prof_varid) ) call check( nf90_put_att(nc_write_id, PH_prof_varid, "description","perturbation geopotential at lamont") ) call check( nf90_put_att(nc_write_id, PH_prof_varid, "unit", "m2 s-2") ) call check( nf90_def_var(nc_write_id, CO2_prof_NAME, NF90_REAL, dimids_prof, CO2_prof_varid) ) call check( nf90_put_att(nc_write_id, CO2_prof_varid, "description", "CO2 concentration at lamont") ) call check( nf90_put_att(nc_write_id, CO2_prof_varid, "unit", "ppmv") ) call check( nf90_def_var(nc_write_id, height_prof_NAME, NF90_REAL, dimids_prof, height_prof_varid) ) call check( nf90_put_att(nc_write_id, height_prof_varid, "description", "total geopotential height at lamont") ) call check( nf90_put_att(nc_write_id, height_prof_varid, "unit", "m AGL") ) call check( nf90_def_var(nc_write_id, T_potential_prof_NAME, NF90_REAL, dimids_prof, T_potential_prof_varid) ) call check( nf90_put_att(nc_write_id, T_potential_prof_varid, "description", "T_potential at lamont") ) call check( nf90_put_att(nc_write_id, T_potential_prof_varid, "unit", "K") ) call check( nf90_def_var(nc_write_id, QVAPOR_prof_NAME, NF90_REAL, dimids_prof, QVAPOR_prof_varid) ) call check( nf90_put_att(nc_write_id, QVAPOR_prof_varid, "description", "QVAPOR at lamont") ) call check( nf90_put_att(nc_write_id, QVAPOR_prof_varid, "unit", "k/kg") ) call check( nf90_def_var(nc_write_id, pressure_prof_NAME, NF90_REAL, dimids_prof, pressure_prof_varid) ) call check( nf90_put_att(nc_write_id, pressure_prof_varid, "description", "total pressure at lamont") ) call check( nf90_put_att(nc_write_id, pressure_prof_varid, "unit", "mb") ) ! End define mode. call check( nf90_enddef(nc_write_id) ) rewind(31) do II = 1,NWRFOUTFiles read(31,'(a)') WRFOUTFileNames(II) FILE_NAME = WRFOUTFileNames(iI) ! Open the file. NF90_NOWRITE tells netCDF we want read-only access to ! the file. call check( nf90_open(FILE_NAME, NF90_NOWRITE, ncid) ) ! Get the varid of the data variable, based on its name. call check( nf90_inq_varid(ncid, "QVAPOR", varid) ) call check( nf90_get_var(ncid, varid, QVAPOR) ) call check( nf90_inq_varid(ncid, "T", varid) ) call check( nf90_get_var(ncid, varid, Tpotential) ) 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_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, "PH", varid) ) call check( nf90_get_var(ncid, varid, PH) ) call check( nf90_inq_varid(ncid, "PHB", varid) ) call check( nf90_get_var(ncid, varid, PHB) ) call check( nf90_inq_varid(ncid, "HGT", varid) ) call check( nf90_get_var(ncid, varid, HGT) ) 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) ) rec_xhu = II do JJ = 1, Nsite I = Isite(JJ)+1 ! add 1 to fotran index J = Jsite(JJ)+1 ! add 1 to fotran index pressure_profile = pressure(I,J,:,1) pressure_profile_edgeT = (pressure_profile(1:NZ-1) + pressure_profile(2:NZ))/2 ! get top edge pressure, collape by 1 pressure_profile_edgeB = pressure_profile_edgeT pressure_profile_edgeB(1) = (PSFC(I,J,1))/100. pressure_profile_edgeB(2:) = pressure_profile_edgeT(1:NZ-2) CO2_profile = CO2_BIO(I,J,:,1) + CO2_ANT(I,J,:,1) - CO2_BCK(I,J,:,1) CO2_BCK_profile = CO2_BCK(I,J,:,1) CO2_ANT_profile = CO2_ANT(I,J,:,1) - CO2_BCK(I,J,:,1) CO2_BIO_profile = CO2_BIO(I,J,:,1) - CO2_BCK(I,J,:,1) XCO2(JJ) = 0.0 XCO2_BIO(JJ) = 0.0 XCO2_BCK(JJ) = 0.0 XCO2_ANT(JJ) = 0.0 do ilayer = 1, NZ-1 XCO2(JJ) = XCO2(JJ) + (pressure_profile_edgeB(ilayer)- pressure_profile_edgeT(ilayer) ) * CO2_profile(ilayer) XCO2_BIO(JJ) = XCO2_BIO(JJ) + (pressure_profile_edgeB(ilayer)- pressure_profile_edgeT(ilayer) ) * CO2_BIO_profile(ilayer) XCO2_BCK(JJ) = XCO2_BCK(JJ) + (pressure_profile_edgeB(ilayer)- pressure_profile_edgeT(ilayer) ) * CO2_BCK_profile(ilayer) XCO2_ANT(JJ) = XCO2_ANT(JJ) + (pressure_profile_edgeB(ilayer)- pressure_profile_edgeT(ilayer) ) * CO2_ANT_profile(ilayer) end do ; XCO2(JJ) = XCO2(JJ)/ (pressure_profile_edgeB(1)-pressure_profile_edgeT(NZ-1)) XCO2_BIO(JJ) = XCO2_BIO(JJ)/ (pressure_profile_edgeB(1)-pressure_profile_edgeT(NZ-1)) XCO2_ANT(JJ) = XCO2_ANT(JJ)/ (pressure_profile_edgeB(1)-pressure_profile_edgeT(NZ-1)) XCO2_BCK(JJ) = XCO2_BCK(JJ)/ (pressure_profile_edgeB(1)-pressure_profile_edgeT(NZ-1)) if (JJ.eq.1) then !for lamont, write profiles count2d = (/NZ,1 /) count2d_s = (/NZ+1,1 /) start2d = (/1, rec_xhu /) non_stag_profile = QVAPOR(I,J,:,1)*1000. call check( nf90_put_var(nc_write_id, QVAPOR_prof_varid, non_stag_profile, start = start2d, & count = count2d) ) non_stag_profile = Tpotential(I,J,:,1)+300. call check( nf90_put_var(nc_write_id, T_potential_prof_varid, non_stag_profile, start = start2d, & count = count2d) ) call check( nf90_put_var(nc_write_id, pressure_prof_varid, pressure_profile, start = start2d, & count = count2d) ) call check( nf90_put_var(nc_write_id, CO2_prof_varid, CO2_profile, start = start2d, & count = count2d) ) non_stag_profile = (PHB(I,J,1:NZ,1)+PHB(I,J,2:NZ+1,1)+PH(I,J,1:NZ,1)+PH(I,J,2:NZ+1,1))/2./9.81-HGT(I,J,1) call check( nf90_put_var(nc_write_id, height_prof_varid, non_stag_profile, start = start2d, & count = count2d) ) stag_profile = PHB(I,J,:,1) ! call check( nf90_put_var(nc_write_id, PHB_prof_varid, PHB(I,J,:,1), start = start2d, & ! may have issue for some system call check( nf90_put_var(nc_write_id, PHB_prof_varid, stag_profile, start = start2d, & count = count2d_s) ) stag_profile = PH(I,J,:,1) call check( nf90_put_var(nc_write_id, PH_prof_varid, stag_profile, start = start2d, & count = count2d_s) ) ! write(*,*) "lamont site elevation = ", HGT(I,J,1) end if end do ! loop over each TCCON site WRITE(33,"(A15,4(2F10.4,2E15.4))") trim(FILE_NAME(15:27)), ((XCO2(JJ), XCO2_BCK(JJ), XCO2_BIO(JJ), XCO2_ANT(JJ)),JJ = 1, Nsite) WRITE(* ,"(A15,4(2F10.4,2E15.4))") trim(FILE_NAME(15:27)), ((XCO2(JJ), XCO2_BCK(JJ), XCO2_BIO(JJ), XCO2_ANT(JJ)),JJ = 1, Nsite) !!!!!!!!!!!!!!!!! write XCO2 into NetCDF file WRITE(*,'(A,I4,1x,6A)'), "writing record",rec_xhu,WRFtime(1), " in",output_FILE_NAME_nc, " !!!!!!!!!!!!!" call check( nf90_put_var(nc_write_id, time_varid, WRFtime(1), start =(/1,rec_xhu/),& count = (/NChar,1/))) countXCO2= (/ Nsite,1 /) start_xco2 = (/ 1, rec_xhu /) 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_ANT_varid, XCO2_ANT, start = start_XCO2,& count = countXCO2) ) call check( nf90_put_var(nc_write_id, XCO2_BCK_varid, XCO2_BCK, start = start_XCO2,& count = countXCO2) ) end do ! loop over files call check( nf90_close(nc_write_id) ) 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 simple_xy_rd