Metgrid, PRES, and the intermediate file format

Dataset specific topics involving WPS.

Re: Metgrid, PRES, and the intermediate file format

Postby alsayed » Wed Apr 11, 2012 4:30 am

Hello all,

I've solved the problem. I've written a simple fortran code that computes the File we used in metgrid.exe.
------
program sample_read
! Fortran 90 version.

! This is a simple program to write data in the WPS intermediate
! format. It is included mainly as an aid to understanding said
! format.

! implicit none
include '/usr/local/netcdf-4.1.3-gcc45/include/netcdf.inc'
real*4, dimension(128,64,10,365,5) :: VARESULT
real*4 , dimension(1:64) :: ALAT
real*4 , dimension(1:128) :: ALON
real*4 , dimension(1:365) :: ATIME
real*4 , dimension(10) :: PRESSURET
integer :: itim
character(len=19) , dimension(1:365) :: get_timeid
character(len=4) variable(5),varn(5)
character*32 :: filename
! WRF CODE DECLARATIONS

integer nl
integer, parameter :: IUNIT = 10
integer, parameter :: OUNIT = 11
integer :: ierr,ns
! format WPS=5
integer :: VERSION=5
! The time, in format "YYYY-MM-DD_HH:mm:ss
character(len=24) :: HDATE
! Forecast time (in hours) of the data in the slab
real :: XFCST = 0
! Start location of data. Could be "CENTER" or "SWCORNER". "SWCORNER" is typical.
character(len=8) :: STARTLOC = "SWCORNER"

! A field name.
character(len=9) :: FIELD
! A field name.
character(len=9) , dimension(1:5) :: FIELD_i

! Units describing the field in the slab.
character(len=25) :: UNITS
! Units describing the field in the slab.
character(len=25), dimension(1:5) :: UNITS_i

! Text description of the field in the slab.
character(len=46) :: DESC
! A field name.
character(len=46) , dimension(1:5) :: DESC_i

! Source of data
character(len=32) :: MAP_SOURCE = "Arpege"
! Pressure-level (Pa) of the data. 200100 Pa indicates surface data; 201300 Pa indicates sea-level pressure
real*4 XLVL
! Slab dimension in the X direction
integer :: NX = 128
! Slab dimension in the Y direction
integer :: NY = 64
! Flag denoting the projection, 0: Cylindrical Equidistant (Lat/lon) projection.
integer :: IPROJ = 0
! Starting latitude (degrees north)
real :: STARTLAT = -87.8638
! Starting longitude (degrees east)
real :: STARTLON = 0
! Latitude increment (degrees) for lat/lon grid
real :: DELTALAT = 2.7672
!Longitude increment (degrees) for lat/lon grid
real :: DELTALON = 2.8125
real :: EARTH_RADIUS = 6367470. * .001
! Logical flag use to indicate if Lambert projected data has earth or model rotated winds.
logical :: IS_WIND_EARTH_REL = .FALSE.

! SLAB is an allocatable array, because we do not necessarily know in
! advance the size of the array we need to read.
! real, allocatable, dimension(:,:) :: SLAB
! I decided to change the allocatable array feature for Fortran 90
! I get a segmentation fault when I try to pass the VARESULT to it
! I could also create a dummy slab, get the VARESULT and then pass it
! to the allocatable one. This takes a lot of memory!
! change made by Michel Mesquita - Nov 08

real , dimension(1:128,1:64) :: SLAB


! Reading PSL (PMSL in Metgrid - see Vtables) data

! Call subroutine I created to read the netcdf files
! read_netcdf(ALAT,ALON,ATIME,VARESULT)
print *, 'Reading subroutine sample_read...'
PRESSURET(1)=100000;PRESSURET(2)=92500;PRESSURET(3)=85000;PRESSURET(4)=70000;PRESSURET(5)=60000
PRESSURET(6)=50000;PRESSURET(7)=40000;PRESSURET(8)=30000;PRESSURET(9)=20000;PRESSURET(10)=10000

