Automated front drawing
Evgeny Himmelreich
evgenyh at IMS.GOV.IL
Tue Feb 8 06:13:39 EST 2005
Michael !
Thanks a lot for procedures
==========================
Evgeny Himmelreich
Israel Meteorological Service
P.O.Box 25, Bet-Dagan, ISRAEL
Tel: +972-3-9682157
http://www.ims.gov.il
==========================
----- Original Message -----
From: "Sestak, Dr. Michael" <michael.sestak at FNMOC.NAVY.MIL>
To: <GRADSUSR at LIST.CINECA.IT>
Sent: Monday, February 07, 2005 6:01 PM
Subject: Automated front drawing
> 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