Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: IMiss=-99, Nvar=5, N365=365 Integer, Parameter :: M15=15, M31=31 ! Running mean windows in days Character (Len=5), Parameter :: vname(Nvar)=(/"pev","e","swvl1","swvl2","swvl3"/) Character (Len=25), Parameter :: var_ln(Nvar)=(/"Potential evaporation (m)","lwe_thickness_of_water_evaporation_amount (m)", & "Volumetric soil water layer 1 (m**3 m**-3)","Volumetric soil water layer 2 (m**3 m**-3)", & "Volumetric soil water layer 3 (m**3 m**-3)"/) Real, Allocatable, Dimension (:) :: Lon, Lat, time Integer :: Nlon, Nlat, Ntime Real, Allocatable, Dimension (:,:,:) :: PET, ET, sw1, sw2, sw3 Real, Allocatable, Dimension (:,:,:) :: Pet_15, Et_15, Pet_31, Et_31 Real, Allocatable, Dimension (:,:,:) :: ESI_15, ESI_31 Character (Len=255) :: OutDir, Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2 Character (Len=1) :: Ch Real*8 :: scale_factor, add_offset, missing_value, fill_value !---------------------------------------------------------------------------------------------- Fname="/home/tanpv/Data/ERA5/Thoa_VJU/ev_sw_2001-2020.nc" Call Read_data (Fname) Allocate (Pet_15(Nlon,Nlat,Ntime),Et_15(Nlon,Nlat,Ntime),Pet_31(Nlon,Nlat,Ntime),Et_31(Nlon,Nlat,Ntime)) Call Calc_Pet_Ave (PET,Pet_15,M15); Call Calc_Pet_Ave (PET,Pet_31,M31) Call Calc_Pet_Ave ( ET, Et_15,M15); Call Calc_Pet_Ave ( ET, Et_31,M31) Allocate (ESI_15(Nlon,Nlat,Ntime),ESI_31(Nlon,Nlat,Ntime)) Call Calc_ESI_1 (ESI_15,Pet_15,Et_15,Miss) Call Calc_ESI_1 (ESI_31,Pet_31,Et_31,Miss) OutDir="/home/cordex_sea/Thoa/" Fname=trim(OutDir)//"esi_soctrang_2001-2020.nc" Call Write_data_to_file (Fname) CONTAINS !---------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------- Subroutine GetAtrr (fname, vname, scale_factor, add_offset, missing_value, fill_value) Character (Len=*) :: fname, vname Integer :: ncid, id_t_dim, varid Real*8 :: scale_factor, add_offset, missing_value, fill_value Logical :: status status = nf90_open(fname, nf90_NoWrite, ncid) status = nf90_inq_varid(ncid,trim(vname), varid) status = nf90_get_att(ncid, varid, "scale_factor", scale_factor) status = nf90_get_att(ncid, varid, "add_offset", add_offset) status = nf90_get_att(ncid, varid, "missing_value", missing_value) status = nf90_get_att(ncid, varid, "_FillValue", fill_value) status = nf90_close(ncid) End Subroutine GetAtrr !---------------------------------------------------------------------------------------------- Subroutine Read_data (Fname) Character (Len=*) :: Fname Real, Allocatable, Dimension (:) :: xlon, xlat, xtime Real*8 :: scale_factor, add_offset, missing_value, fill_value Character (Len=5) :: var_nm Integer :: ivar,Nlo,Nla,Nt Do ivar=1,Nvar var_nm=vname(ivar) if (ivar==1) then Call rd_ncdf_3d (Fname,Lon,Lat,time,Nlon,Nlat,Ntime,PET,trim(var_nm)) Call GetAtrr (Fname,var_nm,scale_factor,add_offset,missing_value,fill_value) Where (PET==missing_value) PET=Miss Where (PET==fill_value) PET=Miss Where (PET/=Miss) PET = PET*scale_factor + add_offset elseif (ivar==2) then Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,ET,trim(var_nm)) Call GetAtrr (Fname,var_nm,scale_factor,add_offset,missing_value,fill_value) Where (ET==missing_value) ET=Miss Where (ET==fill_value) ET=Miss Where (ET/=Miss) ET = ET*scale_factor + add_offset Deallocate (xlon,xlat,xtime) elseif (ivar==3) then Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,sw1,trim(var_nm)) Call GetAtrr (Fname,var_nm,scale_factor,add_offset,missing_value,fill_value) Where (sw1==missing_value) sw1=Miss Where (sw1==fill_value) sw1=Miss Where (sw1/=Miss) sw1 = sw1*scale_factor + add_offset Deallocate (xlon,xlat,xtime) elseif (ivar==4) then Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,sw2,trim(var_nm)) Call GetAtrr (Fname,var_nm,scale_factor,add_offset,missing_value,fill_value) Where (sw2==missing_value) sw2=Miss Where (sw2==fill_value) sw2=Miss Where (sw2/=Miss) sw2 = sw2*scale_factor + add_offset Deallocate (xlon,xlat,xtime) else Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,sw3,trim(var_nm)) Call GetAtrr (Fname,var_nm,scale_factor,add_offset,missing_value,fill_value) Where (sw3==missing_value) sw3=Miss Where (sw3==fill_value) sw3=Miss Where (sw3/=Miss) sw3 = sw3*scale_factor + add_offset Deallocate (xlon,xlat,xtime) endif Enddo End Subroutine Read_data !---------------------------------------------------------------------------------------------- !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 Calc_ESI_1 (ESI,Pet,Et,Miss) Integer :: ilon,ilat,it,iyr,iday,i1,i2 Integer, Parameter :: Nyrs=20 ! Ntime/N365 Real :: Miss, Mean(N365), Std(N365) Real :: X1d(Nyrs), REt(Nlon,Nlat,Ntime), Pet(Nlon,Nlat,Ntime), Et(Nlon,Nlat,Ntime) Real :: ESI(Nlon,Nlat,Ntime) REt=Miss; ESI=Miss Do ilon=1,Nlon Do ilat=1,Nlat Do it=1,Ntime if (Et(ilon,ilat,it)/=Miss.AND.Pet(ilon,ilat,it)/=Miss) then REt(ilon,ilat,it)=Et(ilon,ilat,it)/Pet(ilon,ilat,it) endif Enddo Do iday=1,N365 iyr=0 Do it=1,Ntime i1=mod(it,N365); if (i1==0) i1=N365 if (i1==iday) then iyr=iyr+1; X1d(iyr)=REt(ilon,ilat,it) endif Enddo Mean(iday)=Ave(X1d,iyr,Miss); Std(iday)=Stdev(X1d,iyr,Miss) Enddo Do it=1,Ntime iday=mod(it,N365); if (iday==0) iday=N365 if (REt(ilon,ilat,it)/=Miss.AND.Mean(iday)/=Miss.AND.Std(iday)/=Miss) then ESI(ilon,ilat,it)=(REt(ilon,ilat,it)-Mean(iday))/Std(iday) endif Enddo Enddo Enddo End Subroutine Calc_ESI_1 !---------------------------------------------------------------------------------------------- Subroutine Calc_ESI (ESI,Pet,Et,Miss) Integer :: ilon,ilat,it Real :: Miss, Mean, Std Real :: X1d(Ntime), REt(Nlon,Nlat,Ntime), Pet(Nlon,Nlat,Ntime), Et(Nlon,Nlat,Ntime) Real :: ESI(Nlon,Nlat,Ntime) REt=Miss; ESI=Miss Do ilon=1,Nlon Do ilat=1,Nlat Do it=1,Ntime if (Et(ilon,ilat,it)/=Miss.AND.Pet(ilon,ilat,it)/=Miss) then REt(ilon,ilat,it)=Et(ilon,ilat,it)/Pet(ilon,ilat,it) endif Enddo X1d(1:Ntime)=REt(ilon,ilat,1:Ntime) Mean=Ave(X1d,Ntime,Miss); Std=Stdev(X1d,Ntime,Miss) Do it=1,Ntime if (REt(ilon,ilat,it)/=Miss.AND.Mean/=Miss.AND.Std/=Miss) then ESI(ilon,ilat,it)=(REt(ilon,ilat,it)-Mean)/Std endif Enddo Enddo Enddo End Subroutine Calc_ESI !---------------------------------------------------------------------------------------------- Subroutine Calc_Pet_Ave (X,Y,m) Integer :: ilon,ilat,it,N,m, i,j,k Real :: X(Nlon,Nlat,Ntime), Y(Nlon,Nlat,Ntime), Ax(Ntime), Bx(Ntime) Do ilon=1,Nlon Do ilat=1,Nlat Ax(1:Ntime)=X(ilon,ilat,1:Ntime) Call Moving_Ave (Ax,Ntime,m,Bx,Miss) Y(ilon,ilat,1:Ntime)=Bx(1:Ntime) Enddo Enddo End Subroutine Calc_Pet_Ave !---------------------------------------------------------------------------------------------- Subroutine Moving_Ave (X,N,m,Y,Miss) Integer :: N,m, i,j,k, j1, M2 ! m must be an odd number Real :: X(N), Ax(m), Y(N), tmp, Miss M2=M/2; Y=Miss; Do j=M2+1,N-M2 k=0 Do i=j-M2,j+M2 k=k+1; Ax(k)=X(i) Enddo j1=j+M2; Y(j1)=Ave(Ax,M,Miss) ! Gan vao phan tu thu m (o cuoi) !Y(j)=Ave(Ax,M,Miss) ! Gan vao phan tu thu m/2+1 (o giua) Enddo End Subroutine Moving_Ave !---------------------------------------------------------------------------------------------- Subroutine Write_data_to_file (Fname) Character (Len=*) :: Fname Integer :: it,ifile,ireg,Ndays_60y Real, Allocatable :: X3d(:,:,:),xday(:),X2d(:,:),xreg(:) Character (Len=50) :: var_ln, t_unit, lon_nm, lat_nm, z_nm, t_nm Character (Len=10) :: var_unit !Allocate (xday(Npentad),xreg(Nreg)) !Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo Do it=1,Ntime; time(it)=it-1; Enddo var_unit="-"; t_unit="days since 2001-01-01 "; lon_nm="lon"; lat_nm="lat"; t_nm="time" var_ln="ESI with 15-day window running mean of PET and ET" Call wr_ncdf_3d (Fname,Lon,Lat,time,Nlon,Nlat,Ntime,ESI_15,"esi15",t_unit,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="ESI with 31-day window running mean of PET and ET" Call wr_ncdf_3d_apd (Fname,ESI_31,"esi31",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,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/Stdev.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