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