Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: IMiss=-99, Nmem=11 Character (Len=5), Parameter :: memdir(Nmem)=(/"mem1","mem2","mem3","mem4","mem5","mem6","mem7","mem8","mem9","mem10","cf"/) Integer, Allocatable :: time(:) Integer :: Ntime, i, LenSt, sys Character (Len=255) :: Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2 Character (Len=255), Allocatable :: Files(:) Character (Len=12) :: date0, timeres, idate0, idate1 Character (Len=1) :: Ch Integer :: yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2 Real*8 :: factor, offset, Missing Real, Allocatable :: Prec(:,:,:), lon(:), lat(:), xtime(:) Real, Allocatable, Dimension (:,:,:,:) :: Mem1,Mem2,Mem3,Mem4,Mem5,Mem6,Mem7,Mem8,Mem9,Mem10,Mem11 ! (nlon,nlat,ntime,nfcs) Integer, Allocatable, Dimension (:) :: date_m1,date_m2,date_m3,date_m4,date_m5,date_m6,date_m7,date_m8,date_m9,date_m10,date_m11 ! Fcs date Integer :: Nlon, Nlat, Nfcs, date_fcs, Nlon0, Nlat0 ! Nfcs = So lan du bao, date_fcs = yyyymmdd cua ngay DB Integer :: imem, ifcs, iy, im, id DirIn="/home/tanpv/Progs/S2S/Data_EC/" !Fname="/home/tanpv/Progs/S2S/Data_EC/cf/raw_tp_cf_2002-06-18.nc" !Call getarg(1, Fname) !Call GetTime(Fname,time,Ntime,yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2,factor, offset, Missing) !Call Read_Fname("Fname_list_1.txt", Nfcs) !Call GetTime (Fname, Ntime, date_fcs, factor,offset,Missing) !j=Ntime/4; print*,"New Nt=",j !j=0 !Do i=1,Ntime, 4 ! if (i>1.AND.i<=Ntime) then ! j=j+1 ! print*, j,i-4,i ! endif !Enddo !i=date_fcs/10000; j=(date_fcs-i*10000)/100; ifcs=mod(date_fcs,100) !print*, i,j,ifcs !stop Do imem=4,4 !Nmem !imem=1,Nmem write (Str,"(A,I2.2,A)") "Fname_list_",imem,".txt" Call Read_Fname (Str,Nfcs) ! Return to Nfcs and Files(Nfcs) Do ifcs=1,Nfcs Fname=trim(DirIn)//trim(memdir(imem))//"/"//trim(Files(ifcs)) Call GetTime (Fname, Nt, date_fcs, factor,offset,Missing) ! print*,date_fcs Call Read_Prec (Fname) ! Return to Ntime (Leadtime) and Prec(Nlon,Nlat,Ntime) if (ifcs==1) then Call Init_mem (imem,Nlon,Nlat,Ntime,Nfcs) iy=date_fcs/10000; im=(date_fcs-i*10000)/100; id=mod(date_fcs,100) write(time_unit,"(A,I4.4,A,I2.2,A,I2.2,A)") "hours since ",iy,"-",im,"-",id," 00 " endif Call Save_to_Arr (imem,ifcs) ! Save to Mem1, Mem2,...,Mem11; date_m1, date_m2,...,date_m11 !print*,"Mem=",imem," ifcs=", ifcs," date_fcs=",date_fcs Enddo !print*,">>>Nlon,Nlat,Ntime,Nfcs=",Nlon,Nlat,Ntime,Nfcs write(Fname1,"(A,I2.2,A)") trim(DirIn)//"Daily/"//"Prec_daily_s2s_ecmwf_",imem,".nc" Call Write_mem_to_file (Fname1) Enddo !Call GetTime (Fname, Ntime, date_fcs, factor,offset,Missing) !print*, Ntime !Do i=1,Ntime, 4 ! print*, i !Enddo ! stop !print*,xtime !print*,yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2 !print*,date_fcs, factor, offset !Call Read_Prec (Fname) !Where (Prec==Missing) Prec=Miss ! unpacked_value = packed_value * scale_factor + add_offset !Where (Prec/=Miss) Prec=Prec * factor + offset !print*,Prec(:,:,5)-Prec(:,:,1) CONTAINS !---------------------------------------------------------------------------------------------- Subroutine Write_mem_to_file (Fname1) ! ,X4d,X1d,Nlon,Nlat,Nt,Nfcs) Character (Len=*) :: Fname1 Real, Allocatable :: xtime(:), X4d(:,:,:,:) !, fcs_date(:), days(:) Integer, Allocatable :: fcs_date(:), days(:) Integer :: it Character (Len=50) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm Allocate (xtime(Ntime)) ! Ntime = Nt/4 Num of leadtime days Allocate (X4d(Nlon,Nlat,Ntime,Nfcs),fcs_date(Nfcs),days(Ntime)) !print*,"Nlon,Nlat,Ntime,Nfcs=",Nlon,Nlat,Ntime,Nfcs Do it=1,Ntime; xtime(it)=it-1; Enddo ! Leadtime=0,1,...,Ntime-1 Do it=1,Ntime; days(it)=it; Enddo Select Case (imem) Case (1) X4d(:,:,:,:) = Mem1(:,:,:,:); fcs_date(:) = date_m1(:) Case (2) X4d(:,:,:,:) = Mem2(:,:,:,:); fcs_date(:) = date_m2(:) Case (3) X4d(:,:,:,:) = Mem3(:,:,:,:); fcs_date(:) = date_m3(:) Case (4) X4d(:,:,:,:) = Mem4(:,:,:,:); fcs_date(:) = date_m4(:) Case (5) X4d(:,:,:,:) = Mem5(:,:,:,:); fcs_date(:) = date_m5(:) Case (6) X4d(:,:,:,:) = Mem6(:,:,:,:); fcs_date(:) = date_m6(:) Case (7) X4d(:,:,:,:) = Mem7(:,:,:,:); fcs_date(:) = date_m7(:) Case (8) X4d(:,:,:,:) = Mem8(:,:,:,:); fcs_date(:) = date_m8(:) Case (9) X4d(:,:,:,:) = Mem9(:,:,:,:); fcs_date(:) = date_m9(:) Case (10) X4d(:,:,:,:) = Mem10(:,:,:,:); fcs_date(:) = date_m10(:) Case (11) X4d(:,:,:,:) = Mem11(:,:,:,:); fcs_date(:) = date_m11(:) End Select !print*,fcs_date t_unit=trim(time_unit); var_ln="Daily precipitation fcs"; var_unit="mm/day"; lon_nm="lon"; lat_nm="lat"; z_nm="day"; t_nm="fcs_date" Call wr_ncdf_4d (Fname1,lon,lat,days,fcs_date,nlon,nlat,Ntime,Nfcs,X4d,"pr", & t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) !lon_nm="day" !Call wr_ncdf_1d (Fname2, days, Nfcs, fcs_date, "f_date", "Forecast date", "yymmdd", lon_nm, Miss) End Subroutine Write_mem_to_file !---------------------------------------------------------------------------------------------- Subroutine Save_to_Arr (imem,ifcs) Integer :: imem,ifcs Select Case (imem) Case (1) Mem1(:,:,:,ifcs) = Prec(:,:,:); date_m1(ifcs)=date_fcs Case (2) Mem2(:,:,:,ifcs) = Prec(:,:,:); date_m2(ifcs)=date_fcs Case (3) Mem3(:,:,:,ifcs) = Prec(:,:,:); date_m3(ifcs)=date_fcs Case (4) Mem4(:,:,:,ifcs) = Prec(:,:,:); date_m4(ifcs)=date_fcs Case (5) Mem5(:,:,:,ifcs) = Prec(:,:,:); date_m5(ifcs)=date_fcs Case (6) Mem6(:,:,:,ifcs) = Prec(:,:,:); date_m6(ifcs)=date_fcs Case (7) Mem7(:,:,:,ifcs) = Prec(:,:,:); date_m7(ifcs)=date_fcs Case (8) Mem8(:,:,:,ifcs) = Prec(:,:,:); date_m8(ifcs)=date_fcs Case (9) Mem9(:,:,:,ifcs) = Prec(:,:,:); date_m9(ifcs)=date_fcs Case (10) Mem10(:,:,:,ifcs) = Prec(:,:,:); date_m10(ifcs)=date_fcs Case (11) Mem11(:,:,:,ifcs) = Prec(:,:,:); date_m11(ifcs)=date_fcs End Select End Subroutine Save_to_Arr !---------------------------------------------------------------------------------------------- Subroutine Init_mem (imem,Nlon,Nlat,Ntime,Nfcs) Integer :: imem,Nlon,Nlat,Ntime,Nfcs Select Case (imem) Case (1) if (.not.Allocated(Mem1)) Allocate (Mem1(Nlon,Nlat,Ntime,Nfcs)); Mem1=Miss if (.not.Allocated(date_m1)) Allocate (date_m1(Nfcs)); date_m1=IMiss Case (2) if (.not.Allocated(Mem2)) Allocate (Mem2(Nlon,Nlat,Ntime,Nfcs)); Mem2=Miss if (.not.Allocated(date_m2)) Allocate (date_m2(Nfcs)); date_m2=IMiss Case (3) if (.not.Allocated(Mem3)) Allocate (Mem3(Nlon,Nlat,Ntime,Nfcs)); Mem3=Miss if (.not.Allocated(date_m3)) Allocate (date_m3(Nfcs)); date_m3=IMiss Case (4) if (.not.Allocated(Mem4)) Allocate (Mem4(Nlon,Nlat,Ntime,Nfcs)); Mem4=Miss if (.not.Allocated(date_m4)) Allocate (date_m4(Nfcs)); date_m4=IMiss Case (5) if (.not.Allocated(Mem5)) Allocate (Mem5(Nlon,Nlat,Ntime,Nfcs)); Mem5=Miss if (.not.Allocated(date_m5)) Allocate (date_m5(Nfcs)); date_m5=IMiss Case (6) if (.not.Allocated(Mem6)) Allocate (Mem6(Nlon,Nlat,Ntime,Nfcs)); Mem6=Miss if (.not.Allocated(date_m6)) Allocate (date_m6(Nfcs)); date_m6=IMiss Case (7) if (.not.Allocated(Mem7)) Allocate (Mem7(Nlon,Nlat,Ntime,Nfcs)); Mem7=Miss if (.not.Allocated(date_m7)) Allocate (date_m7(Nfcs)); date_m7=IMiss Case (8) if (.not.Allocated(Mem8)) Allocate (Mem8(Nlon,Nlat,Ntime,Nfcs)); Mem8=Miss if (.not.Allocated(date_m8)) Allocate (date_m8(Nfcs)); date_m8=IMiss Case (9) if (.not.Allocated(Mem9)) Allocate (Mem9(Nlon,Nlat,Ntime,Nfcs)); Mem9=Miss if (.not.Allocated(date_m9)) Allocate (date_m9(Nfcs)); date_m9=IMiss Case (10) if (.not.Allocated(Mem10)) Allocate (Mem10(Nlon,Nlat,Ntime,Nfcs)); Mem10=Miss if (.not.Allocated(date_m10)) Allocate (date_m10(Nfcs)); date_m10=IMiss Case (11) if (.not.Allocated(Mem11)) Allocate (Mem11(Nlon,Nlat,Ntime,Nfcs)); Mem11=Miss if (.not.Allocated(date_m11)) Allocate (date_m11(Nfcs)); date_m11=IMiss End Select End Subroutine Init_mem !---------------------------------------------------------------------------------------------- Subroutine Read_Fname(Fname, Nfile) Character (Len=*) :: Fname Character (Len=100) :: Str Integer :: Nfile, i, IO Open (1,file=Fname,status="old") IO=0; i=0 Do while (IO==0) i=i+1 Read (1,"(A)",IOStat=IO) Str if (IO/=0) then; i=i-1; exit; endif !print*,i," ",trim(Str) Enddo Nfile=i if (Allocated(Files)) Deallocate(Files) Allocate (Files(Nfile)) Rewind(1) Do i=1,Nfile Read (1,"(A)") Files(i) Enddo !print*,Nfile End Subroutine Read_Fname !---------------------------------------------------------------------------------------------- !Subroutine GetTime (fname, time, nt, yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2, scale_factor, add_offset,missing_value) Subroutine GetTime (fname, nt, date_fcs, scale_factor, add_offset,missing_value) Character (Len=*) :: fname Integer, Intent (Out) :: nt, date_fcs Integer :: ncid, id_t_dim, varid, xtype, len, attnum Integer, Allocatable :: time(:) Logical :: status Character (Len=13) :: date0, timeres, idate0, idate1 Character (Len=1) :: Ch Character (Len=100) :: t_unit Integer :: sys,yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2, LenSt Real*8 :: scale_factor, add_offset, missing_value status = nf90_open(fname, nf90_NoWrite, ncid) ! Time dimension: status = nf90_inq_dimid(ncid,"time", id_t_dim) status = nf90_inquire_dimension(ncid, id_t_dim,len=Nt) Allocate (Time(Nt)) status = nf90_inq_varid(ncid,"time", varid) status = nf90_get_var(ncid, varid, Time) write(idate0,"(A,I10)") " ",time(1); write(idate1,"(A,I10)") " ",time(Nt) Deallocate (Time) status = nf90_inquire_attribute(ncid, varid, "time", xtype, len, attnum) status = nf90_get_att(ncid, varid, "units", t_unit) LenSt=Len_trim(t_unit) Do i=1,LenSt Ch=t_unit(i:i) if (Ch=="1".OR.Ch=="2") exit Enddo date0=t_unit(i:i+12) timeres=" "//t_unit(1:4) sys = system ('rm -rf days.txt') !sys = system ('date -u -d "'//trim(date0)//' '//trim(idate0)//trim(timeres)//' " +%Y | cat > days.txt') !sys = system ('date -u -d "'//trim(date0)//' '//trim(idate0)//trim(timeres)//' " +%m | cat >> days.txt') !sys = system ('date -u -d "'//trim(date0)//' '//trim(idate0)//trim(timeres)//' " +%d | cat >> days.txt') !sys = system ('date -u -d "'//trim(date0)//' '//trim(idate0)//trim(timeres)//' " +%H | cat >> days.txt') !sys = system ('date -u -d "'//trim(date0)//' '//trim(idate1)//trim(timeres)//' " +%Y | cat >> days.txt') !sys = system ('date -u -d "'//trim(date0)//' '//trim(idate1)//trim(timeres)//' " +%m | cat >> days.txt') !sys = system ('date -u -d "'//trim(date0)//' '//trim(idate1)//trim(timeres)//' " +%d | cat >> days.txt') !sys = system ('date -u -d "'//trim(date0)//' '//trim(idate1)//trim(timeres)//' " +%H | cat >> days.txt') sys = system ('date -u -d "'//trim(date0)//' '//trim(idate0)//trim(timeres)//' " +%Y%m%d | cat > days.txt') Open (1,file="days.txt") Read (1,*) date_fcs !Read (1,*) yy1; Read (1,*) mm1; Read (1,*) dd1; Read (1,*) hh1; !Read (1,*) yy2; Read (1,*) mm2; Read (1,*) dd2; Read (1,*) hh2; Close(1) ! Lon dimension: status = nf90_inq_varid(ncid,"tp", varid) status = nf90_get_att(ncid, varid, "scale_factor", scale_factor) status = nf90_get_att(ncid, varid, "add_offset", add_offset) status = nf90_get_att(ncid, varid, "missing_value", missing_value) status = nf90_close(ncid) End Subroutine GetTime !---------------------------------------------------------------------------------------------- Subroutine Read_Prec (Fname) Character (Len=*) :: Fname Real, Allocatable :: tmp(:,:,:) Integer :: Nt, i, j, ilon,ilat if (Allocated(lon)) Deallocate (lon) if (Allocated(lat)) Deallocate (lat) if (Allocated(xtime)) Deallocate (xtime) if (Allocated(Prec)) Deallocate (Prec) Call rd_ncdf_3d (Fname,lon,lat,xtime,nlon,nlat,Nt, Prec,"tp") Where (Prec==Missing) Prec=Miss Where (Prec/=Miss) Prec=Prec * factor + offset ! unpacked_value = packed_value * scale_factor + add_offset Ntime = Nt/4 ! Number of days, 4 obs per day Allocate (tmp(nlon,nlat,ntime)) tmp=Miss j=0 Do i=1,Nt,4 if (i>1.AND.i<=Nt) then j=j+1 Do ilon=1,nlon Do ilat=1,nlat if (Prec(ilon,ilat,i)/=Miss.AND.Prec(ilon,ilat,i-4)/=Miss) then tmp(ilon,ilat,j)=Prec(ilon,ilat,i)-Prec(ilon,ilat,i-4) if (tmp(ilon,ilat,j)<0) tmp(ilon,ilat,j)=0 endif !print"(3I5,4F10.1)",j,i-4,i,Prec(ilon,ilat,i-4),Prec(ilon,ilat,i),tmp(ilon,ilat,j) !if (tmp(ilon,ilat,j)<0) print"(3I5,4F10.3)",j,i-4,i,Prec(ilon,ilat,i-4),Prec(ilon,ilat,i),tmp(ilon,ilat,j) Enddo Enddo !pause endif Enddo Deallocate(Prec) Allocate (Prec(nlon,nlat,ntime)) Prec(:,:,:) = tmp(:,:,:) Deallocate (tmp) End Subroutine Read_Prec !---------------------------------------------------------------------------------- Subroutine wr_ncdf_4d (filenm,lon,lat,z,t,nlon,nlat,nz,nt,var,var_nm, & t_unit, var_long_name, var_unit, lon_nm, lat_nm, z_nm, t_nm, flag) character (len = *), intent (in) :: filenm, var_nm, t_unit, var_long_name, var_unit character (len = *) :: lon_nm, lat_nm, z_nm, t_nm integer, parameter :: NDIMS = 4 integer :: ncid, var_id, dims_id(NDIMS) integer :: lon_id, lat_id, z_id, t_id integer, intent (in) :: nlon, nlat, nz, nt real, dimension (nlon,nlat,nz,nt), intent (in) :: var real, intent (in) :: lon(nlon), lat(nlat), flag !, z(nz), t(nt), flag integer, intent (in) :: z(nz), t(nt) print*," Writing netCDF file: ",filenm,"..." call hulk ( nf90_create (filenm, nf90_write, ncid) ) call hulk ( nf90_def_dim (ncid, lon_nm, nlon, dims_id(1)) ) call hulk ( nf90_def_dim (ncid, lat_nm, nlat, dims_id(2)) ) call hulk ( nf90_def_dim (ncid, z_nm, nz, dims_id(3)) ) call hulk ( nf90_def_dim (ncid, t_nm, nt, dims_id(4)) ) call hulk ( nf90_def_var (ncid,lon_nm,nf90_float,dims_id(1),lon_id) ) call hulk ( nf90_def_var (ncid,lat_nm,nf90_float,dims_id(2),lat_id) ) call hulk ( nf90_def_var (ncid,z_nm,NF90_INT,dims_id(3),z_id) ) call hulk ( nf90_def_var (ncid,t_nm,NF90_INT,dims_id(4),t_id) ) call hulk ( nf90_def_var (ncid,var_nm,nf90_float,dims_id(1:4),var_id) ) call hulk ( nf90_put_att (ncid,lon_id,"long_name","Longitude") ) call hulk ( nf90_put_att (ncid,lon_id,"units","degree_east") ) call hulk ( nf90_put_att (ncid,lon_id,"axis","X") ) call hulk ( nf90_put_att (ncid,lat_id,"long_name","Latitude") ) call hulk ( nf90_put_att (ncid,lat_id,"units","degree_north") ) call hulk ( nf90_put_att (ncid,lat_id,"axis","Y") ) call hulk ( nf90_put_att (ncid,z_id, "long_name","z") ) call hulk ( nf90_put_att (ncid,z_id, "units","Lev") ) call hulk ( nf90_put_att (ncid,z_id, "axis","Z") ) call hulk ( nf90_put_att (ncid,t_id,"long_name","Time") ) call hulk ( nf90_put_att (ncid,t_id,"units",trim(t_unit)//" 00:00:00") ) call hulk ( nf90_put_att (ncid,t_id,"calendar","proleptic_gregorian") ) call hulk ( nf90_put_att (ncid,t_id,"axis","T") ) call hulk ( nf90_put_att (ncid,var_id,"long_name",trim(var_long_name))) call hulk ( nf90_put_att (ncid,var_id,"units", trim(var_unit))) call hulk ( nf90_put_att (ncid,var_id,"missing_value",flag) ) call hulk ( nf90_enddef (ncid) ) call hulk ( nf90_put_var (ncid,lon_id,lon) ) call hulk ( nf90_put_var (ncid,lat_id,lat) ) call hulk ( nf90_put_var (ncid,z_id, z) ) call hulk ( nf90_put_var (ncid,t_id, t) ) call hulk ( nf90_put_var (ncid,var_id,var) ) call hulk ( nf90_close (ncid) ) print*," Completed writing file ",filenm End Subroutine wr_ncdf_4d !---------------------------------------------------------------------------------------------- ! <-~~~=| Nf90 checking errors ... SUBROUTINE hulk (stat) use netcdf integer :: stat if ( stat .ne. nf90_noerr) then print*, trim(nf90_strerror(stat)) stop endif END SUBROUTINE hulk !---------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" !Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" END