!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MKARPSVAR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, & 24,5
iorder,iprtopt,intropt,iloc,jloc,x_ext,y_ext, &
hgtext,zext,x2d,y2d,za,varext, &
zsnd,varsnd,arpsbar,arpsprt, &
dxfld,dyfld,rdxfld,rdyfld, &
slopey,alphay,betay, &
tem_ext)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Take data from external file on pressure levels and interpolate in
! the horizontal and vertical to ARPS height levels and horizontal
! locations. Then, form the ARPS mean and perturbation quantities.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! 5/01/1994.
!
! MODIFICATION HISTORY:
!
! 11/21/1994 (KB)
! Added full documentation.
!
! 2000/08/16 (Gene Bassett)
! Added intropt.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx_ext Number of grid points in the x-direction (east/west)
! for the external grid
! ny_ext Number of grid points in the y-direction (north/south)
! for the external grid
! nz_ext Number of grid points in the vertical
! for the external grid
!
! nx Number of grid points in the x-direction (east/west)
! for the ARPS grid
! ny Number of grid points in the y-direction (north/south)
! for the ARPS grid
! nz Number of grid points in the vertical
! for the ARPS grid
!
! iorder order of polynomial for interpolation (1, 2 or 3)
! intropt Option indicating to interpolate perturbation or total
! variables:
! = 1 Interpolate perturbation variables and add to base
! sounding (default);
! = 2 Interploate total variables.
! iprtopt Flag for producing a perturbation variable
! iprtopt = 1 Produce mean and perturbation field
! iprtopt = 0 Produce mean and total field
!
! iloc x-index location of ARPS grid point in the external array
! jloc y-index location of ARPS grid point in the external array
!
! x2d x coordinate of ARPS grid point in external coordinate
! y2d x coordinate of ARPS grid point in external coordinate
!
! zext Array of heights of the external grid
! za Array of heights of ARPS physical heights
! zsnd 1-D array of height representing a mean sounding
! over domain
! varsnd 1-D array of variable representing a mean sounding
! over domain
!
! OUTPUT:
!
! arpsbar 3-D array of mean field on ARPS levels
! arpsprt 3-D array of perturbation (iprtopt=1) or
! total (iprtopt=0) variable at ARPS grid locations
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'bndry.inc'
INCLUDE 'mp.inc'
!
!-----------------------------------------------------------------------
!
! Input variables
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext,ny_ext,nz_ext ! external grid dimensions
INTEGER :: nx,ny,nz ! ARPS grid dimensions
INTEGER :: lvlprof ! levels in mean sounding
INTEGER :: iorder ! interpolating polynomial order
INTEGER :: intropt ! interpolation option
INTEGER :: iprtopt ! perturbation generation flag
INTEGER :: iloc(nx,ny) ! external x-index of ARPS grid point
INTEGER :: jloc(nx,ny) ! external y-index of ARPS grid point
REAL :: x_ext(nx_ext) ! external x-coord
REAL :: y_ext(ny_ext) ! external y-coord
REAL :: hgtext(nx_ext,ny_ext,nz_ext) ! heights of external levels
REAL :: zext(nx,ny,nz_ext) ! heights of external levels
! interpolated to ARPS grid locs
REAL :: x2d(nx,ny)
REAL :: y2d(nx,ny)
REAL :: za(nx,ny,nz) ! ARPS physical heights
REAL :: varext(nx_ext,ny_ext,nz_ext) ! variable to convert
REAL :: zsnd(lvlprof) ! 1-D array of level heights
REAL :: varsnd(lvlprof) ! 1-D array of level-means
!
!-----------------------------------------------------------------------
!
! Output variables
!
!-----------------------------------------------------------------------
!
REAL :: arpsbar( nx, ny, nz) ! 3-D array of level-means
REAL :: arpsprt( nx, ny, nz) ! Output array, perturnbation variable
! or total variable (see iprtopt)
!
!-----------------------------------------------------------------------
!
! Temporary work arrays
!
!-----------------------------------------------------------------------
!
REAL :: dxfld(nx_ext)
REAL :: dyfld(ny_ext)
REAL :: rdxfld(nx_ext)
REAL :: rdyfld(ny_ext)
REAL :: slopey(nx_ext,ny_ext,nz_ext)
REAL :: alphay(nx_ext,ny_ext,nz_ext)
REAL :: betay(nx_ext,ny_ext,nz_ext)
REAL :: tem_ext(nx_ext,ny_ext,nz_ext)
REAL :: tem1(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,ia,ja,ka,kl
REAL :: wlow
REAL :: topprt,botprt,arpstop,arpsbot
REAL :: pntint2d
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Subtract horizontal mean from external fields.
!
!-----------------------------------------------------------------------
!
IF (intropt == 1) THEN
DO j=1,ny_ext
DO i=1,nx_ext
DO k=1,nz_ext
DO kl=2,lvlprof-1
IF(zsnd(kl) > hgtext(i,j,k)) EXIT
END DO
! 11 CONTINUE
wlow=(zsnd(kl)-hgtext(i,j,k))/ &
(zsnd(kl)-zsnd(kl-1))
tem_ext(i,j,k)=varext(i,j,k)- &
((1.-wlow)*varsnd(kl) + wlow*varsnd(kl-1))
END DO
END DO
END DO
ELSE
tem_ext = varext
ENDIF
!
!-----------------------------------------------------------------------
!
! Compute derivative terms
!
!-----------------------------------------------------------------------
!
CALL setdrvy
(nx_ext,ny_ext,nz_ext, &
1,nx_ext,1,ny_ext,1,nz_ext, &
dyfld,rdyfld,tem_ext, &
slopey,alphay,betay)
!
!-----------------------------------------------------------------------
!
! Loop through all ARPS grid points
!
!-----------------------------------------------------------------------
!
DO ja=1,ny
DO ia=1,nx
!
!-----------------------------------------------------------------------
!
! Interpolate from the mean sounding to get arpsbar
!
!-----------------------------------------------------------------------
!
DO ka=1,nz
DO kl=2,lvlprof-1
IF(zsnd(kl) > za(ia,ja,ka)) EXIT
END DO
! 51 CONTINUE
wlow=(zsnd(kl)-za(ia,ja,ka))/ &
(zsnd(kl)-zsnd(kl-1))
arpsbar(ia,ja,ka)= &
(1.-wlow)*varsnd(kl) + wlow*varsnd(kl-1)
END DO
!
!-----------------------------------------------------------------------
!
! Find vertical location
! and interpolate in vertical between two horizontal
! interpolations on the external grid surfaces.
!
! Extrapolation is done by assuming the perturbation
! from mean is constant.
!
!-----------------------------------------------------------------------
!
botprt=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,1), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,1),alphay(1,1,1),betay(1,1,1))
topprt=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,nz_ext), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,nz_ext), &
alphay(1,1,nz_ext),betay(1,1,nz_ext))
DO ka=1,nz
IF(za(ia,ja,ka) < zext(ia,ja,1)) THEN
arpsprt(ia,ja,ka)=botprt
ELSE IF(za(ia,ja,ka) > zext(ia,ja,nz_ext)) THEN
arpsprt(ia,ja,ka)=topprt
ELSE
DO kl=2,nz_ext-1
IF(zext(ia,ja,kl) > za(ia,ja,ka)) EXIT
END DO
! 351 CONTINUE
wlow=(zext(ia,ja,kl)- za(ia,ja,ka))/ &
(zext(ia,ja,kl)-zext(ia,ja,kl-1))
arpstop=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,kl), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,kl), &
alphay(1,1,kl),betay(1,1,kl))
arpsbot=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,kl-1), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,kl-1), &
alphay(1,1,kl-1),betay(1,1,kl-1))
arpsprt(ia,ja,ka)=(1.-wlow)*arpstop+wlow*arpsbot
END IF
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! If total quantity, rather than perturbation quantity desired,
! add the level mean to the interpolated perturbation variable
! at each grid point.
!
!-----------------------------------------------------------------------
!
IF(iprtopt == 0 .and. intropt == 1) THEN
arpsprt = arpsprt + arpsbar
ELSE IF(iprtopt == 1 .and. intropt == 2) THEN
arpsprt = arpsprt - arpsbar
END IF
CALL mpsendrecv2dew(arpsbar, nx, ny, nz, ebc, wbc, 1, tem1)
CALL mpsendrecv2dew(arpsprt, nx, ny, nz, ebc, wbc, 1, tem1)
CALL mpsendrecv2dns(arpsbar, nx, ny, nz, nbc, sbc, 1, tem1)
CALL mpsendrecv2dns(arpsprt, nx, ny, nz, nbc, sbc, 1, tem1)
!
RETURN
END SUBROUTINE mkarpsvar
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MKARPSVAR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE mkarpsvar1(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, & 4,2
iorder,iprtopt,intropt,iloc,jloc,x_ext,y_ext, &
hgtext,zext,x2d,y2d,za,varext, &
zsnd,varsnd,arpsbar,arpsprt, &
trn_ext,var_h0,h0, &
dxfld,dyfld,rdxfld,rdyfld, &
slopey,alphay,betay, &
tem_ext)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Take data from external file on pressure levels and interpolate in
! the horizontal and vertical to ARPS height levels and horizontal
! locations. Then, form the ARPS mean and perturbation quantities.
! Near surface fields are used in low level interpolation.
! (Written from Keith Brewster's original mkarpsvar code)
!
!-----------------------------------------------------------------------
!
! AUTHOR: Fanyou Kong (Modified based on Keith Brewster's mkarpsvar)
! 10/25/2003
!
! MODIFICATION HISTORY:
!
! 05/20/2004 (Yunheng Wang)
! Removed the unused parameter extsfcopt.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx_ext Number of grid points in the x-direction (east/west)
! for the external grid
! ny_ext Number of grid points in the y-direction (north/south)
! for the external grid
! nz_ext Number of grid points in the vertical
! for the external grid
!
! nx Number of grid points in the x-direction (east/west)
! for the ARPS grid
! ny Number of grid points in the y-direction (north/south)
! for the ARPS grid
! nz Number of grid points in the vertical
! for the ARPS grid
!
! iorder order of polynomial for interpolation (1, 2 or 3)
! intropt Option indicating to interpolate perturbation or total
! variables:
! = 1 Interpolate perturbation variables and add to base
! sounding (default);
! = 2 Interploate total variables.
! iprtopt Flag for producing a perturbation variable
! iprtopt = 1 Produce mean and perturbation field
! iprtopt = 0 Produce mean and total field
!
! iloc x-index location of ARPS grid point in the external array
! jloc y-index location of ARPS grid point in the external array
!
! x2d x coordinate of ARPS grid point in external coordinate
! y2d x coordinate of ARPS grid point in external coordinate
!
! zext Array of heights of the external grid
! za Array of heights of ARPS physical heights
! zsnd 1-D array of height representing a mean sounding
! over domain
! varsnd 1-D array of variable representing a mean sounding
! over domain
! trn_ext External grid terrain height
! var_h0 Near surface variable at h0 height above ground
! h0 Near surface variable height (m)
!
! OUTPUT:
!
! arpsbar 3-D array of mean field on ARPS levels
! arpsprt 3-D array of perturbation (iprtopt=1) or
! total (iprtopt=0) variable at ARPS grid locations
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Input variables
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext,ny_ext,nz_ext ! external grid dimensions
INTEGER :: nx,ny,nz ! ARPS grid dimensions
INTEGER :: lvlprof ! levels in mean sounding
INTEGER :: iorder ! interpolating polynomial order
INTEGER :: intropt ! interpolation option
INTEGER :: iprtopt ! perturbation generation flag
INTEGER :: iloc(nx,ny) ! external x-index of ARPS grid point
INTEGER :: jloc(nx,ny) ! external y-index of ARPS grid point
REAL :: x_ext(nx_ext) ! external x-coord
REAL :: y_ext(ny_ext) ! external y-coord
REAL :: hgtext(nx_ext,ny_ext,nz_ext) ! heights of external levels
REAL :: zext(nx,ny,nz_ext) ! heights of external levels
! interpolated to ARPS grid locs
REAL :: x2d(nx,ny)
REAL :: y2d(nx,ny)
REAL :: za(nx,ny,nz) ! ARPS physical heights
REAL :: varext(nx_ext,ny_ext,nz_ext) ! variable to convert
REAL :: zsnd(lvlprof) ! 1-D array of level heights
REAL :: varsnd(lvlprof) ! 1-D array of level-means
REAL :: trn_ext(nx_ext,ny_ext) ! external terrain height (m)
REAL :: var_h0(nx_ext,ny_ext) ! near surface values
REAL :: h0 ! near surface varaiable height (m)
REAL :: tem_var_h0(nx_ext,ny_ext)
!
!-----------------------------------------------------------------------
!
! Output variables
!
!-----------------------------------------------------------------------
!
REAL :: arpsbar( nx, ny, nz) ! 3-D array of level-means
REAL :: arpsprt( nx, ny, nz) ! Output array, perturnbation variable
! or total variable (see iprtopt)
!
!-----------------------------------------------------------------------
!
! Temporary work arrays
!
!-----------------------------------------------------------------------
!
REAL :: dxfld(nx_ext)
REAL :: dyfld(ny_ext)
REAL :: rdxfld(nx_ext)
REAL :: rdyfld(ny_ext)
REAL :: slopey(nx_ext,ny_ext,nz_ext)
REAL :: alphay(nx_ext,ny_ext,nz_ext)
REAL :: betay(nx_ext,ny_ext,nz_ext)
REAL :: tem_ext(nx_ext,ny_ext,nz_ext)
REAL :: tem_1(nx_ext,ny_ext)
REAL :: tem_2(nx_ext,ny_ext)
REAL :: tem_3(nx_ext,ny_ext)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,ia,ja,ka,kl,ks
REAL :: wlow,zsfc
REAL :: topprt,botprt,arpstop,arpsbot
REAL :: pntint2d
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Subtract horizontal mean from external fields.
!
!-----------------------------------------------------------------------
!
IF (intropt == 1) THEN
DO j=1,ny_ext
DO i=1,nx_ext
DO k=1,nz_ext
DO kl=2,lvlprof-1
IF(zsnd(kl) > hgtext(i,j,k)) EXIT
END DO
! 11 CONTINUE
wlow=(zsnd(kl)-hgtext(i,j,k))/ &
(zsnd(kl)-zsnd(kl-1))
tem_ext(i,j,k)=varext(i,j,k)- &
((1.-wlow)*varsnd(kl) + wlow*varsnd(kl-1))
END DO
! processing near surface variable
zsfc=trn_ext(i,j) + h0
DO kl=2,lvlprof-1
IF(zsnd(kl) > zsfc) EXIT
END DO
wlow=(zsnd(kl)-zsfc)/(zsnd(kl)-zsnd(kl-1))
tem_var_h0(i,j)=var_h0(i,j)- &
((1.-wlow)*varsnd(kl) + wlow*varsnd(kl-1))
END DO
END DO
ELSE
tem_ext = varext
tem_var_h0=var_h0
ENDIF
!
!-----------------------------------------------------------------------
!
! Compute derivative terms
!
!-----------------------------------------------------------------------
!
CALL setdrvy(nx_ext,ny_ext,nz_ext, &
1,nx_ext,1,ny_ext,1,nz_ext, &
dyfld,rdyfld,tem_ext, &
slopey,alphay,betay)
! processing near surface variable
CALL setdrvy(nx_ext,ny_ext,1, &
1,nx_ext,1,ny_ext,1,1, &
dyfld,rdyfld,tem_var_h0, &
tem_1,tem_2,tem_3)
!
!-----------------------------------------------------------------------
!
! Loop through all ARPS grid points
!
!-----------------------------------------------------------------------
!
DO ja=1,ny
DO ia=1,nx
!
!-----------------------------------------------------------------------
!
! Interpolate from the mean sounding to get arpsbar
!
!-----------------------------------------------------------------------
!
DO ka=1,nz
DO kl=2,lvlprof-1
IF(zsnd(kl) > za(ia,ja,ka)) EXIT
END DO
! 51 CONTINUE
wlow=(zsnd(kl)-za(ia,ja,ka))/ &
(zsnd(kl)-zsnd(kl-1))
arpsbar(ia,ja,ka)= &
(1.-wlow)*varsnd(kl) + wlow*varsnd(kl-1)
END DO
!
!-----------------------------------------------------------------------
!
! Find vertical location
! and interpolate in vertical between two horizontal
! interpolations on the external grid surfaces.
!
! Extrapolation is done by assuming the perturbation
! from mean is constant.
!
!-----------------------------------------------------------------------
!
! botprt=pntint2d(nx_ext,ny_ext, &
! 2,nx_ext-1,2,ny_ext-1, &
! iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
! iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,1), &
! dxfld,dyfld,rdxfld,rdyfld, &
! slopey(1,1,1),alphay(1,1,1),betay(1,1,1))
topprt=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,nz_ext), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,nz_ext), &
alphay(1,1,nz_ext),betay(1,1,nz_ext))
! find first external layer index (ks) above zsfc
zsfc=0.5*(za(ia,ja,1)+za(ia,ja,2)) + h0
DO ks=1,nz_ext
IF(zext(ia,ja,ks) > zsfc) EXIT ! ks determined
END DO
! vertical loop (arps domain)
DO ka=2,nz
! IF(za(ia,ja,ka) < zext(ia,ja,1)) THEN
! arpsprt(ia,ja,ka)=botprt
! ELSE IF(za(ia,ja,ka) > zext(ia,ja,nz_ext)) THEN
IF(za(ia,ja,ka) > zext(ia,ja,nz_ext)) THEN
arpsprt(ia,ja,ka)=topprt
ELSE
DO kl=1,nz_ext-1
IF(zext(ia,ja,kl) > za(ia,ja,ka)) EXIT
END DO
! 351 CONTINUE
arpstop=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,kl), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,kl), &
alphay(1,1,kl),betay(1,1,kl))
IF(za(ia,ja,ka) >= zext(ia,ja,ks)) THEN ! using external layer
wlow=(zext(ia,ja,kl)- za(ia,ja,ka))/ &
(zext(ia,ja,kl)-zext(ia,ja,kl-1))
arpsbot=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,kl-1), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,kl-1), &
alphay(1,1,kl-1),betay(1,1,kl-1))
ELSE ! using near surface value
wlow=(zext(ia,ja,kl)- za(ia,ja,ka))/ &
(zext(ia,ja,kl)-zsfc)
! The line below reflects rare situation when za(ka)<zsfc (occurring for u,v)
IF ( za(ia,ja,ka) <= zsfc ) wlow=1.0
arpsbot=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_var_h0, &
dxfld,dyfld,rdxfld,rdyfld, &
tem_1,tem_2,tem_3)
END IF
arpsprt(ia,ja,ka)=(1.-wlow)*arpstop+wlow*arpsbot
END IF
END DO
! Extropolating to k=1 level
wlow=(za(ia,ja,3)-za(ia,ja,1))/(za(ia,ja,3)-za(ia,ja,2))
arpsprt(ia,ja,1)=wlow*arpsprt(ia,ja,2)+ &
(1.-wlow)*arpsprt(ia,ja,3)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! If total quantity, rather than perturbation quantity desired,
! add the level mean to the interpolated perturbation variable
! at each grid point.
!
!-----------------------------------------------------------------------
!
IF(iprtopt == 0 .and. intropt == 1) THEN
arpsprt = arpsprt + arpsbar
ELSE IF(iprtopt == 1 .and. intropt == 2) THEN
arpsprt = arpsprt - arpsbar
END IF
!
RETURN
END SUBROUTINE mkarpsvar1
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MKARPSVLZ ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE mkarpsvlz(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, & 2,1
iorder,iprtopt,intropt,iloc,jloc,x_ext,y_ext, &
hgtext,zext,x2d,y2d,za,varext, &
zsnd,vlnsnd,arpsbar,arpsprt, &
dxfld,dyfld,rdxfld,rdyfld, &
slopey,alphay,betay, &
tem_ext)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Take data from external file on pressure levels and interpolate in
! the horizontal and vertical to ARPS height levels and horizontal
! locations. Then, form the ARPS mean and perturbation quantities.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! 5/01/1994.
!
! MODIFICATION HISTORY:
!
! 11/21/1994 (KB)
! Added full documentation.
!
! 2000/08/16 (Gene Bassett)
! Added intropt.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx_ext Number of grid points in the x-direction (east/west)
! for the external grid
! ny_ext Number of grid points in the y-direction (north/south)
! for the external grid
! nz_ext Number of grid points in the vertical
! for the external grid
!
! nx Number of grid points in the x-direction (east/west)
! for the ARPS grid
! ny Number of grid points in the y-direction (north/south)
! for the ARPS grid
! nz Number of grid points in the vertical
! for the ARPS grid
!
! iorder order of polynomial for interpolation (1, 2 or 3)
! intropt Option indicating to interpolate perturbation or total
! variables:
! = 1 Interpolate perturbation variables and add to base
! sounding (default);
! = 2 Interploate total variables.
! iprtopt Flag for producing a perturbation variable
! iprtopt = 1 Produce mean and perturbation field
! iprtopt = 0 Produce mean and total field
!
! iloc x-index location of ARPS grid point in the external array
! jloc y-index location of ARPS grid point in the external array
!
! x2d x coordinate of ARPS grid point in external coordinate
! y2d x coordinate of ARPS grid point in external coordinate
!
! zext Array of heights of the external grid
! za Array of heights of ARPS physical heights
! zsnd 1-D array of height representing a mean sounding
! over domain
! vlnsnd 1-D array of variable representing the log of the
! mean sounding of the variable.
!
! OUTPUT:
!
! arpsbar 3-D array of mean field on ARPS levels
! arpsprt 3-D array of perturbation (iprtopt=1) or
! total (iprtopt=0) variable at ARPS grid locations
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Input variables
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext,ny_ext,nz_ext ! external grid dimensions
INTEGER :: nx,ny,nz ! ARPS grid dimensions
INTEGER :: lvlprof ! levels in mean sounding
INTEGER :: iorder ! interpolating polynomial order
INTEGER :: intropt ! interpolation option
INTEGER :: iprtopt ! perturbation generation flag
INTEGER :: iloc(nx,ny) ! external x-index of ARPS grid point
INTEGER :: jloc(nx,ny) ! external y-index of ARPS grid point
REAL :: x_ext(nx_ext) ! external x-coord
REAL :: y_ext(ny_ext) ! external y-coord
REAL :: hgtext(nx_ext,ny_ext,nz_ext) ! heights of external levels
REAL :: zext(nx,ny,nz_ext) ! heights of external levels
! interpolated to ARPS grid locs
REAL :: x2d(nx,ny)
REAL :: y2d(nx,ny)
REAL :: za(nx,ny,nz) ! ARPS physical heights
REAL :: varext(nx_ext,ny_ext,nz_ext) ! variable to convert
REAL :: zsnd(lvlprof) ! 1-D array of level heights
REAL :: vlnsnd(lvlprof) ! 1-D array of level-means
!
!-----------------------------------------------------------------------
!
! Output variables
!
!-----------------------------------------------------------------------
!
REAL :: arpsbar( nx, ny, nz) ! 3-D array of level-means
REAL :: arpsprt( nx, ny, nz) ! Output array, perturnbation variable
! or total variable (see iprtopt)
!
!-----------------------------------------------------------------------
!
! Temporary work arrays
!
!-----------------------------------------------------------------------
!
REAL :: dxfld(nx_ext)
REAL :: dyfld(ny_ext)
REAL :: rdxfld(nx_ext)
REAL :: rdyfld(ny_ext)
REAL :: slopey(nx_ext,ny_ext,nz_ext)
REAL :: alphay(nx_ext,ny_ext,nz_ext)
REAL :: betay(nx_ext,ny_ext,nz_ext)
REAL :: tem_ext(nx_ext,ny_ext,nz_ext)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,ia,ja,ka,kl
REAL :: wlow
REAL :: topprt,botprt,arpstop,arpsbot
REAL :: pntint2d
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Subtract horizontal mean from external fields.
!
!-----------------------------------------------------------------------
!
IF (intropt == 1) THEN
DO j=1,ny_ext
DO i=1,nx_ext
DO k=1,nz_ext
DO kl=2,lvlprof-1
IF(zsnd(kl) > hgtext(i,j,k)) EXIT
END DO
! 11 CONTINUE
wlow=(zsnd(kl)-hgtext(i,j,k))/ &
(zsnd(kl)-zsnd(kl-1))
tem_ext(i,j,k)=varext(i,j,k)-EXP( &
((1.-wlow)*vlnsnd(kl) + wlow*vlnsnd(kl-1)))
END DO
END DO
END DO
ELSE ! intropt=2
tem_ext = LOG(varext)
ENDIF
!
!-----------------------------------------------------------------------
!
! Compute derivative terms
!
!-----------------------------------------------------------------------
!
CALL setdrvy(nx_ext,ny_ext,nz_ext, &
1,nx_ext,1,ny_ext,1,nz_ext, &
dyfld,rdyfld,tem_ext, &
slopey,alphay,betay)
!
!-----------------------------------------------------------------------
!
! Loop through all ARPS grid points
!
!-----------------------------------------------------------------------
!
DO ja=1,ny
DO ia=1,nx
!
!-----------------------------------------------------------------------
!
! Interpolate from the mean sounding to get arpsbar
!
!-----------------------------------------------------------------------
!
DO ka=1,nz
DO kl=2,lvlprof-1
IF(zsnd(kl) > za(ia,ja,ka)) EXIT
END DO
! 51 CONTINUE
wlow=(zsnd(kl)-za(ia,ja,ka))/ &
(zsnd(kl)-zsnd(kl-1))
arpsbar(ia,ja,ka)=EXP( &
(1.-wlow)*vlnsnd(kl) + wlow*vlnsnd(kl-1))
END DO
!
!-----------------------------------------------------------------------
!
! Find vertical location
! and interpolate in vertical between two horizontal
! interpolations on the external grid surfaces.
!
! Extrapolation is done by assuming the perturbation
! from mean is constant.
!
!-----------------------------------------------------------------------
!
botprt=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,1), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,1),alphay(1,1,1),betay(1,1,1))
topprt=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,nz_ext), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,nz_ext), &
alphay(1,1,nz_ext),betay(1,1,nz_ext))
DO ka=1,nz
IF(za(ia,ja,ka) < zext(ia,ja,1)) THEN
arpsprt(ia,ja,ka)=botprt
ELSE IF(za(ia,ja,ka) > zext(ia,ja,nz_ext)) THEN
arpsprt(ia,ja,ka)=topprt
ELSE
DO kl=2,nz_ext-1
IF(zext(ia,ja,kl) > za(ia,ja,ka)) EXIT
END DO
! 351 CONTINUE
wlow=(zext(ia,ja,kl)- za(ia,ja,ka))/ &
(zext(ia,ja,kl)-zext(ia,ja,kl-1))
arpstop=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,kl), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,kl), &
alphay(1,1,kl),betay(1,1,kl))
arpsbot=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),tem_ext(1,1,kl-1), &
dxfld,dyfld,rdxfld,rdyfld, &
slopey(1,1,kl-1), &
alphay(1,1,kl-1),betay(1,1,kl-1))
arpsprt(ia,ja,ka)=(1.-wlow)*arpstop+wlow*arpsbot
END IF
END DO
END DO
END DO
IF(intropt == 2) THEN
arpsprt = EXP(arpsprt)
ENDIF
!
!-----------------------------------------------------------------------
!
! If total quantity, rather than perturbation quantity desired,
! add the level mean to the interpolated perturbation variable
! at each grid point.
!
!-----------------------------------------------------------------------
!
IF(iprtopt == 0 .and. intropt == 1) THEN
arpsprt = arpsprt + arpsbar
ELSE IF(iprtopt == 1 .and. intropt == 2) THEN
arpsprt = arpsprt - arpsbar
END IF
!
RETURN
END SUBROUTINE mkarpsvlz
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MKARPS2D ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE mkarps2d (nx_ext,ny_ext,nx,ny, & 34,1
iorder,iloc,jloc,x_ext,y_ext, &
x2d,y2d,varext,arpsvar, &
dxfld,dyfld,rdxfld,rdyfld, &
slopey,alphay,betay, &
tem_ext)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Take 2D surface data from external file and interpolate in
! the horizontal to form the ARPS 2D data set
!
!-----------------------------------------------------------------------
!
! AUTHOR: Fanyou Kong
! 6/10/1997.
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx_ext Number of grid points in the x-direction (east/west)
! for the external grid
! ny_ext Number of grid points in the y-direction (north/south)
! for the external grid
!
! nx Number of grid points in the x-direction (east/west)
! for the ARPS grid
! ny Number of grid points in the y-direction (north/south)
! for the ARPS grid
!
! iorder order of polynomial for interpolation (1, 2 or 3)
!
! iloc x-index location of ARPS grid point in the external array
! jloc y-index location of ARPS grid point in the external array
!
! x2d x coordinate of ARPS grid point in external coordinate
! y2d x coordinate of ARPS grid point in external coordinate
!
! varext external 2D data
!
! OUTPUT:
!
! arpsvar 2D array of variable at ARPS grid locations
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Input variables
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext,ny_ext ! external grid dimensions
INTEGER :: nx,ny ! ARPS grid dimensions
INTEGER :: iorder ! interpolating polynomial order
INTEGER :: iloc(nx,ny) ! external x-index of ARPS grid point
INTEGER :: jloc(nx,ny) ! external y-index of ARPS grid point
REAL :: x_ext(nx_ext) ! external x-coord
REAL :: y_ext(ny_ext) ! external y-coord
! interpolated to ARPS grid locs
REAL :: x2d(nx,ny)
REAL :: y2d(nx,ny)
REAL :: varext(nx_ext,ny_ext) ! 2D variable to convert
!
!-----------------------------------------------------------------------
!
! Output variables
!
!-----------------------------------------------------------------------
!
REAL :: arpsvar( nx, ny) ! 2-D array of ARPS grid data
!
!-----------------------------------------------------------------------
!
! Temporary work arrays
!
!-----------------------------------------------------------------------
!
REAL :: dxfld(nx_ext)
REAL :: dyfld(ny_ext)
REAL :: rdxfld(nx_ext)
REAL :: rdyfld(ny_ext)
REAL :: slopey(nx_ext,ny_ext)
REAL :: alphay(nx_ext,ny_ext)
REAL :: betay(nx_ext,ny_ext)
REAL :: tem_ext(nx_ext,ny_ext)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: ia,ja
REAL :: arpsdata
REAL :: pntint2d
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Compute derivative terms
!
!-----------------------------------------------------------------------
!
CALL setdrvy(nx_ext,ny_ext,1, &
1,nx_ext,1,ny_ext,1,1, &
dyfld,rdyfld,varext, &
slopey,alphay,betay)
!
!-----------------------------------------------------------------------
!
! Loop through all ARPS grid points
!
!-----------------------------------------------------------------------
!
DO ja=1,ny
DO ia=1,nx
!
!-----------------------------------------------------------------------
!
! Horizontal interpolation
!
!-----------------------------------------------------------------------
!
arpsdata=pntint2d(nx_ext,ny_ext, &
2,nx_ext-1,2,ny_ext-1, &
iorder,x_ext,y_ext,x2d(ia,ja),y2d(ia,ja), &
iloc(ia,ja),jloc(ia,ja),varext, &
dxfld,dyfld,rdxfld,rdyfld, &
slopey,alphay,betay)
arpsvar(ia,ja)=arpsdata
END DO
END DO
RETURN
END SUBROUTINE mkarps2d
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INITSOILEXT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE initsoilext(nx,ny,nzsoil,nzsoil_ext, nstyp,zpsoil, &
zpsoil_ext,tsoil,tsoil_ext,qsoil,qsoil_ext)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Take soil data from external file and interpolate in
! the vertical to form the ARPS 2D data set profile
!
!-----------------------------------------------------------------------
!
! AUTHOR: Jerry Brotzge
! 05/15/2002
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! for the ARPS grid
! ny Number of grid points in the y-direction (north/south)
! for the ARPS grid
! nzsoil Number of grid points in the soil
!
! nzsoil_ext Number of grid points in the soil external file
!
! nstyp Number of soil types
!
! zpsoil Physical depth of soil layers (m)
!
! zpsoil_ext Physical depth of external model soil layers (m)
!
! tsoil Soil temperature (K)
!
! qsoil Soil moisture (m**3/m**3)
!
! OUTPUT:
!
! tsoil Soil temperature profile (K)
!
! qsoil Soil moisture profile (m**3/m**3)
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!
!-----------------------------------------------------------------------
!
! Input/output variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,kk
INTEGER :: nx,ny ! ARPS grid dimensions
INTEGER :: nzsoil ! ARPS soil levels
INTEGER :: nzsoil_ext ! External file soil levels
INTEGER :: nstyp ! Number of soil types
!
REAL :: zpsoil(nx,ny,nzsoil) ! Physical depth of ARPS soil layers
REAL :: zpsoil_ext(nx,ny,nzsoil_ext) ! Physical depth of ext. file soil layers
REAL :: tsoil(nx,ny,nzsoil) ! Soil temperature (K)
REAL :: qsoil(nx,ny,nzsoil) ! Soil moisture (m**3/m**3)
REAL :: tsoil_ext(nx,ny,nzsoil_ext) ! External file soil temperature (K)
REAL :: qsoil_ext(nx,ny,nzsoil_ext) ! External file soil moisture (m**3/m**3)
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
REAL :: dampdepth
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'grid.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Vertical interpolation
!
!-----------------------------------------------------------------------
!----------------------------------------------------------------
! Set top boundary condition at level = 1
!-------------------------------------------------------------
DO j=1,ny
DO i=1,nx
tsoil (i,j,1) = tsoil_ext(i,j,1)
qsoil (i,j,1) = qsoil_ext(i,j,1)
DO k=1,nzsoil
zpsoil(i,j,k) = zpsoil(i,j,k) + zrefsfc
END DO
DO kk=1,nzsoil_ext
zpsoil_ext(i,j,kk) = zpsoil_ext(i,j,kk) + zrefsfc
END DO
END DO
END DO
!----------------------------------------------------------------------
! Initialize soil temp and moisture profiles
!----------------------------------------------------------------------
! Note: All soil depths are negative (downward)
DO k=2,nzsoil
DO j=1,ny
DO i=1,nx
dampdepth = -0.15 - zpsoil(i,j,k) !Typical damping depth = 15 cm
DO kk=2,nzsoil_ext
!----------------------------------------------------------------------
! Linear fit between initialized levels below damping depth (15 cm)
!----------------------------------------------------------------------
IF (zpsoil(i,j,k) >= zpsoil_ext(i,j,kk) .AND. &
zpsoil(i,j,k) <= zpsoil_ext(i,j,kk-1)) THEN
tsoil(i,j,k) = tsoil_ext(i,j,kk)+((tsoil_ext(i,j,kk)- &
tsoil_ext(i,j,kk-1))* (zpsoil(i,j,k)/zpsoil_ext(i,j,kk)) )
qsoil(i,j,k) = qsoil_ext(i,j,kk)+((qsoil_ext(i,j,kk)- &
qsoil_ext(i,j,kk-1))* (zpsoil(i,j,k)/zpsoil_ext(i,j,kk)) )
!----------------------------------------------------------------------
! Exponential fit to initialization above damping depth (15 cm)
!----------------------------------------------------------------------
ELSE IF (zpsoil(i,j,k) > zpsoil_ext(i,j,kk)) THEN
IF (zpsoil_ext(i,j,kk) < dampdepth) THEN
IF (zpsoil(i,j,k) >= dampdepth) THEN
tsoil(i,j,k)=tsoil(i,j,k-1)+(tsoil_ext(i,j,kk)-tsoil(i,j,k-1))* &
EXP( - (zpsoil(i,j,k)/dampdepth) )
qsoil(i,j,k)=qsoil(i,j,k-1)+(qsoil_ext(i,j,kk)-qsoil(i,j,k-1))* &
EXP( - (zpsoil(i,j,k)/dampdepth) )
ELSE IF (zpsoil(i,j,k) < dampdepth) THEN
tsoil(i,j,k) = tsoil(i,j,k-1)+((tsoil_ext(i,j,kk)- &
tsoil(i,j,k-1))* (zpsoil(i,j,k)/zpsoil_ext(i,j,kk)) )
qsoil(i,j,k) = qsoil(i,j,k-1)+((qsoil_ext(i,j,kk)- &
qsoil(i,j,k-1))* (zpsoil(i,j,k)/zpsoil_ext(i,j,kk)) )
END IF
ELSE IF (zpsoil_ext(i,j,kk) >= dampdepth) THEN
tsoil(i,j,k)=tsoil(i,j,k-1)+(tsoil_ext(i,j,kk)-tsoil(i,j,k-1))* &
EXP( - (zpsoil(i,j,k)/zpsoil_ext(i,j,kk)) )
qsoil(i,j,k)=qsoil(i,j,k-1)+(qsoil_ext(i,j,kk)-qsoil(i,j,k-1))* &
EXP( - (zpsoil(i,j,k)/zpsoil_ext(i,j,kk)) )
END IF
!----------------------------------------------------------------------
! Set constant below lowest initialized level
!----------------------------------------------------------------------
ELSE IF (zpsoil(i,j,k) < zpsoil_ext(i,j,kk)) THEN
tsoil(i,j,k) = tsoil_ext(i,j,kk)
qsoil(i,j,k) = qsoil_ext(i,j,kk)
END IF
END DO
END DO
END DO
END DO
RETURN
END SUBROUTINE initsoilext