PROGRAM arpsensbc,44
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Generate one set of perturbation external BOUNDARY files from two
!  sets of ARPS external boundary files and write it out for the use in
!  ENSEMBLE forecast.
!  The idea is based on the SLAF (Scaled Lagged Average Forecast).
!   The procedure is as follows:
!    1, read file a;
!    2, read in file b;
!    3, find the difference bwtween a and b, store it in a. a=a-b
!    4, read in the base file c or use b as the control file (iread=0).
!    5, generate the perturbation  c=c+a/n (b is used as c in code)
!    6, n is input as the variable iorder; it can be ... -2,-1,+1,+2 ...
!
! AUTHOR: Dingchen Hou
! History: Apr. 30, 1998: developed initially for SAMEX98.
!          Sep. 02, 1999: modified to change the input file to namelist
!          format.
!          Feb-Apr, 2002: (F. KONG) Major modifications to:
!                - add random perturbation option (with irandom = 1,
!                  and iseed provided in arpsensbc.input)
!                - read in domain config directly from data
!                This version allows to produce BC for both SLAF and
!                BMG
!-----------------------------------------------------------------------
!
!  DATA ARRAYS READ IN:
!  (and another set with b repalcing a)
!
!   au        x component of velocity (m/s)
!   av        y component of velocity (m/s)
!   aww        vertical component of velocity in Cartesian
!             coordinates (m/s).
!
!   apt     potential temperature (K)
!   ap      pressure (Pascal)
!
!   aqv       water vapor mixing ratio (kg/kg)
!   aqc       Cloud water mixing ratio (kg/kg)
!   aqr       Rainwater mixing ratio (kg/kg)
!   aqi       Cloud ice mixing ratio (kg/kg)
!   aqs       Snow mixing ratio (kg/kg)
!   aqh       Hail mixing ratio (kg/kg)
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'exbc.inc'
  INCLUDE 'grid.inc'

!-----------------------------------------------------------------------
!
! Dimension nx, ny, nz 
!
!-----------------------------------------------------------------------
  INTEGER :: nch,idummy,hdmpfm
  INTEGER :: nx, ny, nz
!
!-----------------------------------------------------------------------

  REAL, ALLOCATABLE :: au   (:,:,:)    !  u-velocity (m/s)
  REAL, ALLOCATABLE :: av   (:,:,:)    !  v-velocity (m/s)
  REAL, ALLOCATABLE :: aww  (:,:,:)    !  w-velocity (m/s)
  REAL, ALLOCATABLE :: apt  (:,:,:)    !  potential  temperature (K)
  REAL, ALLOCATABLE :: ap   (:,:,:)    !  pressure (Pascal)
  REAL, ALLOCATABLE :: aqv  (:,:,:)    !  water vapor specific humidity
  REAL, ALLOCATABLE :: aqc  (:,:,:)    ! Cloud water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: aqr  (:,:,:)    ! Rain water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: aqi  (:,:,:)    ! Cloud ice mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: aqs  (:,:,:)    ! Snow mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: aqh  (:,:,:)    ! Hail mixing ratio (kg/kg)

  REAL, ALLOCATABLE :: bu   (:,:,:)    !  u-velocity (m/s)
  REAL, ALLOCATABLE :: bv   (:,:,:)    !  v-velocity (m/s)
  REAL, ALLOCATABLE :: bww  (:,:,:)    !  w-velocity (m/s)
  REAL, ALLOCATABLE :: bpt  (:,:,:)    !  potential  temperature (K)
  REAL, ALLOCATABLE :: bp   (:,:,:)    !  pressure (Pascal)
  REAL, ALLOCATABLE :: bqv  (:,:,:)    !  water vapor specific humidity
  REAL, ALLOCATABLE :: bqc  (:,:,:)    ! Cloud water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: bqr  (:,:,:)    ! Rain water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: bqi  (:,:,:)    ! Cloud ice mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: bqs  (:,:,:)    ! Snow mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: bqh  (:,:,:)    ! Hail mixing ratio (kg/kg)
!-----------------------------------------------------------------------
!
!  Work Arrays
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE  :: tem1(:,:,:)

!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: filename1(10),filename2(10),                    &
               filename3(10),filnmout(10)
  INTEGER :: lenfil1(10),lenfil2(10),lenfil3(10),lenfilo(10)
  INTEGER :: iorder,irandom,iseed,iread,ierr
  INTEGER :: ireturn
  INTEGER :: i,j,k
  CHARACTER (LEN=15) :: ctime
  INTEGER :: is,isets, istat
  INTEGER :: sd_id,clipxy, clipz
