Use NetCDF

      Integer, Parameter :: Nmon=12, Nday=31, Yr1=1980, Yr2=2020
      Real, Parameter :: Miss=-99.0, Thre1=35.0

      Character (Len=255) :: Fname, Data_Dir
      Integer :: Nst, Nyear 

      Real, Allocatable :: Tx(:,:,:,:)
      Real, Allocatable :: Lon(:), Lat(:), Reg(:), H(:)
      Character (Len=20), Allocatable :: Stname(:)
      Real, Allocatable :: xyear(:),xmon(:),StnID(:),xday(:)
      Real, Allocatable :: Hdm(:,:,:), Hdy(:,:)  ! Hot day in mon, year
      Real, Allocatable :: Data_avail_days(:,:,:)

      Data_Dir="/home/tanpv/OBS_Stn/Daily/NC/"
      Fname=trim(Data_Dir)//"Tx_Daily_VN.nc"
      Call Read_Data (Fname,"Tx")
      Call Calc_Num_of_Hot_Day
      Call Write_Hdy_to_NC ("NumberOfHotDay.nc")
      CONTAINS
      !------------------------------------------------------------------------
      Subroutine Write_Hdy_to_NC (Fname)
         Character (Len=*) :: Fname
         !Character (Len=20) :: Stn(Nst)
         Character (Len=100) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm

         t_unit="months since 1980-01-01 "
         var_ln="Number of hot day in each month"; var_unit="days/month"
         lon_nm="stn"; lat_nm="year"; t_nm="mon"
         !Do ist=1,Nst; Stn(ist) = Stname(ist); Enddo

         Call wr_ncdf_3d (Fname,StnID,xyear,xmon,Nst,Nyear,Nmon,Hdm, &
              "hdm",t_unit, var_ln, var_unit, lon_nm,lat_nm,t_nm,Miss)
         var_ln="Number of hot day in each year"; var_unit="days/year"
         Call wr_ncdf_2d_apd (Fname,Hdy,"hdy",Nst,Nyear, var_ln, &
              var_unit,lon_nm,lat_nm,Miss)
         var_ln="Longtitude of Stations"; var_unit="Degress in East"
         Call wr_ncdf_1d_apd (Fname,Lon,"lon",Nst,var_ln, var_unit, lon_nm, Miss)
         var_ln="Latitude of Stations"; var_unit="Degress in North"
         Call wr_ncdf_1d_apd (Fname,Lat,"lat",Nst,var_ln, var_unit, lon_nm, Miss)
         var_ln="Sub-region of Stations"; var_unit="-"
         Call wr_ncdf_1d_apd (Fname,Reg,"reg",Nst,var_ln, var_unit, lon_nm, Miss)
         var_ln='Name of Stations'
         Call wr_ncdf_1d_apd_char (Fname,Stname,"stname",Nst, var_ln,"-", lon_nm,Miss)
         ! wr_ncdf_1d_apd_char (filenm,var,var_nm, nlon, var_long_name, var_unit, lon_nm, flag)

         ! Write infor data to text file:
         Open (3,file="NumOfDay_In_mon_Data_Available.txt", &
                       Status="Unknown")
         Do ist=1,Nst
            write (3,"(I3,A5,I2,2A)") ist," Reg ",Int(Reg(ist))," Stn: ",trim(Stname(ist))
            write (3,"(A4,12(A4,I2.2))") "Yrs",("Mon_",imon,imon=1,Nmon)
            Do iyr=1,Nyear !Yr1,Yr2
               write (3,"(I4.4,12I6)") iyr+Yr1-1, (Int(Data_avail_days(ist,iyr,imon)),imon=1,Nmon)
            Enddo
         Enddo
         Close(3)
      End Subroutine Write_Hdy_to_NC
      !------------------------------------------------------------------------
      Subroutine Calc_Num_of_Hot_Day
         Integer :: iday,imon,iyr, ist, Ndys
         Integer :: ic1, ic2, mon,yrs, ic3
         ! Hdm(:,:,:), Hdy(:,:)

         Allocate (Hdm(Nst,Nyear,Nmon),Hdy(Nst,Nyear))
         Allocate (Data_avail_days(Nst,Nyear,Nmon))
         !Allocate (Hdm(Nst,Yr1:Yr2,Nmon),Hdy(Nst,Yr1:Yr2))
         !Allocate (Data_avail_days(Nst,Yr1:Yr2,Nmon))
         Do ist=1,Nst
            Do iyr=1,Nyear ! Yr1,Yr2 ! 1,Nyear
               ic1=0
               Do imon=1,Nmon
                  ic2=0; ic3=0
                  mon=imon; yrs=iyr; Ndys=DayOfMon(mon,yrs)
                  Do iday=1,Ndys
                     if (Tx(ist,iyr,imon,iday)/=Miss) ic3=ic3+1
                     if (Tx(ist,iyr,imon,iday)/=Miss.AND.Tx(ist,iyr,imon,iday)>=Thre1) then
                        ic2=ic2+1; ic1=ic1+1
                     endif
                  Enddo
                  Hdm(ist,iyr,imon)=ic2
                  Data_avail_days(ist,iyr,imon)=ic3
               Enddo
               Hdy(ist,iyr)=ic1
            Enddo
         Enddo
      End Subroutine Calc_Num_of_Hot_Day
      !------------------------------------------------------------------------
      Subroutine Read_Data (Fname,vname)
         Character (Len=*) :: Fname,vname
         !Real, Allocatable :: StnID(:),xday(:)
         Integer :: Nm,Nd

         Call rd_ncdf_4d (Fname,StnID,xyear,xmon,xday,Nst,Nyear,Nm,Nd,Tx,trim(vname))
         print*,xyear;pause
         Deallocate(StnID)
         Call rd_ncdf_1d (Fname,StnID,Nst,Lon,"Lon")
         Deallocate(StnID)
         Call rd_ncdf_1d (Fname,StnID,Nst,Lat,"Lat")
         Deallocate(StnID)
         Call rd_ncdf_1d (Fname,StnID,Nst,Reg,"reg")
         Deallocate(StnID)
         Call rd_ncdf_1d (Fname,StnID,Nst,H,"h")
         Deallocate(StnID)
         Call rd_ncdf_1d_Char (Fname,StnID,Nst,Stname,"Stname")
         !print*, Stname(:)
         !Deallocate(StnID)
      End Subroutine Read_Data
      !------------------------------------------------------------------------
      Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/DayOfMon.f90"
      END