Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: IMiss=-99, Nvar=5, N365=365 Integer, Parameter :: Nyear=20 Integer, Parameter :: Mwd=15 ! Running mean window in days Integer, Parameter :: Swd=7 ! Change window in day 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 (:,:) :: Mask Real, Allocatable, Dimension (:,:,:) :: PET, ET, sw1, sw2, sw3 Real, Allocatable, Dimension (:,:,:) :: Pet_wd, Et_wd ! , Pet_31, Et_31 Real, Allocatable, Dimension (:,:,:) :: ESI_wd ! , ESI_31 Real, Allocatable, Dimension (:,:,:) :: dESI_wd ! Change in ESI over Swd window steps Real, Allocatable, Dimension (:,:,:) :: zESI ! standardized ESI change anomalies Real, Allocatable, Dimension (:,:,:) :: fd_evt ! FD event (0/1) Real, Allocatable, Dimension (:,:,:) :: Ann_Freq, Ann_Dur ! Nlon,Nlat,Nyear Real, Allocatable, Dimension (:,:) :: Mean_Freq, Mean_Dur ! Nlon,Nlat 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 Integer :: No, iday, ievt !---------------------------------------------------------------------------------------------- Real, Allocatable :: X1d(:), M(:), Spells(:), Start(:), Ends(:) Real :: N_Spell, Imax Integer :: Nper, Period, Nper_day ! Real :: X1(N), M1(N), IMAX1, Spells(N), N_Spell_1, I_Start(N), I_End(N) !---------------------------------------------------------------------------------------------- Fname="esi_zesi_fd_soctrang_2001-2020.nc" Call rd_ncdf_3d (Fname,Lon,Lat,time,Nlon,Nlat,Ntime,fd_evt,"fd") Deallocate (Lon,Lat,time) Fname="esi_soctrang_2001-2020.nc" Call rd_ncdf_3d (Fname,Lon,Lat,time,Nlon,Nlat,Ntime,ESI_wd,"esi") Allocate (Mask(Nlon,Nlat)) Mask(:,:) = ESI_wd(:,:,3990) Where (Mask/=Miss) Mask=1 ! Landpoints Period = 5 ! Giai doan Period nam Nper = Nyear/Period ! So giai doan Nper_day = Period*N365 ! So ngay trong 1 giai doan Allocate (X1d(Nper_day),M(Nper_day),Spells(Nper_day),Start(Nper_day),Ends(Nper_day)) Allocate (Ann_Freq(Nlon,Nlat,Nper)) Allocate (Ann_Dur (Nlon,Nlat,Nper)) Allocate (Mean_Freq(Nlon,Nlat), Mean_Dur(Nlon,Nlat)) Ann_Freq=Miss; Ann_Dur=Miss; Mean_Freq=Miss; Mean_Dur=Miss Call Calc_Chars Fname="FD_Chars_Soc_Trang.nc" Call Write_to_file (Fname) CONTAINS !---------------------------------------------------------------------------------------------- Include "Sub_Spell.f90" !---------------------------------------------------------------------------------------------- Subroutine Calc_Chars Integer :: i,ilon,ilat,it, Nspl Real, Allocatable :: xtmp(:) Allocate (xtmp(Nper)) Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then Do iyr=1,Nper ! 1,Nyear !it1=(iyr-1)*N365+1; it2=it1+N365-1 !X1d(1:N365) = fd_evt(ilon,ilat,it1:it2) it1=(iyr-1)*Nper_day+1; it2=it1+Nper_day-1 X1d(1:Nper_day) = fd_evt(ilon,ilat,it1:it2) !Call NSpell_and_Leng (X1d, N365, M, Imax, Spells, N_Spell, Start, Ends, Miss) Call NSpell_and_Leng (X1d,Nper_day,M, Imax, Spells, N_Spell,Start, Ends, Miss) Nspl=Int(N_Spell) Do i=1,Nspl ! if (Spells(i)<=20) then ! print*,"Nspl=",Nspl,Spells(1:Nspl); ! print*, ilon,ilat,Start(i),Ends(i) ! print"(50I2)",Int(X1d(1:Nper_day)) ! pause ! endif Enddo Ann_Freq(ilon,ilat,iyr) = N_Spell Ann_Dur (ilon,ilat,iyr) = Ave (Spells,Nspl,Miss) ! print*,Int(N_Spell),Spells(1:Int(N_Spell)); pause Enddo xtmp(1:Nper) = Ann_Freq(ilon,ilat,1:Nper) Mean_Freq(ilon,ilat) = Ave(xtmp,Nper,Miss) xtmp(1:Nper) = Ann_Dur (ilon,ilat,1:Nper) Mean_Dur (ilon,ilat) = Ave(xtmp,Nper,Miss) endif Enddo Enddo End Subroutine Calc_Chars !---------------------------------------------------------------------------------------------- Subroutine Write_to_file (Fname) Character (Len=*) :: Fname Integer :: it,ifile,ireg,Ndays_60y Real, Allocatable :: X3d(:,:,:),xtime(:),X2d(:,:),xreg(:) Character (Len=150) :: var_ln, t_unit, lon_nm, lat_nm, z_nm, t_nm Character (Len=10) :: var_unit Allocate (xtime(Nper)) Do it=1,Nper; xtime(it)=it-1; Enddo var_unit="-"; t_unit="years since 2001-01-01 "; lon_nm="lon"; lat_nm="lat"; t_nm="time" write(var_ln,"(A)") "Annual number of FD events " Call wr_ncdf_3d (Fname,Lon,Lat,xtime,Nlon,Nlat,Nper,Ann_Freq,"freq",t_unit,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) write(var_ln,"(A,I2.2,A)") "Annual mean of FD duration " Call wr_ncdf_3d_apd (Fname,Ann_Dur,"dur",Nlon,Nlat,Nper,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) write(var_ln,"(A,I2.2,A)") "Long term mean of frequency" Call wr_ncdf_2d_apd (Fname,Mean_Freq,"mean_fre",Nlon,Nlat, var_ln,var_unit,lon_nm,lat_nm,Miss) write(var_ln,"(A,I2.2,A)") "Long term mean of duration" Call wr_ncdf_2d_apd (Fname,Mean_Dur,"mean_dur",Nlon,Nlat,var_ln,var_unit,lon_nm,lat_nm,Miss) End Subroutine Write_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