!
!------------------------------------------------------------------------
!
! NAMELIST blocks for input parameters 
!
!------------------------------------------------------------------------

!  NAMELIST /grid_dims/ nx, ny, nz
!  NAMELIST /ptbcpara/ iorder,iread,dx,dy,dz,ctrlat,ctrlon
  NAMELIST /ptbcpara/ iorder,irandom,iseed,iread,ngbrz,exbcdmp,exbcfmt
  NAMELIST /input_data/ isets,filename1,filename2,filename3
  NAMELIST /outpt_data/ filnmout

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  READ (5,ptbcpara,END=999)
  WRITE(6,'(/a)') 'Namelist block ptbcpara sucessfully read in.'
  WRITE(6,ptbcpara)

  READ (5,input_data,END=999)
  WRITE(6,'(/a)') 'Namelist block input_data sucessfully read in.'
!  WRITE(6,input_data)

  READ (5,outpt_data,END=999)
  WRITE(6,'(/a)') 'Namelist block output_data sucessfully read in.'
!  WRITE(6,outpt_data)

  DO is=1,isets
    WRITE(6,'(/a,i5)')' The data set  No.', is
    lenfil1(is) = LEN_trim(filename1(is))
    lenfil2(is) = LEN_trim(filename2(is))
    lenfil3(is) = LEN_trim(filename3(is))
    lenfilo(is) = LEN_trim(filnmout(is))
    WRITE(6,'(/a,1x,a)') ' data file 1 ', filename1(is)(1:lenfil1(is))
    WRITE(6,'(/a,1x,a)') ' data file 2 ', filename2(is)(1:lenfil2(is))
    WRITE(6,'(/a,1x,a)') ' data file 3 ', filename3(is)(1:lenfil3(is))
    WRITE(6,'(/a,1x,a)') ' output file ', filnmout(is)(1:lenfilo(is))
  END DO
  GO TO 1001
  999 CONTINUE
  PRINT *,'ERROR in reading input file, Program Terminated!'
  STOP
  1001 CONTINUE
!
!-----------------------------------------------------------------------
!
!  Read in the dimensions etc.
!
!-----------------------------------------------------------------------

  IF (exbcfmt == 0) THEN

!  lenfil1(1)=len_trim(filename1(1))
  CALL getunit( nch )
  CALL asnctl ('NEWLOCAL', 1, ierr)
  CALL asnfile(trim(filename1(1)), '-F f77 -N ieee', ierr)

  OPEN(nch,FILE=trim(filename1(1)),STATUS='old',                   &
        FORM='unformatted',IOSTAT=istat)
  IF ( istat /= 0 ) THEN
    PRINT *, 'Problem to open BC file, STOP'
    STOP
  END IF
  READ(nch,ERR=99) nx,ny,nz,dx,dy,dz,ctrlat,ctrlon
  READ(nch,ERR=99) ctime
  READ(nch,ERR=99) ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd,                &
                  qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd,               &
                  qhbcrd,idummy,idummy,idummy,idummy,               &
                  idummy,idummy,idummy,idummy,idummy,               &
                  idummy,idummy,idummy,idummy,idummy,               &
                  idummy,idummy,idummy,idummy,idummy,               &
                  idummy,idummy,idummy,idummy,idummy,               &
                  idummy,idummy,idummy,idummy,idummy

  GO TO 900
  99   CONTINUE
  WRITE (6,*) "READEXBC: ERROR reading data ",        &
        "from file ",trim(filename1(1)(1:lenfil1(1))),"STOP"
  STOP
  900   CONTINUE
  CLOSE (nch)
  CALL retunit( nch )

  ELSE            ! HDF4 format

    CALL hdfopen(trim(filename1(1)), 1, sd_id)
    IF (sd_id < 0) THEN
      WRITE (6,*) "READEXBC: ERROR opening ",                           &
                 trim(filename1(1))," for reading."
      ierr = 1
