Make a netcdf file by fortran (example by open boundary file for FVCOM)-优快云博客

本文链接：https://blog.youkuaiyun.com/weixin_42222933/article/details/144440479

First, compile and run the following f90 code to generate a netcdf file containing the target variables.

      program main
      use netcdf
      implicit none
      !character(len=*),parameter :: FILE_NAME = "sample_obc.nc"
      !character(len=*),parameter :: FILE_NAME = "sample_obc.nc"
      !character(len=*),parameter :: FILE_NAME = "sample_obc.nc"
      !character(len=*),parameter :: FILE_NAME = "sample_obc.nc"
      character(len=*),parameter :: FILE_NAME = "sample_obc.nc"
      integer :: ncid
!===============================================================================!
      integer,parameter :: dims0=0,dims1=1,dims2=2,dims3=3,nrecs=1,DStrLen=26
      !integer,parameter :: nthree = 3, nobc = 1182, siglev = 101,siglay=100
      !integer,parameter :: nthree = 3, nobc = 1182, siglev = 37,siglay=36
      !integer,parameter :: nthree = 3, nobc = 1148, siglev = 37,siglay=36
      !integer,parameter :: nthree = 3, nobc = 1106, siglev = 37,siglay=36
      !integer,parameter :: nthree = 3, nobc = 357, siglev = 37,siglay=36
      integer,parameter :: nthree = 3, nobc = 2, siglev = 21,siglay=20
      integer :: dimids0(dims0)
      integer :: dimids1a(dims1),dimids1b(dims1),dimids1c(dims1),dimids1d(dims1)
      integer :: dimids2a(dims2),dimids2b(dims2),dimids2c(dims2),dimids2d(dims2)
      integer :: dimids3a(dims3)
!===============================================================================!
!============================================================================!
      integer :: nobc_did,siglev_did,DStrLen_did,rec_did,three_did,siglay_did
!============================================================================!
      integer :: start1(dims1)
      integer :: count1a(dims1),count1b(dims1),count1c(dims1),count1d(dims1)
      integer :: start2(dims2)
      integer :: count2a(dims2),count2b(dims2),count2c(dims2)
      integer :: start3(dims3)
      integer :: count3a(dims3)
!============================================================================!
      integer,dimension(nobc) :: obc_nodes,obc_h,elevation
      integer,dimension(nobc) :: lon,lat
      integer,dimension(nobc,siglay) :: obc_siglay
      integer,dimension(nobc,siglev) :: obc_siglev
      integer,dimension(nobc,siglay) :: obc_temp,obc_salinity
      integer,dimension(nobc,siglay) :: obc_u,obc_v
      integer,dimension(1) :: Itime,Itime2
      real,dimension(1) :: time
      character(len=DStrlen) :: Times
!============================================================================!
!============================================================================!
      integer :: obc_nodes_vid,obc_h_vid,obc_siglay_vid,obc_siglev_vid
      integer :: obc_temp_vid,obc_salinity_vid,elevation_vid
      integer :: obc_u_vid,obc_v_vid
      integer :: time_vid,Itime_vid,Itime2_vid,Times_vid
      integer :: lon_vid,lat_vid
!============================================================================!
!============================================================================!
      integer :: i,j,k,l,m,n,o,p
      time=55562.0
      Itime=55562
      Itime2=0
      Times="1857-11-17T00:00:00.000000"
      obc_nodes=0
      obc_h=0.0
      elevation=0.0
      obc_siglay=0.0
      obc_siglev=0.0
      obc_temp=0.0
      obc_salinity=0.0
      obc_u=0.0
      obc_v=0.0
      lon=0.0
      lat=0.0
!============================================================================!
!==================================================================================!
!=================================create filename==================================!
!==================================================================================!

      call check(nf90_create(FILE_NAME,nf90_clobber,ncid))
          print *,"finished create filename"

!==================================================================================!
!===============================define dimensions==================================!
!==================================================================================!
      call check(nf90_def_dim(ncid,'nobc',nobc,nobc_did))
      call check(nf90_def_dim(ncid,'siglay',siglay,siglay_did))
      call check(nf90_def_dim(ncid,'siglev',siglev,siglev_did))
      call check(nf90_def_dim(ncid,'DateStrLen',DStrLen,DStrLen_did))
      call check(nf90_def_dim(ncid,'time',NF90_UNLIMITED,rec_did))
          print *, "finished define dimensions"

