!     Can be driven by gen_atmos.pl
!     Output: sflux_[air,rad,prc]_1.???.nc

!     Compile with: 
!     ifort -Bstatic -O3 -o gen_atmos gen_atmos.f90 -Vaxlib -I/usr/local/netcdf/include/ -L/usr/local/netcdf/lib -lnetcdf
      program fgennc
      implicit none
      include 'netcdf.inc'

! error status return
      integer  iret
! netCDF id
      integer  ncid
! dimension ids
      integer  nx_grid_dim
      integer  ny_grid_dim
      integer  time_dim
! dimension lengths
!      integer  nx_grid_len
!      integer  ny_grid_len
      integer  time_len
!      parameter (nx_grid_len = 349)
!      parameter (ny_grid_len = 277)
      parameter (time_len = NF_UNLIMITED)
! variable ids
      integer  time_id
      integer  lon_id
      integer  lat_id
      integer  uwind_id
      integer  vwind_id
      integer  prmsl_id
      integer  stmp_id
      integer  spfh_id
      integer  dlwrf_id
      integer  dswrf_id
      integer  prate_id

! rank (number of dimensions) for each variable
      integer  time_rank
      integer  lon_rank
      integer  lat_rank
      integer  uwind_rank
      integer  vwind_rank
      integer  prmsl_rank
      integer  stmp_rank
      integer  spfh_rank
      integer  dlwrf_rank
      integer  dswrf_rank
      integer  prate_rank
      parameter (time_rank = 1)
      parameter (lon_rank = 2)
      parameter (lat_rank = 2)
      parameter (uwind_rank = 3)
      parameter (vwind_rank = 3)
      parameter (prmsl_rank = 3)
      parameter (stmp_rank = 3)
      parameter (spfh_rank = 3)
      parameter (dlwrf_rank = 3)
      parameter (dswrf_rank = 3)
      parameter (prate_rank = 3)

! variable shapes
      integer  time_dims(time_rank)
      integer  lon_dims(lon_rank)
      integer  lat_dims(lat_rank)
      integer  uwind_dims(uwind_rank)
      integer  vwind_dims(vwind_rank)
      integer  prmsl_dims(prmsl_rank)
      integer  stmp_dims(stmp_rank)
      integer  spfh_dims(spfh_rank)
      integer  dlwrf_dims(dlwrf_rank)
      integer  dswrf_dims(dswrf_rank)
      integer  prate_dims(prate_rank)

! data variables
      real, allocatable ::  lon(:,:),lat(:,:),times(:),data1(:,:)
! attribute vectors
      integer  intval(4)

!     My own variables
      character(len=50) fname
      character(len=21) date_in
      character(len=4) year_char
      character(len=2) mon_char
      character(len=2) day_char
      integer num_times !# of time records in each file
      integer ni,nj,i_time,data_start(3), data_count(3),i,j,k,ifiletype

!     Inputs
      open(10,file='gen_atmos.in',status='old')
      read(10,*) ifiletype !1: air; 2: rad; 3: prc
      read(10,'(a)')fname !output name; e.g., 'sflux_air_1.002.nc'
      read(10,*) intval(1:3) !year (e.g. 1990); month (no leading 0s); day (no leading 0s)
      close(10)

      if(ifiletype<1.or.ifiletype>3) stop "Unknown ifiletype"
      intval(4)=0 !starting hour offset from UTC
      write(year_char,'(i4)')intval(1)
      write(mon_char,'(i2.2)')intval(2)
      write(day_char,'(i2.2)')intval(3)
      date_in='days since '//year_char//'-'//mon_char//'-'//day_char !in UTC

      num_times=8 !# of time records in each file
      ni=2 !dimension in longitude
      nj=3 !dimension in lattitude

!     Allocate arrays
      allocate(lon(ni,nj),lat(ni,nj),data1(ni,nj),times(num_times)) 

!     Write times (>=0 but can exceed 1)
      times=(/0., 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875/) !in days; UTC

!     Write lon/lat (degrees_east and degrees_north)
      do i=1,ni
        do j=1,nj
          lon(i,j)=-140+(0.+140.)/(ni-1)*(i-1)
          lat(i,j)=0.+(60.-0)/(nj-1)*(j-1)
        enddo !j
      enddo !i

! enter define mode
      iret = nf_create(trim(fname), NF_CLOBBER, ncid)
      call check_err(iret)
