Purpose: To perform basic extrapolations resulting in continuous fields of NASA MERRA reanalyses for selected model levels where pressure values are greater than the surface pressure (i.e., filling in model grid points "below ground" with defined values). This could be useful where a continuous global field is needed, such as in modeling applications or spherical harmonic analyses. Code samples are included, below. Geopotential heights are extrapolated by applying the hydrostatic equation and equation of state, as are temperature fields using adiabatic lapse rate approximations (moist or dry).
Background reference from the MERRA FAQs (https://gmao.gsfc.nasa.gov/research/merra/faq.php)
Q. 6 "Why are there such large discrepancies at 1000MB and 850MB between MERRA and other Reanalyses?"
A. "The GEOS5 data assimilation system used to produce MERRA does not (or did not at the time of production) extrapolate data to pressure levels greater than the surface pressure. These grid points are marked by undefined values. The result is that area averages that include these points will not be representative compared to other data sets without additional screening. Time averages, such as monthly means, may also have substantial differences at the edges of topography. The lowest model level data and surface data are available so that users can produce their own extrapolation. A page discussing this issue is available. See http://gmao.gsfc.nasa.gov/research/merra/pressure_surface.php. "
Below, note regions corresponding to high model topography where contours end due to undefined values in a sample MERRA Reanalysis plot (for example, portions of Rocky Mountains, Andes Mountains, Tibetan Plateau, Antarctica):
With extrapolation, the values are "filled in", as this figure shows (below):
Here is another example with Geopotential height:
And, after extrapolation:
Below are sample codes and scripts that perform extrapolation on 5 variables from a MERRA NetCDF file.
Software needed: GrADS, Fortran compiler, Unix/Linux shell environment.
Need a MERRA NetCDF or HDF file along with the following:
1.Shell executable extrap.csh: Runs programs, creates output. Calls grads scripts and fortran executable
***********************************
#!/bin/csh
grads -bpc 'fill_merra_9999.gs'
#output tot_merra_filled_9999.prs
#view output with check_filled_9999.ctl
#Next run Fortran77 executable to do extrapolation:
set infile="tot_merra_filled_9999.prs"
set outfile="tot_extrapolated.prs"
#
#Note this code multiplies RH*100
extrap_merra << input
$infile
$outfile
input
#Use this grads control file to check on extrapolated fields:
#check_extrapolated.ctl
echo All done with extrapolation
*********************************************
2. Two GrADS control files are needed:
check_filled_9999.ctl (views output w/-9999 set as undefined):
DSET ^tot_merra_filled_9999.prs
undef -9999
TITLE MERRA Analyses
XDEF 288 linear -179.375 1.25
YDEF 144 linear -89.375 1.25
ZDEF 17 LEVELS 1000. 925. 850. 700. 600. 500. 400. 300. 250. 200. 150. 100. 70. 50. 30. 20. 10.
TDEF 8 LINEAR 00Z10Jan2010 3hr
VARS 5
h 17 99 geop height (m)
t 17 99 temperature (K)
u 17 99 u-wind (m/s)
v 17 99 v-wind (m/s)
rh 17 99 relative humidity (*100 is %)
ENDVARS
check_extrapolated.ctl (views extrapolated output):
DSET ^tot_extrapolated.prs
undef -9999
TITLE Extrapolated MERRA Analyses
XDEF 288 linear -179.375 1.25
YDEF 144 linear -89.375 1.25
ZDEF 17 LEVELS 1000. 925. 850. 700. 600. 500. 400. 300. 250. 200. 150. 100. 70. 50. 30. 20. 10.
TDEF 8 LINEAR 00Z10Jan2010 3hr
VARS 5
h 17 99 geop height (m)
t 17 99 temperature (K)
u 17 99 u-wind (m/s)
v 17 99 v-wind (m/s)
rh 17 99 relative humidity (%)
ENDVARS
3. GrADS script fill_merra_9999.gs:
fill_merra_9999.gs is a GrADS script called in extrap.csh.
It specifies "undefined" w/ -9999 using const() and maskout() functions.
A netcdf file in the script was downloaded from the NASA MERRA server
MERRA300.prod.assim.inst3_3d_asm_Cp.20100110.SUB_1.nc
The file contains 17 levels, 3-hourly output, for one day (8 output times):
The file contains the 5 variables that will be saved to a new file with -9999 as the new undefined value.
Variables are geopotential height, temperature, u-wind, v-wind, and RH.
The resulting binary file is tot_merra_filled_9999.prs, viewed with check_filled_9999.ctl
Here is the code for fill_merra_9999.gs:
*Run in extrap.csh using the command:
*grads -bpc 'fill_merra_9999.gs'
'reinit'
'sdfopen MERRA300.prod.assim.inst3_3d_asm_Cp.20100110.SUB_1.nc'
'set x 1 288'
'set y 1 144'
'set gxout fwrite'
'set fwrite tot_merra_filled_9999.prs'
t=1
while (t <= 8)
'set t 't
say 'Analysis time is 't
*height
z=1
while (z<=17)
'set z 'z
'd const(maskout(h,h),-9999.,-u)'
z=z+1
endwhile
*
*Temperature
z=1
while (z<=17)
'set z 'z
'd const(maskout(t,t),-9999.,-u)'
z=z+1
endwhile
*
*u-wind
z=1
while (z<=17)
'set z 'z
'd const(maskout(u,u+1000),-9999.,-u)'
z=z+1
endwhile
*
*v-wind
z=1
while (z<=17)
'set z 'z
'd const(maskout(v,v+1000),-9999.,-u)'
z=z+1
endwhile
*
*RH
z=1
while (z<=17)
'set z 'z
'd const(maskout(rh,rh),-9999.,-u)'
z=z+1
endwhile
*
t = t + 1
endwhile
'disable fwrite'
'quit'
4. Fortran file extrap_merra.f:
Code to perform extrapolation, called in extrap.csh.
Sample compile (Mac Intel):
ifort -g -mssse3 -assume byterecl -save -extend-source -o extrap_merra extrap_merra.f
Requires input of 5 variables (see above, h,t,u,v,rh)
Replaces -9999. values in MERRA binary (tot_merra_filled_9999.prs) w/ new, extrapolated data
Wind (u,v), Rel Humidity (rh) use value of grid point lowest to terrain
Temperature (t) extrapolated using adiabatic lapse rate (moist or dry), dT/DP=const
Geop Height (h) extrapolated by hydrostatic eqn and Eqn of State: dZ/DP=-RT/Pg
Resulting binary tot_merra_filled_9999.prs is viewed with check_extrapolated.ctl
program extrap_merra
c-9-7-12
clab purpose is to fill in -9999. values in a MERRA binary w/extrapolated data
clab u,v,rh just use value of grid point lowest to terrain
clab T extrapolated using adiabatic lapse rate, dT/DP=const, H extrapolated by hydrostatic eqn/eqn of state:
clab dZ/DP=-RT/Pg
c-9-7-12
parameter(ni=288,nj=144,nk=17,nt=8) ! array dimensions, nt analysis times
dimension var(ni,nj),u(ni,nj,nk),vgrads(ni,nj)
dimension h(ni,nj,nk),t(ni,nj,nk),v(ni,nj,nk),rh(ni,nj,nk)
dimension pout(nk),tt(ni,nj,nk),rhh(ni,nj,nk)
dimension hh(ni,nj,nk),uu(ni,nj,nk),vv(ni,nj,nk)
character*120 infile,outfile
data pout/10.,20.,30.,50.,70.,100.,150.,200.,
+ 250.,300.,400.,500.,600.,700.,850.,925.,1000./ !Reanalysis levels
clab read out data for p levels:
do k=1,nk
pout(k)=pout(k)*100.
write(6,*) pout(k)
enddo
clab binary file from merra whose values "below ground" have been
clab filled in with -9999. using grads script fill_merra_9999.gs
read(*,'(A120)')infile
open(UNIT=1,FILE=infile,
+access='direct',type='old',form='unformatted',recl=288*144*4)
clab output file:
read(*,'(A120)')outfile
open(UNIT=17,FILE=outfile,
+access='direct',form='unformatted',recl=288*144*4)
irec1=1
irec2=1
do ntime=1,nt
do nvar=1,5
write(6,*)nvar
if(nvar.eq.1) then
c Read out all data and flip k-index (to aid in filling in from top-down)
do k=1,nk
read(1,rec=irec1)vgrads !height
do i=1,ni
do j=1,nj
h(i,j,nk+1-k)=vgrads(i,j)
enddo
enddo
irec1=irec1+1
enddo
elseif(nvar.eq.2)then !temperature
do k=1,nk
read(1,rec=irec1)vgrads
do i=1,ni
do j=1,nj
t(i,j,nk+1-k)=vgrads(i,j)
enddo
enddo
irec1=irec1+1
enddo
elseif(nvar.eq.3)then !u-wind
do k=1,nk
read(1,rec=irec1)vgrads
do i=1,ni
do j=1,nj
u(i,j,nk+1-k)=vgrads(i,j)
enddo
enddo
irec1=irec1+1
enddo
elseif(nvar.eq.4)then !v-wind
do k=1,nk
read(1,rec=irec1)vgrads
do i=1,ni
do j=1,nj
v(i,j,nk+1-k)=vgrads(i,j)
enddo
enddo
irec1=irec1+1
enddo
else !RH
do k=1,nk
read(1,rec=irec1)vgrads
do i=1,ni
do j=1,nj
rh(i,j,nk+1-k)=vgrads(i,j)
enddo
enddo
irec1=irec1+1
enddo
endif
enddo
*
c Now extrapolate temperature per lapse rate dt/dp
r=287.0
g=9.81
c Assume dt/dp=const, where const is such that dt/dz=6.5C/1km (moist adiab approx)
c This is similar to dt/dp=6.5C/100 mb for lowest couple of km
gammadry=9.8/10000.0
gammamoi=6.5/10000.0
do k=1,nk
do i=1,ni
do j=1,nj
tt(i,j,k)=t(i,j,k)
if(tt(i,j,k).le.-9999.) tt(i,j,k)=tt(i,j,k-1) + gammamoi*(pout(k)-pout(k-1))
c if(tt(i,j,k).le.-9999.) tt(i,j,k)=tt(i,j,k-1) + gammadry*(pout(k)-pout(k-1))
enddo
enddo
enddo
c Now extrapolate to find geopotential height below ground
c Use dz/dp=-RT/Pg to get geopotential height field z at pressure level p
c working downward from the lowest level above the surface
r=287.0
g=9.81
do k=1,nk
do i=1,ni
do j=1,nj
hh(i,j,k)=h(i,j,k)
if(hh(i,j,k).le.-9999.) hh(i,j,k)=hh(i,j,k-1) - ((r*tt(i,j,k))/(pout(k)*g))*(pout(k)-pout(k-1))
enddo
enddo
enddo
c Now extrapolate to find u-wind below ground
do k=1,nk
do i=1,ni
do j=1,nj
uu(i,j,k)=u(i,j,k)
if(uu(i,j,k).le.-9999.) uu(i,j,k)=uu(i,j,k-1)
enddo
enddo
enddo
c Now extrapolate to find v-wind below ground
do k=1,nk
do i=1,ni
do j=1,nj
vv(i,j,k)=v(i,j,k)
if(vv(i,j,k).le.-9999.) vv(i,j,k)=vv(i,j,k-1)
enddo
enddo
enddo
c Now extrapolate to find RH below ground
do k=1,nk
do i=1,ni
do j=1,nj
rhh(i,j,k)=rh(i,j,k)
if(rhh(i,j,k).le.-9999.) rhh(i,j,k)=rhh(i,j,k-1)
enddo
enddo
enddo
C Now reverse k-index to write out all variables from sfc to top
do k=1,nk !geopotential height
do i=1,ni
do j=1,nj
var(i,j)=hh(i,j,nk+1-k)
enddo
enddo
write(17,rec=irec2)var
irec2=irec2+1
enddo
do k=1,nk !temperature
do i=1,ni
do j=1,nj
var(i,j)=tt(i,j,nk+1-k)
enddo
enddo
write(17,rec=irec2)var
irec2=irec2+1
enddo
do k=1,nk !u-wind
do i=1,ni
do j=1,nj
var(i,j)=uu(i,j,nk+1-k)
enddo
enddo
write(17,rec=irec2)var
irec2=irec2+1
enddo
do k=1,nk !v-wind
do i=1,ni
do j=1,nj
var(i,j)=vv(i,j,nk+1-k)
enddo
enddo
write(17,rec=irec2)var
irec2=irec2+1
enddo
do k=1,nk ! RH and multiply x100 to get percent
do i=1,ni
do j=1,nj
var(i,j)=rhh(i,j,nk+1-k)*100.
enddo
enddo
write(17,rec=irec2)var
irec2=irec2+1
enddo
c End iteration of time:
enddo
end
Contact: Lee Byerle
Re: Atmospheric Reanalyses – Recent Progress and Prospects ...