Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: Miss=-99.0 Integer, Parameter :: Nmon=12, N365=365 !, Yr1=1961, Yr2=2020 ! Distribution function's parameters: ! Real, Allocatable :: beta_spi(:), gamm(:), pzero(:) ! Real :: beta_spei(0:2), logLogisticParams(0:Nmon,0:2), DisPars(0:Nmon,0:2) Real, Allocatable :: Lon(:), Lat(:), xlon(:), xlat(:), time(:), xmod(:), Mask(:,:) Integer :: Nlon, Nlat, it, ilon, ilat !------------------------------------------------------------------------------ ! Real, Allocatable, Dimension(:,:,:) :: Pr, Temp, Tmax, Tmin ! Nlon,Nlat,Ntime !------------------------------------------------------------------------------ Real, Allocatable, Dimension(:,:,:) :: Tx_obs ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:,:) :: Tx_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: Tm_obs ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:,:) :: Tm_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: Pr_obs ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:,:) :: Pr_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: Qp_tx_obs ! Nlon,Nlat,N365 (Percentile for HW Real, Allocatable, Dimension(:,:,:,:) :: Qp_tx_mod ! Nlon,Nlat,Nmodel,N365 !------------------------------------------------------------------------------ Real, Allocatable, Dimension(:,:,:) :: HWevt_obs ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:,:) :: HWevt_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: HWsev_obs ! Nlon,Nlat,Ntime Real, Allocatable, Dimension(:,:,:,:) :: HWsev_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: SPEI_obs ! Nlon,Nlat,Ntime - SPEI with a fixed timescale in days (e.g. 31 days, 91 days) Real, Allocatable, Dimension(:,:,:,:) :: SPEI_mod ! Nlon,Nlat,Nmodel,Ntime Real, Allocatable, Dimension(:,:,:) :: DrEvt_obs ! Nlon,Nlat,Ntime - Drought events with a fixed threshold Real, Allocatable, Dimension(:,:,:,:) :: DrEvt_mod ! Nlon,Nlat,Nmodel,Ntime !------------------------------------------------------------------------------ ! Real, Allocatable, Dimension(:,:,:) :: Qp_tx ! Nlon,Nlat,N365 (Percentile for HW ! Real, parameter :: p=0.9 ! p=P(X>Qp) ! Integer, Parameter :: Mwd=15 ! data window for calculating percentile Real :: p ! Read from Namelist Integer :: Mwd ! Read from Namelist !------------------------------------------------------------------------------ Character (Len=255) :: Fname,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, Timescale, Ntime Character (Len=20), Allocatable :: Modelnm(:) Character (Len=6) :: Tx_nm, Tm_nm, Pr_nm ! Varname in files Logical :: status, Empty, calc_qp ! Integer, Allocatable :: TimeScale(:) Integer :: Nlo,Nla,Nt,Nz !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:,:,:,:) :: X4d Real, allocatable, Dimension (:) :: X1d, Xt, SPEI_out, Rain, T2m, Tx, Tm !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Mask_File,Obs_Dir,Inp_Dir,Out_Dir Namelist /Data_Info/ Yr1, Yr2, Nmodel, Timescale, calc_qp, p, Mwd Namelist /Modelname/ Modelnm, Tx_nm, Tm_nm, Pr_nm !------------------------------------------------------------------------------ Open (1, file="Namelist_CDH_SEA.txt",Status="old") Read (1, File_Dir_Names) Read (1, Data_Info) Allocate (Modelnm(Nmodel)) Read (1, Modelname) Ntime=(Yr2-Yr1+1)*N365 ! Read landmask data for SEA region with resolution of CORDEX SEA 2024 (~0.25 deg) Call rd_ncdf_2d (Mask_File,Lon,Lat,Nlon,Nlat,Mask,"land") ! Land points !------------------------------------------------------------------------ ! Read Tx and Tm of (CHIRTS) and Pr (CHIRPS) daily data used as observations Call Read_CHIRTS_Tx_Tm_Pr ! Read obs data (CHIRTS, and CHIRPS) ! Read Tx, Tm and Pr model data from CORDEX SEA: Call Read_Model_Data ! Read model data !------------------------------------------------------------------------ ! Since the calculation of Qp is time-consuming, it is necessary to store the results in files for future use: Fname_evt=trim(Out_Dir)//"Qp_obs.nc" Fname_sev=trim(Out_Dir)//"Qp_mod.nc" if (calc_qp) then ! Need to calculate Qp again Call Calc_percentile_Tx_obs Call Calc_percentile_Tx_mod Call Write_Qp (Fname_evt,Fname_sev) else ! Read from file Call rd_ncdf_4d (Fname_sev,xlon,xlat,xmod,time,nlo,nla,nz,nt,Qp_tx_mod,"qp") Deallocate (xlon,xlat,xmod,time) Call rd_ncdf_3d (Fname_evt,xlon,xlat,time,nlo,nla,nt,Qp_tx_obs,"qp") Deallocate (xlon,xlat,time) endif !------------------------------------------------------------------------ ! Call Calc_HW_events Fevt_obs=trim(Out_Dir)//"HW_evt_obs.nc" Fevt_mod=trim(Out_Dir)//"HW_evt_mod.nc" Fsev_obs=trim(Out_Dir)//"HW_sev_obs.nc" Fsev_mod=trim(Out_Dir)//"HW_sev_mod.nc" ! Call Write_HWevt_and_HWsev (Fevt_obs,Fevt_mod,Fsev_obs,Fsev_mod) Call Calc_spei_Ha (TimeScale) Fname = trim(Out_Dir)//"SPEI_mod.nc" Fname_evt = trim(Out_Dir)//"SPEI_obs.nc" Call Write_Drought_IDX (Fname,Fname_evt) !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ Subroutine Calc_HW_events Integer :: imod,ilon,ilat,it,iday !-------- Obs: ----------- Allocate (HWevt_obs(Nlon,Nlat,Ntime)); HWevt_obs=Miss Allocate (HWsev_obs(Nlon,Nlat,Ntime)); HWsev_obs=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then Do it=1,Ntime iday=mod(it,N365); if (iday==0) iday=N365 if (Tx_obs(ilon,ilat,it)/=Miss.AND.Qp_tx_obs(ilon,ilat,iday)/=Miss) then if (Tx_obs(ilon,ilat,it)>=Qp_tx_obs(ilon,ilat,iday)) then HWevt_obs(ilon,ilat,it)=1 HWsev_obs(ilon,ilat,it)=Tx_obs(ilon,ilat,it)-Qp_tx_obs(ilon,ilat,iday) else HWevt_obs(ilon,ilat,it)=0.0 HWsev_obs(ilon,ilat,it)=0.0 endif endif Enddo endif Enddo Enddo !-------- Models: ----------- Allocate (HWevt_mod(Nlon,Nlat,Nmodel,Ntime)); HWevt_mod=Miss Allocate (HWsev_mod(Nlon,Nlat,Nmodel,Ntime)); HWsev_mod=Miss Do imod=1,Nmodel Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then Do it=1,Ntime iday=mod(it,N365); if (iday==0) iday=N365 if (Tx_mod(ilon,ilat,imod,it)/=Miss.AND.Qp_tx_mod(ilon,ilat,imod,iday)/=Miss) then if (Tx_mod(ilon,ilat,imod,it)>=Qp_tx_mod(ilon,ilat,imod,iday)) then HWevt_mod(ilon,ilat,imod,it)=1 HWsev_mod(ilon,ilat,imod,it)=Tx_mod(ilon,ilat,imod,it)-Qp_tx_mod(ilon,ilat,imod,iday) else HWevt_mod(ilon,ilat,imod,it)=0.0 HWsev_mod(ilon,ilat,imod,it)=0.0 endif endif Enddo endif Enddo Enddo Enddo End Subroutine Calc_HW_events !------------------------------------------------------------------------------ Subroutine Calc_spei_Ha (Timescale) ! Using Hargreaves method Integer :: TimeScale, ilon,ilat,it,imod, ilen Real :: xlat0 print*,"In Calc_spei_all ..." Allocate (Tx(Ntime), Tm(Ntime), Rain(Ntime), SPEI_out(Ntime)) Allocate (SPEI_obs(Nlon,Nlat,Ntime), SPEI_mod(Nlon,Nlat,Nmodel,Ntime)) SPEI_obs=Miss; SPEI_mod=Miss ! OBS: Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then xlat0 = Lat(ilat) Tx (1:Ntime) = Tx_obs(ilon,ilat,1:Ntime) Tm (1:Ntime) = Tm_obs(ilon,ilat,1:Ntime) Rain(1:Ntime) = Pr_obs(ilon,ilat,1:Ntime) Call Calc_SPEI_Harg_Daily_One_Stn (Tx, Tm, Rain, xlat0, Ntime, Timescale, SPEI_out, Miss) SPEI_obs(ilon,ilat,1:Ntime) = SPEI_out(1:Ntime) endif Enddo Enddo ! Models: Do imod=1,Nmodel Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then xlat0 = Lat(ilat) Tx (1:Ntime) = Tx_mod(ilon,ilat,imod,1:Ntime) Tm (1:Ntime) = Tm_mod(ilon,ilat,imod,1:Ntime) Rain(1:Ntime) = Pr_mod(ilon,ilat,imod,1:Ntime) Call Calc_SPEI_Harg_Daily_One_Stn (Tx, Tm, Rain, xlat0, Ntime, TimeScale, SPEI_out, Miss) SPEI_mod(ilon,ilat,imod,1:Ntime) = SPEI_out(1:Ntime) endif Enddo Enddo Enddo End Subroutine Calc_spei_Ha !------------------------------------------------------------------------------ Subroutine Write_Drought_IDX (Fevt,Fobs) Character (Len=*) :: Fevt,Fobs Character (Len=255) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Str Integer :: it,i,ilen,ilon,ilat,ilead,Nlo,Nla,Nz,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, lat_reg Real :: xmod(Nmodel) ! Real, Allocatable, Dimension(:,:,:) :: SPEI_obs ! Nlon,Nlat,Ntime - SPEI with a fixed timescale in days (e.g. 31 days, 91 days) ! Real, Allocatable, Dimension(:,:,:,:) :: SPEI_mod ! Nlon,Nlat,Nmodel,Ntime if (Allocated(time)) Deallocate (time) Allocate (time(Ntime)) t_unit="days since 1983-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="model"; t_nm="day" Do i=1,Ntime; time(i)=i-1; Enddo Do i=1,Nmodel; xmod(i)=i; Enddo Write(Var_ln,"(A,I2.2,A)") "The model drought index (SPEI) with timescale of ",Timescale, " days. Model=1:CNRM, 2:Can, 3:EC-Earth, 4:MIROC, 5:MPI, 6:NorESM" Call wr_ncdf_4d (Fevt,Lon,Lat,xmod,time, Nlon,Nlat,Nmodel,Ntime,SPEI_mod,"mod", & t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Write(Var_ln,"(A,I2.2,A)") "The obs drought index (SPEI) with timescale of ",Timescale, " days calculated from CHIRTS-daily data" Call wr_ncdf_3d (Fobs,Lon,Lat,time,Nlon,Nlat,Ntime,SPEI_obs,"obs",t_unit, var_ln,"-", lon_nm, lat_nm, t_nm, Miss) !Call wr_ncdf_3d_apd (Fevt,SPEI_obs,"obs",Nlon,Nlat,Ntime,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) !Call wr_ncdf_4d_apd (Fevt,X4d,"obs", Nlon,Nlat,Nz,Ntime, & ! Var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) End Subroutine Write_Drought_IDX !------------------------------------------------------------------------------ Subroutine Write_HWevt_and_HWsev (Fevt_obs,Fevt_mod,Fsev_obs,Fsev_mod) Character (Len=*) :: Fevt_obs,Fevt_mod,Fsev_obs,Fsev_mod Character (Len=255) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Str Integer :: it,i,ilen,ilon,ilat,ilead,Nlo,Nla,Nz,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, lat_reg Real :: xmod(Nmodel) if (Allocated(time)) Deallocate (time) Allocate (time(Ntime)) t_unit="days since 1983-01-01 "; var_unit="-"; lon_nm="lon"; lat_nm="lat"; z_nm="model"; t_nm="day" Do i=1,Ntime; time(i)=i-1; Enddo Do i=1,Nmodel; xmod(i)=i; Enddo Var_ln="The model HW events. Model=1:CNRM, 2:Can, 3:EC-Earth, 4:MIROC, 5:MPI, 6:NorESM" Call wr_ncdf_4d (Fevt_mod,Lon,Lat,xmod,time, Nlon,Nlat,Nmodel,Ntime,HWevt_mod,"evt", & t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="The model HW severity. Model=1:CNRM, 2:Can, 3:EC-Earth, 4:MIROC, 5:MPI, 6:NorESM" Call wr_ncdf_4d (Fsev_mod,Lon,Lat,xmod,time, Nlon,Nlat,Nmodel,Ntime,HWsev_mod,"sev", & t_unit,Var_ln,"C deg",lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="The obs HW events calculated from Tx CHIRTS-daily of Climate Hazards Center https://www.chc.ucsb.edu/data/chirtsdaily" Call wr_ncdf_3d (Fevt_obs,Lon,Lat,time,Nlon,Nlat,Ntime,HWevt_obs,"evt",t_unit, var_ln,var_unit, lon_nm, lat_nm, t_nm, Miss) Var_ln="The obs HW severity calculated from Tx CHIRTS-daily of Climate Hazards Center https://www.chc.ucsb.edu/data/chirtsdaily" Call wr_ncdf_3d (Fsev_obs,Lon,Lat,time,Nlon,Nlat,Ntime,HWsev_obs,"sev",t_unit, var_ln,"C deg", lon_nm, lat_nm, t_nm, Miss) ! Call wr_ncdf_3d_apd (Fevt,HWevt_obs,"evt_obs",Nlon,Nlat,Ntime,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Var_ln="The model HW severity. Model=1:CNRM, 2:Can, 3:EC-Earth, 4:MIROC, 5:MPI, 6:NorESM" ! Call wr_ncdf_4d (Fsev,Lon,Lat,xmod,time, Nlon,Nlat,Nmodel,Ntime,HWsev_mod,"sev_mod", & ! t_unit,Var_ln,"-",lon_nm, lat_nm, z_nm, t_nm,Miss) ! Var_ln="The obs HW severity" ! Call wr_ncdf_3d_apd (Fsev,HWsev_obs,"sev_obs",Nlon,Nlat,Ntime,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) !Call wr_ncdf_4d_apd (Fname,SPEI_ha,"spei_ha",Nlon,Nlat,Nlen,Ntime,Var_ln,"-", & ! lon_nm, lat_nm, z_nm, t_nm,Miss) End Subroutine Write_HWevt_and_HWsev !------------------------------------------------------------------------------ Subroutine Write_Qp (Fobs,Fmod) ! write percentiles of obs and models: Character (Len=*) :: Fobs,Fmod Character (Len=255) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Str Integer :: it,i,ilen,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real :: tmpr, tmpt, tmpet, alpha, lat_reg Real :: xmod(Nmodel) if (Allocated(time)) Deallocate (time) Allocate (time(N365)) t_unit="days since 1983-01-01 "; var_unit="C deg"; lon_nm="lon"; lat_nm="lat"; z_nm="model"; t_nm="day" Do i=1,N365; time(i)=i-1; Enddo Do i=1,Nmodel; xmod(i)=i; Enddo Var_ln="The model percentiles. Model=1:CNRM, 2:Can, 3:EC-Earth, 4:MIROC, 5:MPI, 6:NorESM" Call wr_ncdf_4d (Fmod,Lon,Lat,xmod,time, Nlon,Nlat,Nmodel,N365,Qp_tx_mod,"qp", & t_unit,Var_ln,var_unit,lon_nm, lat_nm, z_nm, t_nm,Miss) Var_ln="The obs percentiles." Call wr_ncdf_3d (Fobs,Lon,Lat,time,Nlon,Nlat,N365,Qp_tx_obs,"qp",t_unit, var_ln,var_unit, lon_nm, lat_nm, t_nm, Miss) End Subroutine Write_Qp !------------------------------------------------------------------------------ Subroutine Calc_percentile_Tx_obs Integer :: ilon,ilat,it Real :: Qp(N365) Allocate(Tx(Ntime),Qp_tx_obs(Nlon,Nlat,N365)); Qp_tx_obs=Miss Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then Tx (1:Ntime) = Tx_obs(ilon,ilat,1:Ntime) Call Daily_per_HW (Tx,Ntime,p,Mwd,Qp,Miss) ! p and Mwd are needed to be defined Qp_tx_obs(ilon,ilat,1:N365)=Qp(1:N365) endif Enddo Enddo Deallocate (Tx) End Subroutine Calc_percentile_Tx_obs !------------------------------------------------------------------------------ Subroutine Calc_percentile_Tx_mod Integer :: imod,ilon,ilat,it Real :: Qp(N365) Allocate(Tx(Ntime),Qp_tx_mod(Nlon,Nlat,Nmodel,N365)); Qp_tx_mod=Miss Do imod=1,Nmodel Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)/=Miss) then Tx (1:Ntime) = Tx_mod(ilon,ilat,imod,1:Ntime) Call Daily_per_HW (Tx,Ntime,p,Mwd,Qp,Miss) ! p and Mwd are needed to be defined Qp_tx_mod(ilon,ilat,imod,1:N365)=Qp(1:N365) endif Enddo Enddo Enddo Deallocate (Tx) End Subroutine Calc_percentile_Tx_mod !------------------------------------------------------------------------------ Subroutine Read_CHIRTS_Tx_Tm_Pr Real, parameter :: FillValue=-9999.0 Integer :: iyr,it,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,icas, it1, it2 Character (Len=255) :: Fname Real, Allocatable :: X3d(:,:,:), xtime(:), xlon(:),xlat(:) !,xtime(:) Allocate (Tx_obs(Nlon,Nlat,Ntime)); Tx_obs=Miss Allocate (Tm_obs(Nlon,Nlat,Ntime)); Tm_obs=Miss Allocate (Pr_obs(Nlon,Nlat,Ntime)); Pr_obs=Miss ! Tx: it1=0 Do iyr=Yr1,Yr2 write(Fname,"(A,I4.4,A)") trim(Obs_Dir)//"Tx/"//"Tmax.",iyr,".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"Tmax") Where (X3d==FillValue) X3d=Miss Do it=1,Nt it1=it1+1 Tx_obs(:,:,it1) = X3d(:,:,it) Enddo Deallocate (xlon,xlat,xtime,X3d) Enddo ! Tm: it1=0 Do iyr=Yr1,Yr2 write(Fname,"(A,I4.4,A)") trim(Obs_Dir)//"Tm/"//"Tmin.",iyr,".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"Tmin") Where (X3d==FillValue) X3d=Miss Do it=1,Nt it1=it1+1 Tm_obs(:,:,it1) = X3d(:,:,it) Enddo Deallocate (xlon,xlat,xtime,X3d) Enddo ! Pr: it1=0 Do iyr=Yr1,Yr2 write(Fname,"(A,I4.4,A)") trim(Obs_Dir)//"Pr/"//"pr.",iyr,".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"precip") Where (X3d==FillValue) X3d=Miss Do it=1,Nt it1=it1+1 Pr_obs(:,:,it1) = X3d(:,:,it) Enddo Deallocate (xlon,xlat,xtime,X3d) Enddo !-------------- Check --------- if (it1/=Ntime) then print*,"Something wrong during reading obs data." stop endif End Subroutine Read_CHIRTS_Tx_Tm_Pr !------------------------------------------------------------------------------ Subroutine Read_Model_Data Integer :: imod,it,ilon,ilat,ilead,Nlo,Nla,Nz,Nt,icas, it1, it2 Character (Len=255) :: Fname Real, Allocatable :: X4d(:,:,:,:), X3d(:,:,:), xtime(:), xlon(:),xlat(:),m2(:) !,xtime(:) Allocate (Tx_mod(Nlon,Nlat,Nmodel,Ntime)); Tx_mod=Miss Allocate (Tm_mod(Nlon,Nlat,Nmodel,Ntime)); Tm_mod=Miss Allocate (Pr_mod(Nlon,Nlat,Nmodel,Ntime)); Pr_mod=Miss Do imod=1,Nmodel ! Tx: Fname=trim(Inp_Dir)//trim(Modelnm(imod))//"/hist/"//trim(Tx_nm)//"_"//trim(Modelnm(imod))//"_historical_daily.nc" Call rd_ncdf_4d (Fname,xlon,xlat,m2,xtime,nlo,nla,nz,nt,X4d,trim(Tx_nm)) !Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,trim(Varnm)) if (Nt/=Ntime) then print*,"Something wrong during reading data file: ",trim(Fname) stop endif Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)>0) then Tx_mod(ilon,ilat,imod,1:Ntime) = X4d(ilon,ilat,1,1:Ntime) endif Enddo Enddo Deallocate (xlon,xlat,m2,xtime,X4d) ! Tm: Fname=trim(Inp_Dir)//trim(Modelnm(imod))//"/hist/"//trim(Tm_nm)//"_"//trim(Modelnm(imod))//"_historical_daily.nc" Call rd_ncdf_4d (Fname,xlon,xlat,m2,xtime,nlo,nla,nz,nt,X4d,trim(Tm_nm)) !Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,trim(Varnm)) if (Nt/=Ntime) then print*,"Something wrong during reading data file: ",trim(Fname) stop endif Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)>0) then Tm_mod(ilon,ilat,imod,1:Ntime) = X4d(ilon,ilat,1,1:Ntime) endif Enddo Enddo Deallocate (xlon,xlat,m2,xtime,X4d) ! Pr: Fname=trim(Inp_Dir)//trim(Modelnm(imod))//"/hist/"//trim(Pr_nm)//"_"//trim(Modelnm(imod))//"_historical_daily.nc" !Call rd_ncdf_4d (Fname,xlon,xlat,m2,xtime,nlo,nla,nz,nt,X4d,trim(Pr_nm)) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,trim(Pr_nm)) if (Nt/=Ntime) then print*,"Something wrong during reading data file: ",trim(Fname) stop endif Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)>0) then !Pr_mod(ilon,ilat,imod,1:Ntime) = X4d(ilon,ilat,1,1:Ntime) Pr_mod(ilon,ilat,imod,1:Ntime) = X3d(ilon,ilat,1:Ntime) endif Enddo Enddo !Deallocate (xlon,xlat,m2,xtime,X4d) Deallocate (xlon,xlat,xtime,X3d) Enddo End Subroutine Read_Model_Data !------------------------------------------------------------------------------ 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/Sub_SPEI_Daily_One_Stn_2024.f90" Include "/home/tanpv/Mylibs/Sub_Daily_Percentile.f90" !-------------------------------------------------------------------- END ! Program