module module_read_NWP use netcdf use hdf5 use module_satalite use wgrib2api implicit none !include 'netcdf.inc' public :: init_NWP, Read_NWP, & allocate_NWP,deallocate_nwp integer,public :: n_lat,n_lon,nlevel real(kind=4),public,allocatable, dimension(:,:,:) :: Tt,qt,cfract,ozonet real(kind=4),public,allocatable, dimension(:,:,:,:) :: cloudt real(kind=4),public,allocatable, dimension(:,:) :: LATP, LONP real(kind=4),public,allocatable, dimension(:,:) :: tst,T2mt, q2mt, u10mt, v10mt real(kind=4),public,allocatable, dimension(:,:) :: PSFCt,seaicet !(Pa) real(kind=4),public,allocatable, dimension(:) :: Pt !(hPa) real(kind=4),public,allocatable, dimension(:,:,:) :: T1,q1,cfrac1,ozone1 real(kind=4),public,allocatable, dimension(:,:,:,:) :: cloud1 real(kind=4),public,allocatable, dimension(:,:) :: ts1,T2m1, q2m1, u10m1, v10m1 real(kind=4),public,allocatable, dimension(:,:) :: PSFC1,seaice1 real(kind=4),public,allocatable, dimension(:) :: P real(kind=4),public,allocatable, dimension(:,:,:) :: T2,q2,cfrac2,ozone2 real(kind=4),public,allocatable, dimension(:,:,:,:) :: cloud2 real(kind=4),public,allocatable, dimension(:,:) :: ts2,T2m2, q2m2, u10m2, v10m2 real(kind=4),public,allocatable, dimension(:,:) :: PSFC2,seaice2 real(kind=4),public,allocatable, dimension(:,:,:) :: T3,q3,cfrac3,ozone3 real(kind=4),public,allocatable, dimension(:,:,:,:) :: cloud3 real(kind=4),public,allocatable, dimension(:,:) :: ts3,T2m3, q2m3, u10m3, v10m3 real(kind=4),public,allocatable, dimension(:,:) :: PSFC3,seaice3 integer,private :: nx,ny,nz integer,private :: ncid,ncerror,dimid,varid,iret character*400,private :: invfile contains !================================================================================================= !================================================================================================= subroutine init_NWP(nwpfilename,satfilename) implicit none character*400 :: nwpfilename character*400 :: satfilename character*13 :: filetime character*50 :: L1name integer :: ids,idl print*,'nwpfilename',nwpfilename ids = index(nwpfilename,'gmf.gra.') idl = index(satfilename,'_MULT_') L1name = satfilename(idl-35:idl-22)//satfilename(idl+10:idl+24)//nwpfilename(ids+8:ids+8+12) invfile = '/fy4_ads_data/stss_v2.0/FILE_GRAPES_inv'//trim(adjustl(L1name)) iret = grb2_mk_inv(trim(nwpfilename),trim(invfile)) !print*,'invfile',invfile iret = grb2_inq(trim(nwpfilename),trim(invfile),':TMP:',':surface:',nx=nx,ny=ny) nz = 40 n_lat = ny n_lon = nx nlevel = nz return end subroutine init_NWP !================================================================================================= !================================================================================================= subroutine allocate_NWP(allocstatus) integer :: allocstatus !nx,ny,nz allocate(Tt (nx, ny, nz)) allocate(ozonet (nx, ny, nz)) allocate(qt (nx, ny, nz)) allocate(psfct (nx, ny)) allocate(seaicet (nx, ny)) allocate(T2mt (nx, ny)) allocate(tst (nx, ny)) allocate(q2mt (nx, ny)) allocate(u10mt (nx, ny)) allocate(cloudt (nx, ny, 6, nz-1)) allocate(cfract (nx, ny, nz-1)) allocate(v10mt (nx, ny)) allocate(Pt (nz)) allocstatus = allocated return end subroutine allocate_NWP !================================================================================================= !================================================================================================= subroutine deallocate_nwp(allocstatus) integer :: allocstatus deallocate(Tt ) deallocate(ozonet ) deallocate(qt ) deallocate(psfct ) deallocate(seaicet ) deallocate(T2mt ) deallocate(tst ) deallocate(q2mt ) deallocate(u10mt ) deallocate(cloudt ) deallocate(cfract ) deallocate(v10mt ) deallocate(Pt ) allocstatus = deallocated return end subroutine deallocate_NWP !================================================================================================= !================================================================================================= Subroutine Read_NWP( & nwpfilename , & sfcfilename , & imonth ) integer(kind=4) :: ncid,sfcid integer :: ncerror,varid, len_file, ndims, dimids(4), natts, xtyp character*10 :: vname, dimname character*300 :: sfcfilename,nwpfilename,O3filename integer(kind=4) :: i,j,k,ilay,ilev, nlt,imonth real(kind=4) :: es real(kind=4) :: q2 !water vapor pressure at 2m real(kind=4) :: scale,offset real(kind=4) :: scaleq,offsetq real(kind=4) :: cc character(len=10),allocatable,dimension(:) :: vartemp real(kind=4),allocatable,dimension(:,:) :: Ttt real(kind=4),allocatable,dimension(:,:) :: qtt real(kind=4),allocatable,dimension(:,:,:,:) :: ozonett real(kind=4),allocatable,dimension(:) :: pressor_ozone real(kind=4),allocatable,dimension(:,:) :: sstt real(kind=4),allocatable,dimension(:,:) :: LATPt,LONPt real(kind=4),allocatable,dimension(:,:) :: sktt real(kind=4),allocatable,dimension(:,:) :: T2mtt real(kind=4),allocatable,dimension(:,:) :: q2mtt real(kind=4),allocatable,dimension(:,:) :: u10mtt real(kind=4),allocatable,dimension(:,:) :: v10mtt real(kind=4),allocatable,dimension(:,:) :: seaicett real(kind=4),allocatable,dimension(:,:) :: psfctt real(kind=4),allocatable,dimension(:,:) :: landmasknwp integer(kind=4) :: n_time integer :: hdferror integer(HID_T) :: file_id integer(HID_T) :: dataset_id,dspace_id integer(HSIZE_T),dimension(4) :: dims_4d integer(HSIZE_T) :: maxdimsd(4) O3filename = '/fy4_ads_data/stss_v2.0/ProduceBackgroundData/instrument_geo/2017_monthly_ozone_profile_GRAPES.h5' n_time = 12 allocate(ozonett (n_lon,n_lat,nlevel,n_time)) allocate(vartemp (nlevel)) vartemp = (/':0.1 mb:', ':0.2 mb:', ':0.5 mb:', ':1 mb:', ':1.5 mb:', ':2 mb:', & ':3 mb:' , ':4 mb:' , ':5 mb:' , ':7 mb:', ':10 mb:' , ':20 mb:', & ':30 mb:', ':50 mb:' , ':70 mb:' , ':100 mb:', ':125 mb:', ':150 mb:', & ':175 mb:', ':200 mb:', ':225 mb:', ':250 mb:', ':275 mb:', ':300 mb:', & ':350 mb:', ':400 mb:', ':450 mb:', ':500 mb:', ':550 mb:', ':600 mb:', & ':650 mb:', ':700 mb:', ':750 mb:', ':800 mb:', ':850 mb:', ':900 mb:', & ':925 mb:', ':950 mb:', ':975 mb:', ':1000 mb:' /) ! vartemp = (/':10 mb:' , ':20 mb:', & ! ':30 mb:', ':50 mb:' , ':70 mb:' , ':100 mb:', ':125 mb:', ':150 mb:', & ! ':175 mb:', ':200 mb:', ':225 mb:', ':250 mb:', ':275 mb:', ':300 mb:', & ! ':350 mb:', ':400 mb:', ':450 mb:', ':500 mb:', ':550 mb:', ':600 mb:', & ! ':650 mb:', ':700 mb:', ':750 mb:', ':800 mb:', ':850 mb:', ':900 mb:', & ! ':925 mb:', ':950 mb:', ':975 mb:', ':1000 mb:' /) !==========================read profile================================= call h5open_f(hdferror) call h5fopen_f(trim(O3filename), H5F_ACC_RDONLY_F, file_id, hdferror) call h5dopen_f(file_id, 'Ozone',dataset_id, hdferror) call h5dget_space_f(dataset_id, dspace_id, hdferror) call h5sget_simple_extent_dims_f(dspace_id, dims_4d, maxdimsd,hdferror) ! print*,dims_4d ! stop call h5dread_f(dataset_id,H5T_NATIVE_REAL,ozonett,dims_4d,hdferror) call h5dclose_f(dataset_id,hdferror) call h5fclose_f(file_id,hdferror) call h5close_f(hdferror) do i = 1,nz iret = grb2_inq(trim(nwpfilename),trim(invfile),':TMP:',trim(vartemp(i)), & data2=Ttt) if(i<=10)then iret = grb2_inq(trim(nwpfilename),trim(invfile),':SPFH:',trim(vartemp(11)), & data2=qtt) else iret = grb2_inq(trim(nwpfilename),trim(invfile),':SPFH:',trim(vartemp(i)), & data2=qtt) endif do j=1,n_lat Tt(:,j,i) = Ttt(:,n_lat-j+1) qt(:,j,i) = qtt(:,n_lat-j+1) enddo enddo print*,'+++++',trim(nwpfilename) Pt = (/ 0.1, 0.2, 0.5, 1.0, 1.5 , 2.0, & 3.0, 4.0, 5.0, 7.0, 10.0 , 20.0, & 30., 50., 70., 100., 125. , 150., & 175., 200., 225., 250., 275. , 300., & 350., 400., 450., 500., 550. , 600., & 650., 700., 750., 800., 850. , 900., & 925., 950., 975., 1000. /) ! Pt = (/ 10 , 20, & ! 30, 50 , 70 , 100, 125 , 150, & ! 175, 200, 225, 250, 275 , 300, & ! 350, 400, 450, 500, 550 , 600, & ! 650, 700, 750, 800, 850 , 900, & ! 925, 950, 975, 1000 /) where(Tt<-900. .or. Tt>500. ) Tt = -9999.9 endwhere where(qt<=0. .or. qt>900.) !qt = -9999.9 qt = 1E-11 !20210414 caobei endwhere where(qt>0) qt = qt*1e6*28.966/18.0 endwhere iret = grb2_inq(trim(nwpfilename),trim(invfile),':TMP:',':surface:', & lat=LATPt,lon=LONPt) !=========================================================================== !==================== read surface ========================================= iret = grb2_inq(trim(nwpfilename),trim(invfile),':TMP:',':surface:', & data2=sktt) where(sktt>900.) sktt = -9999.9 endwhere iret = grb2_inq(trim(nwpfilename),trim(invfile),':PRES:',':surface:', & data2=psfctt) where(psfctt>9.99e10) psfctt = -9999.9 endwhere !print*,'skt',ts(1:10,10,1) iret = grb2_inq(trim(nwpfilename),trim(invfile),':TMP:',':2 m above ground:', & data2=T2mtt) where(T2mtt>900.) T2mtt = -9999.9 endwhere iret = grb2_inq(trim(nwpfilename),trim(invfile),':SPFH:',':2 m above ground:', & data2=q2mtt) where(q2mtt<=0. .or. q2mtt>900.) q2mtt = -9999.9 endwhere where(q2mtt>0) q2mtt = q2mtt*1e6*28.966/18.0 endwhere iret = grb2_inq(trim(nwpfilename),trim(invfile),':UGRD:',':10 m above ground:', & data2=u10mtt) where(u10mtt>900.) u10mtt = -9999.9 endwhere iret = grb2_inq(trim(nwpfilename),trim(invfile),':VGRD:',':10 m above ground:', & data2=v10mtt) where(v10mtt>900.) v10mtt = -9999.9 endwhere do i = 1, n_lat tst(:,i) = sktt(:,n_lat-i+1) psfct(:,i) = psfctt(:,n_lat-i+1) T2mt(:,i) = T2mtt(:,n_lat-i+1) u10mt(:,i) = u10mtt(:,n_lat-i+1) v10mt(:,i) = v10mtt(:,n_lat-i+1) ozonet(:,i,:) = ozonett(:,n_lat-i+1,:,imonth)*1e6*28.966/48.0 !ppmv LONP(:,i) = LONPt(:,n_lat-i+1) LATP(:,i) = LATPt(:,n_lat-i+1) q2mt(:,i) = q2mtt(:,n_lat-i+1) enddo seaicet(:,:) = 0.0 !======================== end of read ============================================= !================================================================================== cloudt(:,:,:,:) = 0.0 cfract(:,:,:) = 0.0 ! print*,'++++',LONP(1:5,1) ! print*,'----',LATP(1,1:5) ! stop !===================================================== print*,'read grapes end' print*,'T=',maxval(Tt),minval(Tt) print*,'q=',maxval(qt),minval(qt) print*,'cloud=',maxval(cloudt),minval(cloudt) print*,'cfrac=',maxval(cfract),minval(cfract) print*,'psfc=',maxval(psfct),minval(psfct) print*,'ts=',maxval(tst),minval(tst) print*,'q2m=',maxval(q2mt),minval(q2mt) print*,'T2m=',maxval(T2mt),minval(T2mt) print*,'u10m=',maxval(u10mt),minval(u10mt) print*,'v10m=',maxval(v10mt),minval(v10mt) print*,'ozone=',maxval(ozonet),minval(ozonet) deallocate(Ttt ) deallocate(qtt ) deallocate(ozonett ) deallocate(LATPt ) deallocate(LONPt ) deallocate(sktt ) deallocate(T2mtt ) deallocate(q2mtt ) deallocate(u10mtt ) deallocate(v10mtt ) deallocate(psfctt ) deallocate(vartemp ) return End Subroutine Read_NWP !================================================================================================= !================================================================================================= end module module_read_NWP