Use NetCDF Use neural type(network) :: Net 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 :: Yr1=2001, Yr2=2020, Nfcs=18, Nyrs=(Yr2-Yr1+1) ! 2 thang x 9 lan x 20 nam Integer, Allocatable :: time(:) Character (Len=5) :: hdate(Nfcs),yvarn !---------------------------------------------------------------------------------------------- Real, Allocatable :: EC(:,:,:,:),OBS(:,:,:,:) ! (Nlon,Nlat,Nspl,NLead) Real, Allocatable :: EC_bc(:,:,:,:) ! (Nlon,Nlat,Nspl,NLead) Real, Allocatable, Dimension(:,:,:,:) :: LH, SH, OLR, TX, TM !---------------------------------------------------------------------------------------------- Real, Allocatable :: X_inp(:,:), Y_inp(:), Yest(:,:), Yhq(:,:) !----------------------------------------------------------------------------- Real, Allocatable :: Regs(:,:), Lon(:), Lat(:), xlon(:), xlat(:) Real, Allocatable :: Ymax(:,:,:), Ymin(:,:,:) ! Nlon,Nlat,Nlead Integer :: Nlon, Nlat, Ntime, Nt !----------------------------------------------------------------------------- Real, Allocatable, Dimension (:,:,:) :: Cor_a,Cor_h,PC_a,PC_h,POD_a,POD_h,SR_a,SR_h,POFD_a,POFD_h,FAR_a,FAR_h Real, Allocatable, Dimension (:,:,:) :: Prb_a, Prb_h !----------------------------------------------------------------------------- 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) Logical :: training_mode double precision, Allocatable :: tmp(:) Integer :: ilr, inode, itest, Ntest, max_node, icase ! max_node: Maximum values among Mxvar, Myvar and number of nodes in Layer_size !----------------------------------------------------------------------------- ! Ann_hai: Integer :: Nsample, Nepoc !----------------------------------------------------------------------------- Character (Len=255) :: Dir_ANN_Save, save_file, Inp_Dir, Out_Dir, Mask_File, Fname Character (Len=255) :: DirInEC,DirInObs,Dir_Cli,Dir_Corrected Namelist /ann_info/ Nlayer,learn_rate,Nepoch,training_mode,yvarn Namelist /layer_info/ Layer_size ! Layer_size(1) will be Mxvar, Layer_size(Nlayer) will be Myvar !----------------------------------------------------------------------------- Open (1, file="Namelist_Dr_fcs_Ann_new.txt",Status="old") Read (1, ann_info) ! Nlayer,learn_rate,Nepoch Allocate (Layer_size(Nlayer)) Read (1, layer_info) ! Close (1) print*,Layer_size !----------------------------------------------------------------------------- Open (1,file="hdate.txt", status="old") Do i=1,Nfcs read (1,*) hdate(i) Enddo Close (1) !----------------------------------------------------------------------------- DirInEC ="/home/cordex_sea/Chi/TayNguyen/Store/" DirInObs ="/home/cordex_sea/Chi/TayNguyen/CHIRPS/" Dir_ANN_Save ="/home/cordex_sea/Chi/TayNguyen/ANN_Save/" Dir_Cli ="/home/cordex_sea/Chi/TayNguyen/Cli_Out/" Dir_Corrected ="/home/cordex_sea/Chi/Res/" !----------------------------------------------------------------------------- Fname="/home/tanpv/Data/S2S/2024/Mask_R6_0.125_S2S.nc" Call rd_ncdf_2d (trim(Fname),lon,lat,Nlon,Nlat,Regs,"r6") ! 7 sub-regs mask !----------------------------------------------------------------------------- Nspl=Nyrs*Nfcs Ntime = Nspl Ncase = Npoint () ! So lan tinh = So diem luoi x N_lead_time !Ncase=2300 Call Read_data ! Read all data of TP, LH, SH, OLR, TX, TM from entire domain (no mask for Central Highland), and CHIRPs Mx=22 ! So bien doc lap Mvars=Mx Mxvar=Mx; Myvar=1 !----------------------------------------------------------------------------- inode = Maxval(Layer_size(1:Nlayer)); max_node = Max(inode,Mxvar,Myvar) Allocate (Ymax(Nlon,Nlat,Nlead), Ymin(Nlon,Nlat,Nlead)) if (training_mode) then Allocate (Mnet_layers(Nlayer)) Allocate (Mnet_weight(Ncase,Nlayer-1,max_node,max_node)) Allocate (Mnet_biases(Ncase,Nlayer-1,max_node)) else Fname=trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//"_Ann_new.txt" Call Read_Ann_info (Fname) endif !----------------------------------------------------------------------------- Allocate (EC_bc(Nlon,Nlat,Nspl,Nlead)) icase = 0 Do ilon=1,Nlon !22,31 !1,Nlon Do ilat=1,Nlat !21,25 !18,27 !1,Nlat if (Regs(ilon,ilat)>0) then ! R6 Do ilead=1,Nlead icase = icase+1 print*,icase,"/",Ncase Call Make_data (ilon,ilat,ilead) Call Calc_ANN_New (training_mode,ilead,ilon,ilat,icase) Enddo endif Enddo Enddo if (training_mode) then Fname=trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//"_Ann_new.txt" Call Write_Ann_info (Fname) endif if (training_mode) then ! Output filename Fname=trim(Dir_Corrected)//"S2S_ecmwf_tp_Ann_new_training.nc" else Fname=trim(Dir_Corrected)//"S2S_ecmwf_tp_Ann_new_testing.nc" endif Call Write_BC_Outputs (Fname) !----------------------------------------------------------------------------- CONTAINS !---------------------------------------------------------------------------------------------- Integer Function Npoint () Integer :: ilon,ilat,ilead, i i=0 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then Do ilead=1,Nlead i=i+1 Enddo endif Enddo Enddo Npoint=i End Function Npoint !---------------------------------------------------------------------------------------------- Subroutine Make_data (ilon,ilat,ilead) Integer :: ilon,ilat,ilead,it,i,j ! OBS, EC, LH, SH, OLR, TX, TM Nsample1=Nspl ! Dung luong mau if (Allocated(X_inp)) Deallocate (X_inp) if (Allocated(Y_inp)) Deallocate (Y_inp) Allocate (X_inp(Nsample1,Mx),Y_inp(Nsample1)) j = 3; ! Vi tri diem luoi cach tam o luoi hien tai Do i=1,Nspl Y_inp(i) = OBS(ilon,ilat,i,ilead) X_inp(i,1) = EC (ilon,ilat,i,ilead) X_inp(i,2) = LH (ilon,ilat,i,ilead) X_inp(i,3) = SH (ilon,ilat,i,ilead) X_inp(i,4) = OLR(ilon,ilat,i,ilead) X_inp(i,5) = TX (ilon,ilat,i,ilead) X_inp(i,6) = TM (ilon,ilat,i,ilead) X_inp(i,7) = EC (ilon-j,ilat,i,ilead) X_inp(i,8) = EC (ilon+j,ilat,i,ilead) X_inp(i,9) = EC (ilon,ilat-j,i,ilead) X_inp(i,10) = EC (ilon,ilat+j,i,ilead) X_inp(i,11) = LH (ilon-j,ilat,i,ilead) X_inp(i,12) = LH (ilon+j,ilat,i,ilead) X_inp(i,13) = LH (ilon,ilat-j,i,ilead) X_inp(i,14) = LH (ilon,ilat+j,i,ilead) X_inp(i,15) = SH (ilon-j,ilat,i,ilead) X_inp(i,16) = SH (ilon+j,ilat,i,ilead) X_inp(i,17) = SH (ilon,ilat-j,i,ilead) X_inp(i,18) = SH (ilon,ilat+j,i,ilead) X_inp(i,19) = OLR(ilon-j,ilat,i,ilead) X_inp(i,20) = OLR(ilon+j,ilat,i,ilead) X_inp(i,21) = OLR(ilon,ilat-j,i,ilead) X_inp(i,22) = OLR(ilon,ilat+j,i,ilead) Enddo Nsample1=Nspl End Subroutine Make_data !-------------------------------------------------------------------------------- Subroutine Calc_ANN_New (training_mode,ilead,ilon,ilat,icase) Integer :: icase,ilead,ilon,ilat,ispl,N_spl, i,j,k,icount Real :: XX(Ntime,Mxvar),YY(Ntime),xtmp(Mxvar), ytmp Double precision :: x1d (Mxvar), y1d (Myvar) Logical :: training_mode,OK Integer, Allocatable :: Check(:) if (Allocated (Check)) Deallocate (Check) Allocate (Check(Ntime)) Do i=1,Ntime YY(i)=Y_inp(i) Do j=1,Mxvar XX(i,j)=X_inp(i,j) Enddo Enddo ! icount=0; Check(:)=0 Do i=1,Ntime ytmp=YY(i); xtmp(1:Mxvar)=XX(i,1:Mxvar) OK=Valid_Ar1(xtmp,Mxvar,Miss).AND.(ytmp/=Miss) if (OK) then icount=icount+1; Check(i)=1 endif Enddo N_spl = icount if (Allocated (Out_Y)) Deallocate (Out_Y) if (Allocated (Inp_X)) Deallocate (Inp_X) Allocate (Inp_X(N_spl,Mxvar),Out_Y(N_spl,Myvar)) ! To be used by subroutine neural_train_network_file_TP ispl=0 Do i=1,Ntime if (Check(i)==1) then ispl=ispl+1 Out_Y(ispl,1) = YY(i) Inp_X(ispl,1:Mxvar) = XX(i,1:Mxvar) endif Enddo !-------------------------------------------------------- ! Normalize: (for subroutine neural_train_network_file_TP) Ymax(ilon,ilat,ilead) = Maxval(Out_Y(1:N_spl,1)) Ymin(ilon,ilat,ilead) = Minval(Out_Y(1:N_spl,1)) Do i=1,N_spl Out_Y(i,1) = (Out_Y(i,1)-Ymin(ilon,ilat,ilead))/(Ymax(ilon,ilat,ilead)-Ymin(ilon,ilat,ilead)) Enddo !-------------------------------------------------------- Layer_size(1) = Mxvar; Layer_size(Nlayer) = Myvar ! Update Num of nodes for input (X) and output (Y) layers call neural_init_network(Net, Layer_size) if (training_mode) then call neural_train_network_file_TP(Net, Nepoch, learn_rate) call neural_save_network_TP (icase, net) else call neural_load_network_TP(icase,Nlayer,net) endif if (Allocated (Yest)) Deallocate (Yest) Allocate (Yest(Ntime,Myvar)); Yest = Miss ispl=0 Do i=1,Ntime if (Check(i)==1) then ispl=ispl+1 x1d(1:Mxvar) = Inp_X(ispl,1:Mxvar) y1d(:) = neural_use_network(Net,x1d) Do j=1,Myvar Yest(i,j) = y1d(j)*(Ymax(ilon,ilat,ilead)-Ymin(ilon,ilat,ilead)) + Ymin(ilon,ilat,ilead) Enddo !print"(A,6F8.4)","Yest,Y_inp: ",Yest(i,:), YY(i) endif Enddo EC_bc(ilon,ilat,1:Nspl,ilead) = Yest(1:Ntime,1) End Subroutine Calc_ANN_New !-------------------------------------------------------------------------------- Subroutine Write_Ann_info (Fname) Character (Len=*) :: Fname Open (1,file=Fname,form="unformatted", status="unknown") write (1) Nlayer, max_node, Ncase write (1) Mnet_layers write (1) Mnet_weight write (1) Mnet_biases End Subroutine Write_Ann_info !-------------------------------------------------------------------------------- Subroutine Read_Ann_info (Fname) Character (Len=*) :: Fname Open (1,file=Fname,form="unformatted", status="old") read (1) Nlayer, max_node, Ncase Allocate (Mnet_layers(Nlayer)) Allocate (Mnet_weight(Ncase,Nlayer-1,max_node,max_node)) Allocate (Mnet_biases(Ncase,Nlayer-1,max_node)) read (1) Mnet_layers read (1) Mnet_weight read (1) Mnet_biases End Subroutine Read_Ann_info !-------------------------------------------------------------------------------- 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 :: X4d(:,:,:,:),xmon(:),lead(:),sampl(:) Integer :: ilon,ilat,imon,ilead,imem,ifcs,im,iobs 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/day"; lon_nm="lon"; lat_nm="lat"; z_nm="fcs_time"; t_nm="week" Allocate (lead(Nlead),sampl(Nspl)) Do ilead=1,Nlead; lead (ilead)=ilead; 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="Observation precipitation (CHIRPS)" Call wr_ncdf_4d (Fname,Lon,Lat,sampl,lead,Nlon,Nlat,Nspl,Nlead,Obs,"obs", & t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) var_ln="S2S precipitation of ecmwf forecast" Call wr_ncdf_4d_apd (Fname,EC,"fcs", Nlon,Nlat,Nspl,Nlead, & var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss) End Subroutine Write_BC_Outputs !---------------------------------------------------------------------------------------------- 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 Allocate (EC(Nlon,Nlat,Nspl,Nlead)) Allocate (LH(Nlon,Nlat,Nspl,Nlead), SH(Nlon,Nlat,Nspl,Nlead), OLR(Nlon,Nlat,Nspl,Nlead)) Allocate (TX(Nlon,Nlat,Nspl,Nlead), TM(Nlon,Nlat,Nspl,Nlead)) ispl=0 Do iyr=Yr1,Yr2 Do ifc=1,Nfcs write(Fname,"(A,I4.4,A)") trim(DirInEC)//"LH_SH_OLR_TP_cf_",iyr,"-"//hdate(ifc)//".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"tp") Where (X3d==ec_miss) X3d=Miss ispl=ispl+1 EC(:,:,ispl,:) = X3d(:,:,:) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"slhf") ! LH Where (X3d==ec_miss) X3d=Miss LH(:,:,ispl,:) = X3d(:,:,:) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"sshf") ! SH Where (X3d==ec_miss) X3d=Miss SH(:,:,ispl,:) = X3d(:,:,:) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"ttr") ! OLR Where (X3d==ec_miss) X3d=Miss OLR(:,:,ispl,:) = X3d(:,:,:) Deallocate (xlon,xlat,xtime,X3d) Enddo Enddo ispl=0 Do iyr=Yr1,Yr2 Do ifc=1,Nfcs write(Fname,"(A,I4.4,A)") trim(DirInEC)//"TX_TM_cf_",iyr,"-"//hdate(ifc)//".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"mx2t6") Where (X3d==ec_miss) X3d=Miss ispl=ispl+1 TX(:,:,ispl,:) = X3d(:,:,:) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"mn2t6") ! LH Where (X3d==ec_miss) X3d=Miss TM(:,:,ispl,:) = X3d(:,:,:) Deallocate (xlon,xlat,xtime,X3d) Enddo Enddo ! Read precipitation obs data: Allocate (OBS(Nlon,Nlat,Nspl,Nlead)) ispl=0 Do iyr=Yr1,Yr2 Do ifc=1,Nfcs write(Fname,"(A,I4.4,A)") trim(DirInObs)//"chirps_",iyr,"-"//hdate(ifc)//".nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"precip") Where (X3d==obs_miss) X3d=Miss ispl=ispl+1 OBS(:,:,ispl,:) = X3d(:,:,:) Deallocate (xlon,xlat,xtime,X3d) Enddo Enddo End Subroutine Read_data !---------------------------------------------------------------------------------------------- !------------------------------------------------------------------------------ Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Function_Valid_Ar1.f90" Include "/home/tanpv/Mylibs/Correlation.f90" Include "/home/tanpv/Mylibs/Sub_julian_day.f90" ! Include "/home/tanpv/Mylibs/Sub_Verification.f90" ! Include "/home/tanpv/Mylibs/Sub_neural_2024.f90" !------------------------------------------------------------------------------ End ! Include "/home/tanpv/Mylibs/Student_Prob.f90"