call read_netcdf_2(ALAT,ALON,ATIME,VARESULT)
!print *, 'VARESULT',VARESULT(:,:,:,:,5)
!print *, 'min est', minval(VARESULT(:,:,:,:,2)),maxval(VARESULT(:,:,:,:,2))
!print *, 'min est', minval(VARESULT(:,:,:,:,3)),maxval(VARESULT(:,:,:,:,3))
!print *, 'min est', minval(VARESULT(:,:,:,:,4)),maxval(VARESULT(:,:,:,:,4))
!print *, 'min est', minval(VARESULT(:,:,:,:,5)),maxval(VARESULT(:,:,:,:,5))
!print *, 'PRESSURET',PRESSURET

! Call the subroutine get_time(get_timeid) in order to get
! info for the WRF variable HDATE

call get_time(get_timeid)

! ------------------
variable(1)='zg_d';varn(1)='zg';FIELD_i(1)='GHT';DESC_i(1)='Height';UNITS_i(1)='m'
variable(2)='ta_d';varn(2)='ta';FIELD_i(2)='TT';DESC_i(2)='Temperature';UNITS_i(2)='K'
variable(3)='ua_d';varn(3)='ua';FIELD_i(3)='UU';DESC_i(3)='U';UNITS_i(3)='m s-1'
variable(4)='va_d';varn(4)='va';FIELD_i(4)='VV';DESC_i(4)='V';UNITS_i(4)='m s-1'
variable(5)='hur_';varn(5)='hur';FIELD_i(5)='RH';DESC_i(5)='Relative Humidity';UNITS_i(5)='%'
! -----------------
do nl=1,10
do ns=1,5
DATALOOP : DO itim = 1,365 ! Time dimension

! Declare some WRF variables
HDATE = get_timeid(itim)
! print *, 'The date is: ', HDATE
FIELD = FIELD_i(ns) !"SKINTEMP" ! As in Metgrid - see Vtables
UNITS = UNITS_i(ns) !"K"
DESC = DESC_i(ns) !"Temperature"

! Open and write a file name
XLVL=PRESSURET(nl)
SLAB = VARESULT(:,:,nl,itim,ns)
! print *, ns,itim
! if ( ( ns .eq. 5 ) .and. ( nl .eq. 1 ) .and. (itim .le. 10 ) ) then
! print *, 'relative humidity',nl,itim,SLAB,SLAB
! endif
! print *, 'SLAB',minval(SLAB),maxval(SLAB)
write (filename, 100) HDATE(1:10), HDATE(12:13)
100 format ('FILE:',a,'_',a)

open (IUNIT, FILE=filename, form='unformatted',access='append')
! open (IUNIT, FILE=filename, form='unformatted',access='sequential',status='unknown')
!!!!!! Start the WRF Int Format coding:
print *, 'writing the intermediate file format...', UNITS, DESC
write (IUNIT, IOSTAT=IERR) VERSION

! WRITE the second record, common to all projections:

write (IUNIT) HDATE, XFCST, MAP_SOURCE, FIELD, UNITS, DESC, XLVL, NX, NY, IPROJ

print*, HDATE//" ", XLVL, FIELD

! WRITE the third record, which depends on the projection:

if (IPROJ == 0) then

! This is the Cylindrical Equidistant (lat/lon) projection:
WRITE (IUNIT) STARTLOC, STARTLAT, STARTLON, DELTALAT, DELTALON, EARTH_RADIUS

! elseif (IPROJ == 1) then

! This is the Mercator projection:
! WRITE (IUNIT) STARTLOC, STARTLAT, STARTLON, DX, DY, TRUELAT1, EARTH_RADIUS
! elseif (IPROJ == 3) then

! This is the Lambert Conformal projection:
! WRITE (IUNIT) STARTLOC, STARTLAT, STARTLON, DX, DY, XLONC, TRUELAT1, TRUELAT2, EARTH_RADIUS


! elseif (IPROJ == 4) then

