      SUBROUTINE growth(per,ter,rer,her)

c**************************************************
c Purpose:   Increase the analysis errors by some amount defined
c            by an empirical formula which represents an error in
c            the forecast.
c            Result is the forecast error which is used in the oi
c            covariance calculations.
c
c Programmer:J. Pfaendtner
c
c Modification(s): 1990/04/10 Replacement of NMC obs errors with NASA
c                  obs errors; corrected computation of alpha and beta.
c                  ...M. Seablom
c
c                  1990/03/08 Allows for increasing and generalizing the
c                  number of levels.
c                  ...D. Lamich
c***********************************************
C**********************************************************************
C  Data parallel implementation (CM Fortran) : M. Makivic (NPAC)      *
C  Date:       4/94                                                   *
C**********************************************************************


      PARAMETER (max_levels = 18)
      DIMENSION her (im, jnp, mlev), rer (im , jnp, mlev), per (im, jnp)
      DIMENSION ter (im, jnp)
      COMMON / gridxx / im, jnp, mlev
      DIMENSION rlim (50), rgro (50)
      DIMENSION hlimt (50), hgrot (50)
      DIMENSION hlimx (50), hgrox (50)
      REAL rgro_table (max_levels)
      REAL hgrot_table (max_levels)
      REAL hlimt_table (max_levels)
      REAL hgrox_table (max_levels)
      REAL hlimx_table (max_levels)
      REAL pi, phi
      LOGICAL first
!hpf$ distribute her(block,block,block)
!hpf$ distribute hgrot(block)
!hpf$ distribute hgrot_table(block)
!hpf$ distribute hgrox(block)
!hpf$ distribute hgrox_table(block)
!hpf$ distribute hlimt(block)
!hpf$ distribute hlimt_table(block)
!hpf$ distribute hlimx(block)
!hpf$ distribute hlimx_table(block)
!hpf$ distribute per(block,block)
!hpf$ distribute rer(block,block,block)
!hpf$ distribute rgro(block)
!hpf$ distribute rgro_table(block)
!hpf$ distribute rlim(block)
!hpf$ distribute ter(block,block)

      real  const(im,jnp),deglat(im,jnp)
      real pgro(im,jnp),plim(im,jnp),hgro(im,jnp),hlim(im,jnp)
!hpf$ distribute const(block,block)
!hpf$ distribute deglat(block,block)
!hpf$ distribute pgro(block,block)
!hpf$ distribute plim(block,block)
!hpf$ distribute hgro(block,block)
!hpf$ distribute hlim(block,block)

      DATA rgro_table / 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 
     &.05, .15, .20, .65, .9, 1.2 /
      DATA hgrot_table / 125.0, 125.0, 125.0, 125.0, 125.0, 110.0, 95.0,
     & 80.0, 60.0, 25.0, 21.0, 17.0, 14.0, 13.0, 12.0, 10.5, 9.0, 9.0 /
      DATA hlimt_table / 147.0, 147.0, 147.0, 147.0, 147.0, 132.0, 117.0
     &, 105.0, 96.0, 76.0, 56.0, 37.0, 29.0, 26.9, 24.5, 20.0, 20.0, 20.
     &0 /
      DATA hgrox_table / 70.0, 70.0, 70.0, 70.0, 70.0, 59.0, 48.0, 37.0,
     & 35.0, 32.0, 30.0, 28.0, 26.5, 23.0, 20.0, 15.0, 12.5, 10.5 /
      DATA hlimx_table / 112.5, 112.5, 112.5, 112.5, 112.5, 109.5, 106.5
     &, 103.0, 98.5, 95.0, 91.5, 88.5, 85.0, 78.0, 71.0, 58.0, 56.5, 55.
     &0 /
      DATA first / .true. /
         level_start = max_levels - mlev + 1

         do l=1,mlev
            rgro(l) = rgro_table(level_start+l-1)
            rlim(l) = rgro(l) * 8.0
            hgrot(l) = hgrot_table(level_start+l-1)
            hlimt(l) = hlimt_table(level_start+l-1)
            hgrox(l) = hgrox_table(level_start+l-1)
            hlimx(l) = hlimx_table(level_start+l-1)
         end do

         pgrot = hgrot(mlev)/8.0
         pgrox = hgrox(mlev)/8.0
         plimt = hlimt(mlev)/8.0
         plimx = hlimx(mlev)/8.0

         pi = 4.0*atan(1.0)
         dphi = 180.0 / ( jnp - 1 )
         philo = 15.
         phihi = 45.
         imjnp = im * jnp
         first = .false.

      print 900

c  Augment analysis errors with empirical growth formula
c      forall(j=1:jnp) deglat(:,j)=-90.0+(j-1)*dphi
      do j=1,jnp 
          deglat(:,j)=-90.0+(j-1)*dphi
      end do
      where(abs(deglat).GT. philo.and.abs(deglat) .LT. phihi)
         plim=0.5*(plimt+plimx)+
     &   0.5*(plimt-plimx)*cos((abs(deglat)-philo)/(phihi-philo)*pi)
         pgro=0.5*(pgrot+pgrox)+
     &   0.5*(pgrot-pgrox)*cos((abs(deglat)-philo)/(phihi-philo)*pi)
      endwhere
      where(abs(deglat).LE. philo)
          pgro=pgrot
          plim=plimt
      end where
      where(abs(deglat) .GE. phihi)
          pgro = pgrox
          plim = plimx
      end where
      const = (plim - pgro) / plim
      per=const*per+pgro

      do l = 1, mlev
        if( rlim(l).ne.0.0 ) then
            const = (rlim(l)-rgro(l)) / rlim(l)
            rer(:,:,l) = const * rer(:,:,l) + rgro(l)
        endif
        where(abs(deglat).GT. philo.and.abs(deglat) .LT. phihi)
          hlim=0.5*(hlimt(l)+hlimx(l))+
     &    0.5*(hlimt(l)-hlimx(l))*
     &    cos((abs(deglat)-philo)/(phihi-philo)*pi)
          hgro=0.5*(hgrot(l)+hgrox(l))+
     &    0.5*(hgrot(l)-hgrox(l))*
     &    cos((abs(deglat)-philo)/(phihi-philo)*pi)
        endwhere
        where(abs(deglat).LE. philo)
          hgro=hgrot(l)
          hlim=hlimt(l)
        end where
        where(abs(deglat) .GE. phihi)
          hgro = hgrox(l)
          hlim = hlimx(l)
        end where
        const = (hlim - hgro) / hlim
        her(:,:,l)=const*her(:,:,l)+hgro

      enddo

      RETURN
900   FORMAT (' GROWTH: Adding 6hr growth to analysis errors')
      END
