PROGRAM t_logistic Use NetCDF USE Logistic_Regression IMPLICIT NONE Real, parameter :: Miss=-99.0 Integer, Parameter :: Nmon=12 Character (Len=255) :: Fname, Inp_Dir, Mask_File, Out_Dir Integer :: Ntime,M, i,j,k, ilead, ndf, ier, IO Real :: aa,bb,cc,dd !Real :: PC, POD, TS, H1 Real,Allocatable :: XX(:,:) Real (dp), Allocatable :: X(:,:), beta(:), se_beta(:) Integer, Allocatable :: LG(:),N(:),s1(:) Real (dp) :: tmp, chisq, devnce, tail_prob Logical :: Status !------------------------------------------------------------------------------ Character (Len=10) :: Dr_type Real, Allocatable :: Regs(:,:), Lon(:), Lat(:) Real, Allocatable, Dimension (:,:,:) :: Dr_obs, Dr_fcs ! (Nlead,Nst,Ntime) Real, Allocatable, Dimension (:,:) :: PC, POD, TS, HSS Real, Allocatable, Dimension (:,:,:) :: betas ! Nlead,Ntime,M+1 Integer, Allocatable, Dimension (:) :: idxLon,idxLat Integer :: Nlon,Nlat,Nlead, ilon,ilat,ist, Nt,it, icnt, icount(6), ic1(6), ii Integer :: Yr1,Yr2,Nmods,Ncase,ilen,idr,Nst Real :: Dr_Thres !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Inp_Dir, Mask_File, Out_Dir !, TVar, RVar Namelist /Data_Info/ Yr1, Yr2, Dr_Thres, Dr_type, Nmods, Ncase !------------------------------------------------------------------------------ Open (1, file="/home/tanpv/Progs/562_VNU/Calc_IDXs/Namelist_Dr_Chars_Discr.txt", Status="Old") ! Tinh voi APH 1961-2015 Read (1, File_Dir_Names) Read (1, Data_Info) Close(1) Ntime=(Yr2-Yr1+1)*Nmon ! Maximum timeseries Nlead=6 ! All leadtimes Call rd_ncdf_2d (Mask_File,Lon,Lat,Nlon,Nlat,Regs,"reg") ilen=1 ! Drought timescale: 1: 1-month, 3: 3-month, 6: 6-month idr =1 ! Drought level: 1: Moderate, 2: Severe, 3:Extreme write (Fname,"(A,I2.2,A)") trim(Out_Dir)//"Dr_4_Logistic_Moderate_Len_",ilen,".dat" Open (1,file=Fname,status="Old",form='unformatted') IO=0; it=0; icount(:)=0; ic1(:)=0 Do while (io==0) ! Repeat in timeseries (Ntime) !if (IO/=0) exit Read (1,IOSTAT=io) ilead if (IO/=0) exit Read (1) Nst, M it=it+1 ii=mod(ilead,6); if (ii==0) ii=6 icount(ii)=icount(ii)+1 !print*,icount if (allocated (idxLon)) Deallocate (idxLon); if (allocated (idxLat)) Deallocate (idxLat) if (allocated (LG)) Deallocate (LG); if (allocated (X)) Deallocate (X) if (allocated (N)) Deallocate (N); if (allocated (s1)) Deallocate (s1) if (allocated (beta)) Deallocate (beta) if (allocated (se_beta)) Deallocate (se_beta) if (allocated (XX)) Deallocate (XX) Allocate (idxLon(Nst),idxLat(Nst)) Allocate (LG(Nst), X(Nst,M), N(Nst), s1(Nst)) Allocate (beta(0:M),se_beta(0:M)) Allocate (XX(Nst,M)) if (it==1) then ! first timestep Allocate (betas(Nlead,Ntime,M+1)); betas=Miss Allocate (Dr_obs(Nlead,Nst,Ntime),Dr_fcs(Nlead,Nst,Ntime)) Allocate (PC(Nlead,Ntime),POD(Nlead,Ntime),TS(Nlead,Ntime),HSS(Nlead,Ntime)) Dr_obs=Miss; Dr_fcs=Miss; PC=Miss; POD=Miss; TS=Miss; HSS=Miss; endif N(:) = 1 Read (1) idxLon Read (1) idxLat Read (1) LG Read (1) XX if (IO/=0) exit ! Enddo !print*,icount; stop !print*,LG(1:Nt) !print*,XX(1:50,1) X(:,:)=XX(:,:) CALL logistic(Nst,X,M,LG,N, chisq, devnce, ndf, beta, se_beta, ier) IF (ier /= 0) THEN ic1(ii)=ic1(ii)+1 ! WRITE(*, *) 'Error number', ier ELSE ! WRITE(*, '(a, f9.3, a, f9.3, a, i4, a)') & ! ' Deviance = ', devnce, ' Chi-squared = ', chisq, & ! ' with ', ndf, ' deg. of freedom' ! tail_prob = 1.0_dp - chi_squared(ndf, chisq) ! IF (tail_prob < 0.01_dp) THEN ! WRITE(*, *) '*** Significantly bad fit at the 1% level ***' ! END IF ! WRITE(*, '(a, g12.4/)') ' Chi-squared tail prob. = ', tail_prob ! WRITE(*, *) ' Coefficient Std.error' ! DO i = 0, M ! WRITE(*, '(" ", I3, " ", 2g13.5)') i, beta(i), se_beta(i) ! END DO betas(ilead,it,1:M+1)=beta(0:M) ! Luu he so Beta aa=0; bb=0;cc=0;dd=0 Do ist=1,Nst tmp=beta(0) Do j=1,M tmp=tmp+beta(j)*x(ist,j) ! F = exp(b0 + b1.X1 + b2.X2 + ... + bk.Xk) Enddo tmp=EXP(tmp); tmp=tmp/(1+tmp); ! p = F/(1 + F) if (tmp>0.5) then; s1(ist)=1; else; s1(ist)=0; endif Dr_obs(ilead,ist,it)=LG(ist) Dr_fcs(ilead,ist,it)=tmp ! print"(2I4,F6.3)",LG(i),s1(i),tmp if (s1(ist)==1.AND.LG(ist)==1) aa=aa+1 if (s1(ist)==1.AND.LG(ist)==0) bb=bb+1 if (s1(ist)==0.AND.LG(ist)==1) cc=cc+1 if (s1(ist)==0.AND.LG(ist)==0) dd=dd+1 Enddo !print*,(aa+dd)/(aa+bb+cc+dd); !pc=aa+bb+cc+dd; pc=(aa+dd)/pc PC(ilead,it) =(aa+dd)/(aa+bb+cc+dd) POD(ilead,it)=aa/(aa+cc); TS(ilead,it) =aa/(aa+bb+cc); HSS(ilead,it)=2.0*(aa*dd-bb*cc)/((aa+cc)*(cc+dd)+(aa+bb)*(bb+dd)) print*,"PC, POD, TS, HSS=",PC(ilead,it), POD(ilead,it), TS(ilead,it), HSS(ilead,it) !print*,POD, TS END IF print*,"it=",it pause Enddo print*,icount(:) print*,ic1(:) write(Fname,"(A,I1,A)") trim(Out_Dir)//"Dr_fcs_Logistic_based_on_dridx_for_"//trim(Dr_type)//"_len_",ilen,".nc" !Call Write_Dr_and_Verification (Fname,Dr_type,ilen) !--------------------------------------------------------------------- CONTAINS !------------------------------------------------------------------------------ Subroutine Write_Dr_and_Verification (Fname,Dr_type,ilen) 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 :: ilead,ilen, it 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 write(Var_ln,"(A,I2.2,A)") "Drought events of the ",ilen,"-month observed "//trim(Dr_type)//" drought" Call wr_ncdf_4d (Fname, Lon,Lat,xlead,xtime,Nlon,Nlat,Nlead,Ntime,Dr_obs,"obs", & t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) write(Var_ln,"(A,I2.2,A)") "Drought events of the ",ilen,"-month forecast "//trim(Dr_type)//" drought" Call wr_ncdf_4d_apd (Fname,Dr_fcs,"fcs",Nlon,Nlat,Nlead,Ntime, & Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="Proportion Correction (PC): Time=1-6 corresponding to LT=1-6" Call wr_ncdf_3d_apd (Fname,PC,"pc",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Probability Of Detection (POD): Time=1-6 corresponding to LT=1-6" Call wr_ncdf_3d_apd (Fname,POD,"pod",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Threat Score (TS), or Critical Success Index (CSI): Time=1-6 corresponding to LT=1-6" Call wr_ncdf_3d_apd (Fname,TS,"ts",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) Var_ln="Heidke Skill Score (HSS): Time=1-6 corresponding to LT=1-6" Call wr_ncdf_3d_apd (Fname,HSS,"hss",Nlon,Nlat,Nlead,Var_ln,var_unit,lon_nm,lat_nm,t_nm, Miss) End Subroutine Write_Dr_and_Verification !------------------------------------------------------------------------------ Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" !------------------------------------------------------------------------------ END PROGRAM t_logistic