This code generates two-dimensional analytic perturbation fields. These perturbations can be used, for example, to initialize idealized numerical model simulations. The methodology is described in Bryan et al. (2007, JAS, p. 1249, section 4a); see also Knievel et al. (2007, MWR, 3808, section 2f). https://www2.mmm.ucar.edu/people/bryan/Code/getp…

getperts.F

    PROGRAM getperts
    use singleton
    implicit none

!-----------------------------------------------------------------------
!
!  getperts:  Version 1.04                  Last modified:  12 October 2008
!
!  Author:  George H. Bryan
!           Mesoscale and Microscale Meteorology Division
!           National Center for Atmospheric Research
!           Boulder, Colorado, USA
!           gbryan@ucar.edu
!
!  Description:
!  This is a program to generate idealized two-dimensional perturbations
!  based on a specified spectrum.  The main output file is called perts.dat:
!  it is a GrADS binary file.  The program also creates two more files ...
!  ... spec1d.dat and spec2d.dat ...... that contain information used 
!  to create the perturbations;  they are for diagnostic purposes only. 
!
!  To compile, run the csh script "compile".  You may have to change 
!  the compiler or compile flags, depending on your computer architecture.
!  The code was created and tested on a Linux machine using the Portland
!  Group compiler (pgf90).
!
!  References:
!
!  Bryan, Rotunno, Fritsch (2007, JAS, p. 1249) (see section 4a)
!  Knievel, Bryan, Hacker (2007, MWR, p. 3808) (see section 2f)
!
!  Disclaimer:  This code is made available WITHOUT WARRANTY.
!
!-----------------------------------------------------------------------
!  Begin:  parameters to be set by user

    integer, parameter :: ni = 192    ! number of grid points in x
    integer, parameter :: nj = 192    ! number of grid points in y

    real, parameter :: dx =  125.0    ! grid spacing (m)

    real, parameter :: maxval = 1.00  ! maximum absolute value of perturbations

    integer, parameter :: spectrum_type = 1
                                      ! 1 = white noise (kappa^0)
                                      ! 2 = kappa^(-5/3)
                                      ! 3 = kappa^(-3)

    real, parameter :: large_eddy_scale = 1000000000.0 
                                      ! Large eddy size (m).
                                      ! Above this scale, the spectrum is
                                      ! white.  If you don't want to use
                                      ! this, make it larger than 
                                      ! max(ni,nj)*dx.

    integer, parameter :: filter_type = 2
                                      ! 0 = no filter
                                      ! 1 = remove all scales smaller than 
                                      !     6 Delta
                                      ! 2 = disk:  keep all scales between
                                      !     L1 and L2
                                      ! 3 = ring:  keep only scale L1

    ! For filter_type 2 and 3 only:
    real, parameter :: L1 =  1000.0   ! smaller scale (m)
    real, parameter :: L2 =  8000.0   ! larger scale (m)

!  End:  parameters to be set by user.
!-----------------------------------------------------------------------
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----------------------------------------------------------------------
!            No need to modify anything below here:
!-----------------------------------------------------------------------
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----------------------------------------------------------------------
!  Declared variables

    integer :: i,j,ik,k,nn,irec,nrec,orec,i1,i2,nsum

    real, dimension(ni,nj) :: kappa
    complex, dimension(ni,nj) :: c

    integer jj,j1,j2,imax,imin,jmin
    real a,rand,avg,beta,rmax,rmin,r1,r2,r3
    real val,vmin,threshold
    character*50 fname
    real*8 sum

    integer, parameter :: nk = min(ni,nj)

    integer, dimension(nk) :: np
    real*8, dimension(nk) :: sp

!-----------------------------------------------------------------------

    real, parameter :: pi = 3.14159265358979323
    real, parameter :: i_c   = ni/2+1    ! center point in i
    real, parameter :: j_c   = nj/2+1    ! center point in j

!-----------------------------------------------------------------------
!  open output files

    open(unit=60,file='spec2d.dat',form='unformatted',   &
         access='direct',recl=4,status='unknown')

    nrec=1

    open(unit=57,file='spec1d.dat',form='unformatted',   &
         access='direct',recl=4,status='unknown')
    orec=1