! Gaussian projection
! WRITE (IUNIT) STARTLOC, STARTLAT, STARTLON, NLATS, DELTALON, EARTH_RADIUS

! elseif (IPROJ == 5) then
! This is the Polar Stereographic projection:
! WRITE (IUNIT) STARTLOC, STARTLAT, STARTLON, DX, DY, XLONC, TRUELAT1, EARTH_RADIUS

endif


WRITE (IUNIT) IS_WIND_EARTH_REL


WRITE (IUNIT) SLAB
close (IUNIT)

! Loop back to read/write the next field.


ENDDO DATALOOP
enddo
enddo
write(*,'(/,"End of read loop. Program finished.")')

print *, 'fin_values',minval(SLAB),maxval(SLAB)
end program sample_read
! ***************************************************************
! * Subroutine to read a text file with the TIME strings *
! * Created by Michel Mesquita - Nov $lat_id$8 *
! * Bjerknes Centre for Climate Research *
! ***************************************************************
subroutine get_time(get_timeid)
implicit none

character (19) , dimension(1:365), intent(inout) :: get_timeid
character (10) :: A
character (8) :: B
integer i

print *, 'reading subroutine get_time...'

do i=1,365
open(unit=5, file='date_file')
read(5,*) A, B
get_timeid(i) = A // ' ' // B ! concatenate the strings
! print *, get_timeid
end do
close(5)


end subroutine get_time

subroutine read_netcdf_2(lat,lon,time,vals)
implicit none

include '/usr/local/netcdf-4.1.3-gcc45/include/netcdf.inc'
integer ncid,ierr,dimidx,dimidy,dimidz,dimidt,varid,varidx,varidy
integer varidz,varidt,startA(1),countA(1),startB(4),countB(4)
integer ns,nv,lenx,leny,lenz,var0,lent,niv,err
integer n1,n2,n3
character(len=4) variable(13),varn(13),field_i(13)
real*4 , dimension(1:64) :: lat
real*4 , dimension(1:128) :: lon
real*4 , dimension(1:365) :: time
character(len=70) cha1
character(len=16) cha2
character(len=88) cha_0
character(len=89) cha
character(len=90) cha_1
! real*4 valsx(128),valsy(64),valsz(10),valst(365),vals(365,10,64,128,5)
real*4 valsx(128),valsy(64),valsz(10),valst(365),vals(128,64,10,365,5)
integer, parameter :: version=5 ! Format version (must =5 for WPS format)
integer, parameter :: nx=128
integer, parameter :: ny=64 ! x- and y-dimensions of 2-d array
integer, parameter :: lev=10
integer, parameter :: iproj=0 ! Code for projection of data in array:
! 0 = cylindrical equidistant
! 1 = Mercator
! 3 = Lambert conformal conic
! 4 = Gaussian (global only!)
! 5 = Polar stereographic
real, parameter :: nlats=0 ! Number of latitudes north of equator
! (for Gaussian grids)
real, parameter :: startlat=-87.8638
real, parameter :: startlon=0 ! Lat/lon of point in array indicated by
! startloc string
real, parameter :: deltalat=2.7672
real, parameter :: deltalon=2.8125 ! Grid spacing, degrees
real, parameter :: dx=1.181168
real, parameter :: dy=1.125966 ! Grid spacing, km
real, parameter :: xlonc=0 ! Standard longitude of projection
real, parameter :: truelat1=0
real, parameter :: truelat2=0 ! True latitudes of projection
real, parameter :: earth_radius=6367470. * .001 ! Earth radius, km
logical, parameter :: is_wind_grid_rel=.false. ! Flag indicating whether winds are
! relative to source grid (TRUE) or
! relative to earth (FALSE)
character (len=8), parameter :: startloc='SWCORNER' ! Which point in array is given by
! startlat/startlon; set either
! to 'SWCORNER' or 'CENTER '
!
!
! A MODIFIER ICI LE VARIABLE
!
!
character(len=8) :: field
real, parameter :: sst_intv = 24. ! input SST interval hours
!
! dynamic output variables
!
real :: xfcst ! Forecast hour of data (hours from ifile1)
character(len=24) :: hdate='1979:01:01_00:00:00' ! Valid date for data YYYY:MM:DD_HH:00:00
!
! other variables
!
character(256) :: ifile1, ifile2, ofile
integer, parameter :: iunit1 = 10, iunit2 = 20, ounit = 100
call getarg(1, ifile1)
call getarg(2, ifile2)
call getarg(3, ofile)
call getarg(4, hdate)
xfcst=24

