Use NetCDF ! Xac dinh so ngay bao hoat dong (so observation) ! Moi TC co the ton tai lau hon 1 obs trong 1 gridbox Integer, Parameter :: Life_Mx=120, NTC_max=2000 Real, Parameter :: Lon1=102, Lon2=140, Lat1=6, Lat2=26 Real, Parameter :: Dlon=2, Dlat=2, Miss=-99 Integer, Parameter :: Nlon=Int(Lon2-Lon1)/Int(Dlon), Nlat=Int(Lat2-Lat1)/Int(Dlat) Type TC_rec Integer :: ID Integer :: Life Integer :: YY(Life_Mx) Integer :: MM(Life_Mx) Integer :: DD(Life_Mx) Integer :: HH(Life_Mx) Real :: Lon(Life_Mx) Real :: Lat(Life_Mx) Integer :: CD(Life_Mx) Real :: Pm(Life_Mx) End Type Integer :: TC1(Nlon,Nlat) Integer, Allocatable :: Track(:,:,:), tmp(:,:,:) Integer, Allocatable :: Track_Grade(:,:,:,:) Real, Allocatable :: XTrack(:,:,:), xlon(:), xlat(:), xtime(:) !-------------------------------------------------------------------- Real, Allocatable, Dimension (:,:) :: Xall, X50, X60, X70, X80, X90, X00, X01 ! For paper: Real, Allocatable, Dimension (:,:) :: X75_94, X95_14 ! Hai giai doan 75-94, 95-14 Real, Allocatable, Dimension (:,:,:) :: X7594Gr, X9514Gr ! !-------------------------------------------------------------------- Character (Len=100) :: Fname, St Integer :: itc, i, io, ig, iid, k1, k2, k3, k4, k5, k6 Integer :: nyy, nmm, ndd, nhh, NTC, Nlf, NTC1,NGrade Integer :: ilon, ilat, Yes_TC, it, iyr, iy, Yr1, Yr2, Nyear Real :: x,y,x1,y1,x2,y2 Type (TC_rec) :: TC(NTC_max) !-------------------------------------------------------------------- ! Reading TC track data: Fname="/home/tanpv/Data/TCs/TC_Track_JMA.txt" Open (1, file=Fname, Status="Old") io=0; itc=0 Do while (io==0) read (1,"(A)",IOStat=IO) St if (IO/=0.OR.Len(trim(ST))==0) exit itc=itc+1 read (St,*) ig, iid, nlf TC(itc)%ID = iid TC(itc)%Life = nlf Do i=1, nlf read (1,"(A)") St read (St,*) k1, k2, k3, k4, k5, k6 !print*, k1, k2, k3, k4, k5, k6 nyy=k1/1000000; nmm=mod(k1,1000000)/10000; ndd=mod(k1,10000)/100; nhh=mod(k1,100) if (nyy>20) then nyy=1900+nyy else nyy=2000+nyy endif TC(itc)%YY(i)=nyy; TC(itc)%MM(i)=nmm; TC(itc)%DD(i)=ndd; TC(itc)%HH(i)=nhh TC(itc)%CD(i)=k3 TC(itc)%lat(i)=real(k4)/10. TC(itc)%lon(i)=real(k5)/10. TC(itc)%Pm(i)=k6 Enddo !print"(A,3I5.4,3I2.2,I6)","itc, TC(itc)%ID, nlf= ",itc, & ! TC(itc)%ID,TC(itc)%YY(1),TC(itc)%MM(1),TC(itc)%DD(1),TC(itc)%HH(1), nlf !print*,TC(itc)%CD(1), TC(itc)%lat(1), TC(itc)%lon(1), TC(itc)%Pm(1) Enddo NTC = itc ! Total Number of TCs !-------------------------------------------------------------------- ! Determining the TC tracks: Yr1=TC(1)%YY(1); Yr2=TC(NTC)%YY(1); Nyear=Yr2-Yr1+1 Allocate (Track(1:Nlon,1:Nlat,1:Nyear),tmp(1:Nlon,1:Nlat,1:Nyear)) Track(:,:,:) =0 ! Set to Zero !Do itc=1, NTC ! Search for each TC ! Nlf=TC(itc)%Life ! A TC's life time (number of observations) ! ! tmp(:,:,:) = 0 ! Total TC count in all years at each gridbox ! Do ilon=1, Nlon ! x1=Lon1+(ilon-1)*Dlon; x2=x1+Dlon ! Do ilat=1, Nlat ! y1=Lat1+(ilat-1)*Dlat; y2=y1+Dlat ! it=0 ! Do iyr=Yr1, Yr2 ! Search for all period ! it=it+1 ! Do i=1,Nlf ! Is this TC located in given gridbox? ! iy=TC(itc)%YY(i) ! This TC(itc) occurred in iy year ! x=TC(itc)%Lon(i); y=TC(itc)%Lat(i); Yes_TC=0 ! Yes_TC=TC_Location (x,y,x1,y1,x2,y2,iy,iyr) ! =0 or 1 if iyr=iy, otherwise =0 !if (Yes_TC>0) exit ! This TC has passed in gridbox in iy year ! tmp(ilon,ilat,it)=tmp(ilon,ilat,it)+Yes_TC ! To be increased if Yes_TC=1 ! Enddo ! Enddo ! Enddo ! Enddo ! Track(:,:,:)=Track(:,:,:)+tmp(:,:,:) ! Update NTC count for TC(itc) !Enddo !-------------------------------------------------------------------- ! For ntc: Call All_NTC Do it=1, Nyear print*,"Year=",it+Yr1-1 Do ilon=1,NLon print"(20I6)",(Track(ilon,ilat,it),ilat=Nlat,1,-1) Enddo !pause Enddo !-------------------------------------------------------------------- Call Characteristics Fname="NTC_Decades_NWP.nc" Call Write_Track_to_NC(Fname) !-------------------------------------------------------------------- ! For number of TCs Track(:,:,:) =0 ! Set to Zero Call All_TC Deallocate (Xall,X50,X60,X70,X80,X90,X00,X01,X75_94,X95_14) Call Characteristics Fname="TC_Decades_NWP.nc" Call Write_Track_to_NC(Fname) !-------------------------------------------------------------------- NGrade=3 ! weak=1, moderate=2, strong=3 Allocate (Track_Grade(1:Nlon,1:Nlat,1:Nyear,NGrade)) Track_Grade(:,:,:,:)=0 Do j=1,NGrade Track(:,:,:)=0 ! Set to Zero Call NTC_Grade (j) Track_Grade(:,:,:,j)=Track(:,:,:) Enddo Fname="NTC_75_14_Grade.nc" Call Char_Grade(fname) !-------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------- Subroutine NTC_Grade (Grade) ! All 6-hrs timeslices ! Read alt NTC with Grade over all tracks and store in Track array Integer :: itc,ilon,ilat,it,i,iyr,Yes_TC Real :: x,y,x1,y1,x2,y2 Integer :: TC_id(NTC), Grade, NTC1 ! "1": weak; "2": moderate; "3": strong; -99: Undefine Integer, Allocatable :: TC_Int(:) NTC1=0; TC_id=0 Do itc=1, NTC Nlf=TC(itc)%Life Allocate (TC_Int(Nlf)) TC_Int(1:Nlf)=TC(itc)%CD(1:Nlf) Yes_TC=TC_Inten(TC_Int,Nlf) if (Yes_TC==Grade) then NTC1=NTC1+1 ! TC co Intensity=Grade TC_id(itc)=1 endif Deallocate (TC_Int) Enddo Do itc=1, NTC ! Search for each TC Nlf=TC(itc)%Life ! A TC's life time (number of observations) if (TC_id(itc)==1) then tmp(:,:,:) = 0 ! Total TC count in all years at each gridbox Do ilon=1, Nlon x1=Lon1+(ilon-1)*Dlon; x2=x1+Dlon Do ilat=1, Nlat y1=Lat1+(ilat-1)*Dlat; y2=y1+Dlat it=0 Do iyr=Yr1, Yr2 ! Search for all period it=it+1 Do i=1,Nlf ! Is this TC located in given gridbox? iy=TC(itc)%YY(i) ! This TC(itc) occurred in iy year x=TC(itc)%Lon(i); y=TC(itc)%Lat(i); Yes_TC=0 Yes_TC=TC_Location (x,y,x1,y1,x2,y2,iy,iyr) ! =0 or 1 if iyr=iy, otherwise =0 !if (Yes_TC>0) exit ! This TC has passed in gridbox in iy year tmp(ilon,ilat,it)=tmp(ilon,ilat,it)+Yes_TC ! To be increased if Yes_TC=1 Enddo Enddo Enddo Enddo Track(:,:,:)=Track(:,:,:)+tmp(:,:,:) ! Update NTC count for TC(itc) endif Enddo End Subroutine NTC_Grade !-------------------------------------------------------------------- Integer Function TC_Inten(CD,Life) ! TC intensity: "1": weak; "2": moderate; "3": strong; -99: Undefine ! 2 : Tropical Depression (TD) ! 3 : Tropical Storm (TS) ! 4 : Severe Tropical Storm (STS) ! 5 : Typhoon (TY) ! 6 : Extra-tropical Cyclone (L) ! 7 : Just entering into the responsible area of Japan Meteorological Agency (JMA) ! 8 : Not used ! 9 : Tropical Cyclone of TS intensity or higher Integer :: Life, i, ITe, Inten Integer :: CD(Life) Inten=-99 Do i=1,Life if (CD(i)==2.OR.CD(i)==3) then ITe=1 else if (CD(i)==4.OR.CD(i)==9) then ITe=2 else if (CD(i)==5) then ITe=3 else ITe=-99 endif if (ITe>Inten) Inten=ITe Enddo TC_Inten=Inten End Function TC_Inten !-------------------------------------------------------------------- Subroutine All_TC ! All TC tracks ! Read all NTC over all tracks and store in Track array Integer :: itc,ilon,ilat,it,i,iyr,Yes_TC Real :: x,y,x1,y1,x2,y2 Do itc=1, NTC ! Search for each TC Nlf=TC(itc)%Life ! A TC's life time (number of observations) tmp(:,:,:) = 0 ! Total TC count in all years at each gridbox Do ilon=1, Nlon x1=Lon1+(ilon-1)*Dlon; x2=x1+Dlon Do ilat=1, Nlat y1=Lat1+(ilat-1)*Dlat; y2=y1+Dlat it=0 Do iyr=Yr1, Yr2 ! Search for all period it=it+1 Do i=1,Nlf ! Is this TC located in given gridbox? iy=TC(itc)%YY(i) ! This TC(itc) occurred in iy year x=TC(itc)%Lon(i); y=TC(itc)%Lat(i); Yes_TC=0 Yes_TC=TC_Location (x,y,x1,y1,x2,y2,iy,iyr) ! =0 or 1 if iyr=iy, otherwise =0 if (Yes_TC>0) exit ! This TC has passed in gridbox in iy year Enddo tmp(ilon,ilat,it)=tmp(ilon,ilat,it)+Yes_TC ! To be increased if Yes_TC=1 !Enddo Enddo Enddo Enddo Track(:,:,:)=Track(:,:,:)+tmp(:,:,:) ! Update NTC count for TC(itc) Enddo End Subroutine All_TC !-------------------------------------------------------------------- Subroutine All_NTC ! All 6-hrs timeslices ! Read all NTC over all tracks and store in Track array Integer :: itc,ilon,ilat,it,i,iyr,Yes_TC Real :: x,y,x1,y1,x2,y2 Do itc=1, NTC ! Search for each TC Nlf=TC(itc)%Life ! A TC's life time (number of observations) tmp(:,:,:) = 0 ! Total TC count in all years at each gridbox Do ilon=1, Nlon x1=Lon1+(ilon-1)*Dlon; x2=x1+Dlon Do ilat=1, Nlat y1=Lat1+(ilat-1)*Dlat; y2=y1+Dlat it=0 Do iyr=Yr1, Yr2 ! Search for all period it=it+1 Do i=1,Nlf ! Is this TC located in given gridbox? iy=TC(itc)%YY(i) ! This TC(itc) occurred in iy year x=TC(itc)%Lon(i); y=TC(itc)%Lat(i); Yes_TC=0 Yes_TC=TC_Location (x,y,x1,y1,x2,y2,iy,iyr) ! =0 or 1 if iyr=iy, otherwise =0 !if (Yes_TC>0) exit ! This TC has passed in gridbox in iy year tmp(ilon,ilat,it)=tmp(ilon,ilat,it)+Yes_TC ! To be increased if Yes_TC=1 Enddo Enddo Enddo Enddo Track(:,:,:)=Track(:,:,:)+tmp(:,:,:) ! Update NTC count for TC(itc) Enddo End Subroutine All_NTC !-------------------------------------------------------------------- Subroutine Characteristics Integer :: ilon,ilat,it,N,N20 Real :: XX(Nyear) N=10 ! A decade N20=20 ! 20 nam Allocate (Xall(Nlon,Nlat), X50(Nlon,Nlat), X60(Nlon,Nlat)) Allocate (X70(Nlon,Nlat), X80(Nlon,Nlat), X90(Nlon,Nlat)) Allocate (X00(Nlon,Nlat), X01(Nlon,Nlat)) Allocate (X75_94(Nlon,Nlat), X95_14(Nlon,Nlat)) Do ilon=1, Nlon Do ilat=1, Nlat XX(1:Nyear)=Track(ilon,ilat,1:Nyear); Xall(ilon,ilat)=Ave(XX,Nyear,Miss) XX(1:N)=Track(ilon,ilat,1:N); X50(ilon,ilat)=Ave(XX,N,Miss) XX(1:N)=Track(ilon,ilat,N+1:2*N); X60(ilon,ilat)=Ave(XX,N,Miss) XX(1:N)=Track(ilon,ilat,2*N+1:3*N); X70(ilon,ilat)=Ave(XX,N,Miss) XX(1:N)=Track(ilon,ilat,3*N+1:4*N); X80(ilon,ilat)=Ave(XX,N,Miss) XX(1:N)=Track(ilon,ilat,4*N+1:5*N); X90(ilon,ilat)=Ave(XX,N,Miss) XX(1:N)=Track(ilon,ilat,5*N+1:6*N); X00(ilon,ilat)=Ave(XX,N,Miss) XX(1:N-1)=Track(ilon,ilat,6*N+1:Nyear); X01(ilon,ilat)=Ave(XX,N-1,Miss) XX(1:N20)=Track(ilon,ilat,25:N20+24); X75_94(ilon,ilat)=Ave(XX,N20,Miss) XX(1:N20)=Track(ilon,ilat,25+N20:24+2*N20); X95_14(ilon,ilat)=Ave(XX,N20,Miss) Enddo Enddo !Deallocate (Xall,X50,X60,X70,X80,X90,X00,X01,X75_94,X95_14) End Subroutine Characteristics !-------------------------------------------------------------------- Subroutine Char_Grade(fname) Integer :: ilon,ilat,igr,N,N20 Real :: XX(Nyear),xgrade(NGrade) Character (Len=*) :: Fname Character (Len=50) :: var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm N20=20 ! 20 nam Allocate (X7594Gr(Nlon,Nlat,NGrade), X9514Gr(Nlon,Nlat,NGrade)) Do igr=1,NGrade Do ilon=1, Nlon Do ilat=1, Nlat XX(1:N20)=Track_Grade(ilon,ilat,25:N20+24,igr); X7594Gr(ilon,ilat,igr)=Ave(XX,N20,Miss) XX(1:N20)=Track_Grade(ilon,ilat,25+N20:24+2*N20,igr); X9514Gr(ilon,ilat,igr)=Ave(XX,N20,Miss) Enddo Enddo Enddo ! Write to file: Allocate (xlon(Nlon), xlat(Nlat)) Do ilon=1, Nlon; xlon(ilon)=(Lon1-Dlon/2)+Dlon*ilon; Enddo Do ilat=1, Nlat; xlat(ilat)=(Lat1-Dlat/2)+Dlat*ilat; Enddo xgrade(1)=1; xgrade(2)=2; xgrade(3)=3 var_nm="m7594"; var_ln="Mean of TC days for period 1975-1994"// & " for TC grades of 1 (weak), 2 (moderate), 3 (strong) TC" var_unit="-"; lon_nm="lon"; lat_nm="lat"; t_nm="grade"; t_unit=" years since 1975-01-01 00:00:00 " Call wr_ncdf_3d (Fname,xlon,xlat,xgrade,Nlon,Nlat,NGrade,X7594Gr, & var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm,Miss) var_nm="m9514"; var_ln="Mean of TC days for period 1995-2014"// & " for TC grades of 1 (weak), 2 (moderate), 3 (strong) TC" var_unit="-"; lon_nm="lon"; lat_nm="lat"; t_nm="grade"; t_unit=" years since 1975-01-01 00:00:00 " Call wr_ncdf_3d_apd (Fname,X9514Gr,var_nm,nlon,nlat,NGrade, & var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) End Subroutine Char_Grade !-------------------------------------------------------------------- Subroutine Write_Track_to_NC(Fname) Character (Len=*) :: Fname Integer :: ilon, ilat, it Character (Len=50) :: var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm Allocate (XTrack(Nlon,Nlat,Nyear), xlon(Nlon), xlat(Nlat), xtime(Nyear)) Do ilon=1, Nlon; xlon(ilon)=(Lon1-Dlon/2)+Dlon*ilon; Enddo Do ilat=1, Nlat; xlat(ilat)=(Lat1-Dlat/2)+Dlat*ilat; Enddo Do it = 1, Nyear; xtime(it)=Yr1+it-1; Enddo XTrack(1:Nlon,1:Nlat,1:Nyear) = Track(1:Nlon,1:Nlat,1:Nyear) var_nm="tc";var_ln="TC days"; var_unit="-"; lon_nm="lon" lat_nm="lat"; t_nm="time"; t_unit=" years since 1951-01-01 00:00:00 " Call wr_ncdf_3d (Fname,xlon,xlat,xtime,Nlon,Nlat,Nyear,XTrack, & var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm,Miss) var_nm="mean"; var_ln="Mean of TC days over 1951-2019" Call wr_ncdf_2d_apd (Fname,Xall,var_nm,Nlon,Nlat, var_ln, var_unit, lon_nm, lat_nm, Miss) var_nm="m50"; var_ln="Mean of TC days over 1951-1960" Call wr_ncdf_2d_apd (Fname,X50,var_nm,Nlon,Nlat, var_ln, var_unit, lon_nm, lat_nm, Miss) var_nm="m60"; var_ln="Mean of TC days over 1961-1970" Call wr_ncdf_2d_apd (Fname,X60,var_nm,Nlon,Nlat, var_ln, var_unit, lon_nm, lat_nm, Miss) var_nm="m70"; var_ln="Mean of TC days over 1971-1980" Call wr_ncdf_2d_apd (Fname,X70,var_nm,Nlon,Nlat, var_ln, var_unit, lon_nm, lat_nm, Miss) var_nm="m80"; var_ln="Mean of TC days over 1981-1990" Call wr_ncdf_2d_apd (Fname,X80,var_nm,Nlon,Nlat, var_ln, var_unit, lon_nm, lat_nm, Miss) var_nm="m90"; var_ln="Mean of TC days over 1991-2000" Call wr_ncdf_2d_apd (Fname,X90,var_nm,Nlon,Nlat, var_ln, var_unit, lon_nm, lat_nm, Miss) var_nm="m00"; var_ln="Mean of TC days over 2001-2010" Call wr_ncdf_2d_apd (Fname,X00,var_nm,Nlon,Nlat, var_ln, var_unit, lon_nm, lat_nm, Miss) var_nm="m01"; var_ln="Mean of TC days over 2011-2019" Call wr_ncdf_2d_apd (Fname,X01,var_nm,Nlon,Nlat, var_ln, var_unit, lon_nm, lat_nm, Miss) var_nm="m7594"; var_ln="Mean of TC days over 1975-1994" Call wr_ncdf_2d_apd (Fname,X75_94,var_nm,Nlon,Nlat, var_ln, var_unit, lon_nm, lat_nm, Miss) var_nm="m9514"; var_ln="Mean of TC days over 1995-2014" Call wr_ncdf_2d_apd (Fname,X95_14,var_nm,Nlon,Nlat, var_ln, var_unit, lon_nm, lat_nm, Miss) Deallocate (XTrack,xlon,xlat,xtime) End Subroutine Write_Track_to_NC !-------------------------------------------------------------------- Integer Function TC_Location (x,y,x1,y1,x2,y2,iyr,Yrs) Real :: x,y,x1,y1,x2,y2 Integer :: iyr,Yrs if (iyr==Yrs) then if (x>x1.and.x<=x2.and.y>y1.and.y<=y2) then !print"(6F10.3)",x1,x,x2,y1,y,y2 TC_Location=1 else TC_Location=0 endif else TC_Location=0 endif End Function TC_Location !-------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Ave.f90" END