!-----------------------------------------------------------------------
! Get kappa

    print *
    print *,'  dx,Lx,Ly = ',dx,ni*dx,nj*dx
    print *

    do j=1,nj
    do i=1,ni
      kappa(i,j)=sqrt( ( 2.0*pi*(float(i)-i_c)/(ni*dx))**2   &
                      +( 2.0*pi*(float(j)-j_c)/(nj*dx))**2 )
    enddo
    enddo

    IF( large_eddy_scale .lt. (max(ni,nj)*dx) )THEN

      do j=1,nj
      do i=1,ni
        kappa(i,j)=max(kappa(i,j),2.0*pi/large_eddy_scale)
      enddo
      enddo

    ENDIF

    do j=1,nj
    do i=1,ni
      write(60,rec=nrec) kappa(i,j)
      nrec=nrec+1
    enddo
    enddo

!-----------------------------------------------------------------------
!  Specify the spectrum.  (This is the primary component of the code.)

    print *,'  specifying spectrum'

    sum = 0.0d0

    rmax = -99999999.9
    rmin =  99999999.9

    call random_seed

    do j=1,nj
    do i=1,ni

      call random_number(r1)
      call random_number(r2)

      IF(spectrum_type.eq.1)THEN    ! white (kappa^0)

        c(i,j)=exp( -csqrt(cmplx(-1.0))*2.0*pi*(float(ni)*r1) )  &
              *exp( -csqrt(cmplx(-1.0))*2.0*pi*(float(nj)*r2) )

      ELSEIF(spectrum_type.eq.2)THEN    ! kappa^(-5/3)

        c(i,j)=exp( -csqrt(cmplx(-1.0))*2.0*pi*(float(ni)*r1) )  &
              *exp( -csqrt(cmplx(-1.0))*2.0*pi*(float(nj)*r2) )  &
              *( kappa(i,j)**(-4.0/3.0) )

      ELSEIF(spectrum_type.eq.3)THEN    ! kappa^(-3)

        c(i,j)=exp( -csqrt(cmplx(-1.0))*2.0*pi*(float(ni)*r1) )  &
              *exp( -csqrt(cmplx(-1.0))*2.0*pi*(float(nj)*r2) )  &
              *( kappa(i,j)**(-2.0) )

      ELSE

        print *,'  unknown spectrum_type'
        stop

      ENDIF

    enddo
    enddo

    c(int(i_c),int(j_c)) = 0.0   ! This ensures that the mean is zero.

    nsum=0

!-----------------------------------------------------------------------
!  redefine kappa  ...  ignore large_eddy_scale this time

    do j=1,nj
    do i=1,ni
      kappa(i,j)=sqrt( ( 2.0*pi*(float(i)-i_c)/(ni*dx))**2   &
                      +( 2.0*pi*(float(j)-j_c)/(nj*dx))**2 )
    enddo
    enddo

!-----------------------------------------------------------------------
!  Make sure the array makes sense:  lower-left is conjugate of upper-right.

    do j=1,nj
    do i=1,ni
      if(ni-i+2.ge.1.and.ni-i+2.le.ni.and.   &
         nj-j+2.ge.1.and.nj-j+2.le.nj)then
        c(i,j)=conjg(c(ni-i+2,nj-j+2))
      endif
    enddo
    enddo

!-----------------------------------------------------------------------
!  Optional:  filter the spectrum.
!---------------------------------------------------------
! remove all information at scales smaller than 6-delta.

  IF(filter_type.eq.1)THEN

    print *,'  i_c,j_c  = ',i_c,j_c
    print *,'  i1,j1    = ',(ni/6.0)+0.5,(nj/6.0)+0.5

    do j=1,nj
    do i=1,ni
      IF(  ( (float(i)-i_c)/(float(ni)/6.0+0.5) )**2      &
         + ( (float(j)-j_c)/(float(nj)/6.0+0.5) )**2      &
               .gt. 1.0 )THEN
        c(i,j) = 0.0
      ENDIF
    enddo
    enddo

!---------------------------------------------------------
! Disk:  keep all scales between L1 and L2:

  ELSEIF(filter_type.eq.2)THEN

    do j=1,nj
    do i=1,ni
      IF( kappa(i,j).gt.(2.0*pi/L1) .or.    &
          kappa(i,j).lt.(2.0*pi/L2) )THEN
        c(i,j) = 0.0
      ENDIF
    enddo
    enddo