!      GO TO 900
      STOP
    END IF

    CALL hdfrdc(sd_id,15,"ctime",ctime,istat)

    CALL hdfrdi(sd_id,"nx",nx,istat)
    CALL hdfrdi(sd_id,"ny",ny,istat)
    CALL hdfrdi(sd_id,"nz",nz,istat)
    CALL hdfrdr(sd_id,"dx",dx,istat)
    CALL hdfrdr(sd_id,"dy",dy,istat)
    CALL hdfrdr(sd_id,"dz",dz,istat)
    CALL hdfrdr(sd_id,"dzmin",dzmin,istat)
    CALL hdfrdi(sd_id,"strhopt",strhopt,istat)
    CALL hdfrdr(sd_id,"zrefsfc",zrefsfc,istat)
    CALL hdfrdr(sd_id,"dlayer1",dlayer1,istat)
    CALL hdfrdr(sd_id,"dlayer2",dlayer2,istat)
    CALL hdfrdr(sd_id,"zflat",zflat,istat)
    CALL hdfrdr(sd_id,"strhtune",strhtune,istat)
    CALL hdfrdi(sd_id,"mapproj",mapproj,istat)
    CALL hdfrdr(sd_id,"trulat1",trulat1,istat)
    CALL hdfrdr(sd_id,"trulat2",trulat2,istat)
    CALL hdfrdr(sd_id,"trulon",trulon,istat)
    CALL hdfrdr(sd_id,"sclfct",sclfct,istat)
    CALL hdfrdr(sd_id,"ctrlat",ctrlat,istat)
    CALL hdfrdr(sd_id,"ctrlon",ctrlon,istat)

    CALL hdfrdi(sd_id,"ubcflg",ubcrd,istat)
    CALL hdfrdi(sd_id,"vbcflg",vbcrd,istat)
    CALL hdfrdi(sd_id,"wbcflg",wbcrd,istat)
    CALL hdfrdi(sd_id,"ptbcflg",ptbcrd,istat)
    CALL hdfrdi(sd_id,"prbcflg",prbcrd,istat)
    CALL hdfrdi(sd_id,"qvbcflg",qvbcrd,istat)
    CALL hdfrdi(sd_id,"qcbcflg",qcbcrd,istat)
    CALL hdfrdi(sd_id,"qrbcflg",qrbcrd,istat)
    CALL hdfrdi(sd_id,"qibcflg",qibcrd,istat)
    CALL hdfrdi(sd_id,"qsbcflg",qsbcrd,istat)
    CALL hdfrdi(sd_id,"qhbcflg",qhbcrd,istat)

    CALL hdfrdi(sd_id, 'clipxy', clipxy,istat)
    IF (istat == 0 .AND. clipxy < ngbrz) THEN
      WRITE (6,*) "READEXBC: ERROR, clipxy (ngbrz) in exbc file too small"
      ierr = 1
!      GO TO 900
      STOP
    END IF
    CALL hdfrdi(sd_id, 'clipz', clipz,istat)
    IF (istat == 0 .AND. clipz > rayklow) THEN
      WRITE (6,*) "READEXBC: ERROR, clipz (rayklow) in exbc file too large"
      ierr = 1
!      GO TO 900
      STOP
    END IF
!wdt end block
    ! alternate dump format ...

  END IF

  PRINT *, 'nx,ny,nz:',nx,ny,nz
  PRINT *, 'dx,dy,dz:',dx,dy,dz
  PRINT *, 'ctrlat,ctrlon:',ctrlat,ctrlon
  PRINT *, 'ctime:',ctime

