! Chuong trinh nay dung de hieu chinh mua du bao noi mua (s2s) cua ECMWF ! cho ca hai truong hop training va testing ! Qua trinh training la xay dung ham TP=f(X1,...,Xm) su dung mang ANN mot lop an, mot lop Input, mot lop Output ! Version hien tai chi dung cho 1-2 thang co dinh ma duoc xem la dong nhat ve quy luat hieu chinh ! Neu training thi chuong trinh se tra ve ket qua hieu chinh dong thoi luu lai file trong so va cac tham so cua ANN ! Neu teating thi chuong trinh se doc thong tin trong so va cac tham so cua ANN tu file da luu va tra ve ket qua hieu chinh. ! Cac tham so can chu y: ! 1) Yr1, Yr2: Cac nam duoc tinh ! 2) Lw: Kich thuoc cua so thoi gian se truot (tong hoac trung binh truot) ! 3) So phien du bao trong tap so lieu (vi du, moi thang 9 phien x 2 thang = 18 phien) ! (cac ngay ung voi cac phien du bao duoc luu trong file hdate.txt) 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 :: 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 Integer :: yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2 Real*8 :: factor, offset, Missing Real, Allocatable :: lon(:), lat(:), Regs(:,:) Integer, Allocatable :: lead(:), fcsdate(:,:) !---------------------------------------------------------------------------------------------- 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 Character (Len=5) :: yvarn, Model_nm Real,allocatable, Dimension(:,:) :: Inputs, Outputs, Inputs_X, Outputs_Y Real,allocatable, Dimension(:) :: min_in, max_in, min_out, max_out Real :: mean_error, mean_error_check !------------------------------------------------------------------------------- integer :: Nin, Nhidden, Nout, iin, ihid, iout integer :: i,j,k Real,allocatable, Dimension(:,:) :: weight_ih, weight_ho, moment_ih, moment_ho Real,allocatable, Dimension(:) :: fin, fhid, fout, ftarget, delta_o, delta_h Real, Allocatable :: X_inp(:,:), Y_inp(:), Yest(:,:), Yhq(:,:) ! Training Real :: learning_rate0, learning_rate,momentum, error, tmp, mintmp logical :: debug !------------------------------------------------------------------------------- ! controller_save_info: ! Nin,Nhidden,Nout - Nin=MXvar; Nout=MYvar Integer :: MXvar, MYvar, Mx, Mvars Integer :: Ncase,icase Real,allocatable, Dimension(:,:,:) :: weight_ip ! (Ncase,Nin,Nhidden) - Ncase=So lan tinh, Real,allocatable, Dimension(:,:,:) :: weight_op ! (Ncase,Nhidden,Nout) Real,allocatable, Dimension(:,:) :: min_ip, max_ip, min_op, max_op ! (Ncase) ! For saving in ANN subroutine: Real,allocatable, Dimension(:,:) :: wei_i, wei_o Real,allocatable, Dimension(:) :: min_i, max_i, min_o, max_o !------------------- For ANN ------------------------------------------------------------- !Namelist /ann_info/ Nlayer,learn_rate,Nepoch Namelist /Dir_info/ Dir_home,Dir_Mask,DirInEC,DirInObs,Dir_ANN_Save,Dir_Corrected Namelist /ann_info/ Nhidden,learning_rate,momentum,training_mode,Nepoc,debug Namelist /Model_info/ Yr1, Yr2, Lw, Nfc_per_year !-------------------------------------------------------------------------------- Open (1, file="Namelist_S2S_ecmwf.txt",Status="old") Read (1, Dir_info) Read (1, ann_info) Read (1, Model_info) ! Allocate (Layer_size(Nlayer)) ! Read (1, layer_info) Close (1) Nws=NLead-Lw+1 Nyrs=(Yr2-Yr1+1) ! DirInEC ="/home/tanpv/Progs/S2S_4Giang/Data/EC_4Giang/" ! DirInObs ="/home/tanpv/Progs/S2S_4Giang/Data/CHIRPS_Giang/" ! Dir_ANN_Save ="/home/tanpv/Progs/S2S_4Giang/Data/ANN_Save/" ! Dir_Corrected ="/home/tanpv/Progs/S2S_4Giang/Res/" ! Fname="/home/tanpv/Progs/S2S_4Giang/Data/Mask_Giang.nc" Call rd_ncdf_2d (trim(Dir_Mask),Lon,Lat,Nlon,Nlat,Regs,"vn") Where (IsNaN(Regs)) Regs=Miss Do isea=1,Nseason ! Calculation for each season among 4 seasons Fname=trim(Dir_home)//"hdate_"//trim(Season(isea))//".txt" Open (1,file=trim(Fname), status="old") Nfcs = Nfc_per_year(isea) ! Number of forecast in each year Nspl = Nfcs*Nyrs ! Total forecast over entire period if (Allocated (hdate)) Deallocate (hdate) Allocate (hdate(Nfcs)) ! Days of forecast Do i=1,Nfcs read (1,*) hdate(i) Enddo Close (1) Call Read_data ! Read all data of TP, LH, SH, OLR, TX, TM from entire domain (no mask for Central Highland), and CHIRPs if (training_mode) then MXvar = 22; MYvar = 1 Mvars = MXvar Ncase = Npoint() ! Obtain from Function Npoint if (Allocated (weight_ip)) Deallocate (weight_ip); if (Allocated (weight_op)) Deallocate (weight_op) if (Allocated (min_ip)) Deallocate (min_ip); if (Allocated (max_ip)) Deallocate (max_ip) if (Allocated (min_op)) Deallocate (min_op); if (Allocated (max_op)) Deallocate (max_op) Allocate (weight_ip(1:Ncase,0:MXvar,1:Nhidden),weight_op(Ncase,0:Nhidden,MYvar)) Allocate (min_ip(Ncase,MXvar),max_ip(Ncase,MXvar),min_op(Ncase,MYvar),max_op(Ncase,MYvar)) else !Fname=trim(Dir_ANN_Save)//"Sav_"//trim(Season(isea))//".txt" !Ncase = Npoint() write (Fname,"(A,I4.4,A,I4.4,A)") trim(Dir_ANN_Save)//"Sav_"//trim(Season(isea))//"_",yr1,"-",Yr2,".txt" Call Read_Ann_info (Fname) Mvars = MXvar endif if (Allocated (wei_i)) Deallocate (wei_i); if (Allocated (wei_o)) Deallocate (wei_o) if (Allocated (min_i)) Deallocate (min_i); if (Allocated (max_i)) Deallocate (max_i) if (Allocated (min_o)) Deallocate (min_o); if (Allocated (max_o)) Deallocate (max_o) if (Allocated (EC_bc)) Deallocate (EC_bc) Allocate (wei_i(0:MXvar,1:Nhidden), wei_o(0:Nhidden,MYvar)) ! To be used in ANN Sub Allocate (min_i(MXvar),max_i(MXvar),min_o(MYvar),max_o(MYvar)) ! To be used in ANN Sub Allocate (EC_bc(Nlon,Nlat,Nspl,Nlead)) EC_bc = Miss; icase = 0 Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)>0) then ! for land gridpoint only Do ilead=1,Nlead icase = icase+1 Call Make_data (ilon,ilat,ilead) if (.not.training_mode) Call Load_Ann_info (icase) Call ANN_2025(Inputs_X,Outputs_Y,Nsample1,Mvars,1,Miss) !ANN_2023 (Xinp,Youtp,Nspl,MXvar,MYvar,Miss) if (training_mode) Call Save_Ann_info (icase) !if (training_mode) then ! Call Save_Ann_info (icase) !else ! Call Load_Ann_info (icase) !endif EC_bc(ilon,ilat,1:Nspl,ilead) = Yest(1:Nspl,1) ! Ann outputs Enddo endif Enddo Enddo Call Calc_weekly_data !Fname=trim(Dir_Corrected)//"S2S_"//trim(Season(isea))//"_corrected.nc" if (training_mode) then write (Fname,"(A,I4.4,A,I4.4,A)") trim(Dir_Corrected)//"S2S_"//trim(Season(isea))//"_",yr1,"-",Yr2,"_corrected.nc" else write (Fname,"(A,I4.4,A,I4.4,A)") trim(Dir_Corrected)//"S2S_"//trim(Season(isea))//"_",yr1,"-",Yr2,"_corrected_test.nc" endif Call Write_BC_Outputs (Fname) if (training_mode) then !Fname=trim(Dir_ANN_Save)//"Sav_"//trim(Season(isea))//".txt" write (Fname,"(A,I4.4,A,I4.4,A)") trim(Dir_ANN_Save)//"Sav_"//trim(Season(isea))//"_",yr1,"-",Yr2,".txt" Call Write_Ann_info (Fname) endif Enddo CONTAINS !---------------------------------------------------------------------------------------------- Subroutine Save_Ann_info (icase) Integer :: icase ! Nin=MXvar; Nout=MYvar wei_i, wei_o weight_ip(icase,0:MXvar,1:Nhidden) = wei_i(0:MXvar,1:Nhidden) weight_op(icase,0:Nhidden,1:MYvar) = wei_o(0:Nhidden,1:MYvar) min_ip(icase,1:MXvar) = min_i(1:MXvar); max_ip(icase,1:MXvar) = max_i(1:MXvar) min_op(icase,1:MYvar) = min_o(1:MYvar); max_op(icase,1:MYvar) = max_o(1:MYvar) End Subroutine Save_Ann_info !---------------------------------------------------------------------------------------------- Subroutine Load_Ann_info (icase) Integer :: icase ! Nin=MXvar; Nout=MYvar wei_i, wei_o wei_i(0:MXvar,1:Nhidden) = weight_ip(icase,0:MXvar,1:Nhidden) wei_o(0:Nhidden,1:MYvar) = weight_op(icase,0:Nhidden,1:MYvar) min_i(1:MXvar) = min_ip(icase,1:MXvar); max_i(1:MXvar) = max_ip(icase,1:MXvar) min_o(1:MYvar) = min_op(icase,1:MYvar); max_o(1:MYvar) = max_op(icase,1:MYvar) End Subroutine Load_Ann_info !---------------------------------------------------------------------------------------------- Subroutine Read_Ann_info (Fname) Character (Len=*) :: Fname Open (1,file=trim(Fname),form="unformatted",status="old") read(1) MXvar,Nhidden,MYvar,Ncase if (Allocated (weight_ip)) Deallocate (weight_ip); if (Allocated (weight_op)) Deallocate (weight_op) if (Allocated (min_ip)) Deallocate (min_ip); if (Allocated (max_ip)) Deallocate (max_ip) if (Allocated (min_op)) Deallocate (min_op); if (Allocated (max_op)) Deallocate (max_op) Allocate (weight_ip(1:Ncase,0:MXvar,1:Nhidden),weight_op(Ncase,0:Nhidden,MYvar)) Allocate (min_ip(Ncase,MXvar),max_ip(Ncase,MXvar),min_op(Ncase,MYvar),max_op(Ncase,MYvar)) read(1) weight_ip read(1) weight_op read(1) min_ip read(1) max_ip read(1) min_op read(1) max_op Close (1) End Subroutine Read_Ann_info !---------------------------------------------------------------------------------------------- Subroutine Write_Ann_info (Fname) Character (Len=*) :: Fname Open (1,file=trim(Fname),form="unformatted",status="unknown") write(1) MXvar,Nhidden,MYvar,Ncase write(1) weight_ip write(1) weight_op write(1) min_ip write(1) max_ip write(1) min_op write(1) max_op Close (1) End Subroutine Write_Ann_info !---------------------------------------------------------------------------------------------- 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 Mx=Mvars ! So bien doc lap 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 if (Allocated(Inputs_X)) Deallocate (Inputs_X) if (Allocated(Outputs_Y)) Deallocate (Outputs_Y) Allocate (Inputs_X(Mx,Nsample1), Outputs_Y(1,Nsample1)) i=0 Do it=1,Nspl i=i+1 Outputs_Y(1,i)=Y_inp(it) Do j=1,Mx Inputs_X(j,i)=X_inp(it,j) Enddo Enddo End Subroutine Make_data !---------------------------------------------------------------------------------------------- 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 Calc_Clim_tp Integer :: ilon,ilat,ilead Real :: X(Nspl) ! Real, Allocatable, Dimension(:,:,:) :: TP_cli, Obs_cli ! (Nlon,Nlat,NLead) - Clim over all forecasts ! Allocate (EC(Nlon,Nlat,Nspl,Nlead)) ! Allocate (LH(Nlon,Nlat,Nspl,Nlead), SH(Nlon,Nlat,Nspl,Nlead), OLR(Nlon,Nlat,Nspl,Nlead)) Allocate (TP_cli(Nlon,Nlat,NLead),Obs_cli(Nlon,Nlat,NLead)) TP_cli=Miss; Obs_cli=Miss Do ilon=1,Nlon Do ilat=1,Nlat Do ilead=1,Nlead X(1:Nspl) = EC (ilon,ilat,1:Nspl,ilead) TP_cli (ilon,ilat,ilead) = Ave(X,Nspl,Miss) X(1:Nspl) = OBS(ilon,ilat,1:Nspl,ilead) Obs_cli(ilon,ilat,ilead) = Ave(X,Nspl,Miss) Enddo Enddo Enddo End Subroutine Calc_Clim_tp !---------------------------------------------------------------------------------------------- Subroutine Areal_Mean (X,Nlon,Nlat,NLead,Y,ireg) 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 (Int(Regs(ilon,ilat))==ireg.AND.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); if (Allocated (LH)) Deallocate (LH) if (Allocated (SH)) Deallocate (SH); if (Allocated (OLR)) Deallocate (OLR) if (Allocated (TX)) Deallocate (TX); if (Allocated (TM)) Deallocate (TM) 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: if (Allocated (OBS)) Deallocate (OBS) 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 !---------------------------------------------------------------------------------------------- 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 !---------------------------------------------------------------------------