Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: Miss=-99.0 Integer, Parameter :: Nmon=12, Nyrs=40 !Integer, Parameter :: Max_evt=1000 ! Maximum number of events over all period Yr1-Yr2 Integer, Parameter :: Nidx=2 ! "spi","spei" Integer, Parameter :: Nlen=3, Dr_lev=3 Integer, Parameter :: Len_SPI(Nlen) =(/3,6,12/) Real, parameter :: Dr_thre(Dr_lev) =(/-0.5,-1.0,-1.5/) ! Dry (or Mild), Moderate and Severe Character (Len=3), Parameter :: Dr_lev_nm(Dr_lev)=(/"mil","mod","sev"/) Character (Len=5) :: idx_nm(Nidx)=(/"spi","spei"/) Real, Allocatable :: Lon(:), Lat(:), Reg(:), xlon(:), xlat(:), time(:), xmod(:), Mask(:,:) Character (Len=20), Allocatable :: stname(:) Integer, Parameter :: ENSO_yrs(Nyrs)=(/0,0,1,-1,-1,0,1,1,-1,0,0,1,0,0,1,-1,0,1,-1,-1,-1,0,1,0,1,-1,1,-1,-1,1,-1,-1,0,0,1,1,-1,-1,1,0/) Integer :: Nst, Nlon, Nlat, it, Max_evt !------------------------------------------------------------------------------ ! For HW SEA: Real, Allocatable, Dimension(:,:,:) :: Frq_obs, Dm_obs, Dx_obs ! Nst,Nlen,Nyear Real, Allocatable, Dimension(:,:,:) :: Sev_obs, HWday_obs Real, Allocatable, Dimension(:,:,:) :: Svx_obs ! For over all period: Real, Allocatable, Dimension(:,:) :: Frq_obs_all, Dm_obs_all, Dx_obs_all ! Nst,Nlen Real, Allocatable, Dimension(:,:) :: Sev_obs_all, HWday_obs_all Real, Allocatable, Dimension(:,:) :: Svx_obs_all ! Trend of Chars: Real, Allocatable, Dimension(:,:) :: Tr_Frq, Tr_Dm, Tr_Dx, Tr_Sev, Tr_Svx, Tr_Nday Real, Allocatable, Dimension(:,:) :: pv_Frq, pv_Dm, pv_Dx, pv_Sev, pv_Svx, pv_Nday ! ENSO: Real, Allocatable, Dimension(:,:) :: EL_Frq, EL_Dm, EL_Dx, EL_Sev, EL_Svx, EL_Nday Real, Allocatable, Dimension(:,:) :: LA_Frq, LA_Dm, LA_Dx, LA_Sev, LA_Svx, LA_Nday Real, Allocatable, Dimension(:,:) :: NO_Frq, NO_Dm, NO_Dx, NO_Sev, NO_Svx, NO_Nday !------------------------------------------------------------------------------ Real, Allocatable, Dimension(:,:,:) :: HWevt_obs ! Nst,Ntime,Nlen Real, Allocatable, Dimension(:,:,:) :: HWsev_obs ! Nst,Ntime,Nlen !------------------------------------------------------------------------------ Real, Allocatable, Dimension(:,:,:) :: SPEI_obs ! Nst,Ntime,Nlen - SPEI with a fixed timescale in days (e.g. 31 days, 91 days) !------------------------------------------------------------------------------ Character (Len=255) :: Fname,Fname1,Fname2,Fname_evt,Fname_sev,Fname_hw,Fname_qp Character (Len=255) :: Obs_Dir, Inp_Dir, Out_Dir, Mask_File Character (Len=255) :: Fevt_obs,Fevt_mod,Fsev_obs,Fsev_mod Integer :: Yr1, Yr2, Nmodel, Ntime, Nyear, Mwd, Timescale Real :: p Character (Len=6) :: Tx_nm, Tm_nm, Pr_nm Logical :: status, Empty, calc_qp Real :: Drou_Thres !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:,:,:,:) :: X4d Real, allocatable, Dimension (:) :: X1d, Xt, SPEI_out, Rain, T2m, Tx, Tm !------------------------------------------------------------------------------ Integer :: Nlo,Nla,Nt,Nz,imod Integer :: i,j,k,iidx,ilen !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Inp_Dir,Out_Dir Namelist /Data_Info/ Yr1, Yr2 !------------------------------------------------------------------------------ Open (1, file="Namelist_Chars_CDH.txt",Status="old") Read (1, File_Dir_Names) Read (1, Data_Info) Close(1) Ntime = (Yr2-Yr1)*Nmon ! Bo nam dau tien Nyear = Yr2-Yr1 Max_evt = Ntime ! Ntime=1981-2020=40=Nyrs !------------------------------------------------------------------------ Do iidx=1,Nidx Fname = trim(Inp_Dir)//"spi_spei_stn_1980-2020.nc" Call Read_Drough_IDX (trim(Fname),idx_nm(iidx)) Do ilev=1,Dr_lev ! Drought levels Call Calc_Dr_Evnts_obs (Dr_thre(ilev)) Call Calc_Chars_obs (Dr_thre(ilev)) Call Calc_Trend Call Calc_mean_for_ENSO_phases Fname=trim(Out_Dir)//"Chars_Dr_yearly_"//trim(idx_nm(iidx))//"_"//trim(Dr_lev_nm(ilev))//".nc" Call Write_Chars_yearly_to_NC (Fname) Fname=trim(Out_Dir)//"Chars_Dr_Clim_"//trim(idx_nm(iidx))//"_"//trim(Dr_lev_nm(ilev)) Call Write_Chars_Clim_to_GrADS_and_NC (Fname) Fname=trim(Out_Dir)//"Tr_of_Char_"//trim(idx_nm(iidx))//"_"//trim(Dr_lev_nm(ilev)) Call Write_Trend_to_GrADS_and_NC (Fname) Fname=trim(Out_Dir)//"ENSO_Char_"//trim(idx_nm(iidx))//"_"//trim(Dr_lev_nm(ilev)) Call Write_Mean_ENSO_to_GrADS (Fname) Enddo Enddo print*,"Normal terminal of Chars_Dr_Ha_Obs.f90" !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ Subroutine Write_Chars_yearly_to_NC (Fname) Character (Len=*) :: Fname Real, Allocatable :: stnid(:) Real :: xtime(Nyear),xlen(Nlen) Character (Len=255) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm t_unit="years since 1981-01-01 "; var_unit="-"; lon_nm="stn"; lat_nm="timescale"; t_nm="time" Allocate (stnid(Nst)) Do i=1,Nst; stnid(i)=i; Enddo Do i=1,Nyear; xtime(i)=i; Enddo xlen(:) = Len_SPI(:) var_ln="Yearly Drought frequency"; var_unit="spell/year"; Call wr_ncdf_3d (Fname,stnid,xlen,xtime,Nst,Nlen,Nyear,Frq_obs,"frq",t_unit, var_ln,"-", lon_nm, lat_nm, t_nm, Miss) var_ln="Yearly drought mean duration (month)"; var_unit="month/spell"; Call wr_ncdf_3d_apd (Fname,Dm_obs,"dm",Nst,Nlen,Nyear, var_ln, var_unit, lon_nm, lat_nm, t_nm,Miss) var_ln="Yearly drought maximum duration (month)"; var_unit="month/spell"; Call wr_ncdf_3d_apd (Fname,Dx_obs,"dx",Nst,Nlen,Nyear, var_ln, var_unit, lon_nm, lat_nm, t_nm,Miss) var_ln="Yearly drought mean severity"; var_unit="-" Call wr_ncdf_3d_apd (Fname,Sev_obs,"sev",Nst,Nlen,Nyear, var_ln, var_unit, lon_nm, lat_nm, t_nm,Miss) var_ln="Yearly drought maximum severity"; var_unit="-" Call wr_ncdf_3d_apd (Fname,Svx_obs,"svx",Nst,Nlen,Nyear, var_ln, var_unit, lon_nm, lat_nm, t_nm,Miss) var_ln="Yearly drought month (number of months)"; var_unit="month/year" Call wr_ncdf_3d_apd (Fname,HWday_obs,"nday",Nst,Nlen,Nyear, var_ln, var_unit, lon_nm, lat_nm, t_nm,Miss) Call wr_ncdf_1d_apd (Fname,xlon,"lon",Nst,"Longitude","Deg E", lon_nm, Miss) Call wr_ncdf_1d_apd (Fname,xlat,"lat",Nst,"Latitude","Deg N", lon_nm, Miss) Deallocate (stnid) End Subroutine Write_Chars_yearly_to_NC !------------------------------------------------------------------------------ Subroutine Write_Chars_Clim_to_GrADS_and_NC (Fname) Character (Len=*) :: Fname Integer, Parameter :: Nvar=6 Character (Len=4) :: vname(Nvar)=(/"frq","dm","dx","sev","svx","nday"/) Real, Allocatable :: X(:,:,:), stnid(:) Integer :: i,j,k Character (Len=255) :: Fname1,t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Real :: varid(Nvar),xlen(Nlen) Allocate (X(Nst,Nvar,Nlen),stnid(Nst)) X(:,1,:) = Frq_obs_all(:,:) X(:,2,:) = Dm_obs_all (:,:) X(:,3,:) = Dx_obs_all (:,:) X(:,4,:) = Sev_obs_all(:,:) X(:,5,:) = Svx_obs_all(:,:) X(:,6,:) = HWday_obs_all(:,:) Call Write_GrADS_Stn_3D (vname,X,Nst,Nvar,Nlen,xlon,xlat,Stname,Fname) Fname1=trim(Fname)//".nc" t_unit="years since 1983-01-01 "; var_unit="-"; lon_nm="stn"; lat_nm="dr_char"; t_nm="timescale" Do i=1,Nst; stnid(i)=i; Enddo Do i=1,Nvar; varid(i)=i; Enddo xlen(:) = Len_SPI(:) var_ln="Climatology of Drought characteristics for X=1-6: frq,dm,dx,sev,svx,nday" Call wr_ncdf_3d (Fname1,stnid,varid,xlen,Nst,Nvar,Nlen,X,"cha",t_unit, var_ln,"-", lon_nm, lat_nm, t_nm, Miss) Call wr_ncdf_1d_apd (Fname1,xlon,"lon",Nst,"Longitude","Deg E", lon_nm, Miss) Call wr_ncdf_1d_apd (Fname1,xlat,"lat",Nst,"Latitude","Deg N", lon_nm, Miss) Deallocate (X,stnid) End Subroutine Write_Chars_Clim_to_GrADS_and_NC !------------------------------------------------------------------------------ Subroutine Write_Trend_to_GrADS_and_NC (Fname) Character (Len=*) :: Fname Integer, Parameter :: Nvar=6 Character (Len=4) :: vname(Nvar)=(/"frq","dm","dx","sev","svx","nday"/) Real, Allocatable :: X(:,:,:),Y(:,:,:), stnid(:) Integer :: i,j,k Character (Len=255) :: Fname1,t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Real :: varid(Nvar),xlen(Nlen) Allocate (X(Nst,Nvar,Nlen),Y(Nst,Nvar,Nlen),stnid(Nst)) X(:,1,:) = Tr_Frq(:,:); Y(:,1,:) = pv_Frq(:,:); X(:,2,:) = Tr_Dm (:,:); Y(:,2,:) = pv_Dm (:,:); X(:,3,:) = Tr_Dx (:,:); Y(:,3,:) = pv_Dx (:,:); X(:,4,:) = Tr_Sev(:,:); Y(:,4,:) = pv_Sev(:,:); X(:,5,:) = Tr_Svx(:,:); Y(:,5,:) = pv_Svx(:,:); X(:,6,:) = Tr_Nday(:,:); Y(:,6,:) = pv_Nday(:,:); Call Write_GrADS_Stn_3D (vname,X,Nst,Nvar,Nlen,xlon,xlat,Stname,Fname) Fname1=trim(Fname)//"_pv" Call Write_GrADS_Stn_3D (vname,Y,Nst,Nvar,Nlen,xlon,xlat,Stname,Fname1) Fname1=trim(Fname)//".nc" t_unit="years since 1983-01-01 "; var_unit="-"; lon_nm="stn"; lat_nm="dr_char"; t_nm="timescale" Do i=1,Nst; stnid(i)=i; Enddo Do i=1,Nvar; varid(i)=i; Enddo xlen(:) = Len_SPI(:) var_ln="Changing trend of Drought characteristics for X=1-6: frq,dm,dx,sev,svx,nday" Call wr_ncdf_3d (Fname1,stnid,varid,xlen,Nst,Nvar,Nlen,X,"trend",t_unit, var_ln,"-", lon_nm, lat_nm, t_nm, Miss) var_ln="p-values of Changing trend of Drought characteristics for X=1-6: frq,dm,dx,sev,svx,nday" Call wr_ncdf_3d_apd (Fname1,Y,"pv",Nst,Nvar,Nlen,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) Call wr_ncdf_1d_apd (Fname1,xlon,"lon",Nst,"Longitude","Deg E", lon_nm, Miss) Call wr_ncdf_1d_apd (Fname1,xlat,"lat",Nst,"Latitude","Deg N", lon_nm, Miss) Deallocate (X,Y,stnid) End Subroutine Write_Trend_to_GrADS_and_NC !------------------------------------------------------------------------------ Subroutine Write_Mean_ENSO_to_GrADS (Fname) Character (Len=*) :: Fname Integer, Parameter :: Nvar=6 Character (Len=4) :: vname(Nvar)=(/"frq","dm","dx","sev","svx","nday"/) Real, Allocatable :: X(:,:,:),Y(:,:,:),Z(:,:,:), stnid(:) Integer :: i,j,k Character (Len=255) :: Fname1,t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Real :: varid(Nvar),xlen(Nlen) Allocate (X(Nst,Nvar,Nlen),Y(Nst,Nvar,Nlen),Z(Nst,Nvar,Nlen),stnid(Nst)) X(:,1,:) = EL_Frq(:,:); Y(:,1,:) = LA_Frq(:,:); Z(:,1,:) = NO_Frq(:,:); X(:,2,:) = EL_Dm (:,:); Y(:,2,:) = LA_Dm (:,:); Z(:,2,:) = NO_Dm (:,:); X(:,3,:) = EL_Dx (:,:); Y(:,3,:) = LA_Dx (:,:); Z(:,3,:) = NO_Dx (:,:); X(:,4,:) = EL_Sev(:,:); Y(:,4,:) = LA_Sev(:,:); Z(:,4,:) = NO_Sev(:,:); X(:,5,:) = EL_Svx(:,:); Y(:,5,:) = LA_Svx(:,:); Z(:,5,:) = NO_Svx(:,:); X(:,6,:) = EL_Nday(:,:); Y(:,6,:) = LA_Nday(:,:);Z(:,6,:) = NO_Nday(:,:) Fname1=trim(Fname)//"_EL" Call Write_GrADS_Stn_3D (vname,X,Nst,Nvar,Nlen,xlon,xlat,Stname,Fname1) Fname1=trim(Fname)//"_LA" Call Write_GrADS_Stn_3D (vname,Y,Nst,Nvar,Nlen,xlon,xlat,Stname,Fname1) Fname1=trim(Fname)//"_NO" Call Write_GrADS_Stn_3D (vname,Z,Nst,Nvar,Nlen,xlon,xlat,Stname,Fname1) End Subroutine Write_Mean_ENSO_to_GrADS !------------------------------------------------------------------------------ Subroutine Calc_mean_for_ENSO_phases Integer :: ist,ilen,it,iel,ila,ino !Real :: El(Nyear),La(Nyear),No(Nyear) if (Allocated (EL_Frq)) Deallocate (EL_Frq); if (Allocated (EL_Dm)) Deallocate (EL_Dm); if (Allocated (EL_Dx)) Deallocate (EL_Dx); if (Allocated (EL_Sev)) Deallocate (EL_Sev); if (Allocated (EL_Svx)) Deallocate (EL_Svx); if (Allocated (EL_Nday)) Deallocate (EL_Nday); if (Allocated (LA_Frq)) Deallocate (LA_Frq); if (Allocated (LA_Dm)) Deallocate (LA_Dm); if (Allocated (LA_Dx)) Deallocate (LA_Dx); if (Allocated (LA_Sev)) Deallocate (LA_Sev); if (Allocated (LA_Svx)) Deallocate (LA_Svx); if (Allocated (LA_Nday)) Deallocate (LA_Nday); if (Allocated (NO_Frq)) Deallocate (NO_Frq); if (Allocated (NO_Dm)) Deallocate (NO_Dm); if (Allocated (NO_Dx)) Deallocate (NO_Dx); if (Allocated (NO_Sev)) Deallocate (NO_Sev); if (Allocated (NO_Svx)) Deallocate (NO_Svx); if (Allocated (NO_Nday)) Deallocate (NO_Nday); Allocate (EL_Frq(Nst,Nlen), EL_Dm(Nst,Nlen), EL_Dx(Nst,Nlen), EL_Sev(Nst,Nlen), EL_Svx(Nst,Nlen), EL_Nday(Nst,Nlen)) Allocate (LA_Frq(Nst,Nlen), LA_Dm(Nst,Nlen), LA_Dx(Nst,Nlen), LA_Sev(Nst,Nlen), LA_Svx(Nst,Nlen), LA_Nday(Nst,Nlen)) Allocate (NO_Frq(Nst,Nlen), NO_Dm(Nst,Nlen), NO_Dx(Nst,Nlen), NO_Sev(Nst,Nlen), NO_Svx(Nst,Nlen), NO_Nday(Nst,Nlen)) EL_Frq=Miss; EL_Dm=Miss; EL_Dx=Miss; EL_Sev=Miss; EL_Svx=Miss; EL_Nday=Miss; LA_Frq=Miss; LA_Dm=Miss; LA_Dx=Miss; LA_Sev=Miss; LA_Svx=Miss; LA_Nday=Miss; NO_Frq=Miss; NO_Dm=Miss; NO_Dx=Miss; NO_Sev=Miss; NO_Svx=Miss; NO_Nday=Miss; Call Mean_ENSO (Frq_obs,EL_Frq,LA_Frq,NO_Frq) Call Mean_ENSO (Dm_obs, EL_Dm, LA_Dm, NO_Dm) Call Mean_ENSO (Dx_obs, EL_Dx, LA_Dx, NO_Dx) Call Mean_ENSO (Sev_obs,EL_Sev,LA_Sev,NO_Sev) Call Mean_ENSO (Svx_obs,EL_Svx,LA_Svx,NO_Svx) Call Mean_ENSO (HWday_obs,EL_Nday,LA_Nday,NO_Nday) !iel=0; ila=0; ino=0 !Do ist=1,Nst ! Do ilen=1,Nlen ! Do it=1,Nyear ! if (ENSO_yrs(it)==1) then ! El Nino ! iel=iel+1; El(iel)=Frq_obs(ist,ilen,it) ! else (ENSO_yrs(it)==-1) then ! La Nina ! ila=ila+1; La(ila)=Frq_obs(ist,ilen,it) ! else ! ino=ino+1; No(ino)=Frq_obs(ist,ilen,it) ! endif ! Enddo ! EL_Frq( ! Enddo !Enddo ! Real, Allocatable, Dimension(:,:,:) :: Frq_obs, Dm_obs, Dx_obs ! Nst,Nlen,Nyear ! Real, Allocatable, Dimension(:,:,:) :: Sev_obs, HWday_obs ! Real, Allocatable, Dimension(:,:,:) :: Svx_obs ! Real, Allocatable, Dimension(:,:) :: EL_Frq, EL_Dm, EL_Dx, EL_Sev, EL_Svx, EL_Nday ! Real, Allocatable, Dimension(:,:) :: LA_Frq, LA_Dm, LA_Dx, LA_Sev, LA_Svx, LA_Nday ! Real, Allocatable, Dimension(:,:) :: NO_Frq, NO_Dm, NO_Dx, NO_Sev, NO_Svx, NO_Nday End Subroutine Calc_mean_for_ENSO_phases !------------------------------------------------------------------------------ Subroutine Mean_ENSO (X,E,L,N) Integer :: ist,ilen,it,iel,ila,ino Real, intent(in), Allocatable :: X(:,:,:) !Real, intent(out), Allocatable :: E(:,:),L(:,:),N(:,:) Real, intent(out) :: E(Nst,Nlen),L(Nst,Nlen),N(Nst,Nlen) Real :: El(Nyear),La(Nyear),No(Nyear) Do ist=1,Nst Do ilen=1,Nlen iel=0; ila=0; ino=0 Do it=1,Nyear if (ENSO_yrs(it)==1) then ! El Nino iel=iel+1; El(iel)=X(ist,ilen,it) else if (ENSO_yrs(it)==-1) then ! La Nina ila=ila+1; La(ila)=X(ist,ilen,it) else ino=ino+1; No(ino)=X(ist,ilen,it) endif Enddo E(ist,ilen)=Ave(El,iel,Miss) L(ist,ilen)=Ave(La,ila,Miss) N(ist,ilen)=Ave(No,ino,Miss) Enddo Enddo End Subroutine Mean_ENSO !------------------------------------------------------------------------------ Subroutine Calc_Trend Integer :: ist,ilen Real :: SenSlope,prob Real :: X(Nyear) if (Allocated (Tr_Frq)) Deallocate (Tr_Frq); if (Allocated (Tr_Dm)) Deallocate (Tr_Dm); if (Allocated (Tr_Dx)) Deallocate (Tr_Dx); if (Allocated (Tr_Sev)) Deallocate (Tr_Sev); if (Allocated (Tr_Svx)) Deallocate (Tr_Svx); if (Allocated (Tr_Nday)) Deallocate (Tr_Nday); if (Allocated (pv_Frq)) Deallocate (pv_Frq); if (Allocated (pv_Dm)) Deallocate (pv_Dm); if (Allocated (pv_Dx)) Deallocate (pv_Dx); if (Allocated (pv_Sev)) Deallocate (pv_Sev); if (Allocated (pv_Svx)) Deallocate (pv_Svx); if (Allocated (pv_Nday)) Deallocate (pv_Nday); Allocate (Tr_Frq(Nst,Nlen), Tr_Dm(Nst,Nlen), Tr_Dx(Nst,Nlen)); Allocate (Tr_Sev(Nst,Nlen), Tr_Svx(Nst,Nlen), Tr_Nday(Nst,Nlen)) Tr_Frq=Miss; Tr_Dm=Miss; Tr_Dx=Miss; Tr_Sev=Miss; Tr_Svx=Miss; Tr_Nday=Miss; Allocate (pv_Frq(Nst,Nlen), pv_Dm(Nst,Nlen), pv_Dx(Nst,Nlen)) Allocate (pv_Sev(Nst,Nlen), pv_Svx(Nst,Nlen), pv_Nday(Nst,Nlen)) pv_Frq=Miss; pv_Dm=Miss; pv_Dx=Miss; pv_Sev=Miss; pv_Svx=Miss; pv_Nday=Miss; Do ist=1,Nst Do ilen=1,Nlen X(1:Nyear) = Frq_obs(ist,ilen,1:Nyear) !print*,X;pause Call MK_Sensl(X,Nyear,SenSlope,prob,Miss) Tr_Frq(ist,ilen) = SenSlope; pv_Frq(ist,ilen) = prob !print*,SenSlope,prob X(1:Nyear) = Dm_obs (ist,ilen,1:Nyear) Call MK_Sensl(X,Nyear,SenSlope,prob,Miss) Tr_Dm (ist,ilen) = SenSlope; pv_Dm (ist,ilen) = prob X(1:Nyear) = Dx_obs (ist,ilen,1:Nyear) Call MK_Sensl(X,Nyear,SenSlope,prob,Miss) Tr_Dx (ist,ilen) = SenSlope; pv_Dx (ist,ilen) = prob X(1:Nyear) = Sev_obs(ist,ilen,1:Nyear) Call MK_Sensl(X,Nyear,SenSlope,prob,Miss) Tr_Sev(ist,ilen) = SenSlope; pv_Sev(ist,ilen) = prob X(1:Nyear) = Svx_obs(ist,ilen,1:Nyear) Call MK_Sensl(X,Nyear,SenSlope,prob,Miss) Tr_Svx(ist,ilen) = SenSlope; pv_Svx(ist,ilen) = prob X(1:Nyear) = HWday_obs(ist,ilen,1:Nyear) Call MK_Sensl(X,Nyear,SenSlope,prob,Miss) Tr_Nday(ist,ilen) = SenSlope; pv_Nday(ist,ilen) = prob Enddo Enddo print*,"Calling Calc_Trend ... Done" End Subroutine Calc_Trend !------------------------------------------------------------------------------ Subroutine Calc_Chars_obs (Drou_Thres) Real :: Drou_Thres Integer :: it,it1,it2,ist,ilen Real, Allocatable :: X1d(:), Y1d(:), M1(:), Spells(:), I_Start(:), I_End(:) Integer :: Nspls,Nspl3,istar,iend,Max_evt Real :: IMAX1, N_Spell_1, duration Real, dimension (Nyear) :: Frq, Dm, Dx, Nd, Sev, Svx Integer :: i,j,iyr Max_evt = Ntime Allocate (X1d(Ntime), Y1d(Ntime)) Allocate (Spells(Max_evt), I_Start(Max_evt), I_End(Max_evt)) Allocate (M1(Max_evt)) !--------------------------------------------------------------------- if (Allocated (Frq_obs)) Deallocate (Frq_obs); if (Allocated (Dm_obs)) Deallocate (Dm_obs); if (Allocated (Dx_obs)) Deallocate (Dx_obs); if (Allocated (Sev_obs)) Deallocate (Sev_obs); if (Allocated (HWday_obs)) Deallocate (HWday_obs); if (Allocated (Svx_obs)) Deallocate (Svx_obs); if (Allocated (Frq_obs_all)) Deallocate (Frq_obs_all); if (Allocated (Dm_obs_all)) Deallocate (Dm_obs_all); if (Allocated (Dx_obs_all)) Deallocate (Dx_obs_all); if (Allocated (Sev_obs_all)) Deallocate (Sev_obs_all); if (Allocated (Svx_obs_all)) Deallocate (Svx_obs_all); if (Allocated (HWday_obs_all)) Deallocate (HWday_obs_all); !--------------------------------------------------------------------- Allocate (Frq_obs(Nst,Nlen,Nyear),Dm_obs(Nst,Nlen,Nyear),Dx_obs(Nst,Nlen,Nyear)) Allocate (Sev_obs(Nst,Nlen,Nyear),HWday_obs(Nst,Nlen,Nyear)) Allocate (Svx_obs(Nst,Nlen,Nyear)) Frq_obs=Miss; Dm_obs=Miss; Dx_obs=Miss; Sev_obs=Miss; HWday_obs=Miss; Svx_obs=Miss; !--------------------------------------------------------------------- Allocate (Frq_obs_all(Nst,Nlen),Dm_obs_all(Nst,Nlen),Dx_obs_all(Nst,Nlen)) Allocate (Sev_obs_all(Nst,Nlen),Svx_obs_all(Nst,Nlen),HWday_obs_all(Nst,Nlen)) Frq_obs_all=Miss; Dm_obs_all=Miss; Dx_obs_all=Miss; Sev_obs_all=Miss; HWday_obs_all=Miss; Svx_obs_all=Miss; !--------------------------------------------------------------------- print*, "Calling Calc_Chars_obs...Nst,Nlen,Nyear: ",Nst,Nlen,Nyear ! Obs: HWevt_obs(Nst,Ntime,Nlen) Do ist=1,Nst Do ilen=1,Nlen !if (Mask(ilon,ilat)/=Miss) then X1d(1:Ntime) = HWevt_obs(ist,1:Ntime,ilen) ! Timeseries of Event days Y1d(1:Ntime) = SPEI_obs (ist,1:Ntime,ilen) ! Timeseries of severity Call NSpell_and_Leng (X1d, Ntime, M1, Imax1, Spells, N_Spell_1, I_Start, I_End, Miss) Nspls = N_Spell_1 ! So dot Frq=0; Dm=0; Dx=0; Nd=0; Sev=0; Svx=0 ! Freq, Mean Dur, Max Dur, NHWday, mean Sev Max Sev Do i=1,Nspls duration = Spells(i) ! I_End(i) - I_Start(i) + 1 if (duration >= 1.) then ! Dot han >=2 thang istar = Int(I_Start(i)) ! Thu tu Ngay bat mot HW iend = Int(I_End(i)) ! Thu tu Ngay ket thuc mot HW iyr = istar/Nmon+1; if (mod(istar,Nmon)==0) iyr=iyr-1 ! Thu tu Nam co HW Frq(iyr) = Frq(iyr) + 1 ! Number of HW event in year iyr Dm (iyr) = Dm (iyr) + duration ! Total duration of HW in year iyr Nd (iyr) = Dm (iyr) ! Number of HW days if (Dx(iyr)= |Drou_Thres| Enddo Sev(iyr) = Sev(iyr) + tmp ! Total severity in year iyr if (Svx(iyr)0) then ! Trung binh nam Dm (iyr) = Dm (iyr)/Frq(iyr) Sev(iyr) = Sev(iyr)/Frq(iyr) endif Enddo Frq_obs(ist,ilen,1:Nyear) = Frq(1:Nyear) ! Freq (evt) Dm_obs (ist,ilen,1:Nyear) = Dm (1:Nyear) ! Mean duration (day) Dx_obs (ist,ilen,1:Nyear) = Dx (1:Nyear) ! Max duration (day) HWday_obs(ist,ilen,1:Nyear) = Nd(1:Nyear) ! Num of HW days (day) Sev_obs(ist,ilen,1:Nyear) = Sev(1:Nyear) ! Mean sev (C) Svx_obs(ist,ilen,1:Nyear) = Svx(1:Nyear) ! Max sev (C) ! Clim mean: Frq_obs_all(ist,ilen) = Ave(Frq,Nyear,Miss)! Ave_No_Zero(Frq,Nyear,Miss) ! Clim mean Dm_obs_all (ist,ilen) = Ave_No_Zero(Dm, Nyear,Miss) Dx_obs_all (ist,ilen) = Ave_No_Zero(Dx, Nyear,Miss) Sev_obs_all(ist,ilen) = Ave_No_Zero(Sev,Nyear,Miss) HWday_obs_all(ist,ilen)= Ave_No_Zero(Nd, Nyear,Miss) Svx_obs_all(ist,ilen) = Ave_No_Zero(Svx,Nyear,Miss) !endif Enddo Enddo print*, "Calling Calc_Chars_obs Done" End Subroutine Calc_Chars_obs !------------------------------------------------------------------------------ Subroutine Calc_Dr_Evnts_obs (Drou_Thres) Real :: Drou_Thres Integer :: ist,ilen,it,iday,iday1 Real :: sev ! SPEI_obs(Nst,Ntime,Nlen) if (Allocated (HWevt_obs)) Deallocate (HWevt_obs) ! Used as Dr_events Allocate (HWevt_obs(Nst,Ntime,Nlen)) HWevt_obs=Miss; print*,"Calling Calc_Dr_Evnts_obs ... Nst,Nlen,Ntime: ",Nst,Nlen,Ntime ! Obs: Do ist=1,Nst Do ilen=1,Nlen Do it=1,Ntime sev = SPEI_obs (ist,it,ilen) ! HW severity = Tx - Qp if (sev/=Miss) then if (sev<=Drou_Thres) then ! simultaneously occur both Drought and HW HWevt_obs(ist,it,ilen) = 1 ! HW YES else HWevt_obs(ist,it,ilen) = 0 endif endif Enddo Enddo Enddo End Subroutine Calc_Dr_Evnts_obs !------------------------------------------------------------------------------ Subroutine Read_Drough_IDX (Fname,dridx) Integer :: ist, ilon, ilat, it, it1, it2, Nt, imon, imon1, iyr Integer :: nlat1, nlat0, nlon, nlat, Nle Character (Len=*) :: Fname,dridx ! dridx="spi","spei" Real :: tmp, xmean(12), ic(12) Real, Allocatable :: Stn(:), time(:), Xlen(:), X1d(:), Y1d(:), X3d(:,:,:) !--------------------------------------------------------------------------------- if (Allocated (SPEI_obs)) Deallocate (SPEI_obs) if (Allocated (xlon)) Deallocate (xlon) if (Allocated (xlat)) Deallocate (xlat) if (Allocated (Reg)) Deallocate (Reg) if (Allocated (Stname)) Deallocate (Stname) !--------------------------------------------------------------------------------- it1=Nmon+1 Call rd_ncdf_3d (trim(Fname), Stn, time, xlen, Nst, Nt, nle, X3d, trim(dridx)) ! dridx="spi","spei" Where (IsNaN(X3d)) X3d=Miss Allocate (SPEI_obs(Nst,Ntime,Nlen)); SPEI_obs(1:Nst,1:Ntime,1:Nlen) = X3d(1:Nst,it1:Nt,1:Nlen) Deallocate (Stn, time, xlen, X3d) Call rd_ncdf_1d (Fname, Stn, Nst, xlon, "lon") Deallocate (Stn) Call rd_ncdf_1d (Fname, Stn, Nst, xlat, "lat") Deallocate (Stn) Call rd_ncdf_1d (Fname, Stn, Nst, Reg, "reg") Deallocate (Stn) Call rd_ncdf_1d_Char (Fname, Stn, Nst, Stname, "stname") Deallocate (Stn) End Subroutine Read_Drough_IDX !------------------------------------------------------------------------------ Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Ave.f90" !Include "/home/tanpv/Mylibs/Sum.f90" Include "/home/tanpv/Mylibs/Sub_NSpell_and_Leng.f90" Include "/home/tanpv/Mylibs/Sub_Trend_MK.f90" Include "/home/tanpv/Mylibs/Ave_No_Zero.f90" Include "Sub_GrADS.f90" !Include "/home/tanpv/Mylibs/Sub_SPEI_Daily_One_Stn_2024.f90" !Include "/home/tanpv/Mylibs/Sub_Daily_Percentile.f90" !-------------------------------------------------------------------- END ! Program