!-----------------------------------------------------------------------
!
!  Allocate the variables and initialize the them to zero
!
!-----------------------------------------------------------------------

   ALLOCATE(au(nx,ny,nz),STAT=istat)
  au=0
  ALLOCATE(av(nx,ny,nz),STAT=istat)
  av=0
  ALLOCATE(aww(nx,ny,nz),STAT=istat)
  aww=0
  ALLOCATE(apt(nx,ny,nz),STAT=istat)
  apt=0
  ALLOCATE(ap(nx,ny,nz),STAT=istat)
  ap=0
  ALLOCATE(aqv(nx,ny,nz),STAT=istat)
  aqv=0
  ALLOCATE(aqc(nx,ny,nz),STAT=istat)
  aqc=0
  ALLOCATE(aqr(nx,ny,nz),STAT=istat)
  aqr=0
  ALLOCATE(aqi(nx,ny,nz),STAT=istat)
  aqi=0
  ALLOCATE(aqs(nx,ny,nz),STAT=istat)
  aqs=0
  ALLOCATE(aqh(nx,ny,nz),STAT=istat)
  aqh=0
  
  ALLOCATE(bu(nx,ny,nz),STAT=istat)
  bu=0
  ALLOCATE(bv(nx,ny,nz),STAT=istat)
  bv=0
  ALLOCATE(bww(nx,ny,nz),STAT=istat)
  bww=0
  ALLOCATE(bpt(nx,ny,nz),STAT=istat)
  bpt=0
  ALLOCATE(bp(nx,ny,nz),STAT=istat)
  bp=0
  ALLOCATE(bqv(nx,ny,nz),STAT=istat)
  bqv=0
  ALLOCATE(bqc(nx,ny,nz),STAT=istat)
  bqc=0
  ALLOCATE(bqr(nx,ny,nz),STAT=istat)
  bqr=0
  ALLOCATE(bqi(nx,ny,nz),STAT=istat)
  bqi=0
  ALLOCATE(bqs(nx,ny,nz),STAT=istat)
  bqs=0
  ALLOCATE(bqh(nx,ny,nz),STAT=istat)
  bqh=0

  ALLOCATE(tem1(nx,ny,nz), STAT=istat)
  tem1=0
!
!-----------------------------------------------------------------------
!
!  Read all input data arrays
!
!-----------------------------------------------------------------------
!
  is=0
  101  CONTINUE
  is=is+1

  ierr=0
  PRINT *,filename1(is)(1:lenfil1(is)),lenfil1(is)
!  CALL readexbc(nx,ny,nz,dx,dy,dz,ctrlat,ctrlon,                        &
!             filename1(is)(1:lenfil1(is)),lenfil1(is), ctime,           &
!             ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd,                           &
!             qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd,qhbcrd,                 &
!             au,av,aww,apt,ap,aqv,aqc,aqr,aqi,aqs,aqh, ierr)
  CALL readexbc(nx,ny,nz,                                               &
             filename1(is)(1:lenfil1(is)),lenfil1(is), ctime,           &
             au,av,aww,apt,ap,aqv,aqc,aqr,aqi,aqs,aqh, ierr)

  IF (ierr /= 0) THEN
    PRINT *,'ierr1=', ierr,' file',filename1(is)(1:lenfil1(is)),        &
                          '  not read, STOP!'
    STOP
  END IF
!
!-----------------------------------------------------------------------
!
!  Read in the verification data.
!
!-----------------------------------------------------------------------
!
  ierr=0
!  CALL readexbc(nx,ny,nz,dx,dy,dz,ctrlat,ctrlon,                        &
!             filename2(is)(1:lenfil2(is)),lenfil2(is), ctime,           &
!             ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd,                           &
!             qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd,qhbcrd,                 &
!             bu,bv,bww,bpt,bp,bqv,bqc,bqr,bqi,bqs,bqh, ierr)
  CALL readexbc(nx,ny,nz,                                               &
             filename2(is)(1:lenfil2(is)),lenfil2(is), ctime,           &
             bu,bv,bww,bpt,bp,bqv,bqc,bqr,bqi,bqs,bqh, ierr)
  IF (ierr /= 0) THEN
    PRINT *,'ierr2=', ierr,' file',filename2(is)(1:lenfil2(is)),        &
                          '  not read, STOP!'
    STOP
  END IF

  IF (irandom == 0) THEN
