Use NetCDF Real, Parameter :: Miss=-99. Integer, Parameter :: Yr1=1976, Yr2=2005, Nmon=12, Nmem=7, Nchar=4 Character (Len=10), parameter :: Chname(Nchar)=(/"evt","dur_m","sev_m","inten_m"/) Character (Len=100), parameter :: fmem(Nmem)=(/"CNRM","CSIRO","ECEA","GFDL","HadG","MPI","ENS"/) Character (Len=4) :: BC ! = "ORG" or "QM" Real, Parameter :: pd0=0.30, pd1=0.20, pd2=0.10, pd3=0.05, pd4=0.02 ! Probability to calculate percentiles ! Mo ta thong tin SL dau vao: !--------------------------------------------------------------------------------- Real, allocatable :: xlon(:), xlat(:), xtime(:), Lon(:), Lat(:) Integer :: Nlon, Nlat, Ntime Logical :: Status Real, allocatable :: Mask(:,:), Mask_all(:,:) ! (Nlon,Nlat) !--------------------------------------------------------------------------------- Real, allocatable :: Char_ORG_ENS(:,:,:) ! Ensemble mean or Ensemble median over all members Real, allocatable :: Geo_ORG_ENS(:) ! As above Real, allocatable :: Char_QM_ENS(:,:,:) ! Ensemble mean or Ensemble median over all members Real, allocatable :: Geo_QM_ENS(:) ! As above Real, allocatable :: Char_APH(:,:,:) ! (Nlon,Nlat,Nchar) Real, allocatable :: Char_ORG(:,:,:,:) ! (Nlon,Nlat,Nchar,Nmem) - Number of Events Real, allocatable :: Char_QM (:,:,:,:) ! (Nlon,Nlat,Nchar,Nmem) - Number of Events Real, allocatable :: AVs (:,:,:,:) ! (Nlon,Nlat,Nchar,Nmem) Real, allocatable :: Df_Or_APH(:,:,:,:) ! (Nlon,Nlat,Nchar,Nmem) Real, allocatable :: Df_QM_APH(:,:,:,:) ! (Nlon,Nlat,Nchar,Nmem) Real, allocatable :: Geo_APH(:) ! Ntime Real, allocatable :: Geo_ORG(:,:) ! (Ntime,Nmem) Real, allocatable :: Geo_QM(:,:) ! (Ntime,Nmem) Real, allocatable :: Df_Geo_ORG(:,:) ! (Ntime,Nmem) Real, allocatable :: Df_Geo_QM(:,:) ! (Ntime,Nmem) Real, allocatable :: AV_Geo(:,:) ! (Ntime,Nmem) ! For scatterplot and boxplot: Real, allocatable :: Char_APH_1(:,:,:) ! (Nchar,Landpoint) <<< Chuyen thanh Nmem+1 Real, allocatable :: Char_ORG_1(:,:,:) ! (Nchar,Nmem+1,Landpoint) - Number of Events Real, allocatable :: Char_QM_1 (:,:,:) ! (Nchar,Nmem+1,Landpoint) - Number of Events Real :: Positive_AV(Nmem+1,Nchar+1) Integer :: Nst, iidx, Nyear, Landpoint, Landpts Integer, Parameter :: NDr_idx=3 ! PDSI, SPEI12, SPI12 Character (Len=5), Parameter :: IDX_nm(NDr_idx)= & (/"pdsi","spi12","spe12"/) !(/"PDSI","SPI12","SPE12"/) Character (Len=100) :: Fname, Dir, Fch1, Fch2, Fch3, Fch4 !--------------------------------------------------------------------------------- Call Read_Chars Call Write_all_char_to_NC("NC1/All_Cha_landpoints.nc") Call Calc_AVs Call Write_to_NC("NC/AVs.nc","NC/Del_Geo.nc","Positive_AV.txt") CONTAINS !--------------------------------------------------------------------------------- Subroutine Read_Chars Integer :: ich, ilon, ilat, it, Nt, Nx, imon, imon1, iyr, ic Integer :: Nyear, Nmon, nlat1, nlat0 Real, Allocatable :: X3d(:,:,:), X2d(:,:),xtime(:),idx(:) Real :: Ax(Nmem), a,b,c Fname="/work/users/hmo/tanpv/CORDEX_SEACLID/Land_mask_sea.nc" Call rd_ncdf_2d (trim(Fname), xlon,xlat, nlon,nlat,Mask,"land") ! Landpoint ic=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Mask(ilon,ilat)>0) ic=ic+1; Enddo Enddo Landpoint=ic; Deallocate (xlon, xlat) ! APH Chars: Fname="NC/Chars_APH_BL.nc" Do ich=1,Nchar Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, X3d, trim(Chname(ich))) if (ich==1) then Allocate (Char_APH(nlon,nlat,Nchar),Lon(Nlon),Lat(Nlat)) Lon=xlon; Lat=xlat endif Do ilon=1,Nlon Do ilat=1,Nlat Char_APH(ilon,ilat,ich)=X3d(ilon,ilat,1) ! SPEI Enddo Enddo Deallocate (xlon, xlat, xtime,X3d) Enddo ! ORG Chars: BC="ORIG"; Allocate (Char_ORG(nlon,nlat,Nchar,Nmem)) Do imem=1,Nmem Fname="NC/Chars_"//trim(BC)//"_"//trim(fmem(imem))//".nc" Do ich=1,Nchar Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, X3d, trim(Chname(ich))) Do ilon=1,Nlon Do ilat=1,Nlat Char_ORG(ilon,ilat,ich,imem)=X3d(ilon,ilat,1) Enddo Enddo Deallocate (xlon, xlat, xtime,X3d) Enddo Enddo ! QM Chars: BC="QM"; Allocate (Char_QM(nlon,nlat,Nchar,Nmem)) Do imem=1,Nmem Fname="NC/Chars_"//trim(BC)//"_"//trim(fmem(imem))//".nc" Do ich=1,Nchar Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, X3d, trim(Chname(ich))) Do ilon=1,Nlon Do ilat=1,Nlat Char_QM(ilon,ilat,ich,imem)=X3d(ilon,ilat,1) Enddo Enddo Deallocate (xlon, xlat, xtime,X3d) Enddo Enddo ! Geos: Fname="NC/Geo_APH_BL.nc" Call rd_ncdf_2d (Fname,xtime,idx,Nt,Nx,X2d,"geo") Ntime=Nt Allocate (Geo_APH(Ntime)) Geo_APH(:)=X2d(:,1) Deallocate (xtime,idx,X2d) BC="ORIG"; Allocate (Geo_ORG(Ntime,Nmem+1)) Do imem=1,Nmem Fname="NC/Geo_"//trim(BC)//"_"//trim(fmem(imem))//".nc" Call rd_ncdf_2d (Fname,xtime,idx,Nt,Nx,X2d,"geo") Geo_ORG(:,imem)=X2d(:,1); Deallocate (xtime,idx,X2d) Enddo BC="QM"; Allocate (Geo_QM(Ntime,Nmem+1)) Do imem=1,Nmem Fname="NC/Geo_"//trim(BC)//"_"//trim(fmem(imem))//".nc" Call rd_ncdf_2d (Fname,xtime,idx,Nt,Nx,X2d,"geo") Geo_QM(:,imem)=X2d(:,1); Deallocate (xtime,idx,X2d) Enddo ! Characteristic-based Ensembles: Allocate (Char_ORG_ENS(Nlon,Nlat,Nchar),Char_QM_ENS(Nlon,Nlat,Nchar)) Allocate (Geo_ORG_ENS(Ntime),Geo_QM_ENS(Ntime)) Char_ORG_ENS=Miss;Char_QM_ENS=Miss; Geo_ORG_ENS=Miss; Geo_QM_ENS=Miss Do ich=1,Nchar Do ilon=1,Nlon Do ilat=1,Nlat Ax(1:Nmem)=Char_ORG(ilon,ilat,ich,1:Nmem) Char_ORG_ENS(ilon,ilat,ich)=Ave(Ax,Nmem,Miss) Ax(1:Nmem)=Char_QM(ilon,ilat,ich,1:Nmem) Char_QM_ENS(ilon,ilat,ich)=Ave(Ax,Nmem,Miss) Enddo Enddo Enddo Do it=1,Ntime Ax(1:Nmem)=Geo_ORG(it,1:Nmem); Geo_ORG_ENS(it)=Ave(Ax,Nmem,Miss) Ax(1:Nmem)=Geo_QM (it,1:Nmem); Geo_QM_ENS (it)=Ave(Ax,Nmem,Miss) Geo_ORG(it,Nmem+1) = Geo_ORG_ENS(it) Geo_QM (it,Nmem+1) = Geo_QM_ENS (it) !print*,"it,Geo_QM_ENS=",it,Geo_QM_ENS (it) !print*,Ax(1:Nmem) Enddo !pause ! For scatterplot and boxplot: Allocate (Char_ORG_1(Nchar,Nmem,Landpoint),Char_QM_1(Nchar,Nmem,Landpoint)) Allocate (Char_APH_1(Nchar,Nmem,Landpoint)) Char_APH_1=Miss; Char_ORG_1=Miss; Char_QM_1=Miss Do ich=1,Nchar Do imem=1,Nmem ic=0 Do ilon=1,Nlon Do ilat=1,Nlat if (imem==1) a=Char_APH(ilon,ilat,ich) b=Char_ORG(ilon,ilat,ich,imem); c=Char_QM(ilon,ilat,ich,imem) if (a/=Miss.AND.b/=Miss.AND.c/=Miss) then ic=ic+1 Char_APH_1(ich,imem,ic)=a ! Same values for all imem Char_ORG_1(ich,imem,ic)=b Char_QM_1 (ich,imem,ic)=c endif Enddo Enddo Enddo Enddo End Subroutine Read_Chars !--------------------------------------------------------------------------------- Subroutine Calc_AVs Integer :: it,ilon,ilat,imem,ich,i,j,k,n1,n2 Real :: a, b, c, a1, a2, d ! Them thanh phan Nem+1 la Ensemble theo Chars: Allocate (AVs(Nlon,Nlat,Nchar,Nmem+1)); AVs=Miss Allocate (Df_Or_APH(Nlon,Nlat,Nchar,Nmem+1),Df_QM_APH(Nlon,Nlat,Nchar,Nmem+1)) Allocate (Df_Geo_ORG(Ntime,Nmem+1),Df_Geo_QM(Ntime,Nmem+1)) Allocate (AV_Geo(Ntime,Nmem+1)) Df_Or_APH=Miss; Df_QM_APH=Miss; Df_Geo_ORG=Miss; Df_Geo_QM=Miss Do imem=1,Nmem Do ich=1,Nchar Do ilon=1,Nlon Do ilat=1,Nlat a = Char_APH(ilon,ilat,ich) b = Char_ORG(ilon,ilat,ich,imem) c = Char_QM (ilon,ilat,ich,imem) if (a/=Miss.AND.b/=Miss.AND.c/=Miss) then !Df_Or_APH(ilon,ilat,ich,imem)=b-a !Df_QM_APH(ilon,ilat,ich,imem)=c-a Df_Or_APH(ilon,ilat,ich,imem)=(b-a)/a Df_QM_APH(ilon,ilat,ich,imem)=(c-a)/a a1=ABS(b-a); a2=ABS(c-a) d = Max(a1,a2) if (d/=0) AVs(ilon,ilat,ich,imem)=(a1-a2)/d ! print"(10F8.3)",a,b,c, Df_Or_APH(ilon,ilat,ich,imem),Df_QM_APH(ilon,ilat,ich,imem); pause endif Enddo Enddo Enddo Do it=1,Ntime a = Geo_APH(it); b = Geo_ORG(it,imem); c = Geo_QM(it,imem) if (a/=Miss.AND.b/=Miss.AND.c/=Miss) then Df_Geo_ORG(it,imem)=b-a Df_Geo_QM (it,imem)=c-a !Df_Geo_ORG(it,imem)=(b-a)/a !Df_Geo_QM (it,imem)=(c-a)/a a1=ABS(b-a); a2=ABS(c-a); d = Max(a1,a2) if (d/=0) AV_Geo(it,imem)=(a1-a2)/d endif Enddo Enddo ! Voi Nmem+1: Do ich=1,Nchar Do ilon=1,Nlon Do ilat=1,Nlat a = Char_APH(ilon,ilat,ich) b = Char_ORG_ENS(ilon,ilat,ich) c = Char_QM_ENS(ilon,ilat,ich) if (a/=Miss.AND.b/=Miss.AND.c/=Miss) then !Df_Or_APH(ilon,ilat,ich,Nmem+1)=b-a !Df_QM_APH(ilon,ilat,ich,Nmem+1)=c-a Df_Or_APH(ilon,ilat,ich,Nmem+1)=(b-a)/a Df_QM_APH(ilon,ilat,ich,Nmem+1)=(c-a)/a a1=ABS(b-a); a2=ABS(c-a); d = Max(a1,a2) if (d/=0) AVs(ilon,ilat,ich,imem)=(a1-a2)/d endif Enddo Enddo Enddo Do it=1,Ntime a = Geo_APH(it); b = Geo_ORG_ENS(it); c = Geo_QM_ENS(it) if (a/=Miss.AND.b/=Miss.AND.c/=Miss) then Df_Geo_ORG(it,Nmem+1)=b-a; Df_Geo_QM (it,Nmem+1)=c-a !Df_Geo_ORG(it,Nmem+1)=(b-a)/a; Df_Geo_QM (it,Nmem+1)=(c-a)/a a1=ABS(b-a); a2=ABS(c-a); d = Max(a1,a2) if (d/=0) AV_Geo(it,Nmem+1)=(a1-a2)/d endif print*,"it,Df_Geo_ORG,Df_Geo_QM ",it,Df_Geo_ORG(it,Nmem+1),Df_Geo_QM (it,Nmem+1) Enddo ! pause ! Positive_AV(Nmem+1,Nchar+1) Do imem=1,Nmem+1 Do ich=1,Nchar n1=0; i=0 Do ilon=1,Nlon Do ilat=1,Nlat if (AVs(ilon,ilat,ich,imem)/=Miss) then n1=n1+1 if (AVs(ilon,ilat,ich,imem)>0) i=i+1 endif Enddo Enddo Positive_AV(imem,ich)=real(i)/real(n1)*100.0 Enddo n2=0; j=0 Do it=1,Ntime if (AV_Geo(it,imem)/=Miss) then n2=n2+1 if (AV_Geo(it,imem)>0) j=j+1 endif Enddo Positive_AV(imem,Nchar+1)=real(j)/real(n2)*100.0 Enddo print*,Positive_AV End Subroutine Calc_AVs !--------------------------------------------------------------------------------- Subroutine Write_all_char_to_NC(Fname) ! Write all characteristics at gridpoints to file for scatter plot: Character (Len=*) :: Fname Character (Len=255) :: var_n, var_ln, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm Real :: x_ch(Nchar), xmem(Nmem), X1(Nchar,Nmem),X2(Nchar,Nmem),X3(Nchar,Nmem) Real, Allocatable, Dimension (:,:,:) :: Y1, Y2, Y3 Real, Allocatable :: x_landpts(:) Integer :: i, ich, imem, ilpt, Nlpt Nlpt=0 Do i=1,Landpoint X1(:,:)=Char_APH_1(:,:,i); X2(:,:)=Char_ORG_1(:,:,i); X3(:,:)=Char_QM_1(:,:,i) !print*,x1;print* !print*,x2;print* !print*,x3;pause !if (Ar2_No_Miss(X1,Nchar,Nmem,Miss).AND.Ar2_No_Miss(X2,Nchar,Nmem,Miss).AND.Ar2_No_Miss(X3,Nchar,Nmem,Miss)) then Nlpt=Nlpt+1 !endif Enddo Allocate (Y1(Nchar,Nmem,Nlpt),Y2(Nchar,Nmem,Nlpt),Y3(Nchar,Nmem,Nlpt),x_landpts(Nlpt)) ilpt=0 Do i=1,Landpoint X1(:,:)=Char_APH_1(:,:,i); X2(:,:)=Char_ORG_1(:,:,i); X3(:,:)=Char_QM_1(:,:,i) !if (Ar2_No_Miss(X1,Nchar,Nmem,Miss).AND.Ar2_No_Miss(X2,Nchar,Nmem,Miss).AND.Ar2_No_Miss(X3,Nchar,Nmem,Miss)) then ilpt=ilpt+1 Y1(:,:,ilpt)=Char_APH_1(:,:,i); Y2(:,:,ilpt)=Char_ORG_1(:,:,i); Y3(:,:,ilpt)=Char_QM_1(:,:,i) !endif Enddo lon_nm="cha"; lat_nm="mem"; t_nm="point"; t_unit="months since 1977-01-01 "; var_unit="-" x_ch(1)=1; x_ch(2)=2; x_ch(3)=3; x_ch(4)=4 Do i=1,Nmem; xmem(i)=i; Enddo; Do i=1,Nlpt; x_landpts(i)=i; Enddo; var_ln="APH Drought characteristics values at all land gridpoints: Cha=1:evt; =2:Dur; =3:Sev; =4:Inten_m. Mem=1..6: CNRM,CSIRO,ECEA,GFDL,HadG,MPI" Call wr_ncdf_3d (Fname,x_ch,xmem,x_landpts,Nchar,Nmem,Nlpt,Y1,"aph",t_unit,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="ORG Drought characteristics values at all land gridpoints: Cha=1:evt; =2:Dur; =3:Sev; =4:Inten_m. Mem=1..6: CNRM,CSIRO,ECEA,GFDL,HadG,MPI" Call wr_ncdf_3d_apd (Fname,Y2,"org",Nchar,Nmem,Nlpt,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="QM Drought characteristics values at all land gridpoints: Cha=1:evt; =2:Dur; =3:Sev; =4:Inten_m. Mem=1..6: CNRM,CSIRO,ECEA,GFDL,HadG,MPI" Call wr_ncdf_3d_apd (Fname,Y3,"qm",Nchar,Nmem,Nlpt,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) Deallocate (Y1, Y2, Y3, x_landpts) End Subroutine Write_all_char_to_NC !--------------------------------------------------------------------------------- Logical Function Ar2_No_Miss(X,N,M,Miss) Integer :: N,M,i,j Real :: X(N,M), Miss Logical :: NoMiss NoMiss=.true. Do i=1,N Do j=1,M if (X(i,j)==Miss) then; NoMiss=.false.; exit; endif Enddo Enddo Ar2_No_Miss=NoMiss End Function Ar2_No_Miss !--------------------------------------------------------------------------------- Subroutine Write_to_NC(Fname,Fgeo,txtf) Character (Len=*) :: Fname,Fgeo,txtf Character (Len=255) :: var_n, var_ln, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm Character (Len=10) :: mem_name Integer :: Nmem1 ! =Nmem+1 Real :: x_ch(Nchar), xmem(Nmem+1), Xaph(Ntime,Nmem+1) Nmem1=Nmem+1 ! Thanh phan Nmem+1 la Ensemble mean over Chars lon_nm="lon"; lat_nm="lat"; t_nm="mem"; z_nm="char" t_unit="years since 1977-01-01 " x_ch(1)=1; x_ch(2)=2; x_ch(3)=3; x_ch(4)=4 Do i=1,Nmem1; xmem(i)=i; Enddo ! Added values and differences between models and APH: var_ln="Added values of QM: Z=1:evt; =2:Dur; =3:Sev; =4:Inten_X. t=1..7: CNRM,CSIRO,ECEA,GFDL,HadG,MPI,Ens" var_unit="-" Call wr_ncdf_4d (Fname,Lon,Lat,x_ch,xmem,Nlon,Nlat,Nchar,Nmem1,AVs,"av", & t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) var_ln="Diff between ORG and APH: Z=1:evt; =2:Dur; =3:Sev; =4:Inten_X. t=1..7: CNRM,CSIRO,ECEA,GFDL,HadG,MPI,Ens" Call wr_ncdf_4d_apd (Fname,Df_Or_APH,"df_or_a",Nlon,Nlat,Nchar,Nmem1, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) var_ln="Diff between QM and APH: Z=1:evt; =2:Dur; =3:Sev; =4:Inten_X. t=1..7: CNRM,CSIRO,ECEA,GFDL,HadG,MPI,Ens" Call wr_ncdf_4d_apd (Fname,Df_QM_APH,"df_qm_a",Nlon,Nlat,Nchar,Nmem1, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) ! GE: Allocate (xtime(Ntime)) Do i=1,Ntime; xtime(i)=i; Enddo var_ln="Differences of geographic extent between ORIG and APH" Call wr_ncdf_2d (Fgeo,xtime,xmem,Ntime,Nmem1,Df_Geo_ORG,"ge_org", var_ln, "%", "time", "model", Miss) var_ln="Differences of geographic extent between QM and APH" Call wr_ncdf_2d_apd (Fgeo,Df_Geo_QM,"ge_qm",Ntime,Nmem1, var_ln, "%","time", "model", Miss) var_ln="Added values of QM for geographic extent" Call wr_ncdf_2d_apd (Fgeo,AV_Geo,"av",Ntime,Nmem1, var_ln, "-","time", "model", Miss) var_ln="ORIG geographic extent" Call wr_ncdf_2d_apd (Fgeo,Geo_ORG,"org",Ntime,Nmem1, var_ln, "-","time", "model", Miss) var_ln="QM geographic extent" Call wr_ncdf_2d_apd (Fgeo,Geo_QM,"qm",Ntime,Nmem1, var_ln, "-","time", "model", Miss) Do i=1,Nmem1; Xaph(1:Ntime,i)=Geo_APH(1:Ntime); Enddo var_ln="APH geographic extent" Call wr_ncdf_2d_apd (Fgeo,Xaph,"aph",Ntime,Nmem1, var_ln, "-","time", "model", Miss) Open (3,file=txtf,status="Unknown") write (3,"(10A10)") " ",(trim(Chname(ich)),ich=1,Nchar),"Geo_ext" Do imem=1,Nmem+1 if (imem<=Nmem) then; mem_name=trim(fmem(imem)); else; mem_name="ENS_Ch"; endif write (3,"(A10,10F10.2)") mem_name, (Positive_AV(imem,ich),ich=1,Nchar+1) Enddo Close(3) End Subroutine Write_to_NC !--------------------------------------------------------------------------------- Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/sub_rd_ncdf_Tan.F90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/sub_wr_ncdf_Tan.F90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/Sub_Percentile.f90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/Sub_NSpell_and_Leng.f90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/Ave.f90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/Stdev.f90" Include "/work/users/hmo/remoclic/Phuong/Dr_SEA_2020/Library/Mylibs/Sum.f90" END Program