RegCM postproc

There are 4 steps for the post-processing of RegCM output:

Domain information and 3 methods for map projection
ATM, SRF, RAD variables
Monthly mean, Diurnal mean
sigma to P converter

1. Domain information and 3 methods for map projection

After you have your first successful RegCM results, you could go to the output directory (PracticeRun/output)

cd PracticeRun/output

( or cd /scratch/RegCM/PracticeRun/output ),

There are several files under PracticeRun/output directory:

OUT_HEAD ! Domain information datasets
OUT_HEAD.CTL ! GrADS descriptor file for OUT_HEAD
ATM.1994070100 ! Normal model output, U,V,T,Q,PS, etc
ATM.1994070100.ctl ! GrADS descriptor file for ATM.1994070100
SRF.1994070100 ! Land surface model output (surface variables)
SRF.1994070100.ctl ! GrADS descriptor file for SRF.1994070100
RAD.1994070100 ! Radiation output variables
RAD.1994070100.ctl ! GrADS descriptor file for RAD.1994070100
SAV.1994070600 ! Datasets saved for restart run

To check Domain information, you start up GrADS and enter:

    open OUT_HEAD.CTL
    d xlat
    d xlon

and
    c
    set gxout shaded
    d ht

Now enter:
    c
    set mproj lambert
    set gxout shaded
    d ht
    set gxout contour
    d ht

In OUT_HEAD.CTL, the following line define the map projection
(if you use LAMBERT map projection)
pdef 49 32 lcc 45.39 13.48 24.50 16.00 30.00 60.00 13.48 60000. 60000.
here,
lcc ! LAMBERT projection
49 ! your domain size in west-east direction (jx-2)
32 ! your domain size in south-north direction (iy-2)
45.39 ! latitude of the center point
13.48 ! longitude of the center point
24.50 ! 49/2
16.00 ! 32/2
30.00 ! true latitude 1, 30N
60.00 ! true latitude 2, 60N
13.48 ! longitude of the center point
60000. ! horizontal resolution in meter, 60km
60000. ! horizontal resolution in meter, 60km

please see details of pdef at:
http://grads.iges.org/grads/gadoc/map.html#lambert

The 10 variables in OUT_HEAD are:

ht ! surface elevation
htsd ! surface elevation standard deviation
veg2d ! vegetation type in BATS
landuse ! surface landuse type
xlat ! latitude of cross points
xlong ! longitude of cross points
xmap ! map factors of cross point
dmap ! map factors of dot points
coriol ! coriol force
mask ! land/sea mask

To work better for later, please try the following practice:

cp OUT_HEAD.CTL out_head.ctl

and edit out_head.ctl and make it as the following way:

dset ^OUT_HEAD
title RegCM domain information
options big_endian
undef -1.e34
xdef 49 linear -2.5696 0.6659 ! this line changed
ydef 32 linear 35.7379 0.5399 ! this line changed
zdef 1 levels 1000.00
tdef 1 linear 00z01Jan2001 1mo
vars 11
head 0 99 header information
ht 0 99 surface elevation
htsd 0 99 surface elevation std dev
veg2d 0 99 vegetation type in BATS
landuse 0 99 surface landuse type
xlat 0 99 latitude of cross points
xlong 0 99 longitude of cross points
xmap 0 99 map factors of cross point
dmap 0 99 map factors of dot points
coriol 0 99 coriol force
mask 0 99 land/sea mask
endvars

here you have removed the pdef line and reset xdef and ydef as
normal way, you will practice to use XLAT, XLONG and MASK of RegCM
output to make the map background by starting GrADS and enter:

open out_head.ctl
    set vpage 0.2 10.8 0.2 8.3
    set mproj off
    set xlab off
    set ylab off
    set grads off
    set ccolor rainbow
    set clab on
    set cthick 2
    set cstyle 1
    set gxout contour
    d ht
    set clab off
    set gxout contour
    set grads off
    set ccolor 10
    set cthick 1
    set cstyle 1
    set clevs 1.
    d mask
    set cint 5
    set ccolor 10
    set cthick 2
    set cstyle 5
    d xlat
    set cint 10
    set ccolor 10
    set cthick 2
    set cstyle 5
    d xlong
    set strsiz 0.12
    set string 1 l 4 0
    draw string 0.05 6.6 50N
    draw string 0.05 4.6 45N
    draw string 0.05 2.5 30N
    draw string 1.22 0.6 0
    draw string 4.15 0.6 10E
    draw string 7.2 0.6 20E
    draw string 10.4 0.6 30E
    draw title Topography, RegCM