variable(1)='zg_d';varn(1)='zg';field_i(1)='GHT'
variable(2)='ta_d';varn(2)='ta';field_i(2)='TT'
variable(3)='ua_d';varn(3)='ua';field_i(3)='UU'
variable(4)='va_d';varn(4)='va';field_i(4)='VV'
variable(5)='hur_';varn(5)='hur';field_i(5)='RH'

n1=4
n2=5
n3=5
nv=5
!* Declare functions and variables associated with the netCDF library
! integer NF_OPEN,NF_NOWRITE,NF_INQ_DIMID,NF_INQ_VARID,NF_INQ_DIMLEN
! integer NF_GET_VARA_REAL,NF_CLOSE


!* 1. Open the netCDF file
!* -----------------------
do ns=1,nv
field=field_i(ns)
cha1='/work/crct/mo9378al/ARPEGE/arpege-ieee_grib_netcdf/ARPEGE/sorties/EH1/'
if (ns .le. n1) then
cha=cha1//variable(ns)//'.PLEH1.1979.nc'
elseif ((ns .gt. n1) .and. (ns .le. n2)) then
cha=cha1//variable(ns)//'d.PLEH1.1979.nc'
elseif ((ns .gt. n2) .and. (ns .le. n3)) then
cha=cha1//variable(ns)//'_d.PLEH1.1979.nc'
endif
ierr=NF_OPEN(cha,NF_NOWRITE,ncid)
print *, ierr,ncid,cha
!*
!* 2. Get dimension IDs
!* --------------------
ierr=NF_INQ_DIMID(ncid,'lon',dimidx)
ierr=NF_INQ_DIMID(ncid,'lat',dimidy)
ierr=NF_INQ_DIMID(ncid,'plev',dimidz)
ierr=NF_INQ_DIMID(ncid,'time',dimidt)

!*
!* 3. Get dimension lengths
!* ------------------------
ierr=NF_INQ_DIMLEN(ncid,dimidx,lenx)
ierr=NF_INQ_DIMLEN(ncid,dimidy,leny)
ierr=NF_INQ_DIMLEN(ncid,dimidz,lenz)
ierr=NF_INQ_DIMLEN(ncid,dimidt,lent)
print*,'Dimension lengths',lenx,leny,lenz,lent

!* 4. Get variable IDs
!* -------------------
!* First the variables for the dimensions
ierr=NF_INQ_VARID(ncid,'lon',varidx)
ierr=NF_INQ_VARID(ncid,'lat',varidy)
ierr=NF_INQ_VARID(ncid,'plev',varidz)
ierr=NF_INQ_VARID(ncid,'time',varidt)
!
! IMPORTANT, IMPORTANT
!* ...and for the main data variable
!
ierr=NF_INQ_VARID(ncid,varn(ns),varid)
!*
!* 5. Get contents of dimension variables
!* --------------------------------------
startA(1)=1
countA(1)=lenx
ierr=NF_GET_VARA_REAL(ncid,varidx,startA,countA,valsx)
startA(1)=1
countA(1)=leny
ierr=NF_GET_VARA_REAL(ncid,varidy,startA,countA,valsy)
startA(1)=1
countA(1)=lenz
ierr=NF_GET_VARA_REAL(ncid,varidz,startA,countA,valsz)
startA(1)=1
countA(1)=lent
ierr=NF_GET_VARA_REAL(ncid,varidt,startA,countA,valst)
!*
!* 6. Get values of the main variable
!* ----------------------------------
startB(1)=1
startB(2)=1
startB(3)=1
startB(4)=1
countB(1)=lenx
countB(2)=leny
countB(3)=lenz
countB(4)=lent
! print*,'Dimension lengths',lenx,leny,lenz,lent
! print *,'------------------vals_x est',valsx
! print *, '--------------- vals_y est', valsy
! print *, '--------------- vals_z est', valsz
! print *, '--------------- vals_t est', valst
ierr=NF_GET_VARA_REAL(ncid,varid,startB,countB,vals(:,:,:,:,ns))

