Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: amssng = -99., Miss=amssng !Real, parameter :: Dr_Thres = -1.0 Integer, Parameter :: NMonth = 12, Nmon=NMonth Integer, Parameter :: Nreg = 7, Nmem = 7, Nlead = 6 !Integer, Parameter :: Ncase = 2 ! obs,mod,cli,ann !Character (Len=3), Parameter :: Cases(Ncase)=(/"mod","cli","ann"/) !(/"obs","mod","cli","ann"/) !Character (Len=3), Parameter :: Cases(Ncase)=(/"ann"/) !Character (Len=5), Parameter :: Modeln(Nmods)=(/"RegCM","clWRF"/) !Character (Len=5) :: Modeln(Nmods) Real, Allocatable :: Regs(:,:), Lon(:), Lat(:), xlon(:), xlat(:) Real, Allocatable :: Dr_All_1(:,:,:,:,:,:,:) ! Nlon,Nlat,Ncase,Nmem,Ntime,Nmods,Nlead) Real, Allocatable :: Dr_All_3(:,:,:,:,:,:,:) Real, Allocatable :: Dr_All_6(:,:,:,:,:,:,:) Real, Allocatable, Dimension (:,:,:,:,:) :: IDX_1, IDX_3, IDX_6 ! Nlon,Nlat,Ntime,Nlead,Ncase*Nmem*Nmods) Real, Allocatable, Dimension (:,:,:,:) :: Dr_1, Dr_3, Dr_6 ! Nlon,Nlat,Ntime,Nlead Real, Allocatable :: IDX_Obs_1(:,:,:,:) ! Nlon,Nlat,Ntime,Nlead Real, Allocatable :: IDX_Obs_3(:,:,:,:) Real, Allocatable :: IDX_Obs_6(:,:,:,:) Real, Allocatable :: xmem(:), xlead(:), xtime(:) Real, Allocatable :: SPEI_All(:,:,:,:,:,:) ! Nlon,Nlat,Nmem,Ntime,Ncase,Nlen !------------------------------------------------------------------------------ ! Discr: Real, Allocatable, Dimension (:,:,:,:) :: Dr_Ev_1,Dr_Ev_3,Dr_Ev_6 ! (Nlon,Nlat,Nlead,Ntime) Real, Allocatable, Dimension (:,:,:,:) :: Fc_Ev_1,Fc_Ev_3,Fc_Ev_6 ! Forecast Real, Allocatable, Dimension (:,:,:) :: PC_1, PC_3, PC_6 ! Proportion Correction (Nlon,Nlat,Nlead) Real, Allocatable, Dimension (:,:,:) :: POD_1,POD_3,POD_6 ! Probability Of Detection, Hit Rate Real, Allocatable, Dimension (:,:,:) :: TS_1,TS_3,TS_6 ! Threat Score, or Critical Success Index (CSI) Real, Allocatable, Dimension (:,:,:) :: HSS_1,HSS_3,HSS_6 ! Heidke Skill Score !------------------------------------------------------------------------------ ! ANN: character (len=250) :: input_file, output_file, save_file character (len=250) :: Dir, Dir_ANN_Out, Dir_Cli, Dir_ANN_Save, Dir_Veri, Dir_4Dr logical :: training_mode ! Use for calculating scaling factors Integer :: Nsample, Nepoc, Nsample1 Real,allocatable, Dimension(:,:) :: Inputs, Outputs, Inputs_X, Outputs_Y Real,allocatable, Dimension(:) :: min_in, max_in, min_out, max_out Real :: mean_error, mean_error_check Real, Allocatable :: X_inp(:,:), Y_inp(:), Yest(:,:), Yhq(:,:) integer :: Nin, Nhidden, Nout, iin, ihid, iout integer :: i,j,k Real,allocatable, Dimension(:,:) :: weight_ih, weight_ho, moment_ih, moment_ho Real,allocatable, Dimension(:) :: fin, fhid, fout, ftarget, delta_o, delta_h ! Training Real :: learning_rate0, learning_rate, momentum, error, tmp, mintmp logical :: debug, Run_Cli !------------------------------------------------------------------------------ Character (Len=255) :: Fname, BL_Dir, R_Fname, T_Fname, Out_Dr_Fname Character (Len=255) :: Inp_Dir, Out_Dir, Mask_File Character (Len=3) :: TVar, RVar Character (Len=10) :: Model_nm, Dr_type Logical :: status, Empty Integer :: Yr1, Yr2, NYear, maxyrs, Nmods, Ncase Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, Nt, it, Landpoint Integer :: ilen, imod, ist, iyr, imw, ivar, isys, system, ilon, ilat Integer :: Nt1 Character (Len=3), Allocatable :: Cases(:) Character (Len=5), Allocatable :: Modeln(:) !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:,:,:,:) :: X4d Real, allocatable, Dimension (:) :: X1d,Xt !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Inp_Dir, Mask_File, Out_Dir !, TVar, RVar Namelist /Data_Info/ Yr1, Yr2, Dr_Thres, Dr_type, Nmods, Ncase Namelist /Model_Info/ Modeln, Cases ! Names of Model and Cases !------------------------------------------------------------------------------ Open (1, file="Namelist_Dr_Chars_Discr.txt", Status="Old") ! Tinh voi APH 1961-2015 Read (1, File_Dir_Names) Read (1, Data_Info) Allocate (Modeln(Nmods),Cases(Ncase)) Read (1, Model_Info) Close (1) Nyear = (Yr2 - Yr1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear Nspl = Ncase*Nmem*Nmods Call Read_Regs Dir_ANN_Save="/home/tanpv/Data/562_VNU/0.25_VN/ANN_Save_grd/" !------------------------------------------------------------------------------ !print*,"Beginning: Ntime, Nst=",Ntime, Nst !print*,"Parameters: Nlat,Nlon,Nst,Nlen,Ntime=",Nlat,Nlon,Nst,Nlen,Ntime Call Reinit Call Read_all_data print*,"Reading data done ..." Call Data_organization print*,"Call Data_organization Done" Call Obs_Dr_Events (Thres) print*,"Call Obs_Dr_Events Done" ilen=1; Call Dr_Discr(ilen) print*,"Dr_Discr(1) Done" ilen=3; Call Dr_Discr(ilen) ilen=6; Call Dr_Discr(ilen) !write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Dr_"//trim(Dr_type)//"_Discr.nc" !Call Write_Dr_IDXs_to_file (Fname,Dr_type) CONTAINS !------------------------------------------------------------------------------ Subroutine Reinit Allocate (Dr_All_1(Nlon,Nlat,Ncase,Nmem,Ntime,Nmods,Nlead)) ! Nlon,Nlat,Ncase,Nmem,Ntime,Nmods,Nlead) Allocate (Dr_All_3(Nlon,Nlat,Ncase,Nmem,Ntime,Nmods,Nlead)) Allocate (Dr_All_6(Nlon,Nlat,Ncase,Nmem,Ntime,Nmods,Nlead)) Allocate (IDX_Obs_1(Nlon,Nlat,Nlead,Ntime)) ! Nlon,Nlat,Ntime,Nlead Allocate (IDX_Obs_3(Nlon,Nlat,Nlead,Ntime)) Allocate (IDX_Obs_6(Nlon,Nlat,Nlead,Ntime)) Allocate (IDX_1(Nlon,Nlat,Ntime,Nlead,Nspl)) !, IDX_3, IDX_6 ! Nlon,Nlat,Ntime,Nlead,Ncase*Nmem*Nmods) Allocate (IDX_3(Nlon,Nlat,Ntime,Nlead,Nspl)) Allocate (IDX_6(Nlon,Nlat,Ntime,Nlead,Nspl)) Allocate (Dr_1(Nlon,Nlat,Nlead,Ntime)) ! , Dr_3, Dr_6 ! Nlon,Nlat,Ntime,Nlead --> moderate drought probability (IDX<-1) Allocate (Dr_3(Nlon,Nlat,Nlead,Ntime)) Allocate (Dr_6(Nlon,Nlat,Nlead,Ntime)) Allocate (Dr_Ev_1(Nlon,Nlat,Nlead,Ntime)) Allocate (Dr_Ev_3(Nlon,Nlat,Nlead,Ntime)) Allocate (Dr_Ev_6(Nlon,Nlat,Nlead,Ntime)) Allocate (Fc_Ev_1(Nlon,Nlat,Nlead,Ntime)) Allocate (Fc_Ev_3(Nlon,Nlat,Nlead,Ntime)) Allocate (Fc_Ev_6(Nlon,Nlat,Nlead,Ntime)) Allocate (PC_1(Nlon,Nlat,Nlead)) Allocate (PC_3(Nlon,Nlat,Nlead)) Allocate (PC_6(Nlon,Nlat,Nlead)) Allocate (POD_1(Nlon,Nlat,Nlead)) Allocate (POD_3(Nlon,Nlat,Nlead)) Allocate (POD_6(Nlon,Nlat,Nlead)) Allocate (TS_1(Nlon,Nlat,Nlead)) Allocate (TS_3(Nlon,Nlat,Nlead)) Allocate (TS_6(Nlon,Nlat,Nlead)) Allocate (HSS_1(Nlon,Nlat,Nlead)) Allocate (HSS_3(Nlon,Nlat,Nlead)) Allocate (HSS_6(Nlon,Nlat,Nlead)) ! For ANN: !Allocate (Y_inp(Ntime),X_inp(Ntime,Nspl)) !Allocate (Inputs_X(Nspl,Ntime), Outputs_Y(1,Ntime)) End Subroutine Reinit !------------------------------------------------------------------------------ Subroutine Write_Dr_IDXs_to_file (Fname,Dr_type) Character (Len=*) :: Fname,Dr_type Character (Len=250) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Integer :: it, imem, ilead, ireg, icas, ilen Real :: xlead(Nlead),xtime(Ntime) t_unit="months since 1983-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="lead"; t_nm="mon" Do ilead=1,Nlead; xlead(ilead)=ilead; Enddo Do it=1,Ntime; xtime(it)=it-1; Enddo Var_ln="The 1-month observed "//trim(Dr_type)//" drought events" Call wr_ncdf_4d (Fname, xlon,xlat,xlead,xtime,Nlon,Nlat,Nlead,Ntime,Dr_Ev_1,"obs1", & t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="The 3-month observed "//trim(Dr_type)//" drought events" Call wr_ncdf_4d_apd (Fname,Dr_Ev_3,"obs3",Nlon,Nlat,Nlead,Ntime, & Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="The 6-month observed "//trim(Dr_type)//" drought events" Call wr_ncdf_4d_apd (Fname,Dr_Ev_6,"obs6",Nlon,Nlat,Nlead,Ntime, & Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="The 1-month model-based discriminant analysis "//trim(Dr_type)//" drought events" Call wr_ncdf_4d_apd (Fname,Fc_Ev_1,"dr1",Nlon,Nlat,Nlead,Ntime, & Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="The 3-month model-based discriminant analysis "//trim(Dr_type)//" drought events" Call wr_ncdf_4d_apd (Fname,Fc_Ev_3,"dr3",Nlon,Nlat,Nlead,Ntime, & Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="The 6-month model-based discriminant analysis "//trim(Dr_type)//" drought events" Call wr_ncdf_4d_apd (Fname,Fc_Ev_6,"dr6",Nlon,Nlat,Nlead,Ntime, & Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Proportion Correction (PC) for the case of 1-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,PC_1,"pc1",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Proportion Correction (PC) for the case of 3-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,PC_3,"pc3",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Proportion Correction (PC) for the case of 6-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,PC_6,"pc6",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Probability Of Detection (POD) or Hit Rate (HR) for the case of 1-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,POD_1,"pod1",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Probability Of Detection (POD) or Hit Rate (HR) for the case of 3-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,POD_3,"pod3",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Probability Of Detection (POD) or Hit Rate (HR) for the case of 6-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,POD_6,"pod6",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Threat Score, or Critical Success Index (CSI) for the case of 1-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,TS_1,"ts1",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Threat Score, or Critical Success Index (CSI) for the case of 3-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,TS_3,"ts3",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Threat Score, or Critical Success Index (CSI) for the case of 6-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,TS_6,"ts6",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Heidke Skill Score for the case of 1-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,HSS_1,"hss1",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Heidke Skill Score for the case of 3-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,HSS_3,"hss3",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Heidke Skill Score for the case of 6-month "//trim(Dr_type)//" drought" Call wr_ncdf_3d_apd (Fname,HSS_6,"hss6",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) !Var_ln="Timeseries of Prob. validation of 1-month drought occurrence" !Call wr_ncdf_4d_apd (Fname,PC_1,"pc1",Nlon,Nlat,Nlead,Ntime, & ! Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) !Var_ln="Timeseries of Prob. validation of 3-month drought occurrence" !Call wr_ncdf_4d_apd (Fname,PC_3,"pc3",Nlon,Nlat,Nlead,Ntime, & ! Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) !Var_ln="Timeseries of Prob. validation of 6-month drought occurrence" !Call wr_ncdf_4d_apd (Fname,PC_6,"pc6",Nlon,Nlat,Nlead,Ntime, & ! Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) !Var_ln="Prob. validation of 1-month drought occurrence" !Call wr_ncdf_3d_apd (Fname,PC_All_1,"pca1",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) !Var_ln="Prob. validation of 3-month drought occurrence" !Call wr_ncdf_3d_apd (Fname,PC_All_3,"pca3",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) !Var_ln="Prob. validation of 6-month drought occurrence" !Call wr_ncdf_3d_apd (Fname,PC_All_6,"pca6",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) End Subroutine Write_Dr_IDXs_to_file !------------------------------------------------------------------------------ Subroutine Read_Regs Integer :: ilon,ilat, ireg, ic,j Real :: tmp Call rd_ncdf_2d (Mask_File,xlon,xlat,Nlon,Nlat,Regs,"reg") Nst=0; Do ilon=1,Nlon Do ilat=1,Nlat if (Int(Regs(ilon,ilat))>0) then Nst=Nst+1 endif Enddo Enddo !print*,"Nst=",Nst; pause End Subroutine Read_Regs !------------------------------------------------------------------------------ Subroutine Data_organization Integer :: ilon,ilat,ilead,imod,icas,it,imem,i,j ! Dimension (:,:,:,:,:) :: IDX_1, IDX_2, IDX_3 ! Nlon,Nlat,Ntime,Nlead,Ncase*Nmem*Nmods) IDX_1=Miss; IDX_3=Miss; IDX_6=Miss Do ilead=1,Nlead Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then i=0 Do imod=1,Nmods Do icas=1,Ncase Do imem=1,Nmem i=i+1 IDX_1(ilon,ilat,1:Ntime,ilead,i)=Dr_All_1(ilon,ilat,icas,imem,1:Ntime,imod,ilead) IDX_3(ilon,ilat,1:Ntime,ilead,i)=Dr_All_3(ilon,ilat,icas,imem,1:Ntime,imod,ilead) IDX_6(ilon,ilat,1:Ntime,ilead,i)=Dr_All_6(ilon,ilat,icas,imem,1:Ntime,imod,ilead) Enddo Enddo Enddo !print*,i endif Enddo Enddo Enddo End Subroutine Data_organization !------------------------------------------------------------------------------ Subroutine Dr_Discr_Reg (ireg) Integer :: ireg,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it0, it1, it2 Real, Allocatable :: X0(:,:),X1(:,:),xtmp(:) !Real :: X0(Ntime*Nst,Nspl), X1(Ntime*Nst,Nspl), Xtmp(Nspl) Real, Allocatable :: X(:), F(:,:) Integer, Allocatable :: LG(:), NG(:) Integer :: M,K !----------------------------------------------------------------------- ! He so ham Phan lop: (Nlon,Nlat,Nlead,Nlen,K,M+1) M=Nspl; K=2 Allocate (NG(K)) Allocate (X0(Ntime*Nst,Nspl),X1(Ntime*Nst,Nspl), Xtmp(Nspl)) Do ilead=1,Nlead it1=0; it0=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)==ireg) then !it1=0; it0=0 Do it=1,Ntime if (Dr_Ev_1(ilon,ilat,ilead,it)==1) then Xtmp(1:Nspl)=IDX_1(ilon,ilat,it,ilead,1:Nspl) if (No_missing_Ar1(Xtmp,Nspl,Miss)) then it1=it1+1; X1(it1,1:Nspl)=Xtmp(1:Nspl) endif else if (Dr_Ev_1(ilon,ilat,ilead,it)==0) then Xtmp(1:Nspl)=IDX_1(ilon,ilat,it,ilead,1:Nspl) if (No_missing_Ar1(Xtmp,Nspl,Miss)) then it0=it0+1; X0(it0,1:Nspl)=Xtmp(1:Nspl) endif endif Enddo endif Enddo Enddo NG(1)=it1; NG(2)=it0 Nt=0; Do ig=1,K; Nt=Nt+NG(ig); Enddo ! Nt=Tong so mau (<=Ntime) Allocate (X(Nt*M),LG(Nt*M),F(K,M+1)) it=0 Do j=1,M ! M=Nspl Do i=1,NG(1) it=it+1; X(it)=X1(i,j) Enddo Enddo Do j=1,M ! M=Nspl Do i=1,NG(2) it=it+1; X(it)=X0(i,j) Enddo Enddo Call My_Discr (X,M,K,NG,Nt,F,LG) i1=0; b=0.; c=0. Do i=1,K a=0. Do j=1,NG(i) i1=i1+1; if (LG(i1)==i) a=a+1 !print*,i1,j,LG(i1) Enddo b=b+a; c=c+NG(i) print*, Int(a), NG(i), a/NG(i)*100. Enddo print*,Int(b),Int(c), b/c*100. Deallocate (X,LG,F) pause Enddo End Subroutine Dr_Discr_Reg !------------------------------------------------------------------------------ Subroutine Dr_Discr (ilen) Integer :: ilen,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it0, it1, it2 Real, Allocatable :: Obs_Evt(:,:,:,:),DrIDX(:,:,:,:,:),Fcs(:,:,:,:) Real, Allocatable :: X0(:,:),X1(:,:),XX(:,:),xtmp(:) !Real :: X0(Ntime,Nspl), X1(Ntime,Nspl), Xtmp(Nspl) Real, Allocatable :: X(:), F(:,:), Fx(:), X1d(:), Xdat(:,:) Integer, Allocatable :: LG(:), NG(:), LG1(:) Integer :: M,K, idr Real, Allocatable, Dimension (:,:,:) :: PC, POD, TS, HSS ! Verification scores !----------------------------------------------------------------------- ! He so ham Phan lop: (Nlon,Nlat,Nlead,Nlen,K,M+1) , Fc_Ev_1 if(Allocated (Obs_Evt)) Deallocate (Obs_Evt); if(Allocated (DrIDX)) Deallocate (DrIDX) if(Allocated (Fcs)) Deallocate (Fcs); if(Allocated (X0)) Deallocate (X0) if(Allocated (Xtmp)) Deallocate (Xtmp); if(Allocated (X1)) Deallocate (X1) if(Allocated (PC)) Deallocate (PC); if(Allocated (POD)) Deallocate (POD) if(Allocated (TS)) Deallocate (TS); if(Allocated (HSS)) Deallocate (HSS) if(Allocated (NG)) Deallocate (NG); if(Allocated (Fx)) Deallocate (Fx) Allocate (Obs_Evt(Nlon,Nlat,Nlead,Ntime),DrIDX(Nlon,Nlat,Ntime,Nlead,Nspl)) Allocate (Fcs(Nlon,Nlat,Nlead,Ntime)); Fcs=Miss Allocate (PC (Nlon,Nlat,Nlead)); PC=Miss ! (PC_1(Nlon,Nlat,Nlead)) Allocate (POD(Nlon,Nlat,Nlead)); POD=Miss Allocate (TS (Nlon,Nlat,Nlead)); TS=Miss Allocate (HSS(Nlon,Nlat,Nlead)); HSS=Miss ! Get Data: if (ilen==1) then Obs_Evt=Dr_Ev_1; DrIDX=IDX_1 else if (ilen==3) then Obs_Evt=Dr_Ev_3; DrIDX=IDX_3 else if (ilen==6) then Obs_Evt=Dr_Ev_6; DrIDX=IDX_6 endif M=Nspl; K=2 Allocate (NG(K), X0(Ntime,Nspl), X1(Ntime,Nspl), Xtmp(Nspl),Fx(k)) Do ilead=1,Nlead Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then it1=0; it0=0 Do it=1,Ntime if (Obs_Evt(ilon,ilat,ilead,it)==0) then ! No drought Xtmp(1:Nspl)=DrIDX(ilon,ilat,it,ilead,1:Nspl) if (No_missing_Ar1(Xtmp,Nspl,Miss)) then it0=it0+1; X0(it1,1:Nspl)=Xtmp(1:Nspl) ! X0, it0 endif else if (Obs_Evt(ilon,ilat,ilead,it)==1) then ! Drought Xtmp(1:Nspl)=DrIDX(ilon,ilat,it,ilead,1:Nspl) if (No_missing_Ar1(Xtmp,Nspl,Miss)) then it1=it1+1; X1(it0,1:Nspl)=Xtmp(1:Nspl) ! X1, it1 endif endif Enddo NG(1)=it0; NG(2)=it1 Nt=0; Do ig=1,K; Nt=Nt+NG(ig); Enddo ! Nt=Tong so mau (<=Ntime) if (NG(1)<=0.OR.NG(2)<=0) then print*,"Cannot do for point ",ilon,ilat," LT= ",ilead print*,"NG(1),NG(2)=",NG(1),NG(2) exit endif !Allocate (X(Nt*M),LG(Nt*K),F(K,M+1)) !,LG1(Nt*K),Xdat(Nt,K)) if (Allocated (XX)) Deallocate (XX) if (Allocated (LG)) Deallocate (LG) if (Allocated (F)) Deallocate (F) Allocate (XX(Nt,M),LG(Nt),F(K,M+1)) !it=0 !Do j=1,M ! M=Nspl ! Do i=1,NG(1) ! No drought ! it=it+1; X(it)=X0(i,j) ! Enddo !Enddo !Do j=1,M ! Do i=1,NG(2) ! Drought ! it=it+1; X(it)=X1(i,j) ! Enddo !Enddo !Call My_Discr (X,M,K,NG,Nt,F,LG) it=0 Do i=1,NG(1) it=it+1; XX(it,1:M)=X0(i,1:M) Enddo Do i=1,NG(2) it=it+1; XX(it,1:M)=X1(i,1:M) Enddo Call Discr_2023(K,Nt,M,NG,XX,F,LG) print*,"M,Nt=",M,Nt,NG(1:2) print"(30F8.4)",F(1,1:M+1) ! Re-Forecast using Discriminant functions F: ! a: O=1,F=1; b: O=0,F=1; c: O=1,F=0; d: O=0,F=0 a=0; b=0; c=0; d=0 Do it=1,Ntime Xtmp(1:M)=DrIDX(ilon,ilat,it,ilead,1:M) Do j=1,K tmp=0 Do i1=2,M+1 tmp=tmp+F(j,i1)*Xtmp(i1-1) Enddo Fx(j)=tmp+F(j,1) ! Value of j-th discriminant function Enddo idr=Jmax(Fx,K,Miss); !if (idr==2) idr=0 Fcs(ilon,ilat,ilead,it)=idr-1 print*,it,int(Fcs(ilon,ilat,ilead,it)) if (Obs_Evt(ilon,ilat,ilead,it)==1.AND.Fcs(ilon,ilat,ilead,it)==1) a=a+1 if (Obs_Evt(ilon,ilat,ilead,it)==0.AND.Fcs(ilon,ilat,ilead,it)==1) b=b+1 if (Obs_Evt(ilon,ilat,ilead,it)==1.AND.Fcs(ilon,ilat,ilead,it)==0) c=c+1 if (Obs_Evt(ilon,ilat,ilead,it)==0.AND.Fcs(ilon,ilat,ilead,it)==0) d=d+1 !print*,Int(Obs_Evt(ilon,ilat,ilead,it)),Int(Fcs(ilon,ilat,ilead,it)) Enddo PC (ilon,ilat,ilead)=(a+d)/(a+b+c+d) POD(ilon,ilat,ilead)=a/(a+c) TS (ilon,ilat,ilead)=a/(a+b+c) HSS(ilon,ilat,ilead)=2*(a*d-b*c)/((a+c)*(c+d)+(a+b)*(b+d)) !print*,"FBI=",(a+b)/(a+c) !print*,"PC =",(a+d)/(a+b+c+d) !print*,"POD=",a/(a+c) !print*,"FAR=",b/(a+b) !print*,"POFD=",b/(b+d)," Probability of False Detection" !print*,"TS =",a/(a+b+c), "Threat Score" !ar=(a+b)*(a+c)/(a+b+c+d) !print*,"ETS=",(a-ar)/(a+b+c-ar) !print*,"KSS=",(a*d-b*c)/((a+c)*(b+d)) !print*,"HSS=",2*(a*d-b*c)/((a+c)*(c+d)+(a+b)*(b+d)) i1=0; b=0.; c=0. Do i=1,K a=0.; L=0 Do j=1,NG(i) i1=i1+1; if (LG(i1)==i) a=a+1 !print*,i1,j,LG(i1) Enddo b=b+a; c=c+NG(i) ! print*, Int(a), NG(i), a/NG(i)*100. Enddo print*,Int(b),Int(c), b/c*100. !Deallocate (XX,LG,F) !,LG1,Xdat) !pause endif Enddo Enddo Enddo print*,"Storing results..." if (ilen==1) then Fc_Ev_1=Fcs; PC_1=PC; POD_1=POD; TS_1=TS; HSS_1=HSS else if (ilen==3) then Fc_Ev_3=Fcs; PC_3=PC; POD_3=POD; TS_3=TS; HSS_3=HSS else if (ilen==6) then Fc_Ev_6=Fcs; PC_6=PC; POD_6=POD; TS_6=TS; HSS_6=HSS endif ! Allocate (NG(K), X0(Ntime,Nspl), X1(Ntime,Nspl), Xtmp(Nspl),X1d(M),Fx(k)) ! Allocate (Obs_Evt(Nlon,Nlat,Nlead,Ntime),DrIDX(Nlon,Nlat,Ntime,Nlead,Nspl)) ! Allocate (Fcs(Nlon,Nlat,Nlead,Ntime)); Fcs=Miss ! print*,"Storing results...Done" ! Deallocate (X0,X1,Xtmp,NG) ! Deallocate (Obs_Evt,DrIDX,Fcs,Fx) ! Deallocate (PC,POD,TS,HSS) print*,"Dr_Discr ",ilen," Done" return End Subroutine Dr_Discr !------------------------------------------------------------------------------ Integer Function Jmax (X,N,Miss) ! So thu tu phan tu Max cua X(1..N) Integer :: N,i,j Real :: X(N),tmp,Miss tmp=-999999. Do i=1,N if (X(i)/=Miss.AND.X(i)>tmp) then tmp=X(i); j=i endif Enddo Jmax = j End Function Jmax !------------------------------------------------------------------------------ Subroutine Obs_Dr_Events (Thres) ! Determine observed drought events for each LT, timeslice and timescale based on specific drought threshold Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2 Real :: Thres, a,b,c Dr_Ev_1=Miss; Dr_Ev_3=Miss; Dr_Ev_6=Miss ! (Nlon,Nlat,Nlead,Ntime) Do ilead=1,Nlead Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do it=1,Ntime a=IDX_Obs_1(ilon,ilat,ilead,it) b=IDX_Obs_3(ilon,ilat,ilead,it) c=IDX_Obs_6(ilon,ilat,ilead,it) if (a/=Miss.AND.a<=Thres) then Dr_Ev_1(ilon,ilat,ilead,it)=1 else Dr_Ev_1(ilon,ilat,ilead,it)=0 endif if (b/=Miss.AND.b<=Thres) then Dr_Ev_3(ilon,ilat,ilead,it)=1 else Dr_Ev_3(ilon,ilat,ilead,it)=0 endif if (c/=Miss.AND.c<=Thres) then Dr_Ev_6(ilon,ilat,ilead,it)=1 else Dr_Ev_6(ilon,ilat,ilead,it)=0 endif Enddo endif Enddo Enddo Enddo End Subroutine Obs_Dr_Events !------------------------------------------------------------------------------ Subroutine Read_all_data Integer :: ist,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2 Character (Len=255) :: Fname Real, Allocatable :: X4d(:,:,:,:), time(:), xlon(:),xlat(:),xmem(:) !,xtime(:) ! Dr_Obs_1(:,:,:,:) ! Nlon,Nlat,Ntime,Nlead ! Dr_All_1(:,:,:,:,:,:,:) ! Nlon,Nlat,Ncase,Nmem,Ntime,Nmods,Nlead) Do imod=1,Nmods Do ilead=1,Nlead write (Fname,"(A,I1,A)") trim(Inp_Dir)//"SPEI_"//trim(Modeln(imod))//"_lead_",ilead,".nc" print*,"Reading ",trim(Fname) if (imod==1) then ! .AND.ilead==1) then Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,"obs_01") Do it=1,Nt; IDX_Obs_1(:,:,ilead,it)=X4d(:,:,1,it); Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,"obs_03") Do it=1,Nt; IDX_Obs_3(:,:,ilead,it)=X4d(:,:,1,it); Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,"obs_06") Do it=1,Nt; IDX_Obs_6(:,:,ilead,it)=X4d(:,:,1,it); Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) endif Do icas=1,Ncase Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas))//"_01") Do imem=1,Nmem Do it=1,Nt; Dr_All_1(:,:,icas,imem,it,imod,ilead)=X4d(:,:,imem,it); Enddo Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas))//"_03") Do imem=1,Nmem Do it=1,Nt; Dr_All_3(:,:,icas,imem,it,imod,ilead)=X4d(:,:,imem,it); Enddo Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) Call rd_ncdf_4d (Fname,xlon,xlat,xmem,xtime,Nlo,Nla,Nme,Nt,X4d,trim(Cases(icas))//"_06") Do imem=1,Nmem Do it=1,Nt; Dr_All_6(:,:,icas,imem,it,imod,ilead)=X4d(:,:,imem,it); Enddo Enddo Deallocate (xlon,xlat,xmem,xtime,X4d) Enddo Enddo Enddo Print*,"Reading data ... Done" End Subroutine Read_all_data !------------------------------------------------------------------------------ Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Sub_neural_2023.f90" Include "/home/tanpv/Mylibs/No_missing_Ar1.f90" Include "/home/tanpv/Mylibs/Sub_My_Discr_2023.f90" ! Call Discr_2023(NGROUP,N,NVAR,NG,X,COEF,LG) !Include "/home/tanpv/Mylibs/Sub_SPI.f90" ! Chuong trinh tinh chi so han SPI !Include "/home/tanpv/Mylibs/Sub_SPEI.f90" ! Chuong trinh tinh chi so han SPEI !-------------------------------------------------------------------- END ! Program !Include "Sub_Discr.f90"