Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: IMiss=-99, Nfiles=9, Nreg=7, Pentad=5, N365=365, Nwd=5 ! Nwd=Window data Integer, Parameter :: Ncha=4 ! b(1), number of FD; b(2), mean duration of FD; b(3), mean severity of FD; b(4), mean speed of FD Character (Len=25), Parameter :: Fnames(Nfiles)=(/"VN_tp.nc","VN_t2m.nc","VN_tx.nc","VN_tm.nc","VN_eva.nc","VN_pet.nc","VN_sw1.nc","VN_sw2.nc","VN_sw3.nc"/) Character (Len=10), Parameter :: vname(Nfiles)=(/"tp","t2m","mx2t","mn2t","e","pev","swvl1","swvl2","swvl3"/) !Character (Len=25), Parameter :: Fnames(Nfiles)=(/"tp_daily.nc","t2m_daily.nc","tx_daily.nc","tm_daily.nc","eva_daily.nc", & ! "pet_daily.nc","sw1_daily.nc","sw2_daily.nc","sw3_daily.nc"/) !Character (Len=10), Parameter :: vname(Nfiles)=(/"var228","t2m","mx2t","mn2t","var182","var251","swvl1","swvl2","swvl3"/) Integer, Allocatable :: time(:) Integer :: i, LenSt, sys Character (Len=255) :: OutDir, Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2 Character (Len=255), Allocatable :: Files(:) Character (Len=12) :: date0, timeres, idate0, idate1 Character (Len=1) :: Ch Integer :: yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2 Real*8 :: factor, offset, Missing, fill_value Real, Allocatable :: Prec(:,:,:), lon(:), lat(:), lev(:), xtime(:), Regs(:,:) Real, Allocatable, Dimension (:,:,:) :: sw1 ! (Nlon,Nlat,Ntime) Real, Allocatable, Dimension (:,:) :: sw1_Reg ! (Nreg,Ntime) Real, Allocatable, Dimension (:,:,:) :: psw1 ! (Nlon,Nlat,Ntime) Real, Allocatable, Dimension (:,:,:) :: FDevt ! drt(Nlon,Nlat,Ntime) - 1-drought onset, 2-drought recovery, 0-no drought Real, Allocatable, Dimension (:,:,:) :: FDCha ! (Nlon,Nlat, Ncha) Real, Allocatable, Dimension (:,:) :: psw1_Reg ! (Nreg,Ntime) Real, Allocatable, Dimension (:,:) :: FDevt_Reg ! drt(Nreg,Ntime) Real, Allocatable, Dimension (:,:) :: FDCha_Reg ! (Nreg, Ncha) Real, Allocatable, Dimension (:,:,:,:) :: X4d Integer, Allocatable, Dimension (:) :: date_m1 ! (Nfcs) Integer :: Nlon, Nlat, Nlev, Ntime, Ndays1, Ndays, Npentad, Nlon0, Nlat0 ! Nfcs = So lan du bao, date_fcs = yyyymmdd cua ngay DB Integer :: imem, ifcs, iy, im, id !---------------------------------------------------------------------------------------------- Call Read_Regs ("/home/tanpv/Data/ERA5/Mask_VN_0.25.nc") !VN.grid_0.25.nc") Call Read_sw_percentile_data Call Calculating_FD_Chars OutDir="/home/tanpv/Progs/Flash_Dr/NC/" Call Write_data_to_file (trim(OutDir)//"FD_Char_1961-2020.nc",trim(OutDir)//"FD_Char_1961-2020_Reg.nc") CONTAINS !---------------------------------------------------------------------------------------------- Include "FL_Dr.f90" !---------------------------------------------------------------------------------------------- Subroutine Calculating_FD_Chars Integer :: ilon, ilat, it, ireg Real :: X(Ntime), drt(Ntime), b(Ncha) Allocate (FDevt(Nlon,Nlat,Ntime)) Allocate (FDCha(Nlon,Nlat,Ncha)) Allocate (FDevt_Reg(Nreg,Ntime)) Allocate (FDCha_Reg(Nreg,Ncha)) FDevt=Miss; FDCha=Miss; FDevt_Reg=Miss; FDCha_Reg=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then X(1:Ntime)=psw1(ilon,ilat,1:Ntime) Call flashdrought(Ntime,X,b,drt) FDevt(ilon,ilat,1:Ntime)=drt(1:Ntime) FDCha(ilon,ilat,1:Ncha) =b(1:Ncha) endif Enddo Enddo Do ireg=1,Nreg X(1:Ntime)=psw1_Reg(ireg,1:Ntime) Call flashdrought(Ntime,X,b,drt) FDevt_Reg(ireg,1:Ntime)=drt(1:Ntime) FDCha_Reg(ireg,1:Ncha) =b(1:Ncha) Enddo End Subroutine Calculating_FD_Chars !---------------------------------------------------------------------------------------------- Subroutine Read_Regs (Fname) Character (Len=*) :: Fname Integer :: ilon,ilat,Nlo,Nla, i1, i2,i3,i4,i5,i6,i7,i8,i9 Call rd_ncdf_2d (Fname,Lon,Lat,Nlon,Nlat,Regs,"reg") End Subroutine Read_Regs !---------------------------------------------------------------------------------------------- Subroutine Read_sw_percentile_data Integer :: Nt, it, j, iday, ilon,ilat, ifile, Nobs Integer :: Nlo, Nla Real, Allocatable :: xlon(:),xlat(:),xtime(:),Y1d(:) Real :: tmp Fname="Pentad_sw_1961-2020.nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw1,"sw") Deallocate (xlon,xlat,xtime) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,psw1,"swp") Deallocate (xlon,xlat,xtime) Fname="Pentad_sw_1961-2020_Reg.nc" Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw1_Reg,"sw") Deallocate (xlon,xtime) Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,psw1_Reg,"swp") Deallocate (xlon,xtime) Ntime=Nt End Subroutine Read_sw_percentile_data !---------------------------------------------------------------------------------------------- Subroutine Write_data_to_file (Fname,Fnm1) ! ,X4d,X1d,Nlon,Nlat,Nt,Nfcs) Character (Len=*) :: Fname,Fnm1 ! Daily and Moving Pentad ! Dat(Nlon,Nlat,Nfiles,Ndays) Integer :: it,ifile,ireg,Ndays_60y Real, Allocatable :: X3d(:,:,:),xday(:),X2d(:,:),xreg(:),xcha(:) Character (Len=50) :: t_unit, lon_nm, lat_nm, z_nm, t_nm Character (Len=10) :: var_nm(Nfiles)=(/"pr","t2m","tx","tm","eva","pet","sw1","sw2","sw3"/) Character (Len=10) :: var_unit Character (Len=150) :: var_ln(3)=(/"FD chars: b(1), number of FD; b(2), mean duration of FD; b(3), mean severity of FD; b(4), mean speed of FD", & "Daily volumetric soil water percemtile layer 2 (7 - 28cm)", & "Daily volumetric soil water percemtile layer 3 (28 - 100cm)"/) ! Allocate (FDevt(Nlon,Nlat,Ntime)) ! Allocate (FDCha(Nlon,Nlat,Ncha)) ! Allocate (FDevt_Reg(Nreg,Ntime)) ! Allocate (FDCha_Reg(Nreg,Ncha)) Allocate (xday(Ntime),xreg(Nreg),xcha(Ncha)) Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo Do it=1,Ntime; xday(it)=(it-1)*Pentad; Enddo Do it=1,Ncha; xcha(it)=it; Enddo var_unit="-"; t_unit="days since 1961-01-01 "; lon_nm="lon"; lat_nm="lat"; t_nm="cha" Call wr_ncdf_3d (Fname,Lon,Lat,xcha,Nlon,Nlat,Ncha,FDCha,"fd",t_unit,trim(var_ln(1)),var_unit,lon_nm,lat_nm,t_nm,Miss) ! Call wr_ncdf_3d_apd (Fname,psw2,"sw2",Nlon,Nlat,Ntime,trim(var_ln(2)),var_unit,lon_nm,lat_nm,t_nm,Miss) ! Call wr_ncdf_3d_apd (Fname,psw3,"sw3",Nlon,Nlat,Ntime,trim(var_ln(3)),var_unit,lon_nm,lat_nm,t_nm,Miss) ! Areal mean of daily data: Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo Call wr_ncdf_2d (Fnm1,xreg,xcha,Nreg,Ncha,FDCha_Reg,"fd",trim(var_ln(1)),var_unit,"reg","cha",Miss) ! Call wr_ncdf_2d_apd (Fnm1,psw2_Reg,"sw2",Nreg,Ntime,"Areal mean of "//trim(var_ln(2)),var_unit,"reg","day",Miss) ! Call wr_ncdf_2d_apd (Fnm1,psw3_Reg,"sw3",Nreg,Ntime,"Areal mean of "//trim(var_ln(3)),var_unit,"reg","day",Miss) End Subroutine Write_data_to_file !---------------------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Ave.f90" Include "/home/tanpv/Mylibs/Xmax.f90" Include "/home/tanpv/Mylibs/Xmin.f90" Include "/home/tanpv/Mylibs/Xsum.f90" Include "/home/tanpv/Mylibs/Sub_Percentile.f90" Include "/home/tanpv/Mylibs/Percentile_rank.f90" END