Not bad, right ? We need to to improve the land-sea mask later by using a high resolution land-sea mask (let's say 5 times denser than the one you used).

After you tune the above mapping right, you could copy them to several GrADS script (gs) files: mapoff.gs and mask.gs .

In mapoff.gs file, it contains:

'set mproj off'
'set xlab off'
'set ylab off'
'set grads off'
'set ccolor rainbow'
'set clab on'
'set cthick 2'
'set cstyle 1'

and in mask.gs file, it contains:

'set clab off'
'set gxout contour'
'set grads off'
'set ccolor 10'
'set cthick 1'
'set cstyle 1'
'set clevs 1.'
'd mask'
'set cint 5'
'set ccolor 10'
'set cthick 2'
'set cstyle 5'
'd xlat'
'set cint 10'
'set ccolor 10'
'set cthick 2'
'set cstyle 5'
'd xlong'
'set strsiz 0.12'
'set string 1 l 4 0'
'draw string 0.05 6.6 50N'
'draw string 0.05 4.6 45N'
'draw string 0.05 2.5 30N'
'draw string 1.22 0.6 0'
'draw string 4.15 0.6 10E'
'draw string 7.2 0.6 20E'
'draw string 10.4 0.6 30E'

then start up GrADS, you just need type:

    open out_head.ctl
    set vpage 0.2 10.8 0.2 8.3
    run mapoff.gs
    set gxout contour
    d ht
    run mask.gs
    draw title Topography, RegCM

The above method could be also used for Rotated Mercator Map Projection.

So far, you have learned 3 methods to treat the LAMBERT Map Projection:

method-1: pdef with the default mproj (or 'set mproj latlon')
method-2: pdef and 'set mproj lambert'
method-3: no-pdef and use all model domain information.

2. ATM, SRF, RAD variables

If you have selected to use iotyp=2 (sequential binary output), because
the model output IDATE (date information) for the first writen output and pdef is just worked for iotyp=1 (direct access binary output), you have to use method-3 to view your RegCM output.

Before using GrADS, please:
cp mask.gs mask2.gs

and edit mask2.gs to become:
'set clab off'
'set gxout contour'
'set grads off'
'set ccolor 10'
'set cthick 1'
'set cstyle 1'
'set clevs 1.'
'd mask.2(t=1)' ! this line changed
'set cint 5'
'set ccolor 10'
'set cthick 2'
'set cstyle 5'
'd xlat.2(t=1)' ! this line changed
'set cint 10'
'set ccolor 10'
'set cthick 2'
'set cstyle 5'
'd xlong.2(t=1)' ! this line changed
'set strsiz 0.12'
'set string 1 l 4 0'
'draw string 0.05 6.6 50N'
'draw string 0.05 4.6 45N'
'draw string 0.05 2.5 30N'
'draw string 1.22 0.6 0'
'draw string 4.15 0.6 10E'
'draw string 7.2 0.6 20E'
'draw string 10.4 0.6 30E'

Now start up GrADS, enter:
open ATM.1994070100.ctl
open out_head.ctl

You may want to see what is in ATM.1994070100, so enter:
q file

You will see:
ga-> q file
File 1 : RegCM normal output variables
Descriptor: ATM.1994070100.ctl
Binary: ATM.1994070100
Type = Gridded
Xsize = 49 Ysize = 32 Zsize = 18 Tsize = 21
Number of Variables = 10
u 18 0 westerly wind (m/s)
v 18 0 southerly wind (m/s)
t 18 0 air temperature (degree)
qv 18 0 air specific humidity
qc 18 0 cloud water mixing ratio
psa 0 99 surface pressure (hPa)
rain 0 99 total precipitation
tgb 0 99 groud temperature in BATS
swt 0 99 total soil water in mm H2O
rno 0 99 accumulated infiltration
ga->

There are 5 3D variables (u, v, t, qv, and qc at sigma levels) and 5 2D variables (psa, rain, tgb, swt and rno) at surface.

If you want view surface pressure, enter:
    set vpage 0.2 10.8 0.2 8.3
    run mapoff.gs
    set gxout contour
    d psa
    run mask2.gs
    draw title Surface Pressure of 06Z, Jul 01

Now enter:
    c
    run mapoff.gs
    set gxout contour
    d ave(t,t=1,t=20)
    run mask2.gs
    draw title T at sigma=0.995

Now enter:
    c
    run mapoff.gs
    set gxout vector
    set t 20
    set z 17
    d skip(u,2,2);skip(v,2,2)
    run mask2.gs
    draw title U and V at lev 140hPa

Note: at this moment, all 3D variables are on sigma levels, you will convert them to P-level later.

By the way, in ATM.1994070100.ctl, the line of tdef,

tdef 21 linear 0z01jul1994 6hr

should be:
tdef 20 linear 6z01jul1994 6hr

You could view other fields in ATM.1994070100.
......

You may want to view variables in SRF.1994070100 by enter:
    reinit
    open SRF.1994070100.ctl
    open out_head.ctl
    set vpage 0.2 10.8 0.2 8.3

then if you want view precipitation on July 02, enter:
    run mapoff.gs
    set gxout shaded
    d ave(prca,t=10,t=17)
    set gxout contour
    d ave(prca,t=10,t=17)
    run mask2.gs
    draw title Precipitation of Jul 02

If you want to view the animation of incident solar flux, enter
    c
    run mapoff.gs
    set gxout shaded
    set t 1 40
    d sina
    run mask2.gs
    draw title Incident Solar Flux

Too fast animation, unfortunately. To have a better animation, you need to follow the example file at GrADS webpage:
http://grads.iges.org/grads/gadoc/animation.html

Now, to view air temperature at anemometer level (2m), enter:
    c
    run mapoff.gs
    set gxout contour
    set t 9
    d tanm-273.15
    run mask2.gs
    draw title Air Temperature at 2M

You could view other fields in SRF.1994070100.
......

And you could view other fields in RAD.1994070100.
......

By the way, tdef in SRF.1994070100.ctl and RAD.1994070100.ctl is the right one, the output frequency of
ATM.1994070100         6-hourly
SRF.1994070100         3-hourly
RAD.1994070100         6-hourly.

If your restart run (from 1994070600 to 1994071600) has already finished, you will have 7 files:

ATM.1994070600 ! Normal model output, U,V,T,Q,PS, etc
ATM.1994070600.ctl ! GrADS descriptor file for ATM.1994070600
SRF.1994070600 ! Land surface model output (surface variables)
SRF.1994070600.ctl ! GrADS descriptor file for SRF.1994070600
RAD.1994070600 ! Radiation output variables
RAD.1994070600.ctl ! GrADS descriptor file for RAD.1994070600
SAV.1994071600 ! Datasets saved for restart run

You could view those variables in which you have interest.
......

3. Monthly mean, Diurnal mean

Your work directory now is /scratch/RegCM/PracticeRun/,

    cp ../PostProc/postproc.f .
    cp ../PostProc/Makefile* .
    cp ../postproc.param .

You need to edit postproc.param and reset parameters according to your domain size and output frequency:

        integer nxf,nyf,nz
        parameter (nyf= 34) ! this line changed
        parameter (nxf= 51) ! this line changed
        parameter (nz= 18)
        real dtout,dtbat,dtrad
        parameter (dtout= 6.00)
        parameter (dtbat= 3.00)
        parameter (dtrad= 6.00)

And you could keep the other parameters as the default one.
Once ready,
pgf77 -byteswapio -o postproc postproc.f

Now edit postproc.in :

1994070100, ! idate0 = First date in File (yymmddhh)
1994071600, ! idate1 = Start date for averaging and re-writing
1994080100, ! idate2 = End date for averaging and re-writing
1., ! Number of Days for continual averaging (-1=monthly; 1=daily; 5=5 days)
.true. ! Is there a header file
'output/OUT_HEAD' ! filhead = Header File name
2 ! Number of simulation files listed
'output/ATM.1994070100' ! ATM File: 1
'output/ATM.1994080100' ! ATM File: 2
'output/SRF.1994070100' ! SRF File: 1
'output/SRF.1994080100' ! SRF File: 2
'output/RAD.1994070100' ! RAD File: 1
'output/RAD.1994080100' ! RAD File: 2
.false. ! Write out header?
***** ATM FILES *****
2 ! iotype: 1=I*2 NetCDF; 2=r*4 NetCDF; 3=grads; 4=vis5d
.false. ! Write out all ATM data btwn idate1 & idate2?
'ATM19940701_19940801.NC' ! ATM filename for rewrite
.true. ! Average ATM data btwn idate1 & idate2?
'ATM19940701_19940801AVG.NC' ! ATM average file to be written
.false. ! Avg Diurnal ATM data btwn idate1 & idate2?
'ATM19940701_19940801DIUR.NC' ! ATM diurnal average file to be written
.false. ! Continual avergage ATM data btwn idate1 & idate2?
'ATM19940701_19940801AVGMON.NC' ! ATM continual avg file to be written
***** SRF FILES *****
2 ! iotype: 1=I*2 NetCDF; 2=r*4 NetCDF; 3=grads; 4=vis5d
.false. ! Write out all SRF data btwn idate1 & idate2?
'SRF19940701_19940801.NC' ! SRF filename for rewrite
.true. ! Average SRF data btwn idate1 & idate2?
'SRF19940701_19940801AVG.NC' ! SRF average file to be written
.false. ! Avg Diurnal SRF data btwn idate1 & idate2?
'SRF19940701_19940801DIUR.NC' ! SRF diurnal average file to be written
.false. ! Continual avergage SRF data btwn idate1 & idate2?
'SRF19940701_19940801AVGMON.NC' ! SRF continual avg file to be written
***** RAD FILES *****
2 ! iotype: 1=I*2 NetCDF; 2=r*4 NetCDF; 3=grads; 4=vis5d
.false. ! Write out all RAD data btwn idate1 & idate2?
'RAD19940701_19940801.NC' ! RAD filename for rewrite
.true. ! Average RAD data btwn idate1 & idate2?
'RAD19940701_19940801AVG.NC' ! RAD average file to be written
.false. ! Avg Diurnal RAD data btwn idate1 & idate2?
'RAD19940701_19940801DIUR.NC' ! RAD diurnal average file to be written
.false. ! Continual avergage RAD data btwn idate1 & idate2?
'RAD19940701_19940801AVGMON.NC' ! RAD continual avg file to be written

Suppose you have successfully run both the initial run (from 1994070100 to 1994070600) and the restart run (from 1994070600 to 1994071600)

then you need to set:

1994070100, ! idate0 = First date in File (yymmddhh)
1994070100, ! idate1 = Start date for averaging and re-writing
1994071600, ! idate2 = End date for averaging and re-writing
1., ! Number of Days for continual averaging (-1=monthly; 1=daily; 5=5 days)
.true. ! Is there a header file
'output/OUT_HEAD' ! filhead = Header File name
2 ! Number of simulation files listed
'output/ATM.1994070100' ! ATM File: 1
'output/ATM.1994070600' ! ATM File: 2
'output/SRF.1994070100' ! SRF File: 1
'output/SRF.1994070600' ! SRF File: 2
'output/RAD.1994070100' ! RAD File: 1
'output/RAD.1994070600' ! RAD File: 2
.false. ! Write out header?
***** ATM FILES *****
3 ! iotype: 1=I*2 NetCDF; 2=r*4 NetCDF; 3=grads; 4=vis5d
.false. ! Write out all ATM data btwn idate1 & idate2?
'ATM19940701_19940716.dat' ! ATM filename for rewrite
.true. ! Average ATM data btwn idate1 & idate2?
'ATM19940701_19940716AVG.dat' ! ATM average file to be written
.false. ! Avg Diurnal ATM data btwn idate1 & idate2?
'ATM19940701_19940716DIUR.dat' ! ATM diurnal average file to be written
.false. ! Continual avergage ATM data btwn idate1 & idate2?
'ATM19940701_19940716AVGMON.dat' ! ATM continual avg file to be written
***** SRF FILES *****
3 ! iotype: 1=I*2 NetCDF; 2=r*4 NetCDF; 3=grads; 4=vis5d
.false. ! Write out all SRF data btwn idate1 & idate2?
'SRF19940701_19940716.dat' ! SRF filename for rewrite
.true. ! Average SRF data btwn idate1 & idate2?
'SRF19940701_19940716AVG.dat' ! SRF average file to be written
.false. ! Avg Diurnal SRF data btwn idate1 & idate2?
'SRF19940701_19940716DIUR.dat' ! SRF diurnal average file to be written
.false. ! Continual avergage SRF data btwn idate1 & idate2?
'SRF19940701_19940716AVGMON.dat' ! SRF continual avg file to be written
***** RAD FILES *****
3 ! iotype: 1=I*2 NetCDF; 2=r*4 NetCDF; 3=grads; 4=vis5d
.false. ! Write out all RAD data btwn idate1 & idate2?
'RAD19940701_19940716.dat' ! RAD filename for rewrite
.true. ! Average RAD data btwn idate1 & idate2?
'RAD19940701_19940716AVG.dat' ! RAD average file to be written
.false. ! Avg Diurnal RAD data btwn idate1 & idate2?
'RAD19940701_19940716DIUR.dat' ! RAD diurnal average file to be written
.false. ! Continual avergage RAD data btwn idate1 & idate2?
'RAD19940701_19940716AVGMON.dat' ! RAD continual avg file to be written

and save it, to do the average from 1994070100 to 1994071600, type:
./postproc

Now you may have your average fields in
ATM19940701_19940716AVG.dat
SRF19940701_19940716AVG.dat
RAD19940701_19940716AVG.dat

then you need to copy the sample GrADS descriptor file from
RegCM/Commons/tools/ directory:

    cp ../Commons/tools/ATM19940701_19940801AVG.ctl ATM19940701_19940716AVG.ctl
    cp ../Commons/tools/SRF19940701_19940801AVG.ctl SRF19940701_19940716AVG.ctl
    cp ../Commons/tools/RAD19940701_19940801AVG.ctl RAD19940701_19940716AVG.ctl

and change the dset name in these ctl files accordingly.

Now you could start up GrADS to view your average fields.

To do diurnal average, you need to reedit postproc.in,

1994070100, ! idate0 = First date in File (yymmddhh)
1994070100, ! idate1 = Start date for averaging and re-writing
1994071600, ! idate2 = End date for averaging and re-writing
1., ! Number of Days for continual averaging (-1=monthly; 1=daily; 5=5 days)
.true. ! Is there a header file
'output/OUT_HEAD' ! filhead = Header File name
2 ! Number of simulation files listed
'output/ATM.1994070100' ! ATM File: 1
'output/ATM.1994070600' ! ATM File: 2
'output/SRF.1994070100' ! SRF File: 1
'output/SRF.1994070600' ! SRF File: 2
'output/RAD.1994070100' ! RAD File: 1
'output/RAD.1994070600' ! RAD File: 2
.false. ! Write out header?
***** ATM FILES *****
3 ! iotype: 1=I*2 NetCDF; 2=r*4 NetCDF; 3=grads; 4=vis5d
.false. ! Write out all ATM data btwn idate1 & idate2?
'ATM19940701_19940716.dat' ! ATM filename for rewrite
.false. ! Average ATM data btwn idate1 & idate2?
'ATM19940701_19940716AVG.dat' ! ATM average file to be written
.true. ! Avg Diurnal ATM data btwn idate1 & idate2?
'ATM19940701_19940716DIUR.dat' ! ATM diurnal average file to be written
.false. ! Continual avergage ATM data btwn idate1 & idate2?
'ATM19940701_19940716AVGMON.dat' ! ATM continual avg file to be written
***** SRF FILES *****
3 ! iotype: 1=I*2 NetCDF; 2=r*4 NetCDF; 3=grads; 4=vis5d
.false. ! Write out all SRF data btwn idate1 & idate2?
'SRF19940701_19940716.dat' ! SRF filename for rewrite
.false. ! Average SRF data btwn idate1 & idate2?
'SRF19940701_19940716AVG.dat' ! SRF average file to be written
.true. ! Avg Diurnal SRF data btwn idate1 & idate2?
'SRF19940701_19940716DIUR.dat' ! SRF diurnal average file to be written
.false. ! Continual avergage SRF data btwn idate1 & idate2?
'SRF19940701_19940716AVGMON.dat' ! SRF continual avg file to be written
***** RAD FILES *****
3 ! iotype: 1=I*2 NetCDF; 2=r*4 NetCDF; 3=grads; 4=vis5d
.false. ! Write out all RAD data btwn idate1 & idate2?
'RAD19940701_19940716.dat' ! RAD filename for rewrite
.false. ! Average RAD data btwn idate1 & idate2?
'RAD19940701_19940716AVG.dat' ! RAD average file to be written
.true. ! Avg Diurnal RAD data btwn idate1 & idate2?
'RAD19940701_19940716DIUR.dat' ! RAD diurnal average file to be written
.false. ! Continual avergage RAD data btwn idate1 & idate2?
'RAD19940701_19940716AVGMON.dat' ! RAD continual avg file to be written

and rerun ./postproc, to get the 3 diurnal average files:
ATM19940701_19940716DIUR.dat
SRF19940701_19940716DIUR.dat
RAD19940701_19940716DIUR.dat

The descriptor files for ATM19940701_19940716DIUR.dat is similar
to ATM19940701_19940716AVG.dat, except the tdef is set to 4 instead.

4. SIGMA to P converter

Now quit from your GrADS enviroment, we will use FORTRAN code to
convert varaibales in ATM.1994070100 and ATM.1994070600 from
SIGMA level to P level, and save them to PLEV_VAR.

Your work directory now is /scratch/RegCM/PracticeRun/output/
cp ../../Commons/tools/SIGMAtoP.f SIGMAtoP.f
cp ../../Commons/tools/PLEV_VAR.ctl PLEV_VAR.ctl

edit SIGMAtoP.f, and make sure the parameters setting and ATM files
name are the right one.
jx ! number of points from west to east
ix ! number of points from south to north
kx ! number of vertical levels of SIGMA
np ! number of vertical levels of P
nfile ! how many ATM.yyyymmddhh files you have
inout ! file name of ATM
number ! how many time slices in your ATM.yyyymmddhh file

pgf77 -byteswapio -o SIGMAtoP SIGMAtoP.f
./SIGMAtoP

Hopefully, your job has succesfully finished and you have PLEV_VAR
file now.

Start up GrADS, enter:
open PLEV_VAR.ctl

The pdef is defined in PLEV_VAR.ctl, so you have to use method-1 or
method-2 to make the map projection.

Note: to view wind fields in vector form, you'd better use method-3,
because the wind direction will be absolutely right in method-3.

    set mproj lambert
    set lev 500
    set gxout shaded
    d h
    set gxout contour
    d h
    set gxout vector
    d skip(u,5,5);skip(v,5,5)
    draw title H at 500hPa

Now enter:
    c
    set mproj latlon
    set gxout shaded
    d h
    set gxout contour
    d h
    draw title H at 500hPa

Note: If you use method-1 (pdef with 'set mproj latlon') and method-2, you have to check carefully the shaded aera of you domain if you want to plot wind fields in vector way, and make sure your domain is in rectangle shape.

In a similar way, you could make the SIGMA to P convertor for your Input fields (of ECMWF, or NCEP) under /RegCM/Input/tools by
    cd ../../Input/tools/
    pgf77 -byteswapio -o SIGMAtoP SIGMAtoP.f
    ./SIGMAtoP

You could compare you input and output variables on P-levels.