!
!-----------------------------------------------------------------------
!
!   Find difference = first data - 2nd data set and store in 1st set
!       a=a-b
!  To reduce memory requirements, the difference fields are
!  written to the same arrays as the first data set.
!
!-----------------------------------------------------------------------
!
  CALL difebc(nx,ny,nz,0,tem1,                                          &
                au, av, aww, apt, ap,                                   &
                aqv, aqc, aqr, aqi, aqs, aqh,                           &
                bu, bv, bww, bpt, bp,                                   &
                bqv, bqc, bqr, bqi, bqs, bqh,                           &
                au, av, aww, apt, ap,                                   &
                aqv, aqc, aqr, aqi, aqs, aqh,                           &
                ireturn )

  ELSE IF (irandom == 1) THEN    ! add random perturbation

  do k=1,nz-1
  do i=1,nx
  do j=1,ny-1
  iseed=mod(iseed*7141+54773,259200)
  au(i,j,k) = 2.*0.10*10.0*(iseed/259199.-.5)
  enddo
  enddo
  enddo
  do k=1,nz-1
  do i=1,nx-1
  do j=1,ny
  iseed=mod(iseed*7141+54773,259200)
  av(i,j,k) = 2.*0.10*10.0*(iseed/259199.-.5)
  enddo
  enddo
  enddo
  do k=1,nz-1
  do i=1,nx-1
  do j=1,ny-1
  iseed=mod(iseed*7141+54773,259200)
  apt(i,j,k) = 2.*0.10*11.0*(iseed/259199.-.5)
  enddo
  enddo
  enddo
  do k=1,nz-1
  do i=1,nx-1
  do j=1,ny-1
  iseed=mod(iseed*7141+54773,259200)
  ap(i,j,k) = 2.*0.10*1000.0*(iseed/259199.-.5)
  enddo
  enddo
  enddo
  do k=1,nz-1
  do i=1,nx-1
  do j=1,ny-1
  iseed=mod(iseed*7141+54773,259200)
  aqv(i,j,k) = 2.*0.10*2e-3*(iseed/259199.-.5)
  enddo
  enddo
  enddo

  aww = 0.0
  aqc = 0.0
  aqr = 0.0
  aqi = 0.0
  aqs = 0.0
  aqh = 0.0

  END IF
!
!-----------------------------------------------------------------------
!
!  Read in the base data array if (iread=/=0)
!
!-----------------------------------------------------------------------
!
  IF (iread /= 0) THEN
    ierr=0
!    CALL readexbc(nx,ny,nz,dx,dy,dz,ctrlat,ctrlon,                      &
!               filename3(is)(1:lenfil3(is)),lenfil3(is), ctime,         &
!               ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd,                         &
!               qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd,qhbcrd,               &
!               bu,bv,bww,bpt,bp,bqv,bqc,bqr,bqi,bqs,bqh, ierr)
    CALL readexbc(nx,ny,nz,                                             &
               filename3(is)(1:lenfil3(is)),lenfil3(is), ctime,         &
               bu,bv,bww,bpt,bp,bqv,bqc,bqr,bqi,bqs,bqh, ierr)
    IF (ierr /= 0) THEN
      PRINT *,'ierr3=', ierr,' file',filename3(is)(1:lenfil3(is)),      &
                            '  not read, STOP!'
      STOP
    END IF
  END IF
!-----------------------------------------------------------------------
!  Find the ptb field =   3rd file  +  difference/n
!       b=b+a/n
!  To reduce memory requirements, the perturbation fields are
!  written to the same arrays as the base state fields.
!
!-----------------------------------------------------------------------
!
  IF (iorder /= 0) THEN
    CALL difebc(nx,ny,nz,iorder,tem1,                                   &
                  bu, bv, bww, bpt, bp,                                 &
                  bqv, bqc, bqr, bqi, bqs, bqh,                         &
                  au, av, aww, apt, ap,                                 &
                  aqv, aqc, aqr, aqi, aqs, aqh,                         &
                  bu, bv, bww, bpt, bp,                                 &
                  bqv, bqc, bqr, bqi, bqs, bqh,                         &
                  ireturn )
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          IF (k /= 1) THEN
            IF (bp(i,j,k) >= bp(i,j,k-1)) bp(i,j,k)=bp(i,j,k-1)-1.0
          END IF
          IF (bqv(i,j,k) <= 1.0E-8) THEN
            bqv(i,j,k)=1.0E-8
          END IF
        END DO
      END DO
    END DO
  END IF
!
	print *, filnmout(is)(1:lenfilo(is))
	print *, ctime
	print *, 'ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd:'
	print *, ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd

  CALL writexbc(nx,ny,nz,filnmout(is)(1:lenfilo(is)),lenfilo(is),       &
             ctime,                                                     &
             ubcrd,vbcrd,wbcrd,ptbcrd,prbcrd,                           &
             qvbcrd,qcbcrd,qrbcrd,qibcrd,qsbcrd,qhbcrd,                 &
             bu,bv,bww,bpt,bp,bqv,bqc,bqr,bqi,bqs,bqh)

  IF (is < isets) GO TO 101
  STOP

END PROGRAM arpsensbc
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE DIFEBC                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE difebc(nx,ny,nz,nscale,tem1,                                 & 2,11
           uprt, vprt, wprt, ptprt, pprt,                               &
           qvprt, qc, qr, qi, qs, qh,                                   &
           auprt, avprt, awprt, aptprt, apprt,                          &
           aqvprt, aqc, aqr, aqi, aqs, aqh,                             &
           duprt, dvprt, dwprt, dptprt, dpprt,                          &
           dqvprt, dqc, dqr, dqi, dqs, dqh,                             &
           ireturn)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE: du=u+au/nscale (nscale=/=0) or du=u-au (nscale==0)