!==================================================================================!
!============================define variables(varids)==============================!
!==================================================================================!
      dimids1a = (/ nobc_did /)
      dimids1b = (/ rec_did /)
      dimids2a = (/ nobc_did,rec_did /)
      dimids2b = (/ nobc_did,siglay_did /)
      dimids2c = (/ nobc_did,siglev_did /)
      dimids2d = (/ DStrLen_did,rec_did /)
      dimids3a = (/ nobc_did,siglay_did,rec_did /)

      !call check(nf90_def_var(ncid,"iint",NF90_INT,dimids1a,iint_vid))
      call check(nf90_def_var(ncid,"time",NF90_REAL,dimids1b,time_vid))
      call check(nf90_def_var(ncid,"Itime",NF90_INT,dimids1b,Itime_vid))
      call check(nf90_def_var(ncid,"Itime2",NF90_INT,dimids1b,Itime2_vid))
      call check(nf90_def_var(ncid,"Times",NF90_CHAR,dimids2d,Times_vid))
      call check(nf90_def_var(ncid,"obc_nodes",NF90_INT,dimids1a,obc_nodes_vid))
      call check(nf90_def_var(ncid,"obc_h",NF90_REAL,dimids1a,obc_h_vid))
      call check(nf90_def_var(ncid,"obc_siglay",NF90_REAL,dimids2b,obc_siglay_vid))
      call check(nf90_def_var(ncid,"obc_siglev",NF90_REAL,dimids2c,obc_siglev_vid))
      call check(nf90_def_var(ncid,"obc_temp",NF90_REAL,dimids3a,obc_temp_vid))
      call check(nf90_def_var(ncid,"obc_salinity",NF90_REAL,dimids3a,obc_salinity_vid))
      call check(nf90_def_var(ncid,"obc_u",NF90_REAL,dimids3a,obc_u_vid))
      call check(nf90_def_var(ncid,"obc_v",NF90_REAL,dimids3a,obc_v_vid))
      call check(nf90_def_var(ncid,"elevation",NF90_REAL,dimids2a,elevation_vid))
      call check(nf90_def_var(ncid,"lon",NF90_REAL,dimids1a,lon_vid))
      call check(nf90_def_var(ncid,"lat",NF90_REAL,dimids1a,lat_vid))
          print *, "finished define variables"

!==================================================================================!
!============================input variable attribute==============================!
!==================================================================================!
      call check(nf90_put_att(ncid,lon_vid,"long_name","nodal longitude"))
      call check(nf90_put_att(ncid,lon_vid,"standard_name","longitude"))
      call check(nf90_put_att(ncid,lon_vid,"units","degrees_east"))
!==================================lon=============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,lat_vid,"long_name","nodal latitude"))
      call check(nf90_put_att(ncid,lat_vid,"standard_name","latitude"))
      call check(nf90_put_att(ncid,lat_vid,"units","degrees_north"))
!==================================lat=============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,obc_nodes_vid,"long_name","Open Boundary Node Number"))
      call check(nf90_put_att(ncid,obc_nodes_vid,"grid","obc_grid"))
!==================================obc_nodes=============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,obc_h_vid,"long_name","Open Boundary Depth"))
      call check(nf90_put_att(ncid,obc_h_vid,"grid","obc_grid"))
!==================================obc_h=============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,obc_siglay_vid,"long_name","ocean_sigma/general_coordinate"))
      call check(nf90_put_att(ncid,obc_siglay_vid,"grid","obc_grid"))
!==================================obc_siglay=============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,obc_siglev_vid,"long_name","ocean_sigma/general_coordinate"))
      call check(nf90_put_att(ncid,obc_siglev_vid,"grid","obc_grid"))
!==================================obc_siglev=============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,obc_temp_vid,"long_name","sea_water_temperature"))
      call check(nf90_put_att(ncid,obc_temp_vid,"units","Celcius"))
      call check(nf90_put_att(ncid,obc_temp_vid,"grid","obc_grid"))
!==================================obc_temp=============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,obc_salinity_vid,"long_name","sea_water_salinity"))
      call check(nf90_put_att(ncid,obc_salinity_vid,"units","PSU"))
      call check(nf90_put_att(ncid,obc_salinity_vid,"grid","obc_grid"))
!==================================obc_salinity=============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,obc_u_vid,"long_name","u-conponent of current"))
      call check(nf90_put_att(ncid,obc_u_vid,"units","m/sec"))
      call check(nf90_put_att(ncid,obc_u_vid,"grid","obc_grid"))
!==================================obc_u=============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,obc_v_vid,"long_name","v-conponent of current"))
      call check(nf90_put_att(ncid,obc_v_vid,"units","m/sec"))
      call check(nf90_put_att(ncid,obc_v_vid,"grid","obc_grid"))
!==================================obc_v=============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,elevation_vid,"long_name","Open Boundary Elevation"))
      call check(nf90_put_att(ncid,elevation_vid,"units","meters"))
      call check(nf90_put_att(ncid,elevation_vid,"grid","obc_grid"))
!==================================elevation=============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,time_vid,"long_name","time"))
      call check(nf90_put_att(ncid,time_vid,"units","days since 1858-11-17 00:00:00"))
      call check(nf90_put_att(ncid,time_vid,"format","modified julian day (MJD)"))
      call check(nf90_put_att(ncid,time_vid,"time_zone","UTC"))
!==================================time============================================!
!==================================================================================!
      call check(nf90_put_att(ncid,Itime_vid,"units","days since 1858-11-17 00:00:00"))
      call check(nf90_put_att(ncid,Itime_vid,"format","modified julian day (MJD)"))
      call check(nf90_put_att(ncid,Itime_vid,"time_zone","UTC"))