! define dimensions
      iret = nf_def_dim(ncid, 'nx_grid', ni, nx_grid_dim)
      call check_err(iret)
      iret = nf_def_dim(ncid, 'ny_grid', nj, ny_grid_dim)
      call check_err(iret)
      iret = nf_def_dim(ncid, 'time', NF_UNLIMITED, time_dim)
      call check_err(iret)
! define variables
      time_dims(1) = time_dim
      iret = nf_def_var(ncid, 'time', NF_REAL, time_rank, time_dims, time_id)
      call check_err(iret)
      lon_dims(2) = ny_grid_dim
      lon_dims(1) = nx_grid_dim
      iret = nf_def_var(ncid, 'lon', NF_REAL, lon_rank, lon_dims, lon_id)
      call check_err(iret)
      lat_dims(2) = ny_grid_dim
      lat_dims(1) = nx_grid_dim
      iret = nf_def_var(ncid, 'lat', NF_REAL, lat_rank, lat_dims, lat_id)
      call check_err(iret)
      uwind_dims(3) = time_dim
      uwind_dims(2) = ny_grid_dim
      uwind_dims(1) = nx_grid_dim

      if(ifiletype==1) then
!       air
        iret = nf_def_var(ncid, 'uwind', NF_REAL, uwind_rank, uwind_dims,uwind_id)
        call check_err(iret)
        vwind_dims(3) = time_dim
        vwind_dims(2) = ny_grid_dim
        vwind_dims(1) = nx_grid_dim
        iret = nf_def_var(ncid, 'vwind', NF_REAL, vwind_rank, vwind_dims,vwind_id)
        call check_err(iret)
        prmsl_dims(3) = time_dim
        prmsl_dims(2) = ny_grid_dim
        prmsl_dims(1) = nx_grid_dim
        iret = nf_def_var(ncid, 'prmsl', NF_REAL, prmsl_rank, prmsl_dims,prmsl_id)
        call check_err(iret)
        stmp_dims(3) = time_dim
        stmp_dims(2) = ny_grid_dim
        stmp_dims(1) = nx_grid_dim
        iret = nf_def_var(ncid, 'stmp', NF_REAL, stmp_rank, stmp_dims, stmp_id)
        call check_err(iret)
        spfh_dims(3) = time_dim
        spfh_dims(2) = ny_grid_dim
        spfh_dims(1) = nx_grid_dim
        iret = nf_def_var(ncid, 'spfh', NF_REAL, spfh_rank, spfh_dims, spfh_id)
        call check_err(iret)
      else if(ifiletype==2) then
!       rad
        dlwrf_dims(3) = time_dim
        dlwrf_dims(2) = ny_grid_dim
        dlwrf_dims(1) = nx_grid_dim
        iret = nf_def_var(ncid, 'dlwrf', NF_REAL, dlwrf_rank, dlwrf_dims,dlwrf_id)
        call check_err(iret)
        dswrf_dims(3) = time_dim
        dswrf_dims(2) = ny_grid_dim
        dswrf_dims(1) = nx_grid_dim
        iret = nf_def_var(ncid, 'dswrf', NF_REAL, dswrf_rank, dswrf_dims,dswrf_id)
        call check_err(iret)
      else
!       prcp
        prate_dims(3) = time_dim
        prate_dims(2) = ny_grid_dim
        prate_dims(1) = nx_grid_dim
        iret = nf_def_var(ncid, 'prate', NF_REAL, prate_rank, prate_dims,prate_id)
        call check_err(iret)
      endif

! assign attributes
      iret = nf_put_att_text(ncid, time_id, 'long_name', 4, 'Time')
      call check_err(iret)
      iret = nf_put_att_text(ncid, time_id, 'standard_name', 4, 'time')
      call check_err(iret)
      iret = nf_put_att_text(ncid, time_id, 'units', 21, date_in)