!      (this is done by the subroutine SUBTRBC)
!
!  Subtract the forecast fields from the interpolated verification
!  fields (names beginning with "a") and output to the difference
!  fields (names beginning with "d").  The input and difference
!  fields may share the same storage location.  For this subroutine
!  it is assumed the forecast and corresponding verification
!  data are at the same physical location, however, the physical
!  location may differ between variables. That is uprt and auprt
!  are at the same location, but that may differ from pprt and apprt.
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster  Ou School of Meteorology. April 1992
!
!  MODIFICATION HISTORY:
!   14 May 1992  (KB) changed from arps2.5 to arps3.0
!   03 Aug 1992  (KB) updated to account for changes in arps3.0
!
!   09/07/1995  (KB)
!   Added differencing of surface (soil) fields.
!
!-----------------------------------------------------------------------
!
!  INPUT :
!    nx,ny,nz Array dimensions for forecast field.
!
!    FORECAST FIELDS:
!
!    uprt     perturbation x component of velocity (m/s)
!    vprt     perturbation y component of velocity (m/s)
!    wprt     perturbation vertical component of velocity in Cartesian
!             coordinates (m/s).
!
!    ptprt    perturbation potential temperature (K)
!    pprt     perturbation pressure (Pascal)
!
!    qvprt    perturbation water vapor mixing ratio (kg/kg)
!    qc       Cloud water mixing ratio (kg/kg)
!    qr       Rainwater mixing ratio (kg/kg)
!    qi       Cloud ice mixing ratio (kg/kg)
!    qs       Snow mixing ratio (kg/kg)
!    qh       Hail mixing ratio (kg/kg)
!
!    INTERPOLATED VERIFICATION FIELDS:
!
!    auprt     perturbation x component of velocity (m/s)
!    avprt     perturbation y component of velocity (m/s)
!    awprt     perturbation vertical component of velocity in Cartesian
!              coordinates (m/s).
!
!    aptprt    perturbation potential temperature (K)
!    apprt     perturbation pressure (Pascal)
!
!    aqvprt    perturbation water vapor mixing ratio (kg/kg)
!    aqc       Cloud water mixing ratio (kg/kg)
!    aqr       Rainwater mixing ratio (kg/kg)
!    aqi       Cloud ice mixing ratio (kg/kg)
!    aqs       Snow mixing ratio (kg/kg)
!    aqh       Hail mixing ratio (kg/kg)
!
!  OUTPUT :
!
!    DIFFERENCE FIELDS (may share storage with forecast fields
!                       or interpolated fields in calling program):
!
!    duprt     perturbation x component of velocity (m/s)
!    dvprt     perturbation y component of velocity (m/s)
!    dwprt     perturbation vertical component of velocity in Cartesian
!              coordinates (m/s).
!
!    dptprt    perturbation potential temperature (K)
!    dpprt     perturbation pressure (Pascal)
!
!    dqvprt    perturbation water vapor mixing ratio (kg/kg)
!    dqc       Cloud water mixing ratio (kg/kg)
!    dqr       Rainwater mixing ratio (kg/kg)
!    dqi       Cloud ice mixing ratio (kg/kg)
!    dqs       Snow mixing ratio (kg/kg)
!    dqh       Hail mixing ratio (kg/kg)
!
!
!    tem1      Work array
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz         ! 3 dimensions of array
  INTEGER :: nscale
!
!-----------------------------------------------------------------------
!
!  Model Arrays
!
!-----------------------------------------------------------------------
!
  REAL :: uprt   (nx,ny,nz)    ! Perturbation u-velocity (m/s)
  REAL :: vprt   (nx,ny,nz)    ! Perturbation v-velocity (m/s)
  REAL :: wprt   (nx,ny,nz)    ! Perturbation w-velocity (m/s)
  REAL :: ptprt  (nx,ny,nz)    ! Perturbation potential temperature (K)
  REAL :: pprt   (nx,ny,nz)    ! Perturbation pressure (Pascal)
  REAL :: qvprt  (nx,ny,nz)    ! Perturbation water vapor specific humidity
  REAL :: qc     (nx,ny,nz)    ! Cloud water mixing ratio (kg/kg)
  REAL :: qr     (nx,ny,nz)    ! Rain water mixing ratio (kg/kg)
  REAL :: qi     (nx,ny,nz)    ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs     (nx,ny,nz)    ! Snow mixing ratio (kg/kg)
  REAL :: qh     (nx,ny,nz)    ! Hail mixing ratio (kg/kg)

