Program Dr_IDXs_SEACLID ! Calculating PDSI, SPI, SPEI use netcdf !------------------------------------------------------------------------------ Real, parameter :: amssng = -99., Miss=amssng Integer, Parameter :: NVar = 2, Nmodels=7, Nseason=4, Nq=1000 Integer, Parameter :: NMonth = 12, Nmon=NMonth Integer, Allocatable :: Len_SPI(:) Real, Allocatable :: DatBL(:,:), DatTMP(:,:), Rain(:) ! For calculating SPI, SPEI Real, Allocatable :: beta_spi(:), gamm(:), pzero(:), SPI(:,:,:), SPE(:,:,:) Real, Allocatable, Dimension (:,:) :: mon_ori_pdsi, mon_ori_wplm, mon_ori_zind Real, Allocatable, Dimension (:,:) :: mon_sel_pdsi, mon_sel_wplm, mon_sel_zind Real, Allocatable :: Dat_Mon(:,:,:,:), DatNor(:,:,:) ! For calculating PDSI Real, Allocatable :: AWC(:), xlon(:), xlat(:), Mask(:,:), xlon1(:), xlat1(:) Real, Allocatable :: Mask_Land(:), Lon(:), Lat(:) ! (Nlon1*Nlat1) - All gridpoints !Real, Allocatable, Dimension (:,:) :: Prec, Temp Real, allocatable, dimension (:,:,:) :: Temp_Inp, Rain_Inp Real, Allocatable :: T2m(:,:,:), Pr(:,:,:) !------------------------------------------------------------------------------ Character (Len=255) :: Fname1,Fname, R_Fname, T_Fname, BL_Dir, PDSI_Dir, PDSI_WDir, Data_dir1,Data_dir2,Data_dir3, MaskFname, AWC_Fname, Out_Dir Character (Len=3) :: TVar, RVar Character (Len=50) :: Model_nm(Nmodels) Logical :: status, Empty Integer :: YBas1, YBas2, NYear, maxyrs Integer :: NLen, Nst, Nlon, Nlat, IO, Ntime, Nt, it Integer :: imod, ist, iyr, imw, ivar, isys, system, ilon, ilat !------------------------------------------------------------------------------ Real, Allocatable :: R_ORG(:,:,:,:),T_ORG(:,:,:,:),R_QM(:,:,:,:),T_QM(:,:,:,:) ! (Nlon,Nlat,Nmodels,Ntime) Real, Allocatable :: R_ORG_S(:,:,:),T_ORG_S(:,:,:),R_QM_S(:,:,:),T_QM_S(:,:,:) ! (Nmodels,Nseason,Landpoint*Ntime/4) Real, Allocatable :: R_ORG_q(:,:,:),T_ORG_q(:,:,:),R_QM_q(:,:,:),T_QM_q(:,:,:) ! (Nmodels,Nseason,Nq) Logical :: mms_to_monthly ! True, False Character (Len=4) :: QMBC ! "ORIG", "QM" Integer :: Landpoint,Ntime_S !------------------------------------------------------------------------------ ! For writing to NC file: Real, allocatable, Dimension (:, :, :) :: X1, X2, X3, X4, X5, X6 Real, allocatable, Dimension (:) :: Xt !------------------------------------------------------------------------------ Namelist /File_Dir_Names/ Data_dir1,Data_dir2,Data_dir3, Model_nm, MaskFname, Out_Dir, & TVar, RVar Namelist /Data_Info/ YBas1, YBas2 !------------------------------------------------------------------------------ Open (1, file="Namelist_4season.txt", Status="Old") Read (1, File_Dir_Names) Read (1, Data_Info) Call Proc_land_mask_AWC ! Return xlon, xlat, Nlon, Nlat, Nst, Mask, Mask_Land, Lon, Lat Nyear = (YBas2 - YBas1 + 1); Ntime = Nyear*NMonth; maxyrs = NYear Ntime_S=Landpoint*(Ntime/4) !print*,"Nlon,Nlat,Landpoint,Ntime,Ntime_S: ",Nlon,Nlat,Landpoint,Ntime,Ntime_S; pause Allocate(Temp_Inp(1:Nlon,1:Nlat,1:Ntime)); Allocate(Rain_Inp(1:Nlon,1:Nlat,1:Ntime)); Allocate (R_ORG(Nlon,Nlat,Nmodels,Ntime),T_ORG(Nlon,Nlat,Nmodels,Ntime)) Allocate (R_QM(Nlon,Nlat,Nmodels,Ntime), T_QM(Nlon,Nlat,Nmodels,Ntime)) Allocate (R_ORG_S(Nmodels,Nseason,Ntime_S),T_ORG_S(Nmodels,Nseason,Ntime_S)) Allocate (R_QM_S(Nmodels,Nseason,Ntime_S), T_QM_S(Nmodels,Nseason,Ntime_S)) R_ORG=Miss; T_ORG=Miss; R_QM=Miss; T_QM=Miss; R_ORG_S=Miss; T_ORG_S=Miss; R_QM_S=Miss; T_QM_S=Miss ! Allocate (Rain(Ntime),DatBL(Ntime,Nst),DatTMP(Ntime,Nst)) ! Allocate (beta_spi(Nmonth), gamm(Nmonth), pzero(Nmonth)) ! Allocate (SPI(Nst,Ntime,NLen), SPE(Nst,Ntime,NLen)) ! Allocate (Dat_Mon(NVar,Nlon*Nlat,YBas1:YBas2, Nmonth)) ! Allocate (DatNor (NVar,Nst,Nmonth)) !, Rain_Mon(Nst,Ntime)) QMBC="ORIG"; mms_to_monthly=.True. Do imod=1,Nmodels if (imod0) then Do it=1,Ntime imon=mod(it,Nmon); if (imon==0) imon=Nmon ! print*,"imon: ", imon Select Case (imon) Case (12,1,2) it1=it1+1; R_ORG_S(imod,1,it1)=R_ORG(ilon,ilat,imod,it); T_ORG_S(imod,1,it1)=T_ORG(ilon,ilat,imod,it) R_QM_S (imod,1,it1)=R_QM (ilon,ilat,imod,it); T_QM_S (imod,1,it1)=T_QM (ilon,ilat,imod,it) Case (3,4,5) it2=it2+1; R_ORG_S(imod,2,it2)=R_ORG(ilon,ilat,imod,it); T_ORG_S(imod,2,it2)=T_ORG(ilon,ilat,imod,it) R_QM_S (imod,2,it2)=R_QM (ilon,ilat,imod,it); T_QM_S (imod,2,it2)=T_QM (ilon,ilat,imod,it) Case (6,7,8) it3=it3+1; R_ORG_S(imod,3,it3)=R_ORG(ilon,ilat,imod,it); T_ORG_S(imod,3,it3)=T_ORG(ilon,ilat,imod,it) R_QM_S (imod,3,it3)=R_QM (ilon,ilat,imod,it); T_QM_S (imod,3,it3)=T_QM (ilon,ilat,imod,it) Case (9,10,11) it4=it4+1; R_ORG_S(imod,4,it4)=R_ORG(ilon,ilat,imod,it); T_ORG_S(imod,4,it4)=T_ORG(ilon,ilat,imod,it) R_QM_S (imod,4,it4)=R_QM (ilon,ilat,imod,it); T_QM_S (imod,4,it4)=T_QM (ilon,ilat,imod,it) End Select Enddo endif Enddo Enddo !Enddo print*,it1,it2,it3,it4 Enddo ! Allocate (R_ORG_S(Nmodels,Nseason,Ntime_S),T_ORG_S(Nmodels,Nseason,Ntime_S)) ! Allocate (R_QM_S(Nmodels,Nseason,Ntime_S), T_QM_S(Nmodels,Nseason,Ntime_S)) End Subroutine Make_season_data !------------------------------------------------------------------------------ Subroutine Proc_land_mask_AWC Integer :: ist, ilon, ilat, nlo, nla Real, Allocatable :: xtmp(:,:),lontmp(:),lattmp(:) Call rd_ncdf_2d (trim(MaskFname),xlon,xlat,Nlon,Nlat,Mask, "land") Nst=Nlon*Nlat Allocate (Mask_Land(Nst)) Allocate (lon(Nst), lat(Nst)) !, AWC(Nst)) ist=0; Landpoint=0 Do ilat=1, Nlat Do ilon=1, Nlon ist=ist+1; Mask_Land(ist)=Mask(ilon,ilat); if (Mask(ilon,ilat)>0) Landpoint=Landpoint+1 ! AWC(ist)=xtmp(ilon,ilat) Lon(ist)=xlon(ilon); Lat(ist)=xlat(ilat) Enddo Enddo ! Deallocate (lontmp,lattmp,xtmp) End Subroutine Proc_land_mask_AWC !------------------------------------------------------------------------------ Subroutine Read_Pr_T2m Real, Allocatable :: XTime(:) Real, Allocatable :: xlon(:), xlat(:) Integer :: Nlon1, Nlat1 print*,"==================>In Read_Pr_T2m: " print*,"Reading ",trim(R_Fname)," and ", trim(T_Fname) Call rd_ncdf_3d (trim(R_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, Pr, trim(RVar)) Deallocate (xlon, xlat, XTime) Call rd_ncdf_3d (trim(T_Fname), xlon, xlat, XTime, Nlon1, Nlat1, Nt, T2m, trim(TVar)) print*,"Nlon1, Nlat1, Nt=",Nlon1, Nlat1, Nt Call Proc_Pr_T2m Deallocate (Pr, T2m) End Subroutine Read_Pr_T2m !------------------------------------------------------------------------------ Subroutine Proc_Pr_T2m Integer, parameter :: NMon_Skip = 0 !72 Real :: TMean(Nmonth),RMean(Nmonth),xtas(Nyear),xpr(Nyear),xt75(Nmonth),xp75(Nmonth) Integer :: iyr,i,j,it,ict,icr,imon,imon1 Real, Allocatable :: Xtmp(:,:) print*,"Beginning Proc_Pr_T2m...:", Ntime,Nt,NMon_Skip Where(T2m==-99.9) T2m=Miss; Temp_Inp(:,:,1:Ntime)=T2m(:,:,NMon_Skip+1:Nt); ! SL tu 1976 - 2005 Where(Pr ==-99.9) Pr =Miss; Rain_Inp(:,:,1:Ntime)=Pr(:,:,NMon_Skip+1:Nt) ! Monthly data, mm/s ! Fill up missing data: Do j=1, Nlat Do i=1, Nlon if (Mask(i,j) == 1) then Do imon=1, Nmonth TMean(imon)=0; RMean(imon)=0; ict=0; icr=0 Do it=1, Ntime imon1=mod(it,Nmonth); if (imon1==0) imon1=Nmonth if (imon1==imon) then if (Temp_Inp(i,j,it)/=Miss) then !TMean(imon)=TMean(imon)+Temp_Inp(i,j,it); ict=ict+1 ict=ict+1; xtas(ict)=Temp_Inp(i,j,it) endif if (Rain_Inp(i,j,it)>=0) then !/=Miss) then !RMean(imon)=RMean(imon)+Rain_Inp(i,j,it); icr=icr+1 icr=icr+1; xpr(icr)=Rain_Inp(i,j,it) endif endif Enddo !TMean(imon)=TMean(imon)/Real(ict) !RMean(imon)=RMean(imon)/Real(icr) if (ict>Nyear/3*2.AND.icr>Nyear/3*2) then TMean(imon)=Percentile(xtas,ict,0.5,Miss) RMean(imon)=Percentile(xpr, icr,0.5,Miss) xt75 (imon)=Percentile(xtas,ict,0.75,Miss) xp75 (imon)=Percentile(xpr, icr,0.75,Miss) else TMean(imon)=Miss; RMean(imon)=Miss; xt75 (imon)=Miss; xp75 (imon)=Miss endif Enddo Do it=1, Ntime imon=mod(it,Nmonth); if (imon==0) imon=Nmonth if (Temp_Inp(i,j,it)==Miss) Temp_Inp(i,j,it)=TMean(imon) if (Rain_Inp(i,j,it)==Miss) Rain_Inp(i,j,it)=RMean(imon) if (Temp_Inp(i,j,it)>2.0*xt75(imon)) Temp_Inp(i,j,it)=xt75(imon) if (Rain_Inp(i,j,it)>3.0*xp75(imon)) Rain_Inp(i,j,it)=xp75(imon) Enddo endif Enddo Enddo ! Convert Pr from mm/s to mm/month and K to C: Allocate (Xtmp(1:Nlon,1:Nlat)) it=0 Do iyr=YBas1, YBas2 Do imon=1, Nmonth it=it+1 Xtmp(:,:)=Rain_Inp(:,:,it) !Where (Xtmp/=Miss) Xtmp=Xtmp*DayOfMon(imon,iyr) ! Daily to Monthly if (mms_to_monthly) then Where (Xtmp/=Miss) Xtmp=Xtmp*86400.0*DayOfMon(imon,iyr) ! mm/s to Monthly endif Where (Xtmp/=Miss) Xtmp=Xtmp; Rain_Inp(:,:,it)=Xtmp(:,:) Xtmp(:,:)=Temp_Inp(:,:,it) if (mms_to_monthly) then Where (Xtmp/=Miss) Xtmp=Xtmp-273.15 ! K to C endif Temp_Inp(:,:,it)=Xtmp(:,:) Enddo Enddo ! Assign Miss values to all Sea-points: Do it=1, Ntime Do j=1, Nlat Do i=1, Nlon if (Mask(i,j)==Miss) then Rain_Inp(i,j,it)=Miss; Temp_Inp(i,j,it)=Miss endif Enddo Enddo Enddo Deallocate(Xtmp) End Subroutine Proc_Pr_T2m !------------------------------------------------------------------------------ Subroutine Calc_Q_Q Integer :: imod,isea,it,i,j,k Real :: q(Nq), p Real, Allocatable :: X1d(:) ! Real, Allocatable :: R_ORG_q(:,:,:),T_ORG_q(:,:,:),R_QM_q(:,:,:),T_QM_q(:,:,:) ! (Nmodels,Nseason,Nq) ! x = percentile(Y,N1,p,Miss) print*,"In Calc_Q_Q....... ",Ntime_S; pause Allocate (X1d(Ntime_S)) Allocate (R_ORG_q(Nmodels,Nseason,Nq),T_ORG_q(Nmodels,Nseason,Nq)) Allocate (R_QM_q(Nmodels,Nseason,Nq), T_QM_q(Nmodels,Nseason,Nq)) R_ORG_q=Miss; T_ORG_q=Miss; R_QM_q=Miss; T_QM_q=Miss Do i=1,Nq; q(i)=Real(i)/Real(Nq+1); Enddo Do imod=1,Nmodels Do isea=1,Nseason print*,imod,isea print*,"1" X1d(1:Ntime_S)=R_ORG_S(imod,isea,1:Ntime_S) Do i=1,Nq p = q(i) !print*,p,Ntime_S,X1d(100:110) R_ORG_q(imod,isea,i) = percentile(X1d,Ntime_S,p,Miss) Enddo print*,"2" X1d(1:Ntime_S)=T_ORG_S(imod,isea,1:Ntime_S) Do i=1,Nq p = q(i) T_ORG_q(imod,isea,i) = percentile(X1d,Ntime_S,p,Miss) Enddo print*,"3" X1d(1:Ntime_S)=R_QM_S(imod,isea,1:Ntime_S) Do i=1,Nq p = q(i) R_QM_q(imod,isea,i) = percentile(X1d,Ntime_S,p,Miss) Enddo print*,"4" X1d(1:Ntime_S)=T_QM_S(imod,isea,1:Ntime_S) Do i=1,Nq p = q(i) T_QM_q(imod,isea,i) = percentile(X1d,Ntime_S,p,Miss) Enddo Enddo Enddo End Subroutine Calc_Q_Q !------------------------------------------------------------------------------ Subroutine Write_All_Dr_IDX_to_NetCDF(Fname, Fname1) !, Nlon, Nlat, Ntime) Integer :: ist, jlat, ilon, it, ilen, ix Character (Len=50) :: t_unit, Var_nm, Var_ln, var_unit, lon_nm, lat_nm, t_nm Character (Len=*) :: Fname, Fname1 Real :: X_season(Nseason), X_model(Nmodels) ! Allocate (R_ORG_S(Nmodels,Nseason,Ntime_S),T_ORG_S(Nmodels,Nseason,Ntime_S)) ! Allocate (R_QM_S(Nmodels,Nseason,Ntime_S), T_QM_S(Nmodels,Nseason,Ntime_S)) !if (.Not.Allocated(X1)) Allocate(X1(Nlon, Nlat, Ntime)) !if (.Not.Allocated(X2)) Allocate(X2(Nlon, Nlat, Ntime)) !print*,R_ORG_S(:,:,100:110) !print*,R_QM_S(:,:,100:110) !print*,T_ORG_S(:,:,100:110) !print*,T_QM_S(:,:,100:110) if (.Not.Allocated(Xt)) Allocate(Xt(Ntime_S)) Do it=1, Ntime_S; Xt(it) = it; Enddo Do it=1, Nseason; X_season(it) = it; Enddo Do it=1, Nmodels; X_model(it) = it; Enddo Write(t_unit,"(A,I4.4,A)") "months since ",YBas1,"-01-01" var_unit="mm/month"; lon_nm="Model"; lat_nm="Season"; t_nm="month"; Var_nm="pr_org" Var_ln="Original model Precipitation in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH" Call wr_ncdf_3d (trim(Fname),X_model,X_season,Xt,Nmodels,Nseason,Ntime_S,R_ORG_S,trim(Var_nm),trim(t_unit),trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Var_nm="pr_qm" Var_ln="QM BC model Precipitation in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH" Call wr_ncdf_3d_apd (trim(Fname),R_QM_S,trim(Var_nm),Nmodels,Nseason,Ntime_S,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Var_nm="tas_org" Var_ln="Original model temperature in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH" Call wr_ncdf_3d_apd (trim(Fname),T_ORG_S,trim(Var_nm),Nmodels,Nseason,Ntime_S,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Var_nm="tas_qm" Var_ln="QM BC model temperature in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH" Call wr_ncdf_3d_apd (trim(Fname),T_QM_S,trim(Var_nm),Nmodels,Nseason,Ntime_S,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Deallocate (Xt); Allocate(Xt(Nq)) Do it=1, Nq; Xt(it) = it; Enddo Write(t_unit,"(A,I4.4,A)") "months since ",YBas1,"-01-01" var_unit="mm/month"; lon_nm="Model"; lat_nm="Season"; t_nm="quantile"; Var_nm="pr_org" Var_ln="Original model Precipitation in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH" Call wr_ncdf_3d (trim(Fname1),X_model,X_season,Xt,Nmodels,Nseason,Nq,R_ORG_q,trim(Var_nm),trim(t_unit),trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Var_nm="pr_qm" Var_ln="QM BC model Precipitation in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH" Call wr_ncdf_3d_apd (trim(Fname1),R_QM_q,trim(Var_nm),Nmodels,Nseason,Nq,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Var_nm="tas_org" Var_ln="Original model temperature in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH" Call wr_ncdf_3d_apd (trim(Fname1),T_ORG_q,trim(Var_nm),Nmodels,Nseason,Nq,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) Var_nm="tas_qm" Var_ln="QM BC model temperature in the seasons. Models are: 1-CNRM,2-CSIRO,3-ECEA,4-HadG,5-GFDL,6-MPI,7-APH" Call wr_ncdf_3d_apd (trim(Fname1),T_QM_q,trim(Var_nm),Nmodels,Nseason,Nq,trim(Var_ln),trim(var_unit),trim(lon_nm),trim(lat_nm),trim(t_nm),amssng) print*,"************ Write to file Done !" End Subroutine Write_All_Dr_IDX_to_NetCDF !------------------------------------------------------------------------------ Integer Function DayOfMon(Mon, Year) Integer :: Mon, Year Select Case (Mon) Case (1,3,5,7,8,10,12) DayOfMon = 31 Case (4,6,9,11) DayOfMon = 30 Case (2) if (mod(Year,4)==0.AND.mod(Year,100)/=0) then DayOfMon = 29 else DayOfMon = 28 endif End Select End Function DayOfMon !-------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Sub_RW_Station_nc.f90" Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" 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 Include "Percentile.f90" !-------------------------------------------------------------------- END ! Program