!==================================Itime===========================================!
!==================================================================================!
      call check(nf90_put_att(ncid,Itime2_vid,"units","msec since 00:00:00"))
      call check(nf90_put_att(ncid,Itime2_vid,"time_zone","UTC"))
!==================================Itime2==========================================!
!==================================================================================!
!==================================================================================!
      call check(nf90_put_att(ncid,Times_vid,"time_zone","UTC"))
!==================================Times===========================================!
!==================================================================================!
        call check(nf90_put_att(ncid,nf90_global,"title","FVCOM open boundary"))
        call check(nf90_put_att(ncid,nf90_global,"source","FVCOM_3.0"))
        call check(nf90_put_att(ncid,nf90_global,"modeler","Xue's Group @MTU"))

          print *, "finished input variable attribute"
!==================================================================================!
      call check(nf90_enddef(ncid))

!==================================================================================!
!================================input data========================================!
!==================================================================================!
         call check(nf90_put_var(ncid,lon_vid,lon))
         call check(nf90_put_var(ncid,lat_vid,lat))
         call check(nf90_put_var(ncid,obc_nodes_vid,obc_nodes))
         call check(nf90_put_var(ncid,obc_h_vid,obc_h))
         call check(nf90_put_var(ncid,obc_siglay_vid,obc_siglay))
         call check(nf90_put_var(ncid,obc_siglev_vid,obc_siglev))
      start1 = (/ 1 /)
      count1a = (/ 1 /)
      do i=1,nrecs
      start1=i
         call check(nf90_put_var(ncid,time_vid,time,start=start1,count=count1a))
         call check(nf90_put_var(ncid,Itime_vid,Itime,start=start1,count=count1a))
         call check(nf90_put_var(ncid,Itime2_vid,Itime2,start=start1,count=count1a))
      end do

      start2 = (/ 1, 1 /)
      count2a = (/ nobc, 1 /)
      count2b = (/ DStrLen, 1 /)
      do i=1,nrecs
      start3(2)=i
         call check(nf90_put_var(ncid,elevation_vid,elevation,start=start2,count=count2a))
         call check(nf90_put_var(ncid,Times_vid,Times,start=start2,count=count2b))
      end do

      start3 = (/ 1, 1, 1 /)
      count3a = (/ nobc, siglay, 1 /)
      do i=1,nrecs
      start3(3)=i
         call check(nf90_put_var(ncid,obc_temp_vid,obc_temp,start=start3,count=count3a))
         call check(nf90_put_var(ncid,obc_salinity_vid,obc_salinity,start=start3,count=count3a))
         call check(nf90_put_var(ncid,obc_u_vid,obc_u,start=start3,count=count3a))
         call check(nf90_put_var(ncid,obc_v_vid,obc_v,start=start3,count=count3a))
      end do
            print *,"finished data input"
!============================================================================!
!============================================================================!
!============================================================================!
!============================================================================!
      
      call check(nf90_close(ncid))
      
      print *,"finished creating file ", FILE_NAME


      end program

      subroutine check(status)
      use netcdf
      integer, intent ( in) :: status
           if(status /=nf90_noerr)Then
              write(*,*) trim(nf90_strerror(status))
              stop
           end If
      end subroutine check
                                                                                                                                                                   251,7         Bot

Run the f90 by:

ifort -I/pexue6/local_px/netcdf4.6.1/include makenetcdf_obc.f90 -L/pexue6/local_px/netcdf4.6.1/lib -lnetcdf -lnetcdff  -o makenetcdf_obc

Second, run the following script

obcfile="input/2005_leem_fine_julian_obc.nc"
ncatted -O -h -a history,global,d,, ${obcfile}
ncatted -O -h -a type,global,d,, ${obcfile}
ncatted -O -a type,global,a,c,"FVCOM TIME SERIES ELEVATION FORCING FILE" ${obcfile}

Comments:

ncatted: This is the command used to modify NetCDF attributes.
-O: This option stands for “overwrite.” It tells ncatted to overwrite the existing NetCDF file with the modified version. Without this flag, the command would attempt to create a new file.
-h: This option tells ncatted to preserve the header information in the NetCDF file while modifying the attributes. It ensures that the header is not deleted or overwritten.
-a history,global,d,,: This is the attribute modification specification:

history, type: The name of the attribute to modify.
global: This indicates that the attribute is a global attribute (applies to the whole file).
d: This specifies the type of the modification. In this case, d means delete. So, it deletes the history attribute.
,,: The commas here indicate that no new values are being assigned to the history attribute. In this case, it simply deletes the attribute.

-a type,global,a,c,"FVCOM TIME SERIES ELEVATION FORCING FILE":

This specifies the modification to be made:
type: The name of the attribute to be modified.
global: The attribute is a global attribute, meaning it applies to the entire NetCDF file (not to a specific variable).
a: The operation to perform on the attribute. The a means add or append.
c: The data type of the attribute. In this case, c indicates that the attribute value is a character string.
"FVCOM TIME SERIES ELEVATION FORCING FILE": The value to assign to the type attribute. This is the string that will be set as the value of the type attribute.