!-----------------------------------------------------------------------
!
!  Verification data interpolated to model grid
!
!-----------------------------------------------------------------------
!
  REAL :: auprt  (nx,ny,nz)    ! Perturbation u-velocity (m/s)
  REAL :: avprt  (nx,ny,nz)    ! Perturbation v-velocity (m/s)
  REAL :: awprt  (nx,ny,nz)    ! Perturbation w-velocity (m/s)
  REAL :: aptprt (nx,ny,nz)    ! Perturbation potential temperature (K)
  REAL :: apprt  (nx,ny,nz)    ! Perturbation pressure (Pascal)
  REAL :: aqvprt (nx,ny,nz)    ! Perturbation water vapor specific humidity
  REAL :: aqc    (nx,ny,nz)    ! Cloud water mixing ratio (kg/kg)
  REAL :: aqr    (nx,ny,nz)    ! Rain water mixing ratio (kg/kg)
  REAL :: aqi    (nx,ny,nz)    ! Cloud ice mixing ratio (kg/kg)
  REAL :: aqs    (nx,ny,nz)    ! Snow mixing ratio (kg/kg)
  REAL :: aqh    (nx,ny,nz)    ! Hail mixing ratio (kg/kg)


!-----------------------------------------------------------------------
!
!  Difference arrays
!
!-----------------------------------------------------------------------
!
  REAL :: duprt  (nx,ny,nz)    ! perturbation x component of velocity (m/s)
  REAL :: dvprt  (nx,ny,nz)    ! perturbation y component of velocity (m/s)
  REAL :: dwprt  (nx,ny,nz)    ! perturbation vertical component of
                               ! velocity in Cartesian coordinates (m/s)
  REAL :: dptprt (nx,ny,nz)    ! perturbation potential temperature (K)
  REAL :: dpprt  (nx,ny,nz)    ! perturbation pressure (Pascal)
  REAL :: dqvprt (nx,ny,nz)    ! perturbation water vapor mixing ratio (kg/kg)
  REAL :: dqc    (nx,ny,nz)    ! Cloud water mixing ratio (kg/kg)
  REAL :: dqr    (nx,ny,nz)    ! Rainwater mixing ratio (kg/kg)
  REAL :: dqi    (nx,ny,nz)    ! Cloud ice mixing ratio (kg/kg)
  REAL :: dqs    (nx,ny,nz)    ! Snow mixing ratio (kg/kg)
  REAL :: dqh    (nx,ny,nz)    ! Hail mixing ratio (kg/kg)

  REAL :: tem1   (nx,ny,nz)    ! A work array
  INTEGER :: ireturn, i,j,k
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: is,js,ks,ie,je,ke
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  is=1
  js=1
  ks=1
  ie=nx-1
  je=ny-1
  ke=nz-1

!-----------------------------------------------------------------------
!
!  Scalars
!
!-----------------------------------------------------------------------

  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem1(i,j,k)=0.0
      END DO
    END DO
  END DO

  PRINT *, ' ptprt: '
  CALL subtrbc(nx,ny,nz, ptprt,tem1,aptprt,tem1,dptprt,nscale,          &
             is,js,ks,ie,je,ke)
  PRINT *, ' pprt: '
  CALL subtrbc(nx,ny,nz,  pprt, tem1, apprt, tem1, dpprt,nscale,        &
             is,js,ks,ie,je,ke)
  PRINT *, ' qvprt: '
  CALL subtrbc(nx,ny,nz, qvprt,tem1,aqvprt,tem1,dqvprt,nscale,          &
             is,js,ks,ie,je,ke)
  PRINT *, ' qc: '
  CALL subtrbc(nx,ny,nz,    qc, tem1,   aqc,  tem1,   dqc,nscale,       &
             is,js,ks,ie,je,ke)
  PRINT *, ' qr: '
  CALL subtrbc(nx,ny,nz,    qr, tem1,   aqr,  tem1,   dqr,nscale,       &
             is,js,ks,ie,je,ke)
  PRINT *, ' qi: '
  CALL subtrbc(nx,ny,nz,    qi, tem1,   aqi,  tem1,   dqi,nscale,       &
             is,js,ks,ie,je,ke)
  PRINT *, ' qs: '
  CALL subtrbc(nx,ny,nz,    qs, tem1,   aqs,  tem1,   dqs,nscale,       &
             is,js,ks,ie,je,ke)
  PRINT *, ' qh: '
  CALL subtrbc(nx,ny,nz,    qh, tem1,   aqh,  tem1,   dqh,nscale,       &
             is,js,ks,ie,je,ke)

