Use NetCDF Real, Parameter :: Miss=-99.0 Integer, Parameter :: Nreg=7, Nmon=12, IMiss=-99, NLead=6, Nmem=6 Integer, Parameter :: Yr1=2011, Yr2=2020, Ntime0=(Yr2-Yr1+1)*Nmon, Ntime=Ntime0+NLead Integer, Allocatable :: time(:) Integer :: i, LenSt, sys Character (Len=255) :: Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2 Character (Len=255) :: Dir, DirOut Character (Len=10) :: Model_Name Real, Allocatable :: lon(:), lat(:), Regs(:,:) !---------------------------------------------------------------------------------------------- ! Input data: Real, Allocatable, dimension(:,:,:,:) :: Pr, T2m, Tx, Tm ! (Nlon,Nlat,Nlead,Ntime) Real, Allocatable, dimension(:,:,:,:) :: h500, h850, olr, q2, rh2, slp, u200, u850, w700 Real, Allocatable, dimension(:,:,:) :: Pr_O, T2m_O, Tx_O, Tm_O ! (Nlon,Nlat,Ntime) Real, Allocatable :: xlon(:), xlat(:), xtime(:) Real, Allocatable :: X4d(:,:,:,:), X3d(:,:,:) Integer :: Nlon, Nlat, nlo, nla !---------------------------------------------------------------------------------------------- ! Model Climatology: Real, Allocatable, dimension(:,:,:,:) :: Pr_M_cli, T2m_M_cli, Tx_M_cli, Tm_M_cli ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:) :: Pr_O_cli, T2m_O_cli, Tx_O_cli, Tm_O_cli ! (Nlon,Nlat,Nmon) Real, Allocatable, dimension(:,:,:,:) :: Pr_BS, T2m_BS, Tx_BS, Tm_BS ! (Nlon,Nlat,Nlead,Nmon) !---------------------------------------------------------------------------------------------- ! Differences between Model and Obs: Real, Allocatable, dimension(:,:,:,:) :: Pr_Df, T2_Df, Tx_Df, Tm_Df ! (Nlon,Nlat,Nlead,Ntime) !---------------------------------------------------------------------------------------------- ! Verification Scores: Real, Allocatable, dimension(:,:,:,:) :: Pr_ME, Pr_MAE, Pr_RME, Pr_RMAE, Pr_COR, Pr_pval ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:,:) :: T2_ME, T2_MAE, T2_COR, T2_pval ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:,:) :: Tx_ME, Tx_MAE, Tx_COR, Tx_pval ! (Nlon,Nlat,Nlead,Nmon) Real, Allocatable, dimension(:,:,:,:) :: Tm_ME, Tm_MAE, Tm_COR, Tm_pval ! (Nlon,Nlat,Nlead,Nmon) !---------------------------------------------------------------------------------------------- Real, Allocatable :: RMEreg(:,:,:,:), RMAEreg(:,:,:,:), CORreg(:,:,:,:) ! (Nreg,Mon,LT,Mem) Real, Allocatable :: RMEwreg(:,:,:,:), RMAEwreg(:,:,:,:), CORwreg(:,:,:,:) ! (Nreg,Mon,Nws,Mem) Integer :: ilon,ilat,it,imon, imem Character (Len=1),Allocatable :: IRegs(:,:) Fname="/home/tanpv/Data/562_VNU/Mask_VN_0.1.nc" Call rd_ncdf_2d (Fname,lon,lat,Nlon,Nlat,Regs,"reg") ! 7 sub-regs mask !Allocate (IRegs(Nlon,Nlat)); Where (Regs==Miss); IRegs=" "; elsewhere; IRegs="*"; End where !do ilat=Nlat,1,-1 ! print"(200A1)",(IRegs(ilon,ilat),ilon=1,Nlon) !enddo !pause Call Read_OBS ! Read obs data of Pr, T2m, Tx, Tm ! Read model data: Dir ="/home/tanpv/Data/562_VNU/2011_2020_VN/" DirOut="/work/users/hmo/tanpv/From_HOME/562/Model_Data/" ! Model_Name="clWRF" Model_Name="RegCM" Do imem=1,Nmem Call Read_Fcs (Dir,Model_Name,imem) Call Read_Predictors (Dir,Model_Name,imem) Call Mask_Data print*,"Writing data ..." ! Dir="/work/users/hmo/tanpv/From_HOME/562/Model_Data/" Call Write_Model_Obs (DirOut,Model_Name,imem) Enddo CONTAINS !---------------------------------------------------------------------------------------------- Subroutine Write_Model_Obs (Dir,Model_Name,imem) Integer :: imon, id, imem Character (Len=*) :: Dir,Model_Name Character (Len=100) :: var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm Real, Allocatable :: xtime(:) lon_nm="lon"; lat_nm="lat"; t_nm="mon"; t_unit="months since 1982-01-01 " if (Allocated(X3d)) Deallocate (X3d) Allocate (X3d(Nlon,Nlat,Ntime),xtime(Ntime)) Do imon=1,Ntime; xtime(imon)=imon; Enddo Do ilead=1,NLead !write(Fname,"(4A,I1.1,3A)") trim(Dir),trim(Model_Name),"/","lead",ilead,"_",trim(Model_Name),"_VN.nc" write(Fname,"(2A,I1.1,3A,I2.2,A)") trim(Dir),"lead",ilead,"_",trim(Model_Name),"_",imem,"_VN.nc" print*,trim(Fname) var_ln="Observed Monthly Prec"; var_unit="mm/month" Call wr_ncdf_3d(Fname,lon,lat,xtime,Nlon,Nlat,Ntime,Pr_O,"pr_o",t_unit,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Observed Monthly T2m"; var_unit="C" Call wr_ncdf_3d_apd (Fname,T2m_O,"t2_o",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Observed Monthly Tx"; var_unit="C" Call wr_ncdf_3d_apd (Fname,Tx_O,"tx_o",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Observed Monthly Tm"; var_unit="C" Call wr_ncdf_3d_apd (Fname,Tm_O,"tm_o",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) ! var_ln="Model Monthly Prec"; var_unit="mm/month" X3d(:,:,:)=Pr(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"pr",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Model Monthly T2m"; var_unit="C" X3d(:,:,:)=T2m(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"t2",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Model Monthly Tx"; X3d(:,:,:)=Tx(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"tx",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Model Monthly Tm"; X3d(:,:,:)=Tm(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"tm",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) ! var_ln="Model Geopotential height at 500 hPa (m)"; var_unit="m" X3d(:,:,:)=h500(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"h500",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Model Geopotential height at 850 hPa (m)"; var_unit="m" X3d(:,:,:)=h850(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"h850",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Model TOA OUTGOING LONG WAVE (W m-2)"; var_unit="w/m2" X3d(:,:,:)=olr(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"olr",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Model QV at 2 M (kg kg-1)"; var_unit="kg/kg" X3d(:,:,:)=q2(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"q2",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Model Relative Humidity at 2m (%)"; var_unit="%" X3d(:,:,:)=rh2(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"rh2",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Model Sea Levelp Pressure (hPa)"; var_unit="hPa" X3d(:,:,:)=slp(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"slp",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Model x-wind component at 200 hPa (m s-1)"; var_unit="m/s" X3d(:,:,:)=u200(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"u200",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Model x-wind component at 850 hPa (m s-1)"; var_unit="m/s" X3d(:,:,:)=u850(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"u850",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) var_ln="Model z-wind component at 700 hPa (m s-1)"; var_unit="m/s" X3d(:,:,:)=w700(:,:,ilead,:) Call wr_ncdf_3d_apd (Fname,X3d,"w700",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss) Enddo ! h500, h850, olr, q2, rh2, slp, u200, u850, w700 End Subroutine Write_Model_Obs !---------------------------------------------------------------------------------------------- Subroutine Mask_Data Integer :: ilon,ilat,ilead,it Do ilead=1,NLead Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)==Miss) then h500(ilon,ilat,ilead,:)=Miss; h850(ilon,ilat,ilead,:)=Miss; olr(ilon,ilat,ilead,:)=Miss q2(ilon,ilat,ilead,:)=Miss; rh2(ilon,ilat,ilead,:)=Miss; slp(ilon,ilat,ilead,:)=Miss u200(ilon,ilat,ilead,:)=Miss; u850(ilon,ilat,ilead,:)=Miss; w700(ilon,ilat,ilead,:)=Miss Pr(ilon,ilat,ilead,:)=Miss; T2m(ilon,ilat,ilead,:)=Miss; Tx(ilon,ilat,ilead,:)=Miss; Tm(ilon,ilat,ilead,:)=Miss endif Enddo Enddo Enddo End Subroutine Mask_Data !---------------------------------------------------------------------------------------------- Subroutine Read_Predictors (Dir,Model_Name,imem) Character (Len=*) :: Dir,Model_Name Integer :: imem, Nla, Nlo, Nt, Yr0, it1, it2, ilead, k0 Character (Len=255) Fname ! h500, h850, olr, q2, rh2, slp, u200, u850, w700 if (Allocated(h500)) Deallocate (h500); if (Allocated(h850)) Deallocate (h850); if (Allocated(olr)) Deallocate (olr) if (Allocated(q2)) Deallocate (q2); if (Allocated(rh2)) Deallocate (rh2); if (Allocated(slp)) Deallocate (slp) if (Allocated(u200)) Deallocate (u200); if (Allocated(u850)) Deallocate (u850); if (Allocated(w700)) Deallocate (w700) Allocate (h500(Nlon,Nlat,Nlead,Ntime), h850(Nlon,Nlat,Nlead,Ntime), olr(Nlon,Nlat,Nlead,Ntime)) Allocate (q2(Nlon,Nlat,Nlead,Ntime), rh2(Nlon,Nlat,Nlead,Ntime), slp(Nlon,Nlat,Nlead,Ntime)) Allocate (u200(Nlon,Nlat,Nlead,Ntime), u850(Nlon,Nlat,Nlead,Ntime), w700(Nlon,Nlat,Nlead,Ntime)) h500=Miss; h850=Miss; olr=Miss; q2=Miss; rh2=Miss; slp=Miss; u200=Miss; u850=Miss; w700=Miss; if (Allocated(X3d)) Deallocate (X3d) k0=4 ! So timeslice (thang) khong co hoac khong dung Do ilead=1,NLead write(Fname,"(4A,I1.1,3A,I2.2,A)") trim(Dir),trim(Model_Name),"/","lead",ilead,"_",trim(Model_Name),"_",imem,"_current_pre_VN.nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"h500") if (ilead==1.AND.Nt<117) print*,"===>",trim(Fname) it1=ilead+k0; it2=it1+Nt-1 h500(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"h850") h850(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"olr") olr(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"q2") q2(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"rh2") rh2(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"slp") slp(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"u200") u200(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"u850") u850(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"w700") w700(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Enddo End Subroutine Read_Predictors !---------------------------------------------------------------------------------------------- Subroutine Read_Fcs (Dir,Model_Name,imem) ! Read forecast data Character (Len=*) :: Dir,Model_Name Integer :: imem, Nla, Nlo, Nt, Yr0, it1, it2, ilead, k0 Character (Len=255) Fname !Dir="/home/tanpv/Data/562_VNU/Rfc_VN/" if (Allocated(Pr)) Deallocate (Pr) if (Allocated(T2m)) Deallocate (T2m) if (Allocated(Tx)) Deallocate (Tx) if (Allocated(Tm)) Deallocate (Tm) Allocate (Pr(Nlon,Nlat,Nlead,Ntime), T2m(Nlon,Nlat,Nlead,Ntime)) Allocate (Tx(Nlon,Nlat,Nlead,Ntime), Tm(Nlon,Nlat,Nlead,Ntime)) ! (Nlon,Nlat,Nlead,Ntime) Pr=Miss; T2m=Miss; Tx=Miss; Tm=Miss if (Allocated(X3d)) Deallocate (X3d) k0=4 ! Tru k0 thang dau cua timeseries do khong chay mo hinh Do ilead=1,NLead write(Fname,"(4A,I1.1,3A,I2.2,A)") trim(Dir),trim(Model_Name),"/","lead",ilead,"_",trim(Model_Name),"_",imem,"_current_VN.nc" !print*,trim(Fname) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"pr") !print*,X3d;pause it1=ilead+k0; it2=it1+Nt-1 ! Missing data will be filled at "first k0 values" and "last NLead-ilead values" ! Nt = 117, Ntime = 126 (??) if (ilead==1.AND.Nt<117) print*,"===>",trim(Fname) Pr(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tas") T2m(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tasmax") Tx(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tasmin") Tm(:,:,ilead,it1:it2)=X3d(:,:,1:Nt) Deallocate (xlon,xlat,xtime,X3d) Enddo End Subroutine Read_Fcs !---------------------------------------------------------------------------------------------- Subroutine Read_OBS ! Read observed gridded data of Pr, T2m, Tx, Tm Integer :: j, it, ilon, ilat, Nla, Nlo, Nt, Ntt, Yr0, it1, it2, imon Character (Len=255) Fname Allocate (Pr_O(Nlon,Nlat,Ntime), T2m_O(Nlon,Nlat,Ntime), Tx_O(Nlon,Nlat,Ntime), Tm_O(Nlon,Nlat,Ntime)) Pr_O=Miss; T2m_O=Miss; Tx_O=Miss; Tm_O=Miss Yr0=1980; ! Rain data are available from Jan 1980 it1=(Yr1-Yr0)*Nmon+1 ! Lay tu 1/2011 it2=it1+Ntime-1 ! Lay den 6/2021 ! Tam thoi lay it2 tuong ung voi 12/2020: !it2=it1+Ntime0-1 !print*,it1,it2,it2-it1+1,Yr0+(it1-1)/Nmon; pause Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/R.1980.2018.mm.nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"rain") ! Pr_O(:,:,1:Ntime)=X3d(:,:,it1:it2) ! ***** Tam thoi thay it2 bang Nt vi Nt hien tai ung voi 2018 Ntt=Nt-it1+1; Pr_O(:,:,1:Ntt)=X3d(:,:,it1:Nt) !print*,"OBS-R: Nt,it1,Ntt=",Nt,it1,Ntt Deallocate (xlon,xlat,xtime,X3d) Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/T2m.1961.2018.mm.nc" Yr0=1961; ! Temp data are available from Jan 1961 it1=(Yr1-Yr0)*Nmon+1; it2=it1+Ntime-1 ! Tam thoi lay it2 tuong ung voi 12/2020: !it2=it1+Ntime0-1 !print*,it1,it2,it2-it1+1,Yr0+(it1-1)/Nmon; pause Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"t2m") ! T2m_O(:,:,1:Ntime)=X3d(:,:,it1:it2) ! ***** Tam thoi thay it2 bang Nt vi Nt hien tai ung voi 2018 Ntt=Nt-it1+1; T2m_O(:,:,1:Ntt)=X3d(:,:,it1:Nt) Deallocate (xlon,xlat,xtime,X3d) Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/Tx.1961.2018.mm.nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tx") ! Tx_O(:,:,1:Ntime)=X3d(:,:,it1:it2) ! ***** Tam thoi thay it2 bang Nt vi Nt hien tai ung voi 2018 Ntt=Nt-it1+1; Tx_O(:,:,1:Ntt)=X3d(:,:,it1:Nt) Deallocate (xlon,xlat,xtime,X3d) Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/Tm.1961.2018.mm.nc" Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tm") ! Tm_O(:,:,1:Ntime)=X3d(:,:,it1:it2) ! ***** Tam thoi thay it2 bang Nt vi Nt hien tai ung voi 2018 Ntt=Nt-it1+1; Tm_O(:,:,1:Ntt)=X3d(:,:,it1:Nt) !print*,"OBS-T: Nt,it1,Ntt=",Nt,it1,Ntt; pause Deallocate (xlon,xlat,xtime,X3d) End Subroutine Read_OBS !---------------------------------------------------------------------------------------------- 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/DayOfMon.f90" END !---------------------------------------------------------------------------