!---------------------------------------------------------
! Ring:  keep only scale L1:

  ELSEIF(filter_type.eq.3)THEN

    !--------------------------------
    ! find closest kappa value

    vmin = 99999999.0

    do j=1,nj
    do i=1,ni
      val = abs(kappa(i,j)-(2.0*pi/L1))
      if(val.lt.vmin)then
        vmin = val
        imin = i
        jmin = j
      endif
    enddo
    enddo

    threshold = max( abs(kappa(imin  ,jmin+1)-(2.0*pi/L1)) ,   &
                     abs(kappa(imin  ,jmin-1)-(2.0*pi/L1)) ,   &
                     abs(kappa(imin+1,jmin  )-(2.0*pi/L1)) ,   &
                     abs(kappa(imin-1,jmin  )-(2.0*pi/L1))  )

    do j=1,nj
    do i=1,ni
      IF( abs(kappa(i,j)-(2.0*pi/L1)).gt.threshold )THEN
        c(i,j) = 0.0
      ENDIF
    enddo
    enddo

  ENDIF

!-----
! end optional filtering

!-----------------------------------------------------------------------
!  Shift and perform inverse FFT

    c=cshift(c,ni/2,1)
    c=cshift(c,nj/2,2)

    sum=0.0d0

    do j=1,nj
    do i=1,ni
      sum = sum+real(c(i,j)*conjg(c(i,j)))
    enddo
    enddo

    print *
    print *,'  RHS = ',sum
    print *

    print *,'  calculating fft'

    do j=1,nj
    do i=1,ni
      c(i,j)=c(i,j)*sqrt(float(ni*nj))
    enddo
    enddo

    c=fft(c,inv=.true.)

!-----------------------------------------------------------------------
!  Adjust for user-specified maximum value

    open(unit=71,file='perts.dat',status='unknown',   &
         form='unformatted',access='direct',recl=4)

!  get max/min/avg

    sum=0.0d0
    avg=0.0
    rmax=-99999999.
    rmin=99999999.

    do j=1,nj
    do i=1,ni
      a=real(c(i,j))
      avg=avg+a
      rmax=max(rmax,a)
      rmin=min(rmin,a)
      sum=sum+a**2
    enddo
    enddo

    sum=sum/float(ni*nj)
    avg=avg/(ni*nj)

    print *
    print *,'  LHS = ',sum
    print *

    print *,'  avg,max,min = ',avg,rmax,rmin

!  get adjustment factors

    beta=maxval/max(rmax,abs(rmin))
    rand=0.0
    print *
    print *,'  adjusting:'
    print *,'  beta,rand = ',1.0/beta,rand

    avg=0.0
    rmax=-99999999.
    rmin=99999999.

!  get new max/min/avg

    sum=0.0d0

    irec=1
    do j=1,nj
    do i=1,ni
      a=real(c(i,j))
      a=a*beta + rand
      avg=avg+a
      rmax=max(rmax,a)
      rmin=min(rmin,a)
      write(71,rec=irec) a
      irec=irec+1
      write(60,rec=nrec) a
      nrec=nrec+1
      c(i,j)=cmplx(a)
      sum=sum+a**2
    enddo
    enddo

    avg=avg/(ni*nj)
    sum=sum/float(ni*nj)

    print *
    print *,'  New LHS = ',sum
    print *

    print *,'  avg,max,min = ',avg,rmax,rmin
    print *

    close(unit=71)

!-----------------------------------------------------------------------
!  Go back to spectral space:  peform FFT, shift array.

    print *,'  reverse fft:'

    c=fft(c)

    do j=1,nj
    do i=1,ni
      c(i,j)=c(i,j)/sqrt(float(ni*nj))
    enddo
    enddo

    c=cshift(c,ni/2,1)
    c=cshift(c,nj/2,2)

    sum=0.0d0

    do j=1,nj
    do i=1,ni
      sum=sum+c(i,j)*conjg(c(i,j))
    enddo
    enddo

    print *
    print *,'  New RHS = ',sum
    print *

