Use NetCDF Implicit None Real, Parameter :: Miss=-99.0, Lat_min=-11.0 ! Cat phan Autralia Integer, parameter :: Yr1=1960, Yr2=2019, Nmon=12 Integer, parameter :: Nenso=3, NDr_lev=4 Real, parameter :: D0=-0.5, D1=-1.0, D2=-1.5, D3=-2.0 ! mid, moderate, severe, extreme drought Real, parameter :: El=0.5, La=-0.5 Real, parameter :: W0=0.5, W1=1.0, W2=1.5, W3=2.0 ! mid, moderate, severe, extreme wet Integer :: Nreg, Ntime Logical :: Status !--------------------------------------------------------------------------------- Real, allocatable :: Areal_mean_Dr(:,:) ! (Nreg,Ntime) Real, Allocatable :: ONI(:) Real, Allocatable :: Dr_EN(:,:,:) ! (Nenso,Nreg,Ntime) Real, Allocatable :: Freq_Dr_EN(:,:,:) ! (Nreg,Nenso,NDr_lev) Real, Allocatable :: Wet_EN(:,:,:) ! (Nenso,Nreg,Ntime) Real, Allocatable :: Freq_Wet_EN(:,:,:) ! (Nreg,Nenso,NDr_lev) !--------------------------------------------------------------------------------- Character (Len=255) :: Fname !--------------------------------------------------------------------------------- Ntime=(Yr2-Yr1+1)*Nmon Call Read_Dr_ClimIDX Call ENSO_Drought Fname="/home/tanpv/Progs/2021/NC/Dr_ENSO.txt" Call Freq_Drought_ENSO (Fname) CONTAINS !================================================================================= Subroutine Freq_Drought_ENSO (Fname) Character (Len=*) :: Fname Integer :: it, ireg, ien, ilev Real :: tmp1, nc1, tmp2, nc2, tmp3, nc3 Character (Len=3), Allocatable :: Regn(:) Allocate (Regn(Nreg)) Regn(1)="R1"; Regn(2)="R2"; Regn(3)="R3"; Regn(4)="R4"; Regn(5)="SEA" Allocate (Freq_Dr_EN(Nreg,Nenso,NDr_lev)) Allocate (Freq_Wet_EN(Nreg,Nenso,NDr_lev)) Do ireg=1,Nreg Do ien=1,Nenso Do ilev=1,NDr_lev nc1=0; nc2=0; tmp1=0; tmp2=0 Do it=1,Ntime if (Dr_EN(ien,ireg,it)>0) tmp1=tmp1+1 if (Wet_EN(ien,ireg,it)>0) tmp2=tmp2+1 if (Int(Dr_EN(ien,ireg,it))==ilev) then !.AND. & !Int(Dr_EN(ien,ireg,it))==ilev) then nc1=nc1+1 endif if (Int(Wet_EN(ien,ireg,it))==ilev) then nc2=nc2+1 endif Enddo Freq_Dr_EN(ireg,ien,ilev)=nc1 /tmp1*100.0 !/real(Ntime)*100.0 Freq_Wet_EN(ireg,ien,ilev)=nc2 /tmp2*100.0 !/real(Ntime)*100.0 Enddo print"(A,15F8.1)","Fre_Dr: ",Freq_Dr_EN(ireg,ien,:), SUM(Freq_Dr_EN(ireg,ien,:)) print"(A,15F8.1)","Fre_Wt: ",Freq_Wet_EN(ireg,ien,:), SUM(Freq_Wet_EN(ireg,ien,:)) Enddo Enddo Open (3,File=Fname,status="unknown") write (3,*) "Drought and ENSO" write (3,"(A6,3(A26,6x))") " ", "El Nino", "La Nina", "Neutral" write (3,"(A6,12(A7,1x))") "Reg",("D0","D1","D2","D3",ien=1,Nenso) Do ireg=1,Nreg write (3,"(A3,3x,12F8.2)") Regn(ireg), & ((Freq_Dr_EN(ireg,ien,ilev),ilev=1,NDr_lev),ien=1,Nenso) Enddo write (3,*) "Wet and ENSO" write (3,"(A6,3(A26,6x))") " ", "El Nino", "La Nina", "Neutral" write (3,"(A6,12(A7,1x))") "Reg",("W0","W1","W2","W3",ien=1,Nenso) Do ireg=1,Nreg write (3,"(A3,3x,12F8.2)") Regn(ireg), & ((Freq_Wet_EN(ireg,ien,ilev),ilev=1,NDr_lev),ien=1,Nenso) Enddo Close (3) End Subroutine Freq_Drought_ENSO !--------------------------------------------------------------------------------- Subroutine ENSO_Drought Integer :: it, ireg, ien, i1, i2, i3 Real :: spei Allocate (Dr_EN(Nenso,Nreg,Ntime)) Allocate (Wet_EN(Nenso,Nreg,Ntime)) Dr_EN=Miss; Wet_EN=Miss Do ireg=1,Nreg i1=0; i2=0; i3=0 Do it=1,Ntime spei=Areal_mean_Dr(ireg,it) if (ONI(it)>=El) then !spei=Areal_mean_Dr(ireg,it) Dr_EN(1,ireg,it)=Dr_lev(spei) Wet_EN(1,ireg,it)=Wet_lev(spei) i1=i1+1 else if (ONI(it)<=La) then !spei=Areal_mean_Dr(ireg,it) Dr_EN(2,ireg,it)=Dr_lev(spei) Wet_EN(2,ireg,it)=Wet_lev(spei) i2=i2+1 else !spei=Areal_mean_Dr(ireg,it) Dr_EN(3,ireg,it)=Dr_lev(spei) Wet_EN(3,ireg,it)=Wet_lev(spei) i3=i3+1 endif !print*,Dr_EN(1:Nenso,ireg,it) Enddo print*,"i1, i2, i3, SUM = ", i1, i2, i3, i1+i2+i3 ; pause Enddo End Subroutine ENSO_Drought !--------------------------------------------------------------------------------- Real Function Dr_lev (spei) Real :: spei, tmp if (spei<=D0) tmp=1 if (spei<=D1) tmp=2 if (spei<=D2) tmp=3 if (spei<=D3) tmp=4 Dr_lev=tmp End Function Dr_lev !--------------------------------------------------------------------------------- Real Function Wet_lev (spei) Real :: spei, tmp if (spei>=W0) tmp=1 if (spei>=W1) tmp=2 if (spei>=W2) tmp=3 if (spei>=W3) tmp=4 Wet_lev=tmp End Function Wet_lev !--------------------------------------------------------------------------------- Subroutine Read_Dr_ClimIDX Character (Len=255) :: Fname, Dir Real, Allocatable :: xreg(:), time(:), xscal(:), X3d(:,:,:) Real, Allocatable :: Year(:),Month(:), tmp(:,:) Integer :: it, Nregs, Nt, Nt1, Nt2, nscal, ireg Integer :: Nyr, Nmn, j, i, k, iyr Dir="/work/projects/G09/SEACLID/GPCC_CRU_SEA/" Fname=trim(Dir)//"SPEI_SPI_GPCC_CRU_Reg_59_19_New.nc" Call rd_ncdf_3d (Fname, xreg, xscal, time, Nregs, nscal, Nt,X3d, "spe12") Nreg = Nregs !print*,"Nt,Ntime=", Nt,Ntime; pause Allocate (Areal_mean_Dr(Nreg,Ntime)) Do ireg=1,Nreg Areal_mean_Dr(ireg,1:Ntime)=X3d(ireg,nscal,Nmon+1:Nt) Enddo !Do it=1,Nt; print*,X3d(1:Nregs,nscal,it); enddo; pause !Do it=1,Ntime; print*,Areal_mean_Dr(1:Nregs,it); enddo Deallocate (X3d,xreg,time,xscal) Fname="/work/projects/G09/SEACLID/ClimIDX/Clim_IDX_2020.nc" Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, "oni") Allocate (ONI(Ntime)) j=0 Do iyr=Yr1,Yr2 Do i=1,Nyr if (Year(i)==iyr) then Do k=1,Nmon j=j+1; ONI(j)=tmp(i,k) Enddo endif Enddo Enddo Deallocate (Year,Month,tmp) !Do it=1,Ntime; print "(20F10.3)",Areal_mean_Dr(1:Nregs,it),ONI(it); enddo; pause End Subroutine Read_Dr_ClimIDX !----------------------------------------------------------------------------------! Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" END