Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: IMiss=-99, Nreg=7, Nmon=12, Pentad=5, DayOfYrs=365 ! Period = Do dai cua mot thoi doan (=1: 1 nam; =5: 5 nam) ! Integer, Parameter :: Period=1, N5y=Period*DayOfYrs ! 365 days x 1 years (period) Integer, Parameter :: Period=5, N5y=Period*DayOfYrs ! 1825 ! 1825=365 days x 5 years (period) Integer, Parameter :: Ncha=6 ! b(1): number of FD; b(2): mean duration of FD; b(3): mean severity of FD; ! b(4): max speed of FD; b(5): FD intensification rate; b(6): FD intensity Integer, Parameter :: DayOfMons(Nmon)=(/31,28,31,30,31,30,31,31,30,31,30,31/) Integer, Allocatable :: time(:) Integer :: i, LenSt, sys Character (Len=255) :: OutDir, Fname, t_unit, DirIn, Str, time_unit, Fname1, Fnm1, Fnm2 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,sw2,sw3 ! (Nlon,Nlat,Ntime) Real, Allocatable, Dimension (:,:) :: sw1_Reg,sw2_Reg,sw3_Reg ! (Nreg,Ntime) Real, Allocatable, Dimension (:,:,:,:) :: psw5 ! (Nlon,Nlat,Nyrs,N5y) - Tinh cho tung thoi ky 5 nam Real, Allocatable, Dimension (:,:,:,:) :: FDCha5 ! (Nlon,Nlat,Ncha,Nyrs) Real, Allocatable, Dimension (:,:,:) :: psw5_Reg ! (Nreg,Nyrs,N5y) Real, Allocatable, Dimension (:,:,:) :: FDCha5_Reg ! (Nreg,Ncha,Nyrs) Real, Allocatable, Dimension (:,:,:) :: FDFreq ! (Nlon,Nlat,Nmon) - Frequency of occurrence of FD in month Real, Allocatable, Dimension (:,:) :: FDFreq_Reg ! (Nreg,Nmon) - Frequency of occurrence of FD in month Real, Allocatable, Dimension (:,:,:) :: FDNDay ! (Nlon,Nlat,Nmon) - Num FD day in month Real, Allocatable, Dimension (:,:) :: FDNDay_Reg ! (Nreg,Nmon) - Num FD day in month 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, Nyrs, Ndays, Npentad, Nlon0, Nlat0 ! Nfcs = So lan du bao, date_fcs = yyyymmdd cua ngay DB Integer :: imem, ifcs, iy, im, id !---------------------------------------------------------------------------------------------- OutDir="/home/tanpv/Progs/Flash_Dr/NC/" 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 if (Period==5) then Fname=trim(OutDir)//"FD_Char_1961-2020_5yrs_Daily.nc"; Fname1=trim(OutDir)//"FD_Char_1961-2020_5yrs_Daily_Reg.nc" Fnm1 =trim(OutDir)//"FD_Freq_1961-2020_5yrs_Daily.nc"; Fnm2 =trim(OutDir)//"FD_Freq_1961-2020_5yrs_Daily_Reg.nc" Call Write_data_to_file (Fname,Fnm1,Fname1,Fnm2) ! Call Write_data_to_file (trim(OutDir)//"FD_Char_1961-2020_5yrs_Daily.nc",trim(OutDir)//"FD_Char_1961-2020_5yrs_Daily_Reg.nc") else Fname=trim(OutDir)//"FD_Char_1961-2020_1yrs_Daily.nc"; Fname1=trim(OutDir)//"FD_Char_1961-2020_1yrs_Daily_Reg.nc" Fnm1 =trim(OutDir)//"FD_Freq_1961-2020_1yrs_Daily.nc"; Fnm2 =trim(OutDir)//"FD_Freq_1961-2020_1yrs_Daily_Reg.nc" Call Write_data_to_file (Fname,Fnm1,Fname1,Fnm2) ! Call Write_data_to_file (trim(OutDir)//"FD_Char_1961-2020_1yrs_Daily.nc",trim(OutDir)//"FD_Char_1961-2020_1yrs_Daily_Reg.nc") endif CONTAINS !---------------------------------------------------------------------------------------------- !Include "FL_Dr_1.f90" Include "Detect_event.f90" Include "Sub_Spell.f90" !---------------------------------------------------------------------------------------------- Subroutine FD_Chars (X,YFD,N,b) Real, Parameter :: Th40=40., Th20=20., Speed=5. Integer, Parameter :: Lmin=15 ! Do dai toi thieu Lmin days Integer :: N, it, i1, i2, i, j, k, k1, cnt, ic, kid Real :: X(N), Y40(N), Y20(N), b(Ncha), YFD(N) Real :: M40(N), IMAX40, Spells40(N), N_Spell40, Start40(N), End40(N) Real :: M20(N), IMAX20, Spells20(N), N_Spell20, Start20(N), End20(N) Real :: MFD(N), IMAXFD, SpellsFD(N), N_SpellFD, StartFD(N), EndFD(N) Integer :: N40,N20,I40,I20, ist, ien Real :: Score, Spd, Spdx, Spd1, swmin, osev(n), odur(n), ospd(n) Real :: RDi, RD(n) ! FD intensification rate (RD) = Sum[ (sw(i-1)-sw(i)) ]/Duration ! It's modified by: RD = [(sw(ist-1)-swmin)/sw(ist-1)*100]/[it(sw=swmin)-ist+1] ! (Mukherjee, S., & Mishra, A. K. (2022). A multivariate flash drought indicator for identifying global hotspots and associated climate controls. ! Geophysical Research Letters, 49, e2021GL096804. https://doi. org/10.1029/2021GL096804 Logical :: Less,flag print*,"Begin FD_Chars" Less=.True.; Score=1 ! sw <=40 Call Events (X,Y40,N,Th40,Score,Less,Miss) Score=1 ! sw<=20 Call Events (X,Y20,N,Th20,Score,Less,Miss) Call NSpell_and_Leng (Y40,N,M40,Imax40,Spells40,N_Spell40,Start40,End40,Miss) Call NSpell_and_Leng (Y20,N,M20,Imax20,Spells20,N_Spell20,Start20,End20,Miss) N40=N_Spell40; N20=N_Spell20 if (N40<=0) then print*,"No drought"; return endif YFD=0 ! =1: Flash Drought; =0: No cnt=0 ! So FD Loop1: Do i=1,N40 ! Num of spell with X<40% if (Int(Spells40(i))ien) exit ! Out of spell if (Y20(j)/=1) Cycle Loop2 ! Until meet X<20% if (j>1.AND.X(j)/=Miss.AND.X(j)Spdx) Spdx=Spd if (Spd>=Speed.AND.(j-ist+1)>Pentad) then ! Decline speed must be no less than 5%/pentad and last at least one Pentad if (flag) then; cnt=cnt+1; flag=.false.; endif k=j Loop4: Do k=k+1; if (Y20(k)/=1.OR.k>ien) then; k=k-1; exit; endif ! Until X>20% if (k<=ien.AND.X(k)/=Miss.AND.X(k)Spdx) Spdx=Spd1 Enddo odur(cnt)=k-ist+1 ! Duration in days osev(cnt)=0 Do it=ist,k osev(cnt)=osev(cnt)+(Th40-X(it)) Enddo ospd(cnt)=Spdx ! Maximum speed RDi=0 ! Intensification rate if (kid>ist) then Do it=ist,kid if (it==ist) then RDi=RDi+(Th40-X(it))/Th40*100. else RDi=RDi+(X(it-1)-X(it))/X(it-1)*100. endif Enddo endif if (kid>ist) then RD(cnt)=RDi/Real(kid-ist+1) else RD(cnt)=Miss endif endif Enddo Loop2 if (odur(cnt)1) then; cnt=cnt-1; else; cnt=0; endif; else YFD(ist:k)=1. endif Enddo Loop1 if (cnt==0) then b=0. else ! Trung binh nam: b(1)=cnt/Real(Period); ! Number of FD per period (1 or 5 year) b(2)=Xsum(odur,cnt,Miss)/real(cnt)/Real(Period); ! Mean duration (day/event) b(3)=Xsum(osev,cnt,Miss)/real(cnt)/Real(Period); ! Mean severity (%/event) b(4)=Xsum(ospd,cnt,Miss)/real(cnt)/Real(Period); ! Mean of maximum decline speed (%/Pentad) b(5)=Xsum(RD,cnt,Miss)/real(cnt)/Real(Period) ! Mean intensification rate (%/event) b(6)=b(3)/b(2) ! Mean intensity (%/day) endif print*,"End of FD_Chars: ",b End Subroutine FD_Chars !---------------------------------------------------------------------------------------------- Subroutine FD_Freq (X,N,Y,Yfq) ! N=Nyrs*DayOfYrs: Do dai chuoi tat ca cac ngay trong giai doan NC Integer :: N, imon, iyr, Nyear, it, it1, i1, i2, i, j, k, k1, Nday Real :: X(N), Y(Nmon), Yfq(Nmon), X1(DayOfYrs), Y1(N/DayOfYrs,Nmon), Yw(N/DayOfYrs) Real :: tmp iyr=0; Nyear=N/DayOfYrs Do it1=1,N,DayOfYrs i1=it1; i2=i1+DayOfYrs-1 X1(1:DayOfYrs)=X(i1:i2) iyr=iyr+1 it=0 Do imon=1,Nmon Nday=DayOfMons(imon) i=0 Do k=1,Nday it=it+1; if (X1(it)==1) i=i+1 Enddo Y1(iyr,imon)=i ! So ngay trong thang co FD Enddo Enddo Do imon=1,Nmon Yw(1:Nyear)=Y1(1:Nyear,imon) Y(imon)=Ave(Yw,Nyear,Miss) tmp=0 Do j=1,Nyear if (Y1(j,imon)>Pentad) tmp=tmp+1 Enddo Yfq(imon)=tmp/Real(Nyear) Enddo End Subroutine FD_Freq !---------------------------------------------------------------------------------------------- Subroutine Calculating_FD_Chars Integer :: ilon, ilat, it, it1,it2, ireg, Nt Real :: X(Ntime), drt(Ntime), b(Ncha), Y(Nmon), Y1(Nyrs), Yfrq(Nmon) Real :: XX(Ntime*Nyrs) ! b(1), number of FD; b(2), mean duration of FD; b(3), mean severity of FD; b(4), mean speed of FD Nt=Ntime*Nyrs Allocate (FDCha5(Nlon,Nlat,Ncha,Nyrs)) Allocate (FDCha5_Reg(Nreg,Ncha,Nyrs)) Allocate (FDNDay(Nlon,Nlat,Nmon),FDNDay_Reg(Nreg,Nmon)) Allocate (FDFreq(Nlon,Nlat,Nmon),FDFreq_Reg(Nreg,Nmon)) !Allocate (FDevt(Nlon,Nlat,Ntime),FDevt_Reg(Nreg,Ntime)) FDCha5=Miss; FDCha5_Reg=Miss; FDFreq=Miss; FDFreq_Reg=Miss; FDNDay=Miss; FDNDay_Reg=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do it=1,Nyrs X(1:Ntime)=psw5(ilon,ilat,it,1:Ntime) !Call flashdrought(Ntime,X,b,drt) Call FD_Chars (X,drt,Ntime,b) !b(1) = b(1)/real(Pentad) ! number of FD per year FDCha5(ilon,ilat,1:Ncha,it) =b(1:Ncha) !FDevt (ilon,ilat,1:Ntime)=drt(1:Ntime) ! Call FD_Freq (drt,Ntime,Y) ! Yfrq(it,1:Nmon)=Y(1:Nmon) it1=(it-1)*Ntime+1; it2=it*Ntime; XX(it1:it2)=drt(1:Ntime) Enddo !Do imon=1,Nmon ! Y1(1:Nyrs)=Yfrq(1:Nyrs); FDFreq(ilon,ilat,imon)=Ave(Y1,Nyrs,Miss) !Enddo Call FD_Freq (XX,Nt,Y,Yfrq); FDNDay(ilon,ilat,1:Nmon)=Y(1:Nmon); FDFreq(ilon,ilat,1:Nmon)=Yfrq(1:Nmon) endif Enddo Enddo Do ireg=1,Nreg Do it=1,Nyrs X(1:Ntime)=psw5_Reg(ireg,it,1:Ntime) !Call flashdrought(Ntime,X,b,drt) Call FD_Chars (X,drt,Ntime,b) !b(1) = b(1)/real(Pentad) ! number of FD per year FDCha5_Reg(ireg,1:Ncha,it) =b(1:Ncha) !FDevt_Reg (ireg,1:Ntime)=drt(1:Ntime) ! Call FD_Freq (drt,Ntime,Y) ! Yfrq(it,1:Nmon)=Y(1:Nmon) it1=(it-1)*Ntime+1; it2=it*Ntime; XX(it1:it2)=drt(1:Ntime) Enddo !Do imon=1,Nmon ! Y1(1:Nyrs)=Yfrq(1:Nyrs); FDFreq_Reg(ireg,imon)=Ave(Y1,Nyrs,Miss) !Enddo Call FD_Freq (XX,Nt,Y,Yfrq); FDNDay_Reg(ireg,1:Nmon)=Y(1:Nmon); FDFreq_Reg(ireg,1:Nmon)=Yfrq(1:Nmon); 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, i1,i2 Integer :: Nlo, Nla Real, Allocatable :: xlon(:),xlat(:),xtime(:),Y1d(:) Real :: tmp Fname=trim(OutDir)//"percentile_sw_1961-2020_daily.nc" ! gridpoint Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw1,"sw01") Deallocate (xlon,xlat,xtime) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw2,"sw02") Deallocate (xlon,xlat,xtime) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,sw3,"sw03") Deallocate (xlon,xlat,xtime) ! 3-layer mean: Allocate (psw1(nlo,nla,nt)); psw1=Miss Do ilon=1,Nlon Do ilat=1,Nlat Do ireg=1,Nreg if (Regs(ilon,ilat)>0) then Do it=1,Nt j=0; tmp=0 if (sw1(ilon,ilat,it)/=Miss) then; j=j+1; tmp=tmp+sw1(ilon,ilat,it); endif if (sw2(ilon,ilat,it)/=Miss) then; j=j+1; tmp=tmp+sw2(ilon,ilat,it); endif if (sw3(ilon,ilat,it)/=Miss) then; j=j+1; tmp=tmp+sw3(ilon,ilat,it); endif if (j>0) psw1(ilon,ilat,it)=tmp/real(j) Enddo endif Enddo Enddo Enddo ! Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nt,psw1,"swp") ! Deallocate (xlon,xlat,xtime) Fname=trim(OutDir)//"percentile_sw_1961-2020_daily_Reg.nc" ! Regional mean Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw1_Reg,"sw01") Deallocate (xlon,xtime) Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw2_Reg,"sw02") Deallocate (xlon,xtime) Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,sw3_Reg,"sw03") Deallocate (xlon,xtime) ! 3-layer mean: Allocate (psw1_Reg(Nreg,Nt)); psw1_Reg=Miss Do ireg=1,Nreg Do it=1,Nt j=0; tmp=0 if (sw1_Reg(ireg,it)/=Miss) then; j=j+1; tmp=tmp+sw1_Reg(ireg,it); endif if (sw2_Reg(ireg,it)/=Miss) then; j=j+1; tmp=tmp+sw2_Reg(ireg,it); endif if (sw3_Reg(ireg,it)/=Miss) then; j=j+1; tmp=tmp+sw3_Reg(ireg,it); endif if (j>0) psw1_Reg(ireg,it)=tmp/real(j) Enddo Enddo ! Call rd_ncdf_2d (Fname,xlon,xtime,nlo,nt,psw1_Reg,"swp") ! Deallocate (xlon,xtime) Nyrs=Nt/N5y Ntime=N5y Allocate (psw5(Nlon,Nlat,Nyrs,N5y),psw5_Reg(Nreg,Nyrs,N5y)) Do ilon=1,Nlon Do ilat=1,Nlat Do j=1,Nyrs i1=(j-1)*N5y+1; i2=j*N5y psw5(ilon,ilat,j,1:N5y)=psw1(ilon,ilat,i1:i2) Enddo Enddo Enddo Do ireg=1,Nreg Do j=1,Nyrs i1=(j-1)*N5y+1; i2=j*N5y psw5_Reg(ireg,j,1:N5y)=psw1_Reg(ireg,i1:i2) Enddo Enddo End Subroutine Read_sw_percentile_data !---------------------------------------------------------------------------------------------- Subroutine Write_data_to_file (Fname,Fname1,Fnm1,Fnm2) ! ,X4d,X1d,Nlon,Nlat,Nt,Nfcs) Character (Len=*) :: Fname,Fname1,Fnm1,Fnm2 ! Daily and Moving Pentad ! Dat(Nlon,Nlat,Nfiles,Ndays) Integer :: it,ifile,ireg,Ntime Real, Allocatable :: X3d(:,:,:),xday(:),X2d(:,:),xreg(:),xcha(:) Character (Len=50) :: t_unit, lon_nm, lat_nm, z_nm, t_nm Character (Len=10) :: var_unit Character (Len=150) :: var_ln(3)=(/"FD chars are b(1): num of FD per year (event); b(2): mean duration (day);"// & " b(3): mean severity (%); b(4): mean of maximum decline speed (%/pentad)"// & " b(5): intensification rate (%); b(6): intensity (%/)", & "FD frequency (number of FD days per month", & "Probability of occurrence"/) ! b(1)=cnt/Real(Period); ! Number of FD per period (1 or 5 year) ! b(2)=Xsum(odur,cnt,Miss)/real(cnt)/Real(Period); ! Mean duration (day/event) ! b(3)=Xsum(osev,cnt,Miss)/real(cnt)/Real(Period); ! Mean severity (%/event) ! b(4)=Xsum(ospd,cnt,Miss)/real(cnt)/Real(Period); ! Mean of maximum decline speed (%/Pentad) ! b(5)=Xsum(RD,cnt,Miss)/real(cnt)/Real(Period) ! Mean intensification rate (%/event) ! b(6)=b(3)/b(2) ! Mean intensity (%/day) ! Allocate (FDCha5(Nlon,Nlat,Ncha,Nyrs)) ! Allocate (FDCha5_Reg(Nreg,Ncha,Nyrs)) Ntime=Nyrs Allocate (xday(Ntime),xreg(Nreg),xcha(Ncha)) Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo Do it=1,Ntime; xday(it)=(it-1)*Period; Enddo Do it=1,Ncha; xcha(it)=it; Enddo var_unit="-"; t_unit="years since 1961-01-01 "; lon_nm="lon"; lat_nm="lat"; z_nm="cha"; t_nm="year" Call wr_ncdf_4d (Fname,Lon,Lat,xcha,xday,Nlon,Nlat,Ncha,Nyrs,FDCha5,"fd", & t_unit, trim(var_ln(1)), var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) Deallocate (xday); Allocate (xday(Nmon)) Do it=1,Nmon; xday(it)=(it); Enddo; t_unit="months since 1961-01-01 "; var_unit="days/mon"; t_nm="month" Call wr_ncdf_3d (Fname1,Lon,Lat,xday,Nlon,Nlat,Nmon,FDFreq,"frq",t_unit,trim(var_ln(3)),var_unit,lon_nm,lat_nm,t_nm,Miss) Call wr_ncdf_3d_apd (Fname1,FDNDay,"nday",Nlon,Nlat,Nmon,trim(var_ln(2)),"-",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 lon_nm="reg"; lat_nm="cha"; t_nm="year" Call wr_ncdf_3d (Fnm1,xreg,xcha,xday,Nreg,Ncha,Nyrs,FDCha5_Reg,"fd",t_unit,trim(var_ln(1)),var_unit,lon_nm,lat_nm,t_nm,Miss) var_unit="days/mon";lon_nm="reg"; lat_nm="month" Call wr_ncdf_2d (Fnm2,xreg,xday,Nreg,Nmon,FDFreq_Reg,"frq",trim(var_ln(3)),var_unit,lon_nm,lat_nm,Miss) Call wr_ncdf_2d_apd (Fnm2,FDNDay_Reg,"nday",Nreg,Nmon,trim(var_ln(2)),"-",lon_nm,lat_nm,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