!-----------------------------------------------------------------------
!
!  u wind components
!
!-----------------------------------------------------------------------

  ie=nx
  PRINT *, ' uprt: '
  CALL subtrbc(nx,ny,nz,uprt,tem1,auprt,tem1,duprt,nscale,              &
             is,js,ks,ie,je,ke)

!-----------------------------------------------------------------------
!
!  v wind components
!
!-----------------------------------------------------------------------

  ie=nx-1
  je=ny
  PRINT *, ' vprt: '
  CALL subtrbc(nx,ny,nz,vprt,tem1,avprt,tem1,dvprt,nscale,              &
             is,js,ks,ie,je,ke)

!-----------------------------------------------------------------------
!
!  w wind components
!
!-----------------------------------------------------------------------

  je=ny-1
  ke=nz
  PRINT *, ' wprt: '
  CALL subtrbc(nx,ny,nz,wprt,tem1,awprt,tem1,dwprt,nscale,              &
             is,js,ks,ie,je,ke)
!
  RETURN
END SUBROUTINE difebc
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE SUBTRBC                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE subtrbc(nx,ny,nz, a,abar,b,bbar,c,nscale,                    & 11
           istr,jstr,kstr,iend,jend,kend)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Subtracts 2 three-dimensional arrays, represented by
!  means plus perturbations.
!
!  AUTHOR: Keith Brewster  OU School of Meteorology.  Feb 1992
!
!  MODIFICATION HISTORY:
!   11 Aug 1992  (KB) changed from arps2.5 to arps3.0
!
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    a        perturbation data array
!    abar     mean data array
!    b        perturbation data array to subtract from a
!    bbar     mean data array to subtract from a
!
!  OUTPUT:
!    c        difference array a-b
!             (may share storage in calling program with array a or b)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! 3 dimensions of array
  INTEGER :: nscale

  REAL :: a   (nx,ny,nz)       ! data array
  REAL :: abar(nx,ny,nz)       ! base state of data array a
  REAL :: b   (nx,ny,nz)       ! data array to subtract from a
  REAL :: bbar(nx,ny,nz)       ! base state of data arrya b
  REAL :: c   (nx,ny,nz)       ! difference array a-b
  INTEGER :: istr,jstr,kstr
  INTEGER :: iend,jend,kend
  INTEGER :: i,j,k,imid,jmid,kmid
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  imid=nint(0.5*FLOAT(istr+iend))
  jmid=nint(0.5*FLOAT(jstr+jend))
  kmid=nint(0.5*FLOAT(kstr+kend))
!
!-----------------------------------------------------------------------
!
!  Tell us about a sample input point
!
!-----------------------------------------------------------------------
!
  IF (nscale == 0) THEN
    PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)),   &
                     ' b= ',(b(imid,jmid,kmid)+bbar(imid,jmid,kmid))
  ELSE
    PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)),   &
                     ' b= ',b(imid,jmid,kmid)
  END IF
!
!-----------------------------------------------------------------------
!
!  Subtraction
!
!-----------------------------------------------------------------------
!
  DO k=kstr,kend
    DO j=jstr,jend
      DO i=istr,iend
        IF (nscale == 0) THEN
          c(i,j,k)=a(i,j,k)+abar(i,j,k)-(b(i,j,k)+bbar(i,j,k))
        ELSE
          c(i,j,k)=a(i,j,k)+b(i,j,k)/FLOAT(nscale)
        END IF
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Tell us about a sample output point
!
!-----------------------------------------------------------------------
!
  IF (nscale == 0) THEN
    PRINT *, '         c= ',c(imid,jmid,kmid)
  ELSE
    PRINT *, '         c= ',c(imid,jmid,kmid)+abar(imid,jmid,kmid)
  END IF
  RETURN
END SUBROUTINE subtrbc