!-----------------------------------------------------------------------
!  Writeout final 2d spectrum

    do j=1,nj
    do i=1,ni
      write(60,rec=nrec) real( c(i,j) )
      nrec=nrec+1
    enddo
    enddo

    do j=1,nj
    do i=1,ni
      write(60,rec=nrec) aimag( c(i,j) )
      nrec=nrec+1
    enddo
    enddo

    do j=1,nj
    do i=1,ni
      write(60,rec=nrec) abs( c(i,j) )
      nrec=nrec+1
    enddo
    enddo

!-----------------------------------------------------------------------
!  Calculate and writeout final 1d spectrum

    do k=1,nk
      sp(k)=0.0
      np(k)=0
    enddo

    do j=1,nj
    do i=1,ni
      if(ni.gt.nj)then
        ik=nint( sqrt( ((i-i_c)*float(nj)/float(ni))**2+(j-j_c)**2) )
      elseif(nj.gt.ni)then
        ik=nint( sqrt( (i-i_c)**2+((j-j_c)*float(ni)/float(nj))**2) )
      else
        ik=nint( sqrt( (i-i_c)**2+(j-j_c)**2) )
      endif
      if(ik.gt.0 .and. ik.le.(nk/2))then
        sp(ik)=sp(ik)+(abs(c(i,j)))**2
        np(ik)=np(ik)+1
      endif
    enddo
    enddo

    sum = 0.0d0
    do k=1,nk/2-1
      sum = sum+sp(k)
      write(57,rec=orec) log10( sngl(sp(k)) )
      orec=orec+1
    enddo
    print *
    print *,'  New RHS after 1d calculation = ',sum
    print *

!-----------------------------------------------------------------------
!  Create the GrADS descriptor files.

    open(unit=61,file='perts.ctl',status='unknown')
    write(61,101) 'perts.dat   '
    write(61,103)
    write(61,104)
    write(61,115) ni,0.5*0.001*dx,0.001*dx
    write(61,116) nj,0.5*0.001*dx,0.001*dx
    write(61,117) 1,1.0,1.0
    write(61,108)
    write(61,109) 1
    write(61,110) 'perts     ',1
    write(61,112)
    close(unit=61)

!-----------------------------------------------------------------------

    open(unit=61,file='spec2d.ctl',status='unknown')
    write(61,101) 'spec2d.dat  '
    write(61,103)
    write(61,104)
    write(61,115) ni,0.5*0.001*dx,0.001*dx
    write(61,116) nj,0.5*0.001*dx,0.001*dx
    write(61,117) 1,1.0,1.0
    write(61,108)
    write(61,109) 5
    write(61,110) 'k         ',1
    write(61,110) 'w7        ',1
    write(61,110) 'w10       ',1
    write(61,110) 'w11       ',1
    write(61,110) 'w12       ',1
    write(61,112)
    close(unit=61)

!-----------------------------------------------------------------------

    open(unit=61,file='spec1d.ctl',status='unknown')
    write(61,101) 'spec1d.dat  '
    write(61,103)
    write(61,104)
    write(61,105) nk/2-1
    do k=1,nk/2-1
      write(61,*) log10( 2.0*pi*k/(nk*dx) )
    enddo
    write(61,116) 1,1.0,1.0
    write(61,117) 1,1.0,1.0
    write(61,107)
    write(61,108)
    write(61,109)  1
    write(61,110) 'w1        ',1
    write(61,112)
    close(unit=61)

!-----------------------------------------------------------------------

101 format('dset ^',a12)
102 format('byteswapped')
103 format('title CM1 output')
104 format('undef -99999999.')
105 format('xdef ',i5,' levels')
115 format('xdef ',i5,' linear ',f8.4,' ',f8.4)
106 format('ydef ',i5,' levels')
116 format('ydef ',i5,' linear ',f8.4,' ',f8.4)
107 format('zdef 1 linear 1 1')
117 format('zdef ',i5,' linear ',f8.4,' ',f8.4)
108 format('tdef 2 linear 01Z03JUL2000    1HR')
109 format('vars ',i4)
110 format(a10,' ',i5,' 99 data')
112 format('endvars')

120 format('   ',f8.6)

!-----------------------------------------------------------------------
!  All done.

    STOP 99999
    END program getperts