!'
      call check_err(iret)
      iret = nf_put_att_int(ncid, time_id, 'base_date', NF_INT, 4, intval)
      call check_err(iret)
      iret = nf_put_att_text(ncid, lon_id, 'long_name', 9, 'Longitude')
      call check_err(iret)
      iret = nf_put_att_text(ncid, lon_id, 'standard_name', 9, 'longitude')
      call check_err(iret)
      iret = nf_put_att_text(ncid, lon_id, 'units', 12, 'degrees_east')
      call check_err(iret)
      iret = nf_put_att_text(ncid, lat_id, 'long_name', 8, 'Latitude')
      call check_err(iret)
      iret = nf_put_att_text(ncid, lat_id, 'standard_name', 8, 'latitude')
      call check_err(iret)
      iret = nf_put_att_text(ncid, lat_id, 'units', 13, 'degrees_north')
      call check_err(iret)

      if(ifiletype==1) then !air
        iret = nf_put_att_text(ncid, uwind_id, 'long_name', 39, 'Surface Eastward Air Velocity (10m AGL)')
        call check_err(iret)
        iret = nf_put_att_text(ncid, uwind_id, 'standard_name', 13, 'eastward_wind')
        call check_err(iret)
        iret = nf_put_att_text(ncid, uwind_id, 'units', 3, 'm/s')
        call check_err(iret)
        iret = nf_put_att_text(ncid, vwind_id, 'long_name', 40, 'Surface Northward Air Velocity (10m AGL)')
        call check_err(iret)
        iret = nf_put_att_text(ncid, vwind_id, 'standard_name', 14, 'northward_wind')
        call check_err(iret)
        iret = nf_put_att_text(ncid, vwind_id, 'units', 3, 'm/s')
        call check_err(iret)
        iret = nf_put_att_text(ncid, prmsl_id, 'long_name', 23, 'Pressure reduced to MSL')
        call check_err(iret)
        iret = nf_put_att_text(ncid, prmsl_id, 'standard_name', 25, 'air_pressure_at_sea_level')
        call check_err(iret)
        iret = nf_put_att_text(ncid, prmsl_id, 'units', 2, 'Pa')
        call check_err(iret)
        iret = nf_put_att_text(ncid, stmp_id, 'long_name', 32, 'Surface Air Temperature (2m AGL)')
        call check_err(iret)
        iret = nf_put_att_text(ncid, stmp_id, 'standard_name', 15, 'air_temperature')
        call check_err(iret)
        iret = nf_put_att_text(ncid, stmp_id, 'units', 1, 'K')
        call check_err(iret)
        iret = nf_put_att_text(ncid, spfh_id, 'long_name', 34, 'Surface Specific Humidity (2m AGL)')
        call check_err(iret)
        iret = nf_put_att_text(ncid, spfh_id, 'standard_name', 17, 'specific_humidity')
        call check_err(iret)
        iret = nf_put_att_text(ncid, spfh_id, 'units', 1, '1')
        call check_err(iret)

      else if(ifiletype==2) then !rad

        iret = nf_put_att_text(ncid, dlwrf_id, 'long_name', 33, 'Downward Long Wave Radiation Flux')
        call check_err(iret)
        iret = nf_put_att_text(ncid, dlwrf_id, 'standard_name', 40, 'surface_downwelling_longwave_flux_in_air')
        call check_err(iret)
        iret = nf_put_att_text(ncid, dlwrf_id, 'units', 5, 'W/m^2')
        call check_err(iret)
        iret = nf_put_att_text(ncid, dswrf_id, 'long_name', 34, 'Downward Short Wave Radiation Flux')
        call check_err(iret)
        iret = nf_put_att_text(ncid, dswrf_id, 'standard_name', 41, 'surface_downwelling_shortwave_flux_in_air')
        call check_err(iret)
        iret = nf_put_att_text(ncid, dswrf_id, 'units', 5, 'W/m^2')
        call check_err(iret)

      else !prcp
        iret = nf_put_att_text(ncid, prate_id, 'long_name', 26, 'Surface Precipitation Rate')
        call check_err(iret)
        iret = nf_put_att_text(ncid, prate_id, 'standard_name', 18, 'precipitation_flux')
        call check_err(iret)
        iret = nf_put_att_text(ncid, prate_id, 'units', 8, 'kg/m^2/s')
        call check_err(iret)

      endif

      iret = nf_put_att_text(ncid, NF_GLOBAL, 'Conventions', 6, 'CF-1.0')
      call check_err(iret)
!     leave define mode
      iret = nf_enddef(ncid)
      call check_err(iret)
       
!     Write times
      iret = nf_put_vara_real(ncid, time_id, 1, num_times, times)
      call check_err(iret)

!     Write lon/lat
      data_start(1:2)=1
      data_count(1)=ni
      data_count(2)=nj
      iret = nf_put_vara_real(ncid, lon_id, data_start, data_count, lon)
      iret = nf_put_vara_real(ncid, lat_id, data_start, data_count, lat)
      call check_err(iret)

!     Write record variables for time loop
      do i_time=1,num_times
        data_start(1) = 1
        data_start(2) = 1
        data_start(3) = i_time
        data_count(1) = ni
        data_count(2) = nj
        data_count(3) = 1

        if(ifiletype==1) then !air
