Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: Nmon=12, IMiss=-99, NLead=46 !, Lw=5, Nws=NLead-Lw+1 ! Nws=So tuan; Lw = Do dai tuan Integer, Parameter :: Nfile = 46 ! Integer, Parameter :: Nvar=3 ! Pr, Tx, Tm ! Character (Len=4), Parameter :: yvarn(Nvar)=(/"prec","tmax","tmin"/) ! Integer, Parameter :: Nseason=4 ! 4 seasons ! Character (Len=5), Parameter :: Season(Nseason)=(/"12-02","03-05","06-08","09-11"/) !Integer, Parameter :: Nfc_per_year(Nseason)=(/ 25 , 26 , 27 , 26/) !Integer, Parameter :: Yr1=2003, Yr2=2022, Nyrs=(Yr2-Yr1+1) !, Nspl=Nyrs*Nfcs ! 2 thang x 9 lan x 20 nam ! Integer :: Nfc_per_year(Nseason) ! Obtained from Namelist Integer :: Yr1, Yr2, Nyrs, Nspl, Lw, Nws Integer, Allocatable :: time(:) Integer :: LenSt, sys Character (Len=255) :: Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2 Character (Len=255) :: DirInEC, DirInObs, DirInStn Character (Len=12) :: date0, timeres, idate0, idate1 ! Character (Len=5),Allocatable :: hdate(:) ! (Nfcs) ! forecast dates Character (Len=5) :: hdate(Nfile) Integer :: yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2 Real*8 :: factor, offset, Missing Real, Allocatable :: lon(:), lat(:), Regs(:,:) Integer, Allocatable :: lead(:), fcsdate(:,:) !---------------------------------------------------------------------------------------------- Real :: Pr_Fcs(Nfile,NLead) ! Nfile: Initial date Real, Allocatable :: EC(:,:,:,:),OBS(:,:,:,:) ! (Nlon,Nlat,Nspl,NLead) Real, Allocatable :: EC_bc(:,:,:,:) ! (Nlon,Nlat,Nspl,NLead) Real, Allocatable :: EC_w(:,:,:,:),Obs_w(:,:,:,:) ! (Nlon,Nlat,Nspl,Nws) - Tong/Trung binh tuan Real, Allocatable :: EC_Org_w(:,:,:,:) Real, Allocatable, Dimension(:,:,:) :: TP_cli, Obs_cli ! (Nlon,Nlat,NLead) - Clim over all forecasts Real, Allocatable, Dimension(:,:,:,:) :: LH, SH, OLR, TX, TM !---------------------------------------------------------------------------------------------- Integer :: ireg, ilead Integer :: Nlayer ! Num of layers including Input, output and hidden laysers Integer :: Nepoch ! Num of epochs double precision :: learn_rate ! Learning_Rate Integer, Allocatable :: Layer_size(:) ! (Nlayer) !---------------------------------------------------------------------------------------------- Real, Allocatable :: xlon(:), xlat(:), xtime(:) Real, Allocatable :: X4d(:,:,:,:), X3d(:,:,:) Integer :: Nlon, Nlat, date_fcs, Nlon0, Nlat0 ! Nfcs = So lan du bao, date_fcs = yyyymmdd cua ngay DB Integer :: Nlo, Nla, Ntime ! Nday from 2000-01-01 to 2019-12-31 Integer :: Nfcs ! Number of forecast over all years !---------------------------------------------------------------------------------------------- ! Verification Scores: !Real, Allocatable :: R_ME(:,:,:), R_MAE(:,:,:), COR(:,:,:), Prb(:,:,:) ! (Lon,Lat,LT) !Real, Allocatable :: MEe(:,:,:), MAEe(:,:,:), CORe(:,:,:), Prbe(:,:,:) ! (Lon,Lat,LT) !---------------------------------------------------------------------------------------------- !Real, Allocatable :: RMEreg(:), RMAEreg(:), CORreg(:) ! (LT) Integer :: imem, ifcs, iy, im, id, NLd, iobs !-------------------------------------------------------------------------------- ! ANN: character (len=250) :: input_file, output_file, save_file character (len=250) :: Dir_home, Dir_Mask, Dir_ANN_Out, Dir_Cli, Dir_ANN_Save, Dir_Veri, Dir_Corrected logical :: training_mode ! Use for calculating scaling factors Integer :: Nsample, Nepoc, Nsample1 !---------------------------------------------------------------------------------------------- Open (1,file="hdate_Oct_07.txt",status="old") Do i=1,Nfile Read (1,"(A5)") hdate(i) Enddo Call Read_data Call Heatmap Call Write_Heatmap ("Head_map.nc") CONTAINS !---------------------------------------------------------------------------------------------- Subroutine Heatmap Integer :: i,j,k Real :: X(Nlon,Nlat,NLead), Y(NLead) Pr_Fcs = Miss j=0 Do i=Nfile,1,-1 X(:,:,:) = EC(:,:,:,i) Call Areal_Mean (X,Nlon,Nlat,NLead,Y) j=j+1 Pr_Fcs(j,j:NLead) = Y(1:(NLead-j+1)) ! Pr_Fcs(j,1:NLead) = Y(1:NLead) ! print*,j,NLead,NLead-j+1; pause Enddo End Subroutine Heatmap !---------------------------------------------------------------------------------------------- Subroutine Write_Heatmap (Fname) Character (Len=*) :: Fname Real :: Lead(NLead), Fcs(Nfile) Integer :: i,j,k Do i=1,NLead; Lead(i)=i-1; Enddo Do i=1,Nfile; Fcs (i)=i+2; Enddo Call wr_ncdf_2d (Fname,Lead,Fcs,NLead,Nfile,Pr_Fcs,"tp5day","5-day moving sum of reainfall","mm","leadtime","ini_date",Miss) End Subroutine Write_Heatmap !---------------------------------------------------------------------------------------------- Subroutine Write_BC_Outputs (Fname) Character (Len=*) :: Fname Character (Len=255) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, var_nm Real, Allocatable :: lead(:),sampl(:) Integer :: ilon,ilat,imon,ilead,imem,ifcs,im,iobs, i,j,iyr, sys1,sys2 Integer :: dfc(Nfcs), mfc(Nfcs) ! forecast date and month - hdate(Nfcs) - Yr1=2001, Yr2=2020 Integer :: days_diff, jd_ref, jd_target t_unit="days since 1980-01-01 "; var_unit="mm/week"; lon_nm="lon"; lat_nm="lat"; z_nm="fcs_time"; t_nm="week" if (Allocated (lead)) Deallocate (lead); if (Allocated (sampl)) Deallocate (sampl) Allocate (lead(Nlead),sampl(Nspl)) Do ilead=1,Nlead; lead (ilead)=ilead; Enddo Do ispl =1,Nspl; sampl(ilead)=ispl; Enddo Do i=1,Nfcs read(hdate(i),"(I2)") mfc(i); read(hdate(i),"(3x,I2)") dfc(i) Enddo i=0 Do iyr=Yr1,Yr2 Do j=1,Nfcs i=i+1; jd_ref = julian_day(1980, 1, 1) jd_target = julian_day(iyr, mfc(j), dfc(j)) sampl(i)=jd_target - jd_ref !iyr*10000+mfc(j)*100+dfc(j) Enddo Enddo var_ln="Moving weekly sum with a 1-day step of Observation precipitation (CHIRPS)" Call wr_ncdf_4d (Fname,Lon,Lat,sampl,lead,Nlon,Nlat,Nspl,Nlead,Obs_w,"obs", & t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) var_ln="Moving weekly sum with a 1-day step of S2S corrected precipitation of ecmwf forecast" Call wr_ncdf_4d_apd (Fname,EC_w,"fcs", Nlon,Nlat,Nspl,Nlead, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) if (.not.training_mode) then var_ln="Moving weekly sum with a 1-day step of S2S un-corrected precipitation of ecmwf forecast" Call wr_ncdf_4d_apd (Fname,EC_Org_w,"org", Nlon,Nlat,Nspl,Nlead, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) endif End Subroutine Write_BC_Outputs !---------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------- Subroutine Areal_Mean (X,Nlon,Nlat,NLead,Y) Integer :: Nlon,Nlat,NLead, ireg Real :: X(Nlon,Nlat,NLead), Y(NLead) ! EC(ilon,ilat,1:NLead,ifcs,imem) Integer :: ilon,ilat, ilead Real :: tmp, ntmp Y(:)=Miss Do ilead=1,NLead tmp=0; ntmp=0 Do ilon=1,Nlon Do ilat=1,Nlat if (X(ilon,ilat,ilead)/=Miss) then ntmp=ntmp+1; tmp=tmp+X(ilon,ilat,ilead) endif Enddo Enddo if (ntmp>0) Y(ilead)=tmp/ntmp Enddo End Subroutine Areal_Mean !---------------------------------------------------------------------------------------------- Subroutine Read_data Real, Parameter :: ec_miss=-32767.0, obs_miss=-9999.0 Integer :: imem, ifcs Real, Allocatable :: xlon(:), xlat(:), xtime(:), X3d(:,:,:) Integer :: Nlo,Nla,NLt,i,ifc,ispl,iyr ! Real, Allocatable, Dimension(:,:,:,:) :: LH, SH, OLR, TX, TM ! Read precipitation forecast data: if (Allocated (EC)) Deallocate (EC); !Allocate (EC(Nlon,Nlat,Nlead,Nfile)) Do ifc=1,Nfile !write(Fname,"(A)") "tp_cf_2025-"//hdate(ifc)//".nc" write(Fname,"(A)") "Prec_5days_2025-"//hdate(ifc)//".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"tp") !"tp5day") !"tp") if (ifc==1) then Nlon=Nlo; Nlat=Nla Allocate (EC(Nlon,Nlat,Nlead,Nfile)) endif Where (X3d==ec_miss) X3d=Miss EC(:,:,:,ifc) = X3d(:,:,:) Deallocate (xlon,xlat,xtime,X3d) Enddo End Subroutine Read_data !---------------------------------------------------------------------------------------------- Subroutine Calc_weekly_data ! Moving Sum ! Real, Allocatable :: EC_w(:,:,:,:),Obs_w(:,:,:,:) ! (Nlon,Nlat,Nspl,Nws) - Tong/Trung binh tuan Integer :: ilon,ilat,ilead,ispl Real :: X(Nlead), Y(Nlead) if (Allocated (EC_w)) Deallocate (EC_w); if (Allocated (Obs_w)) Deallocate (Obs_w) Allocate (EC_w(Nlon,Nlat,Nspl,Nlead),Obs_w(Nlon,Nlat,Nspl,Nlead)) EC_w = Miss; Obs_w = Miss if (.not.training_mode) then if (Allocated (EC_Org_w)) Deallocate (EC_Org_w) Allocate (EC_Org_w(Nlon,Nlat,Nspl,Nlead)) EC_Org_w = Miss endif Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do ispl=1,Nspl X(1:Nlead)=EC_bc(ilon,ilat,ispl,1:Nlead) Call Moving_Ave (X,Nlead,Lw,Y,Miss) EC_w(ilon,ilat,ispl,1:Nlead)=Y(1:Nlead) !print"(22F6.1)",EC_w(ilon,ilat,ispl,10:30) X(1:Nlead)=OBS(ilon,ilat,ispl,1:Nlead) Call Moving_Ave (X,Nlead,Lw,Y,Miss) Obs_w(ilon,ilat,ispl,1:Nlead)=Y(1:Nlead) if (.not.training_mode) then X(1:Nlead)=EC(ilon,ilat,ispl,1:Nlead) Call Moving_Ave (X,Nlead,Lw,Y,Miss) EC_Org_w(ilon,ilat,ispl,1:Nlead)=Y(1:Nlead) endif Enddo endif Enddo Enddo End Subroutine Calc_weekly_data !---------------------------------------------------------------------------------------------- Subroutine Moving_Ave (X,N,m,Y,Miss) Integer :: N,m, i,j,k,m2 ! m must be an odd number Real :: X(N), Y(N), A(m), tmp, Miss Y(:) = Miss; m2=m/2 Do j=1,N-m+1 k = j+m-1 A(1:m) = X(j:k) tmp = Ave(A,m,Miss) if (tmp/=Miss) then Y(j+m2)=tmp*m ! Tong truot Gan cho gia tri giua else Y(j+m2)=Miss endif Enddo End Subroutine Moving_Ave !--------------------------------------------------------------------------------- Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Ave.f90" ! Include "/home/tanpv/Mylibs/Sub_neural_2025.f90" Include "/home/tanpv/Mylibs/Sub_julian_day.f90" END !---------------------------------------------------------------------------