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(:,:,:) Real, Allocatable :: XTrack(:,:,:), xlon(:), xlat(:), xtime(:) !-------------------------------------------------------------------- Real, Allocatable, Dimension (:,:) :: Xall, X50, X60, X70, X80, X90, X00, X01 !-------------------------------------------------------------------- Character (Len=100) :: Fname, St Integer :: itc, i, io, ig, iid, k1, k2, k3, k4, k5, k6 Integer :: nyy, nmm, ndd, nhh, NTC, Nlf 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 !-------------------------------------------------------------------- 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_jma.nc" Call Write_Track_to_NC(Fname) !-------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------- Subroutine Characteristics Integer :: ilon,ilat,it,N Real :: XX(Nyear) N=10 ! A decade 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)) 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) Enddo Enddo End Subroutine Characteristics !-------------------------------------------------------------------- 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) 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