!------------------------------------------------------------------------------
!       Calculate uwind: Surface Eastward Air Velocity (10m AGL) in m/s
!       store values in data1(1:ni,1:nj) (at lon(1:ni,1:nj) and lat(1:ni,1:nj)).
!       for invalid value use 9.999e+20
        do i=1,ni
          do j=1,nj
            data1(i,j)=(lon(i,j)+lat(i,j))/10+times(i_time)*10
          enddo !j
        enddo !i

        iret = nf_put_vara_real(ncid,uwind_id,data_start,data_count,data1)

!       Calculate vwind: Surface Northward Air Velocity (10m AGL) in m/s
!       store values in data1(1:ni,1:nj) (at lon(1:ni,1:nj) and lat(1:ni,1:nj)).
!       for invalid value use 9.999e+20
        do i=1,ni
          do j=1,nj
            data1(i,j)=(lon(i,j)-lat(i,j))/10-times(i_time)*10
          enddo !j
        enddo !i

        iret = nf_put_vara_real(ncid,vwind_id,data_start,data_count,data1)

!       Calculate prmsl: Pressure reduced to MSL in Pa
!       store values in data1(1:ni,1:nj) (at lon(1:ni,1:nj) and lat(1:ni,1:nj)).
!       for invalid value use 9.999e+20
        do i=1,ni
          do j=1,nj
            data1(i,j)=1.e5
          enddo !j
        enddo !i

        iret = nf_put_vara_real(ncid,prmsl_id,data_start,data_count,data1)

!       Calculate stmp: Surface Air Temperature (2m AGL) in K
!       store values in data1(1:ni,1:nj) (at lon(1:ni,1:nj) and lat(1:ni,1:nj)).
!       for invalid value use 9.999e+20
        do i=1,ni
          do j=1,nj
            data1(i,j)=298
          enddo !j
        enddo !i

        iret = nf_put_vara_real(ncid,stmp_id,data_start,data_count,data1)

!       Calculate spfh: Surface Specific Humidity (2m AGL) (dimensionless)
!       store values in data1(1:ni,1:nj) (at lon(1:ni,1:nj) and lat(1:ni,1:nj)).
!       for invalid value use 9.999e+20
        do i=1,ni
          do j=1,nj
            data1(i,j)=1.e-2
          enddo !j
        enddo !i

        iret = nf_put_vara_real(ncid,spfh_id,data_start,data_count,data1)

!------------------------------------------------------------------------------
        else if(ifiletype==2) then !rad
!------------------------------------------------------------------------------
!       Calculate dlwrf: Downward Long Wave Radiation Flux in W/m^2
!       store values in data1(1:ni,1:nj) (at lon(1:ni,1:nj) and lat(1:ni,1:nj)).
!       for invalid value use 9.999e+20
        do i=1,ni
          do j=1,nj
            data1(i,j)=abs(lon(i,j)+lat(i,j))/10 !400
          enddo !j
        enddo !i

        iret = nf_put_vara_real(ncid,dlwrf_id,data_start,data_count,data1)

!       Calculate dswrf: Downward Short Wave Radiation Flux in W/m^2
!       store values in data1(1:ni,1:nj) (at lon(1:ni,1:nj) and lat(1:ni,1:nj)).
!       for invalid value use 9.999e+20
        do i=1,ni
          do j=1,nj
            data1(i,j)=abs(lon(i,j)-lat(i,j))/10 !40
          enddo !j
        enddo !i

        iret = nf_put_vara_real(ncid,dswrf_id,data_start,data_count,data1)

!------------------------------------------------------------------------------
        else  !prcp
!------------------------------------------------------------------------------
!       Calculate prate: Surface Precipitation Rate in kg/m^2/s
!       store values in data1(1:ni,1:nj) (at lon(1:ni,1:nj) and lat(1:ni,1:nj)).
!       for invalid value use 9.999e+20
        do i=1,ni
          do j=1,nj
            data1(i,j)=abs(lon(i,j)+lat(i,j))/10 !1.e-6
          enddo !j
        enddo !i

        iret = nf_put_vara_real(ncid,prate_id,data_start,data_count,data1)

!------------------------------------------------------------------------------
        endif
      enddo !i_time=1,num_times
       
      iret = nf_close(ncid)
      call check_err(iret)
      stop
      end
       
      subroutine check_err(iret)
      integer iret
      include 'netcdf.inc'
      if (iret .ne. NF_NOERR) then
      print *, nf_strerror(iret)
      stop
      endif
      end