! print *, 'here_here_here',minval(vals),maxval(vals),vals(:,:,:,:,ns)
!* 7. Close the file up
!* --------------------
err=NF_CLOSE(ncid)
enddo

end

-------------

Mouhamad
alsayed
 
Posts: 23
Joined: Fri Nov 18, 2011 7:30 am

Re: Metgrid, PRES, and the intermediate file format

Postby Louise » Wed Jan 16, 2013 9:59 am

Hi,

Thanks a lot for posting this code. I modified it for my case and it worked really well, after I discovered that it needs to be compiled with the same compiler and compiler options as metgrid itself.
Louise
 
Posts: 4
Joined: Wed Dec 05, 2012 3:51 pm

Re: Metgrid, PRES, and the intermediate file format

Postby hnlim » Fri Mar 29, 2013 3:54 pm

Does anyone know where I can find a list of all variables that I need to write into the intermediate format because my meteorological data are not in grib format?

Thanks
Agnes
hnlim
 
Posts: 74
Joined: Wed Apr 16, 2008 11:51 am

Re: Metgrid, PRES, and the intermediate file format

Postby xserris » Thu Jan 15, 2015 9:46 pm

amercer.msstate, I wanted to try to run a simulation using MERRA but cannot find a Vtable for such. You said earlier that you made a Vtable for using MERRA, could you post/send the Vtable or tell me how to construct my own?
xserris
 
Posts: 1
Joined: Thu Jan 15, 2015 9:38 pm

Re: Metgrid, PRES, and the intermediate file format

Postby UCDavisStudent » Wed Feb 18, 2015 3:44 pm

Hey all,

I am also trying to use MERRA reanalysis data as an input for WRF. Has anyone had any luck doing this?

