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, Ftrend, Ftrend_reg 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 (:,:,:) :: Sen_FDCha5 ! (Nlon,Nlat,Ncha) Real, Allocatable, Dimension (:,:,:) :: Prb_FDCha5 ! (Nlon,Nlat,Ncha) Real, Allocatable, Dimension (:,:) :: Sen_FDCha5_Reg ! (Nreg,Ncha) Real, Allocatable, Dimension (:,:) :: Prb_FDCha5_Reg ! (Nreg,Ncha) !---------------------------------------------------------------------------------------------- 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 !---------------------------------------------------------------------------------------------- ! For calculating MFDI: Real, Allocatable, Dimension (:,:,:) :: FDx, FDz ! (Ngrid,Ncha,Nyrs) - Dac trung tai cac diem luoi Real, Allocatable, Dimension (:,:) :: Mx, Sx ! (Ncha,Nyrs) Real, Allocatable, Dimension (:,:) :: Cv, Mz, Sz, MPI ! (Ngrid,Nyrs) Real, Allocatable, Dimension (:,:,:) :: MFDI ! Nlon,Nlat,Nyrs) - Multi FD Index Real, Allocatable, Dimension (:,:) :: Sen_MFDI, Prob_MFDI Integer :: Ngrid !---------------------------------------------------------------------------------------------- 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 Call Calc_MFDI Call Calc_Trend if (Period==5) then Fname=trim(OutDir)//"FD_Char_1961-2020_5yrs_Daily_New.nc"; Fname1=trim(OutDir)//"FD_Char_1961-2020_5yrs_Daily_Reg_New.nc" Fnm1 =trim(OutDir)//"FD_Freq_1961-2020_5yrs_Daily_New.nc"; Fnm2 =trim(OutDir)//"FD_Freq_1961-2020_5yrs_Daily_Reg_New.nc" Ftrend=trim(OutDir)//"Trend_Char_1961-2020_5yrs_Daily_New.nc"; Ftrend_reg=trim(OutDir)//"Trend_Char_1961-2020_5yrs_Daily_Reg_New.nc" Call Write_data_to_file (Fname,Fnm1,Fname1,Fnm2,Ftrend,Ftrend_reg) Call Write_MFDI_to_File (trim(OutDir)//"MFDI_and_Trend_5yrs_Daily_New.nc") ! 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_New.nc"; Fname1=trim(OutDir)//"FD_Char_1961-2020_1yrs_Daily_Reg_New.nc" Fnm1 =trim(OutDir)//"FD_Freq_1961-2020_1yrs_Daily_New.nc"; Fnm2 =trim(OutDir)//"FD_Freq_1961-2020_1yrs_Daily_Reg_New.nc" Ftrend=trim(OutDir)//"Trend_Char_1961-2020_1yrs_Daily_New.nc"; Ftrend_reg=trim(OutDir)//"Trend_Char_1961-2020_1yrs_Daily_Reg_New.nc" Call Write_data_to_file (Fname,Fnm1,Fname1,Fnm2,Ftrend,Ftrend_reg) Call Write_MFDI_to_File (trim(OutDir)//"MFDI_and_Trend_1yrs_Daily_New.nc") ! 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 Calc_MFDI ! Mukherjee, S., & Mishra, A. K. (2022): Integer :: i,j,ilon,ilat,it,icha,ireg, igrd Real, Allocatable :: X1d(:) igrd=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) igrd=igrd+1 Enddo Enddo Ngrid=igrd Allocate (FDx(Ngrid,Ncha,Nyrs),FDz(Ngrid,Ncha,Nyrs),Mx(Ncha,Nyrs),Sx(Ncha,Nyrs)) Allocate (Cv(Ngrid,Nyrs), Mz(Ngrid,Nyrs), Sz(Ngrid,Nyrs), MPI(Ngrid,Nyrs)) Allocate (X1d(Ngrid)) igrd=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then igrd=igrd+1 Do it=1,Nyrs FDx(igrd,1:Ncha,it)=FDCha5(ilon,ilat,1:Ncha,it) Enddo endif Enddo Enddo ! Tinh Trung binh va Do lech chuan: Do j=1,Ncha Do it=1,Nyrs X1d(1:Ngrid)=FDx(1:Ngrid,j,it) Mx(j,it)=Ave(X1d,Ngrid,Miss); Sx(j,it)=Stdev(X1d,Ngrid,Miss) Enddo Enddo Deallocate (X1d) ! Chuan hoa: Do i=1,Ngrid Do j=1,Ncha Do it=1,Nyrs FDz(i,j,it)=100.0+(FDx(i,j,it)-Mx(j,it))/Sx(j,it)*10.0 Enddo Enddo Enddo ! Ket hop (Aggregation): Allocate (X1d(Ncha)) Do i=1,Ngrid Do it=1,Nyrs X1d(1:Ncha)=FDz(i,1:Ncha,it) Mz(i,it)=Ave(X1d,Ncha,Miss); Sz(i,it)=Stdev(X1d,Ncha,Miss) Cv(i,it)=Sz(i,it)/Mz(i,it) MPI(i,it)=Mz(i,it)+Sz(i,it)*Cv(i,it) Enddo Enddo Deallocate (X1d) Allocate (MFDI(Nlon,Nlat,Nyrs)); MFDI=Miss i=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then i=i+1; MFDI(ilon,ilat,1:Nyrs)=MPI(i,1:Nyrs) endif Enddo Enddo End Subroutine Calc_MFDI !---------------------------------------------------------------------------------------------- Subroutine Write_MFDI_to_File (Fname) Character (Len=*) :: Fname Integer :: it,ifile,ireg,Ntime Real, Allocatable :: X3d(:,:,:),xyrs(:),X2d(:,:),xreg(:),xcha(:) Character (Len=50) :: var_unit,t_unit, lon_nm, lat_nm, z_nm, t_nm Character (Len=255) :: var_ln(3)=(/"Multi Flash Drought Index and its trend", & "Trend (Sen's slope) of the MPI", & "Significant level (p-value)"/) Allocate (xyrs(Nyrs)) Do it=1,Nyrs; xyrs(it)=it*Period; Enddo var_unit="-"; t_unit="years since 1961-01-01 "; lon_nm="lon"; lat_nm="lat"; t_nm="year" Call wr_ncdf_3d (Fname,Lon,Lat,xyrs,Nlon,Nlat,Nyrs,MFDI,"mpi",t_unit,trim(var_ln(1)),var_unit,lon_nm,lat_nm,t_nm,Miss) Call wr_ncdf_2d_apd (Fname,Sen_MFDI,"sen",Nlon,Nlat,trim(var_ln(2)), var_unit, lon_nm, lat_nm, Miss) Call wr_ncdf_2d_apd (Fname,Prob_MFDI,"prob",Nlon,Nlat,trim(var_ln(3)),var_unit, lon_nm, lat_nm, Miss) End Subroutine Write_MFDI_to_File !---------------------------------------------------------------------------------------------- Subroutine Calc_Trend Integer :: ilon,ilat,it,icha,ireg, NumYear, N1 Real :: X1d(Nyrs), SenSlope,prob NumYear=10 ! Do dai toi thieu Allocate (Sen_FDCha5(Nlon,Nlat,Ncha),Prb_FDCha5(Nlon,Nlat,Ncha)) Allocate (Sen_FDCha5_Reg(Nreg,Ncha),Prb_FDCha5_Reg(Nreg,Ncha)) Allocate (Sen_MFDI(Nlon,Nlat),Prob_MFDI(Nlon,Nlat)) Sen_FDCha5=Miss; Prb_FDCha5=Miss; Sen_FDCha5_Reg=Miss; Prb_FDCha5_Reg=Miss Sen_MFDI=Miss; Prob_MFDI=Miss ! At gridpoints: Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do icha=1,Ncha X1d(1:Nyrs)=FDCha5(ilon,ilat,icha,1:Nyrs) Call MK_Sen(X1d,Nyrs,NumYear,Miss,SenSlope,prob,N1) Sen_FDCha5(ilon,ilat,icha)=SenSlope Prb_FDCha5(ilon,ilat,icha)=prob Enddo ! For MFDI: X1d(1:Nyrs)=MFDI(ilon,ilat,1:Nyrs) Call MK_Sen(X1d,Nyrs,NumYear,Miss,SenSlope,prob,N1) Sen_MFDI(ilon,ilat)=SenSlope Prob_MFDI(ilon,ilat)=prob endif Enddo Enddo ! For sub-regs: Do ireg=1,Nreg Do icha=1,Ncha X1d(1:Nyrs)=FDCha5_Reg(ireg,icha,1:Nyrs) Call MK_Sen(X1d,Nyrs,NumYear,Miss,SenSlope,prob,N1) Sen_FDCha5_Reg(ireg,icha)=SenSlope Prb_FDCha5_Reg(ireg,icha)=prob Enddo Enddo End Subroutine Calc_Trend !---------------------------------------------------------------------------------------------- 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 ! Note that, Severity is defined as difference of SW between initial time (ist) and terminal time (iter), which iter is time at which SW begin rising 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) then !.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 ospd(cnt)=Spd ! Onset 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. if (X(it-1)>=X(it)) then ! Chi tinh cac truong hop dang giam RDi=RDi+(X(it-1)-X(it))/X(it-1)*100. else RDi=RDi+0.0 endif endif Enddo endif if (kid>ist) then RD(cnt)=RDi/Real(kid-ist+1) else RD(cnt)=0 ! 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); ! Onset 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) ! Trung binh nhieu nam so ngay FD trong tung thang tmp=0 Do j=1,Nyear if (Y1(j,imon)>Pentad) tmp=tmp+1 ! So nam co so ngay FD trong thang > 5 Enddo Yfq(imon)=tmp/Real(Nyear) ! Tan suat cac thang co FD (cac thang co so ngay FD>5) 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,Ftrend,Ftrend_reg) ! ,X4d,X1d,Nlon,Nlat,Nt,Nfcs) Character (Len=*) :: Fname,Fname1,Fnm1,Fnm2,Ftrend,Ftrend_reg ! 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=255) :: 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): Onset decline speed (%/pentad)"// & " b(5): intensification rate (%); b(6): intensity (%/)", & "FD frequency (number of FD days per month", & "Probability of occurrence"/) ! 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*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: Deallocate (xday); Allocate (xday(Ntime)); Do it=1,Ntime; xday(it)=it*Period; Enddo t_unit="years since 1961-01-01 " 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) Deallocate (xday); Allocate (xday(Nmon)); Do it=1,Nmon; xday(it)=(it); Enddo; t_unit="months since 1961-01-01 " 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) !-------------------------------------------------------------------- ! Writing to Trend files: ! Real, Allocatable, Dimension (:,:,:) :: Sen_FDCha5 ! (Nlon,Nlat,Ncha) ! Real, Allocatable, Dimension (:,:,:) :: Prb_FDCha5 ! (Nlon,Nlat,Ncha) ! Real, Allocatable, Dimension (:,:) :: Sen_FDCha5_Reg ! (Nreg,Ncha) ! Real, Allocatable, Dimension (:,:) :: Prb_FDCha5_Reg ! (Nreg,Ncha) var_unit="per_yr_5yr"; lon_nm="lon"; lat_nm="lat"; t_nm="year" Call wr_ncdf_3d (Ftrend,Lon,Lat,xcha,Nlon,Nlat,Ncha,Sen_FDCha5,"sen",t_unit,trim(var_ln(1)),var_unit,lon_nm,lat_nm,t_nm,Miss) Call wr_ncdf_3d_apd (Ftrend,Prb_FDCha5,"prb",Nlon,Nlat,Ncha,trim(var_ln(1)),var_unit,lon_nm,lat_nm,t_nm,Miss) lon_nm="reg"; lat_nm="cha"; Call wr_ncdf_2d (Ftrend_reg,xreg,xcha,Nreg,Ncha,Sen_FDCha5_Reg,"sen",trim(var_ln(1)),var_unit,lon_nm,lat_nm,Miss) Call wr_ncdf_2d_apd (Ftrend_reg,Prb_FDCha5_Reg,"prb",Nreg,Ncha,trim(var_ln(1)),var_unit,lon_nm,lat_nm,Miss) End Subroutine Write_data_to_file !-------------------------------------------------------------------- Subroutine Correlation(X,Y,N,Corr,prob,Miss) Integer :: N, N1, i, j Real :: X(n), Y(n), Miss, Xtb, Ytb, Sx, Sy, Rxy, tmp, t N1=0; Xtb=0; Ytb=0; Sx=0; Sy=0; Rxy=0 Do i=1,N if (X(i)/=Miss.AND.Y(i)/=Miss) then n1=n1+1 Xtb=Xtb+X(i); Sx=Sx+X(i)*X(i) Ytb=Ytb+Y(i); Sy=Sy+Y(i)*Y(i) Rxy=Rxy+X(i)*Y(i) endif Enddo if (N1>2) then Xtb=Xtb/Real(N1); Ytb=Ytb/Real(N1) Sx=Sx/Real(N1)-Xtb*Xtb; Sy=Sy/Real(N1)-Ytb*Ytb tmp=Rxy/Real(N1) Rxy=(tmp-Xtb*Ytb)/Sqrt(Sx*Sy) t = ABS(Rxy/ ( Sqrt(1-Rxy**2) / Sqrt(Real(N1-2)) ) ) Prob = T_Dist (t, Real(N1-2)) else Rxy=Miss; prob=Miss endif Corr=Rxy End Subroutine Correlation !---------------------------------------------------------------------------------------------- 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" Include "/home/tanpv/Mylibs/Sub_Trends.f90" !---------------------------------------------------------------------------------------------- END !---------------------------------------------------------------------------------------------- !Include "/home/tanpv/Mylibs/Student_Prob.f90" !----------------------------------------------------------------------------------------------