Automated front drawing
Sestak, Dr. Michael
michael.sestak at FNMOC.NAVY.MIL
Mon Feb 7 11:01:29 EST 2005
I thought it would only take a little while to clean up the code for public
consumption, but it is more thoroughly enmeshed in our operational system
than I thought, so this isn't as nice as I would like it to be (at least in
this version, only the needed grids and related factors are passed as
parameters instead of being spread across a bunch of common blocks). Also,
as I said before, like Vladimir Corvus script, this only produces frontal
areas. I also tinkered around with another subroutine to extract the
locations of the proper edges which best represent the actual front
location, but that code seems to have been lost a few months back when one
of our machines had a hard crash. I don't think it was a very complicated
program, but I don't have the time to recreate it here. I think I did
something really simple like assuming most of the grid has ggtheta values
that are about the same, so the average of the grid would be slightly above
this baseline value. So, I collected all the points that would be on this
average contour, then did some sorting to make lists of these points so
those of one frontal zone would be together. The final step would be to use
draw line to have grads produce the front lines. That's pretty simple.
Making them the right color for cold and warm fronts and drawing points and
bubbles would be much more difficult. The articles cited previously (such
as Hewson, T, 1998, Objective fronts, Meteorological Applications,
5(3):37-67, D, W. McCann and J. P Whistler, 2001, Problems and solutions for
drawing fronts objectively, Meteorological Applications 8:193-203) give some
ideas about the graphics production end as well as alternatives for the
basic frontal zone algorithm. The rest is left as an example to the student
(because I clearly do not have time to do any more on this).
The inputs are
pgeom String designating the model geometry. If the pgeom starts
with global, calculations are modified for global wrap-around
The only geometry this algorithm handles are
latitude/longitude grids
im Number of grid points in x (longitude) direction
jm Number of grid points in y (latitude) direction
dlat Grid cell size in latitude (Y) direction
dlon Grid cell size in longitude (X) direction
h1000 1000 mb height
h700 700 mb height
ggt Gridded frontal zone parameter
subroutine ggtheta(pgeom,im,jm,dlat,dlon,h1000,h700,ggt)
C
C...................START PROLOGUE.....................................
C
C SCCS IDENTIFICATION: %W% %G%
C
C CONFIGURATION IDENTIFICATION:
C DAF
C
C MODULE NAME: ggtheta
C
C DESCRIPTION: Perform objective frontal analysis
C
C COPYRIGHT: (c) 1993 FLENUMMETOCCEN
C U.S. GOVERNMENT DOMAIN
C ALL RIGHTS RESERVED
C
C CONTRACT NUMBER AND TITLE:
C GS-09K-90BHD0001
C ADP Technical Services for GSA Pacific Zone
C
C REFERENCES: None
C
C CLASSIFICATION: Unclassified
C
C RESTRICTIONS: None
C
C COMPUTER/OPERATING SYSTEM DEPENDENCIES:
C None
C
C LIBRARIES OF RESIDENCE:
C N/A
C
C USAGE: call getz10
C
C PARAMETERS:
C
C Name Type Description
C ------------ ------- ------------------------------------------------
C imjm integer Size of grid (rows x columns)
C h1000 real Grid of 1000 mb heights
C h700 real Grid of 700 mb heights
C
C
C
C FILES: None
C
C DATABASES: None
C
C NON-FILE INPUT/OUTPUT:
C None
C
C ERROR CONDITIONS:
C None
C
C ADDITIONAL COMMENTS:
C None
C
C....................MAINTENANCE SECTION................................
C
C MODULES CALLED:
C Name Description
C ------------ --------------------------------------------------------
C grdpfg compute the gradient of a field
C
C LOCAL VARIABLES AND STRUCTURES:
C Name Type Description
C ------------ ------- ------------------------------------------------
C g real work array; equivalenced to tmp2 (common //)
C gx real array containing x component of grad(grad(T));
C equivalenced to tmp1 (common //)
C gy real array containing y component of grad(grad(T));
C equivalenced to tmp3 (common //)
C tbar real array containing average temperature of 1000-
C 700 mb layer; equivalenced to tmp1 (common //)
C tx real array containing x component of grad(T);
C equivalenced to tmp4 (common //)
C ty real array containing y component of grad(T);
C equivalenced to tmp5 (common //)
C Donald W. McCann and James P Whistler real array
containing frontal analysis parameter
C (K/m**2); equivalenced to fldout
C (common /fldpfg/)
C ij integer loop index
C fac real factor to convert thickness to average
C temperature for a layer
C
C METHOD:
C grad(grad(T) * grad(T)) * grad(T)
C ggt = - .5* ---------------------------------
C grad(T) * grad(T)
C
C
C COMPILER DEPENDENCIES:
C FORTRAN 90
C
C COMPILER OPTIONS: None
C
C MAKEFILE: /a/ops/met/daf/src/sub/make_daflib
C
C RECORD OF CHANGES:
C
C <<CHANGE NOTICE>> version 1.1 (30 Dec 1993) -- Rennick, M.A. (CSC)
C initial installation
C <<CHANGE NOTICE>> version 1.2 (12 Oct 2001) -- Sestak, M.L.
C modification to use coamps grids
C
C....................END PROLOGUE......................................
C
implicit none
c Parameters
character (len=6) :: pgeom
integer :: im,jm,imjm
real :: dlat,dlon
real, dimension(imjm) :: h1000, h700, ggt
c Local variables
real, allocatable :: g(:), gx(:), gy(:), tbar(:), tx(:), ty(:)
integer :: ij, istat
real :: fac
real :: afac, asold, e0, grav, rad, rgas, zrok, cp
data afac, asold, grav, e0, rad, rgas, zrok
2 / 5412., 6135., 9.81, 6.11, 6.371229e6, 287., 273.15/
c
imjm = im*jm
allocate (g(imjm),gx(imjm),gy(imjm),tbar(imjm),stat=istat)
if (istat.ne.0) then
write(*,*) 'ggtheta: Failed to allocate arrays'
return
end if
allocate (tx(imjm),ty(imjm),stat=istat)
if (istat.ne.0) then
write(*,*) 'ggtheta: Failed to allocate arrays'
return
end if
c compute mean T from 1000-700mb thickness
fac = grav/(rgas*alog(1000.0/700.0))
c
do 110 ij=1,imjm
tbar(ij) = fac*(h700(ij) - h1000(ij))
110 continue
c compute grad(tbar)
call grdpfg(im,jm,dlat,dlon,tbar,g,tx,ty,0)
c compute grad(tbar)/|grad(tbar)|
do 120 ij=1,imjm
if (g(ij) .eq. 0.) then
tx(ij) = 0.
ty(ij) = 0.
else
tx(ij) = tx(ij)/g(ij)
ty(ij) = ty(ij)/g(ij)
end if
120 continue
c compute grad(g)
call grdpfg(im,jm,dlat,dlon,g,g,gx,gy,0)
c compute ggtheta
do 130 ij=1,imjm
ggt(ij) = - (gx(ij)*tx(ij) + gy(ij)*ty(ij))
130 continue
deallocate (g,gx,gy,tbar,stat=istat)
deallocate (tx,ty,stat=istat)
return
end
subroutine grdpfg(im,jm,dlat,dlon,fld,grd,gx,gy,nflag)
C
C...................START PROLOGUE.....................................
C
C SCCS IDENTIFICATION: %W% %G%
C
C CONFIGURATION IDENTIFICATION:
C DAF
C
C MODULE NAME: grdpfg
C
C DESCRIPTION: Compute gradient of a field in spherical coordinates
C
C COPYRIGHT: (c) 1993 FLENUMMETOCCEN
C U.S. GOVERNMENT DOMAIN
C ALL RIGHTS RESERVED
C
C CONTRACT NUMBER AND TITLE:
C GS-09K-90BHD0001
C ADP Technical Services for GSA Pacific Zone
C
C REFERENCES: None
C
C CLASSIFICATION: Unclassified
C
C RESTRICTIONS: None
C
C COMPUTER/OPERATING SYSTEM DEPENDENCIES:
C None
C
C LIBRARIES OF RESIDENCE:
C N/A
C
C USAGE: call grdpfg(fld, grd, gx, gy, nflag)
C
C PARAMETERS:
C Name Type Usage Description
C ------------ ------- ------- ----------------------------------------
C fld real input array containing field from which to
C compute the gradient
C grd real output array containing |grad(fld)|
C gx real output array containing x component of
C grad(fld)
C gy real output array containing y component of
C grad(fld)
C nflag integer input flag to determine if fld should be
C multiplied by cos(phi) before
C differentiating wrt y
C
C COMMON BLOCKS:
C Common blocks are documented where they are defined
C in the code within include files. This module uses
C the following common blocks
C Name Description
C ------------ --------------------------------------------------------
C phypfg physical constants common to several DAF modules
C
C FILES: None
C
C DATABASES: None
C
C NON-FILE INPUT/OUTPUT:
C None
C
C ERROR CONDITIONS:
C None
C
C ADDITIONAL COMMENTS:
C None
C
C....................MAINTENANCE SECTION................................
C
C MODULES CALLED:
C None
C
C LOCAL VARIABLES AND STRUCTURES:
C Name Type Description
C ------------ ------- ------------------------------------------------
C i integer loop index
C im1 integer array index
C ip1 integer array index
C j integer loop index
C cosphi real cos(latitude)
C divfac real array containing value by which to multiply fld
C before computing derivatives
C dlam real interval between grid points in x-direction
C (radians)
C dx real array containing interval between grid points
C in x-direction (meters)
C dphi real interval between grid points in y-direction
C (radians)
C phi real latitude (radians)
C pi real pi
C rdx2 real reciprocal of twice the x-interval (1/meters)
C rdy2 real reciprocal of twice the y-interval (1/meters)
C
C METHOD:
C fld must defined on a lon,lat grid with
C dlam = 360/im
C dphi = 180/(jm-1)
C
C gx = x-component of grad(fld)
C gy = y-component of grad(fld)
C grd = magnitude of grad(fld)
c
C
C INCLUDE FILES:
C Name Description
C ------------ --------------------------------------------------------
C parpfg.h parameter statements common to several daf modules
C phypfg.h common block /phypfg/
C
C COMPILER DEPENDENCIES:
C FORTRAN 90
C
C COMPILER OPTIONS: None
C
C MAKEFILE: /a/ops/met/daf/src/sub/make_daflib
C
C RECORD OF CHANGES:
C
C <<CHANGE NOTICE>> version 1.1 (30 Dec 1993) -- Rennick, M.A. (CSC)
C initial installation
C <<CHANGE NOTICE>> version 1.2 (12 Oct 2001) -- Sestak, M.L.
C modification to use coamps grids
C
C....................END PROLOGUE......................................
C
implicit none
integer :: im,jm,nflag
real :: dlat,dlon
real :: fld(im,jm), grd(im,jm), gx(im,jm), gy(im,jm)
integer :: i, im1, ip1, j
real :: dy, rdy, mlr, rdx
real :: cosphi, divfac(jm), dlam, dx(jm)
real :: dphi, phi, pi, rdx2, rdy2
c
pi = 4.*atan(1.)
c longitude and latitude increments in radians
dphi = dlat*pi/180.
dlam = dlon*pi/180.
c meridianal grid increment in meters
dy = dlam*rad
rdy = 1./dy
rdy2 = 1./(2.*dy)
c zonal grid increments in meters
mlr = minlat * pi / 180.
do j=1,jm
phi = mlr + (j-1)*dphi
cosphi = cos(phi)
dx(j) = rad*cosphi*dlam
divfac(j) = 1. - nflag + nflag*cosphi
end do
c x and y components, interior points, global and regional
do j=2,jm-1
rdx2 = 1./(2.*dx(j))
do i=2,im-1
gx(i,j) = rdx2*(fld(i+1,j) - fld(i-1,j))
gy(i,j) = rdy2*(fld(i,j+1)*divfac(j+1) - fld(i,j-1)
2 *divfac(j-1))
end do
end do
c x and y components, poles global
if (pgeom(1:6)=='global') then
call zilch(gx(1,1),im)
call zilch(gx(1,jm),im)
do i=1,im
gy(i,1) = rdy*(fld(i,2)*divfac(2) - fld(i,1)*divfac(1))
gy(i,jm)=rdy*(fld(i,jm)*divfac(jm)-fld(i,jm-1)*divfac(jm-1))
end do
c x and y components, north and south edges, regional
else
do i=2,im-1
rdx2 = 1./(2.*dx(1))
gx(i,1) = rdx2*(fld(i+1,1) - fld(i-1,1))
gy(i,1) = rdy*(fld(i,2)*divfac(2) - fld(i,1)*divfac(1))
rdx2 = 1./(2.*dx(jm))
gx(i,jm) = rdx2*(fld(i+1,jm) - fld(i-1,jm))
gy(i,jm)=rdy*(fld(i,jm)*divfac(jm)-fld(i,jm-1)*divfac(jm-1))
end do
rdx = 1./dx(1)
gx(1,1) = rdx*(fld(2,1) - fld(1,1))
gx(im,1) = rdx*(fld(im,1) - fld(im-1,1))
gy(1,1) = rdy*(fld(1,2)*divfac(2) - fld(1,1)*divfac(1))
gy(im,1) = rdy*(fld(im,2)*divfac(2) - fld(im,1)*divfac(1))
rdx = 1./dx(jm)
gx(1,jm) = rdx*(fld(2,jm) - fld(1,jm))
gx(im,jm) = rdx*(fld(im,jm) - fld(im-1,jm))
gy(1,jm)=rdy*(fld(1,jm)*divfac(jm)-fld(1,jm-1)*divfac(jm-1))
gy(im,jm)=rdy*(fld(im,jm)*divfac(jm)-fld(im,jm-1)*divfac(jm-1))
end if
c x znd y components, east/west wraparound, global
if (pgeom(1:6)=='global') then
do j=2,jm-1
rdx2 = 1./(2.*dx(j))
gx(1,j) = rdx2*(fld(2,j) - fld(im,j))
gy(1,j) = rdy2*(fld(1,j+1)*divfac(j+1) - fld(1,j-1)
2 *divfac(j-1))
gx(im,j) = rdx2*(fld(1,j) - fld(im-1,j))
gy(im,j) = rdy2*(fld(im,j+1)*divfac(j+1) - fld(im,j-1)
2 *divfac(j-1))
end do
c x znd y components, east/west edges, regional
else
do j=2,jm-1
rdx = 1./dx(j)
gx(1,j) = rdx*(fld(2,j) - fld(1,j))
gy(1,j) = rdy2*(fld(1,j+1)*divfac(j+1) - fld(1,j-1)
2 *divfac(j-1))
gx(im,j) = rdx*(fld(im,j) - fld(im-1,j))
gy(im,j) = rdy2*(fld(im,j+1)*divfac(j+1) - fld(im,j-1)
2 *divfac(j-1))
end do
end if
do j=1,jm
do i=1,im
grd(i,j) = sqrt(gx(i,j)**2 + gy(i,j)**2)
end do
end do
return
end
-----Original Message-----
From: Evgeny Himmelreich [mailto:evgenyh at IMS.GOV.IL]
Sent: Sunday, February 06, 2005 9:23 PM
To: GRADSUSR at LIST.CINECA.IT
Subject: Re: Script to plot fronts from NCEP data
Hi Prodromos !
Check in online archive for 2005
for "Automatic front drawing alorithm"
1. Dr. Micheal Sestak promiced at the beginning of January,
to post fortran routine, for fronts location determination.
2. Vladimir Corvus posted GRADS script for frontal zones determination and
drawing
Best Regards
Evgeny Himelreich
==========================
Evgeny Himmelreich
Israel Meteorological Service
P.O.Box 25, Bet-Dagan, ISRAEL
Tel: +972-3-9682157
http://www.ims.gov.il
==========================
----- Original Message -----
From: "Prodromos Zanis" <pzanis at GEOL.UOA.GR>
To: <GRADSUSR at LIST.CINECA.IT>
Sent: Thursday, February 03, 2005 2:38 PM
Subject: Script to plot fronts from NCEP data
> Dear all
>
> I would like to ask if there is any grads script to plot fronts overlying
on
> surface pressure maps based on NCEP reannalysis data.
>
> Sorry for the trouble
>
> best regards
> Prodromos Zanis
>
>
>
>
> --
> ****************************************************
> Dr. Prodromos Zanis
> Centre for Atmospheric Physics and Climatology
> Academy of Athens
> 3rd September, Athens 15784, Greece
> Tel. +30 210 8832048
> Fax: +30 210 8832048
> e-mail: pzanis at geol.uoa.gr
> Web address: http://users.auth.gr/~zanis/
> *****************************************************
More information about the gradsusr
mailing list