Also which data product from the MERRA data holdings (http://disc.sci.gsfc.nasa.gov/daac-bin/DataHoldings.pl) did you use?

Thank you,
UCDavisStudent
UCDavisStudent
 
Posts: 1
Joined: Wed Feb 18, 2015 3:41 pm

Re: Metgrid, PRES, and the intermediate file format

Postby dstanfel1 » Mon Aug 03, 2015 6:47 pm

Hello,

I am trying to use data from the ESM2M earth system model as initial conditions for the WRF. I have tried to use the sample program code posted here but do not really understand it. When I tried to compile it using gfortran as the compiler I got the following error:

dstanfel:WRF_MODEL ninjawarrior343$ gfortran -v NETCDF_metgrid.f
Driving: gfortran -mmacosx-version-min=10.10.3 -v NETCDF_metgrid.f -l gfortran -shared-libgcc
Using built-in specs.
COLLECT_GCC=gfortran
COLLECT_LTO_WRAPPER=/opt/local/libexec/gcc/x86_64-apple-darwin14/4.9.2/lto-wrapper
Target: x86_64-apple-darwin14
Configured with: /opt/local/var/macports/build/_opt_mports_dports_lang_gcc49/gcc49/work/gcc-4.9.2/configure --prefix=/opt/local --build=x86_64-apple-darwin14 --enable-languages=c,c++,objc,obj-c++,lto,fortran,java --libdir=/opt/local/lib/gcc49 --includedir=/opt/local/include/gcc49 --infodir=/opt/local/share/info --mandir=/opt/local/share/man --datarootdir=/opt/local/share/gcc-4.9 --with-local-prefix=/opt/local --with-system-zlib --disable-nls --program-suffix=-mp-4.9 --with-gxx-include-dir=/opt/local/include/gcc49/c++/ --with-gmp=/opt/local --with-mpfr=/opt/local --with-mpc=/opt/local --with-isl=/opt/local --disable-isl-version-check --with-cloog=/opt/local --disable-cloog-version-check --enable-stage1-checking --disable-multilib --enable-lto --enable-libstdcxx-time --with-as=/opt/local/bin/as --with-ld=/opt/local/bin/ld --with-ar=/opt/local/bin/ar --with-bugurl=https://trac.macports.org/newticket --with-pkgversion='MacPorts gcc49 4.9.2_2'
Thread model: posix
gcc version 4.9.2 (MacPorts gcc49 4.9.2_2)
COLLECT_GCC_OPTIONS='-mmacosx-version-min=10.10.3' '-v' '-shared-libgcc' '-mtune=core2'
/opt/local/libexec/gcc/x86_64-apple-darwin14/4.9.2/f951 NETCDF_metgrid.f -ffixed-form -fPIC -quiet -dumpbase NETCDF_metgrid.f -mmacosx-version-min=10.10.3 -mtune=core2 -auxbase NETCDF_metgrid -version -fintrinsic-modules-path /opt/local/lib/gcc49/gcc/x86_64-apple-darwin14/4.9.2/finclude -o /var/folders/2_/gmrv3l3d7bx42vz02qzfk1nr0000gn/T//ccNdyDKl.s
GNU Fortran (MacPorts gcc49 4.9.2_2) version 4.9.2 (x86_64-apple-darwin14)
compiled by GNU C version 4.9.2, GMP version 6.0.0, MPFR version 3.1.2-p10, MPC version 1.0.3
GGC heuristics: --param ggc-min-expand=100 --param ggc-min-heapsize=131072
GNU Fortran (MacPorts gcc49 4.9.2_2) version 4.9.2 (x86_64-apple-darwin14)
compiled by GNU C version 4.9.2, GMP version 6.0.0, MPFR version 3.1.2-p10, MPC version 1.0.3
GGC heuristics: --param ggc-min-expand=100 --param ggc-min-heapsize=131072

gfortran: internal compiler error: Killed: 9 (program f951)
Abort trap: 6

Has anyone encountered this error before? Also, has anyone been able to convert netcdf data directly into the intermediate file format for my case? Any help would be greatly appreciated in understanding how to do this conversion.

Thanks,

David
dstanfel1
 
Posts: 5
Joined: Thu Jun 25, 2015 12:23 am

Re: Metgrid, PRES, and the intermediate file format

Postby sourabh » Mon Jul 29, 2019 1:26 am

I used the same code but for single time (time step = 1)
and I ma getting error:-

"""reading subroutine get_time...
At line 211 of file test.f90 (unit = 5, file = 'date_file')
Fortran runtime error: End of file """"

Error termination. Backtrace:
#0 0x1514c3bee31a
#1 0x1514c3beeec5
#2 0x1514c3bef68d
#3 0x1514c3d65a33
#4 0x1514c3d5e9c4
#5 0x1514c3d600f9
#6 0x563fe85d6055
#7 0x563fe85d6c90
#8 0x563fe85d7cc4
#9 0x1514c3803b96
#10 0x563fe85d5df9
#11 0xffffffffffffffff


So I am gone through the line mentioned in code

subroutine get_time(get_timeid)
implicit none

character (19) , dimension(1:365), intent(inout) :: get_timeid
character (10) :: A
character (8) :: B
integer i

print *, 'reading subroutine get_time...'

do i=1,365
open(unit=5, file='date_file')
read(5,*) A, B
get_timeid(i) = A // ' ' // B ! concatenate the strings
print *, get_timeid
end do
close(5)

so help me out to get out from this error.
sourabh
 
Posts: 1
Joined: Mon Jul 29, 2019 1:21 am

Previous

Return to Working with Various Datasets

Who is online

Users browsing this forum: No registered users and 3 guests

cron