! FCOPY, ICOPY, etc.: C pre-processor macros for C->N and N->C
! horizontal interpolation. These exist solely to reduce code
! complexity.
! Copy from C array to N array:
#define FCOPY(N,C,i,j,k) \
N(k,i-nits+1)=W1(i,j)*C(II(i,j) ,JJ(i,j) ,k) \
+W2(i,j)*C(II(i,j)+1,JJ(i,j) ,k) \
+W3(i,j)*C(II(i,j)+a,JJ(i,j)-1,k) \
+W4(i,j)*C(II(i,j)+a,JJ(i,j)+1,k)
! Copy from C array to (k,used) indexed N array:
#define uFCOPY(N,C,i,j,k) \
N(k,used)=W1(i,j)*C(II(i,j) ,JJ(i,j) ,k) \
+W2(i,j)*C(II(i,j)+1,JJ(i,j) ,k) \
+W3(i,j)*C(II(i,j)+a,JJ(i,j)-1,k) \
+W4(i,j)*C(II(i,j)+a,JJ(i,j)+1,k)
! ICOPY: Grab from C array, but do not assign:
#define ICOPY(C,i,j,k) \
(W1(i,j)*C(II(i,j) ,JJ(i,j) ,k) \
+ W2(i,j)*C(II(i,j)+1,JJ(i,j) ,k) \
+ W3(i,j)*C(II(i,j)+a,JJ(i,j)-1,k) \
+ W4(i,j)*C(II(i,j)+a,JJ(i,j)+1,k))
! IKJCOPY: grab from IKJ C array, do not assign:
#define IKJCOPY(C,i,k,j) \
(W1(i,j)*C(II(i,j) ,k,JJ(i,j) ) \
+ W2(i,j)*C(II(i,j)+1,k,JJ(i,j) ) \
+ W3(i,j)*C(II(i,j)+a,k,JJ(i,j)-1) \
+ W4(i,j)*C(II(i,j)+a,k,JJ(i,j)+1))
! ICOPY2D: grab from IJ C array, do not assign:
#define ICOPY2D(C,i,j) \
(W1(i,j)*C(II(i,j) ,JJ(i,j) ) \
+ W2(i,j)*C(II(i,j)+1,JJ(i,j) ) \
+ W3(i,j)*C(II(i,j)+a,JJ(i,j)-1) \
+ W4(i,j)*C(II(i,j)+a,JJ(i,j)+1))
! Copying from N array to C array:
#define UPCOPY(C,N,i,j,k,ni,nj) \
C(k,i-istart+1)=( N(ni,nj+2,k) \
+ N(ni-a ,nj+1,k) + N(ni+1-a,nj+1,k) \
+ N(ni-1,nj ,k) + N(ni,nj ,k) + N(ni+1,nj ,k) \
+ N(ni-a ,nj-1,k) + N(ni+1-a,nj-1,k) \
+ N(ni,nj-2,k) \
) / 9
! Average to C points from N points without assignment:
#define NGRAB(N,ni,nj,nk) \
( N(ni,nj+2,nk) \
+ N(ni-a ,nj+1,nk) + N(ni+1-a,nj+1,nk) \
+ N(ni-1,nj ,nk) + N(ni,nj ,nk) + N(ni+1,nj ,nk) \
+ N(ni-a ,nj-1,nk) + N(ni+1-a,nj-1,nk) \
+ N(ni,nj-2,nk) \
) / 9
! Average to C points from N points without assignment on an IKJ grid:
#define NGRABIKJ(N,ni,nk,nj) \
( N(ni,nk,nj+2) \
+ N(ni-a ,nk,nj+1) + N(ni+1-a,nk,nj+1) \
+ N(ni-1,nk,nj ) + N(ni,nk,nj ) + N(ni+1,nk,nj ) \
+ N(ni-a ,nk,nj-1) + N(ni+1-a,nk,nj-1) \
+ N(ni,nk,nj-2) \
) / 9
! Average to C points from N points without assignment, no vertical levels:
#define NGRAB2D(N,ni,nj) \
( N(ni,nj+2) \
+ N(ni-a ,nj+1) + N(ni+1-a,nj+1) \
+ N(ni-1,nj ) + N(ni,nj ) + N(ni+1,nj ) \
+ N(ni-a ,nj-1) + N(ni+1-a,nj-1) \
+ N(ni,nj-2) \
) / 9
! Copying from N array to I array:
#define N2ICOPY(C,N,i,j,k,ni,nj) \
C(i,j,k)=( N(ni,nj+2,k) \
+ N(ni-a ,nj+1,k) + N(ni+1-a,nj+1,k) \
+ N(ni-1,nj ,k) + N(ni,nj ,k) + N(ni+1,nj ,k) \
+ N(ni-a ,nj-1,k) + N(ni+1-a,nj-1,k) \
+ N(ni,nj-2,k) \
) / 9
module module_interp_nmm 22
use module_model_constants
, only: g, R_D, p608
implicit none
private
public :: interp_T_PD_Q
public :: nmm_interp_pd, nmm_keep_pd, nmm_method_linear
public :: c2b_fulldom, c2n_fulldom, n2c_fulldom
public :: c2b_mass, c2n_mass, n2c_mass
public :: c2n_massikj, n2c_massikj
public :: c2b_copy3d, c2n_copy3d, n2c_copy3d
public :: c2b_copy2d, c2n_copy2d, n2c_copy2d
public :: c2b_near2d, c2n_near2d, n2c_near2d
public :: c2b_inear2d, c2n_inear2d, n2c_inear2d
public :: c2n_near3dikj, c2n_sst, c2n_near3d
! unimplemented: public :: nmm_method_quadratic, nmm_method_spline_log, nmm_method_copy, nmm_method_spline
! integer, parameter :: nmm_method_copy = 5
! integer, parameter :: nmm_method_spline_log = 4
! integer, parameter :: nmm_method_spline = 3
! integer, parameter :: nmm_method_quadratic = 2
integer, parameter :: nmm_method_linear = 1
integer, parameter :: nmm_interp_pd = 1
integer, parameter :: nmm_keep_pd = 0
! Constants from the original base_state_parent:
REAL, PARAMETER :: LAPSR=6.5E-3, GI=1./G,D608=0.608
REAL, PARAMETER :: COEF3=287.05*GI*LAPSR, COEF2=-1./COEF3
REAL, PARAMETER :: TRG=2.0*R_D*GI,LAPSI=1.0/LAPSR
REAL, PARAMETER :: P_REF=103000.
integer, parameter :: EConst=0, ECopy=1, EExtrap=2 ! MUST match module_dm
contains
! ********************************************************************
! subs *_NEAR2D -- horizontally interpolate via nearest neighbor
! method, but do not vertically interpolate. Only handles 2D H grid
! arrays.
! ********************************************************************
subroutine c2b_near2d (inear,jnear,cfield, & 1
fbxs,fbxe,fbys,fbye, &
cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte)
implicit none
integer, parameter :: bdyw = 1
integer, intent(in):: &
cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte
integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
real, intent(in) :: cfield(cims:cime,cjms:cjme)
! Output field boundary info:
real,dimension(nims:nime,bdyw) :: fbys,fbye
real,dimension(njms:njme,bdyw) :: fbxs,fbxe
integer :: j,i,a,nx,nz,ck,k,j1,j2,add
real :: weight
nx=min(nide-1,nite)-max(nids,nits)+1
j1=max(njts-1,njds)
if(mod(j1,2)/=1) j1=j1+1
j2=min(njte+1,njde-1)
if(mod(j2,2)/=1) j2=j2-1
if(nits==1) then
i=1
do j=j1,j2,2
fbxs(j,1)=cfield(inear(i,j),jnear(i,j))
enddo
endif
if(nite>=nide-1) then
i=nide-1
do j=j1,j2,2
fbxe(j,1)=cfield(inear(i,j),jnear(i,j))
enddo
endif
if(njts==1) then
j=1
do i=max(nits-1,nids),min(nite+1,nide-1)
fbys(i,1)=cfield(inear(i,j),jnear(i,j))
enddo
endif
if(njte>=njde-1) then
j=njde-1
do i=max(nits-1,nids),min(nite+1,nide-1)
fbye(i,1)=cfield(inear(i,j),jnear(i,j))
enddo
endif
end subroutine c2b_near2d
subroutine n2c_near2d (& 1
cfield,nfield,ipos,jpos, &
cids, cide, cjds, cjde, &
cims, cime, cjms, cjme, &
cits, cite, cjts, cjte, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte)
implicit none
integer, intent(in) :: &
cids, cide, cjds, cjde, &
cims, cime, cjms, cjme, &
cits, cite, cjts, cjte, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte, &
ipos,jpos
real, intent(out) :: cfield(cims:cime,cjms:cjme)
real, intent(in) :: nfield(nims:nime,njms:njme)
integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j
jstart=MAX(jpos+1,cjts)
jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
bigj: do j=jstart,jend
a=mod(j,2)
istart=MAX(ipos+a,cits)
iend=MIN(ipos+(nide-nids)/nri-1,cite)
nj = (j-jpos)*nrj + 1
iloop: do i=istart,iend
ni = (i-ipos)*nri + 2 - a
cfield(i,j)=nfield(ni,nj)
enddo iloop
enddo bigj
end subroutine n2c_near2d
subroutine c2n_near2d (inear,jnear, & 1
cfield,nfield,imask, &
cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte)
implicit none
integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
integer, intent(in):: cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte
integer, intent(in), dimension(nims:nime,njms:njme) :: imask
real, intent(in) :: cfield(cims:cime,cjms:cjme)
real, intent(out) :: nfield(nims:nime,njms:njme)
integer :: j,i,a,nx
nx=min(nide-1,nite)-max(nids,nits)+1
bigj: do j=max(njds,njts),min(njde-1,njte)
interploop: do i=max(nids,nits),min(nide-1,nite)
if(imask(i,j)/=0) cycle interploop
nfield(i,j)=cfield(inear(i,j),jnear(i,j))
enddo interploop
end do bigj
end subroutine c2n_near2d
subroutine c2n_near3d (inear,jnear, & 1
cfield,nfield,imask, &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
implicit none
integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte
integer, intent(in), dimension(nims:nime,njms:njme) :: imask
real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
real, intent(out) :: nfield(nims:nime,njms:njme,nkms:nkme)
integer :: j,i,a,nx,k
nx=min(nide-1,nite)-max(nids,nits)+1
bigk: do k=nkds,nkde
bigj: do j=max(njds,njts),min(njde-1,njte)
interploop: do i=max(nids,nits),min(nide-1,nite)
if(imask(i,j)/=0) cycle interploop
nfield(i,j,k)=cfield(inear(i,j),jnear(i,j),k)
enddo interploop
end do bigj
enddo bigk
end subroutine c2n_near3d
subroutine c2n_near3dikj (inear,jnear, & 1
cfield,nfield,imask, &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
implicit none
integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte
integer, intent(in), dimension(nims:nime,njms:njme) :: imask
real, intent(in) :: cfield(cims:cime,ckms:ckme,cjms:cjme)
real, intent(out) :: nfield(nims:nime,nkms:nkme,njms:njme)
integer :: j,i,a,nx,k
nx=min(nide-1,nite)-max(nids,nits)+1
bigj: do j=max(njds,njts),min(njde-1,njte)
bigk: do k=nkds,nkde
interploop: do i=max(nids,nits),min(nide-1,nite)
if(imask(i,j)/=0) cycle interploop
nfield(i,k,j)=cfield(inear(i,j),k,jnear(i,j))
enddo interploop
enddo bigk
end do bigj
end subroutine c2n_near3dikj
subroutine c2n_sst (inear,jnear, & 1
cfield,nfield, &
cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte)
implicit none
integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
integer, intent(in):: cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte
real, intent(in) :: cfield(cims:cime,cjms:cjme)
real, intent(out) :: nfield(nims:nime,njms:njme)
integer :: j,i,a,nx
nx=min(nide-1,nite)-max(nids,nits)+1
bigj: do j=max(njds,njts),min(njde-1,njte)
interploop: do i=max(nids,nits),min(nide-1,nite)
nfield(i,j)=cfield(inear(i,j),jnear(i,j))
enddo interploop
end do bigj
end subroutine c2n_sst
! ********************************************************************
! subs *_INEAR2D -- horizontally interpolate integers via nearest
! neighbor method, but do not vertically interpolate. Only handles
! 2D H grid integer arrays.
! ********************************************************************
subroutine c2b_inear2d (inear,jnear,cfield, & 1
fbxs,fbxe,fbys,fbye, &
cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte)
implicit none
integer, parameter :: bdyw = 1
integer, intent(in):: &
cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte
integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
integer, intent(in) :: cfield(cims:cime,cjms:cjme)
! Output field boundary info:
integer,dimension(nims:nime,bdyw) :: fbys,fbye
integer,dimension(njms:njme,bdyw) :: fbxs,fbxe
integer :: j,i,a,nx,nz,ck,k,j1,j2,add
nx=min(nide-1,nite)-max(nids,nits)+1
j1=max(njts-1,njds)
if(mod(j1,2)/=1) j1=j1+1
j2=min(njte+1,njde-1)
if(mod(j2,2)/=1) j2=j2-1
if(nits==1) then
i=1
do j=j1,j2,2
fbxs(j,1)=cfield(inear(i,j),jnear(i,j))
enddo
endif
if(nite>=nide-1) then
i=nide-1
do j=j1,j2,2
fbxe(j,1)=cfield(inear(i,j),jnear(i,j))
enddo
endif
if(njts==1) then
j=1
do i=max(nits-1,nids),min(nite+1,nide-1)
fbys(i,1)=cfield(inear(i,j),jnear(i,j))
enddo
endif
if(njte>=njde-1) then
j=njde-1
do i=max(nits-1,nids),min(nite+1,nide-1)
fbye(i,1)=cfield(inear(i,j),jnear(i,j))
enddo
endif
end subroutine c2b_inear2d
subroutine n2c_inear2d (& 1
cfield,nfield,ipos,jpos, &
cids, cide, cjds, cjde, &
cims, cime, cjms, cjme, &
cits, cite, cjts, cjte, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte)
implicit none
integer, intent(in) :: &
cids, cide, cjds, cjde, &
cims, cime, cjms, cjme, &
cits, cite, cjts, cjte, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte, &
ipos,jpos
integer, intent(out) :: cfield(cims:cime,cjms:cjme)
integer, intent(in) :: nfield(nims:nime,njms:njme)
integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j
jstart=MAX(jpos+1,cjts)
jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
bigj: do j=jstart,jend
a=mod(j,2)
istart=MAX(ipos+a,cits)
iend=MIN(ipos+(nide-nids)/nri-1,cite)
nj = (j-jpos)*nrj + 1
iloop: do i=istart,iend
ni = (i-ipos)*nri + 2 - a
cfield(i,j)=nfield(ni,nj)
enddo iloop
enddo bigj
end subroutine n2c_inear2d
subroutine c2n_inear2d (inear,jnear, & 1
cfield,nfield,imask, &
cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte)
implicit none
integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
integer, intent(in):: cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte
integer, intent(in), dimension(nims:nime,njms:njme) :: imask
integer, intent(in) :: cfield(cims:cime,cjms:cjme)
integer, intent(out) :: nfield(nims:nime,njms:njme)
integer :: j,i,a,nx
nx=min(nide-1,nite)-max(nids,nits)+1
bigj: do j=max(njds,njts),min(njde-1,njte)
interploop: do i=max(nids,nits),min(nide-1,nite)
if(imask(i,j)/=0) cycle interploop
nfield(i,j)=cfield(inear(i,j),jnear(i,j))
enddo interploop
end do bigj
end subroutine c2n_inear2d
! ********************************************************************
! subs *_COPY2D -- horizontally interpolate but do not vertically
! interpolate. Only handles 2D arrays (H or V grid).
! ********************************************************************
subroutine c2b_copy2d (II,JJ,W1,W2,W3,W4,cfield, & 3
fbxs,fbxe,fbys,fbye, &
cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte, hgrid)
implicit none
integer, parameter :: bdyw = 1
integer, intent(in):: &
cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte
logical, intent(in) :: hgrid
real, intent(in), dimension(nims:nime,njms:njme) :: &
W1,W2,W3,W4
integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
real, intent(in) :: cfield(cims:cime,cjms:cjme)
! Output field boundary info:
real,dimension(nims:nime,bdyw) :: fbys,fbye
real,dimension(njms:njme,bdyw) :: fbxs,fbxe
integer :: j,i,a,nx,nz,ck,k,j1,j2,add
real :: weight
nx=min(nide-1,nite)-max(nids,nits)+1
if(hgrid) then
a=1
else
a=0
endif
j1=max(njts-1,njds)
if(mod(j1,2)/=a) j1=j1+1
j2=min(njte+1,njde-1)
if(mod(j2,2)/=a) j2=j2-1
if(nits==1) then
i=1
do j=j1,j2,2
if(hgrid) then
a=1-mod(JJ(i,j),2)
else
a=mod(JJ(i,j),2)
endif
fbxs(j,1)=ICOPY2D(cfield,i,j)
enddo
endif
if(nite>=nide-1) then
i=nide-1
do j=j1,j2,2
if(hgrid) then
a=1-mod(JJ(i,j),2)
else
a=mod(JJ(i,j),2)
endif
fbxe(j,1)=ICOPY2D(cfield,i,j)
enddo
endif
if(njts==1) then
j=1
do i=max(nits-1,nids),min(nite+1,nide-1)
if(hgrid) then
a=1-mod(JJ(i,j),2)
else
a=mod(JJ(i,j),2)
endif
fbys(i,1)=ICOPY2D(cfield,i,j)
enddo
endif
if(njte>=njde-1) then
j=njde-1
do i=max(nits-1,nids),min(nite+1,nide-1)
if(hgrid) then
a=1-mod(JJ(i,j),2)
else
a=mod(JJ(i,j),2)
endif
fbye(i,1)=ICOPY2D(cfield,i,j)
enddo
endif
end subroutine c2b_copy2d
subroutine n2c_copy2d (& 4
cfield,nfield,ipos,jpos, &
cids, cide, cjds, cjde, &
cims, cime, cjms, cjme, &
cits, cite, cjts, cjte, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte, &
hgrid)
implicit none
integer, intent(in) :: &
cids, cide, cjds, cjde, &
cims, cime, cjms, cjme, &
cits, cite, cjts, cjte, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte, &
ipos,jpos
real, intent(out) :: cfield(cims:cime,cjms:cjme)
real, intent(in) :: nfield(nims:nime,njms:njme)
logical, intent(in) :: hgrid
integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j
jstart=MAX(jpos+1,cjts)
jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
bigj: do j=jstart,jend
a=mod(j,2)
if(.not.hgrid) then
a=1-a
endif
istart=MAX(ipos+a,cits)
iend=MIN(ipos+(nide-nids)/nri-1,cite)
nj = (j-jpos)*nrj + 1
iloop: do i=istart,iend
ni = (i-ipos)*nri + 2 - a
cfield(i,j)=NGRAB2D(nfield,ni,nj)
enddo iloop
enddo bigj
end subroutine n2c_copy2d
subroutine c2n_copy2d (II,JJ,W1,W2,W3,W4, & 3
cfield,nfield,imask, &
cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte, hgrid)
implicit none
logical, intent(in) :: hgrid
real, intent(in), dimension(nims:nime,njms:njme) :: &
W1,W2,W3,W4
integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
integer, intent(in):: cims, cime, cjms, cjme, &
nids, nide, njds, njde, &
nims, nime, njms, njme, &
nits, nite, njts, njte
integer, intent(in), dimension(nims:nime,njms:njme) :: imask
real, intent(in) :: cfield(cims:cime,cjms:cjme)
real, intent(out) :: nfield(nims:nime,njms:njme)
integer :: j,i,a,nx
nx=min(nide-1,nite)-max(nids,nits)+1
bigj: do j=max(njds,njts),min(njde-1,njte)
interploop: do i=max(nids,nits),min(nide-1,nite)
if(imask(i,j)/=0) cycle interploop
if(hgrid) then
a=1-mod(JJ(i,j),2)
else
a=mod(JJ(i,j),2)
endif
nfield(i,j)=ICOPY2D(cfield,i,j)
enddo interploop
end do bigj
end subroutine c2n_copy2d
! ********************************************************************
! subs *_COPY3D -- horizontally interpolate but do not vertically
! interpolate
! ********************************************************************
subroutine c2b_copy3d (II,JJ,W1,W2,W3,W4,cfield, & 2
fbxs,fbxe,fbys,fbye, &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, hgrid)
implicit none
integer, parameter :: bdyw = 1
logical, intent(in) :: hgrid
integer, intent(in):: &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte
real, intent(in), dimension(nims:nime,njms:njme) :: &
W1,W2,W3,W4
integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
! Output field boundary info:
real,dimension(nims:nime,nkms:nkme,bdyw) :: fbys,fbye
real,dimension(njms:njme,nkms:nkme,bdyw) :: fbxs,fbxe
integer :: j,i,a,nx,nz,ck,k,j1,j2
real :: weight
nx=min(nide-1,nite)-max(nids,nits)+1
nz=nkde-nkds+1
j1=max(njts-1,njds)
if(mod(j1,2)/=0) j1=j1+1
j2=min(njte+1,njde-1)
if(mod(j2,2)/=0) j2=j2-1
if(nits==1) then
i=1
do j=j1,j2,2
if(hgrid) then
a=1-mod(JJ(i,j),2)
else
a=mod(JJ(i,j),2)
endif
kloop1: do k=nkds,min(nkde,nkte)
fbxs(j,k,1)=ICOPY(cfield,i,j,k)
enddo kloop1
enddo
endif
if(nite>=nide-1) then
i=nide-1
do j=j1,j2,2
if(hgrid) then
a=1-mod(JJ(i,j),2)
else
a=mod(JJ(i,j),2)
endif
kloop2: do k=nkds,min(nkde,nkte)
fbxe(j,k,1)=ICOPY(cfield,i,j,k)
enddo kloop2
enddo
endif
if(njts==1) then
j=1
do i=max(nits-1,nids),min(nite+1,nide-1)
if(hgrid) then
a=1-mod(JJ(i,j),2)
else
a=mod(JJ(i,j),2)
endif
kloop3: do k=nkts,min(nkde,nkte)
fbys(i,k,1)=ICOPY(cfield,i,j,k)
enddo kloop3
enddo
endif
if(njte>=njde-1) then
j=njde-1
do i=max(nits-1,nids),min(nite+1,nide-1)
if(hgrid) then
a=1-mod(JJ(i,j),2)
else
a=mod(JJ(i,j),2)
endif
kloop4: do k=nkts,min(nkde,nkte)
fbye(i,k,1)=ICOPY(cfield,i,j,k)
enddo kloop4
enddo
endif
end subroutine c2b_copy3d
subroutine n2c_copy3d (& 2
cfield,nfield,ipos,jpos, &
cids, cide, cjds, cjde, ckds, ckde, &
cims, cime, cjms, cjme, ckms, ckme, &
cits, cite, cjts, cjte, ckts, ckte, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, &
hgrid)
implicit none
logical, intent(in) :: hgrid
integer, intent(in) :: &
cids, cide, cjds, cjde, ckds, ckde, &
cims, cime, cjms, cjme, ckms, ckme, &
cits, cite, cjts, cjte, ckts, ckte, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, &
ipos,jpos
real, intent(out) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
real, intent(in) :: nfield(nims:nime,njms:njme,nkms:nkme)
integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
integer :: nx,nz,jstart,jend,istart,iend,ni,nj,a,i,j,k,nk
real :: weight
nx=min(nide-2,nite)-max(nids+1,nits)+1
nz=nkde-nkds+1
jstart=MAX(jpos+1,cjts)
jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
bigj: do j=jstart,jend
if(hgrid) then
a=mod(j,2)
else
a=1-mod(j,2)
endif
istart=MAX(ipos+a,cits)
iend=MIN(ipos+(nide-nids)/nri-1,cite)
nj = (j-jpos)*nrj + 1
iloop: do i=istart,iend
do k=nkds,nkde
ni = (i-ipos)*nri + 2 - a
cfield(i,j,k)=NGRAB(nfield,ni,nj,k)
enddo
enddo iloop
enddo bigj
end subroutine n2c_copy3d
subroutine c2n_copy3d (II,JJ,W1,W2,W3,W4, & 2
cfield,nfield,imask, &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, &
hgrid)
implicit none
real, intent(in), dimension(nims:nime,njms:njme) :: &
W1,W2,W3,W4
integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
logical, intent(in) :: hgrid
integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte
integer, intent(in), dimension(nims:nime,njms:njme) :: imask
real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
real, intent(out) :: nfield(nims:nime,njms:njme,nkms:nkme)
integer :: j,i,a,nx,nz,k
nx=min(nide-1,nite)-max(nids,nits)+1
nz=nkde-nkds+1
bigj: do j=max(njds,njts),min(njde-1,njte)
interploop: do i=max(nids,nits),min(nide-1,nite)
if(imask(i,j)/=0) cycle interploop
kinterploop: do k=nkds,nkde
if(hgrid) then
a=1-mod(JJ(i,j),2)
else
a=mod(JJ(i,j),2)
endif
nfield(i,j,k) = ICOPY(cfield,i,j,k)
enddo kinterploop
enddo interploop
end do bigj
end subroutine c2n_copy3d
! ********************************************************************
! subs *_MASS -- horizontally and vertically interpolate. Vertical
! interpolation uses the iinfo and winfo arrays (and equivalent
! boundary arrays) generated in the *_fulldom subroutines.
! ********************************************************************
subroutine c2b_mass (II,JJ,W1,W2,W3,W4,cfield, & 1
ibxs,ibxe,ibys,ibye, &
wbxs,wbxe,wbys,wbye, &
fbxs,fbxe,fbys,fbye, &
emethod,evalue, &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
implicit none
integer, parameter :: bdyw = 1
integer, intent(in):: &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte
integer, intent(in) :: emethod ! extrapolation method
real, intent(in) :: evalue
real, intent(in), dimension(nims:nime,njms:njme) :: &
W1,W2,W3,W4
integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
! Output field boundary info:
real,dimension(nims:nime,nkms:nkme,bdyw) :: fbys,fbye
real,dimension(njms:njme,nkms:nkme,bdyw) :: fbxs,fbxe
! Weights and indices:
real,dimension(nims:nime,nkms:nkme,bdyw) :: wbys,wbye
real,dimension(njms:njme,nkms:nkme,bdyw) :: wbxs,wbxe
integer,dimension(nims:nime,nkms:nkme,bdyw) :: ibys,ibye
integer,dimension(njms:njme,nkms:nkme,bdyw) :: ibxs,ibxe
integer :: j,i,a,nx,nz,ck,k,j1,j2
real :: weight
nx=min(nide-1,nite)-max(nids,nits)+1
nz=nkde-nkds+1
j1=max(njts-1,njds)
if(mod(j1,2)/=1) j1=j1+1
j2=min(njte+1,njde-1)
if(mod(j2,2)/=1) j2=j2-1
if(nits==1) then
i=1
do j=j1,j2,2
a=1-mod(JJ(i,j),2)
kloop1: do k=nkds,min(nkde,nkte)
weight=wbxs(j,k,1)
ck=ibxs(j,k,1)
if(ck>1) then
fbxs(j,k,1) = &
ICOPY(cfield,i,j,ck) * weight + &
ICOPY(cfield,i,j,ck-1) * (1.0-weight)
else
! Extrapolate
if(emethod==EConst) then
fbxs(j,k,1)=evalue
else if(emethod==ECopy) then
fbxs(j,k,1)=ICOPY(cfield,i,j,1)
else ! assume 2: linear extrap
fbxs(j,k,1)=evalue * weight + &
ICOPY(cfield,i,j,1) * (1.0-weight)
endif
endif
enddo kloop1
enddo
endif
if(nite>=nide-1) then
i=nide-1
do j=j1,j2,2
a=1-mod(JJ(i,j),2)
kloop2: do k=nkds,min(nkde,nkte)
weight=wbxe(j,k,1)
ck=ibxe(j,k,1)
if(ck>1) then
fbxe(j,k,1) = &
ICOPY(cfield,i,j,ck) * weight + &
ICOPY(cfield,i,j,ck-1) * (1.0-weight)
else
! Extrapolate
if(emethod==EConst) then
fbxe(j,k,1)=evalue
else if(emethod==ECopy) then
fbxe(j,k,1)=ICOPY(cfield,i,j,1)
else ! assume 2: linear extrap
fbxe(j,k,1)=evalue * weight + &
ICOPY(cfield,i,j,1) * (1.0-weight)
endif
endif
enddo kloop2
enddo
endif
if(njts==1) then
j=1
do i=max(nits-1,nids),min(nite+1,nide-1)
a=1-mod(JJ(i,j),2)
kloop3: do k=nkts,min(nkde,nkte)
weight=wbys(i,k,1)
ck=ibys(i,k,1)
if(ck>1) then
fbys(i,k,1) = &
ICOPY(cfield,i,j,ck) * weight + &
ICOPY(cfield,i,j,ck-1) * (1.0-weight)
else
! Extrapolate
if(emethod==EConst) then
fbys(i,k,1)=evalue
else if(emethod==ECopy) then
fbys(i,k,1)=ICOPY(cfield,i,j,1)
else ! assume 2: linear extrap
fbys(i,k,1)=evalue*weight + &
ICOPY(cfield,i,j,1)*(1.0-weight)
endif
endif
enddo kloop3
enddo
endif
if(njte>=njde-1) then
j=njde-1
do i=max(nits-1,nids),min(nite+1,nide-1)
a=1-mod(JJ(i,j),2)
kloop4: do k=nkts,min(nkde,nkte)
weight=wbye(i,k,1)
ck=ibye(i,k,1)
if(ck>1) then
fbye(i,k,1) = &
ICOPY(cfield,i,j,ck) * weight + &
ICOPY(cfield,i,j,ck-1) * (1.0-weight)
else
if(emethod==EConst) then
fbye(i,k,1)=evalue
else if(emethod==ECopy) then
fbye(i,k,1)=ICOPY(cfield,i,j,1)
else ! assume 2: linear extrap
fbye(i,k,1)=evalue*weight + &
ICOPY(cfield,i,j,1)*(1.0-weight)
endif
endif
enddo kloop4
enddo
endif
end subroutine c2b_mass
subroutine n2c_mass (& 1
cfield,nfield,iinfo,winfo,ipos,jpos, &
emethod, evalue, &
cids, cide, cjds, cjde, ckds, ckde, &
cims, cime, cjms, cjme, ckms, ckme, &
cits, cite, cjts, cjte, ckts, ckte, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
implicit none
integer, intent(in) :: &
cids, cide, cjds, cjde, ckds, ckde, &
cims, cime, cjms, cjme, ckms, ckme, &
cits, cite, cjts, cjte, ckts, ckte, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, &
ipos,jpos, emethod
real, intent(in) :: evalue
real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: winfo
integer, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: iinfo
real, intent(out) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
real, intent(in) :: nfield(nims:nime,njms:njme,nkms:nkme)
integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j,k,nk
real :: weight
nz=nkde-nkds+1
jstart=MAX(jpos+1,cjts)
jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
bigj: do j=jstart,jend
a=mod(j,2)
istart=MAX(ipos+a,cits)
iend=MIN(ipos+(nide-nids)/nri-1,cite)
nj = (j-jpos)*nrj + 1
iloop: do i=istart,iend
ni = (i-ipos)*nri + 2 - a
kinterploop: do k=nkds,nkde
weight=winfo(i,j,k)
nk=iinfo(i,j,k)
if(nk>1) then
cfield(i,j,k) = &
NGRAB(nfield,ni,nj,nk) * weight + &
! pjj/cray - source line limit in Cray compiler
NGRAB(nfield,ni,nj,nk-1) &
* (1.0-weight)
else
! Extrapolate
if(emethod==EConst) then
cfield(i,j,k)=evalue
elseif(emethod==ECopy) then
cfield(i,j,k)=NGRAB(nfield,ni,nj,1)
else ! Assume 2: linear extrap
cfield(i,j,k)=evalue*weight + &
NGRAB(nfield,ni,nj,1)*(1.0-weight)
endif
endif
end do kinterploop
enddo iloop
enddo bigj
end subroutine n2c_mass
subroutine n2c_massikj (& 1
cfield,nfield,iinfo,winfo,ipos,jpos, &
emethod, evalue, &
cids, cide, cjds, cjde, ckds, ckde, &
cims, cime, cjms, cjme, ckms, ckme, &
cits, cite, cjts, cjte, ckts, ckte, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
implicit none
integer, intent(in) :: &
cids, cide, cjds, cjde, ckds, ckde, &
cims, cime, cjms, cjme, ckms, ckme, &
cits, cite, cjts, cjte, ckts, ckte, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, &
ipos,jpos, emethod
real, intent(in) :: evalue
real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: winfo
integer, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: iinfo
real, intent(out) :: cfield(cims:cime,ckms:ckme,cjms:cjme)
real, intent(in) :: nfield(nims:nime,nkms:nkme,njms:njme)
integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j,k,nk
real :: weight
nz=nkde-nkds+1
jstart=MAX(jpos+1,cjts)
jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
bigj: do j=jstart,jend
a=mod(j,2)
istart=MAX(ipos+a,cits)
iend=MIN(ipos+(nide-nids)/nri-1,cite)
nj = (j-jpos)*nrj + 1
kinterploop: do k=nkds,nkde
iloop: do i=istart,iend
ni = (i-ipos)*nri + 2 - a
weight=winfo(i,j,k)
nk=iinfo(i,j,k)
if(nk>1) then
cfield(i,k,j) = &
NGRABIKJ(nfield,ni,nk,nj) * weight + &
! pjj/cray - source line limit in Cray compiler
NGRABIKJ(nfield,ni,nk-1,nj) &
* (1.0-weight)
else
! Extrapolate
if(emethod==EConst) then
cfield(i,k,j)=evalue
elseif(emethod==ECopy) then
cfield(i,k,j)=NGRABIKJ(nfield,ni,1,nj)
else ! Assume 2: linear extrap
cfield(i,k,j)=evalue*weight + &
NGRABIKJ(nfield,ni,1,nj)*(1.0-weight)
endif
endif
enddo iloop
end do kinterploop
enddo bigj
end subroutine n2c_massikj
subroutine c2n_mass (II,JJ,W1,W2,W3,W4, & 1
cfield,nfield,iinfo,winfo,imask, &
emethod, evalue, &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
implicit none
real, intent(in), dimension(nims:nime,njms:njme) :: &
W1,W2,W3,W4
integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, emethod
real, intent(in) :: evalue
real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: winfo
integer, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: iinfo
integer, intent(in), dimension(nims:nime,njms:njme) :: imask
real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
real, intent(out) :: nfield(nims:nime,njms:njme,nkms:nkme)
integer :: j,i,a,nx,nz,ck,k
real :: weight
nx=min(nide-1,nite)-max(nids,nits)+1
nz=nkde-nkds+1
bigj: do j=max(njds,njts),min(njde-1,njte)
interploop: do i=max(nids,nits),min(nide-1,nite)
if(imask(i,j)/=0) cycle interploop
kinterploop: do k=nkds,nkde
a=1-mod(JJ(i,j),2)
weight=winfo(i,j,k)
ck=iinfo(i,j,k)
if(ck>1) then
nfield(i,j,k) = &
ICOPY(cfield,i,j,ck) * weight + &
ICOPY(cfield,i,j,ck-1) * (1.0-weight)
else
! Extrapolate
if(emethod==EConst) then
nfield(i,j,k)=evalue
elseif(emethod==ECopy) then
nfield(i,j,k)=ICOPY(cfield,i,j,1)
else ! Assume 2: linear extrap
nfield(i,j,k)=evalue*weight + &
ICOPY(cfield,i,j,1)*(1.0-weight)
endif
endif
enddo kinterploop
enddo interploop
end do bigj
end subroutine c2n_mass
subroutine c2n_massikj (II,JJ,W1,W2,W3,W4, & 1
cfield,nfield,iinfo,winfo,imask, &
emethod, evalue, &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
implicit none
real, intent(in), dimension(nims:nime,njms:njme) :: &
W1,W2,W3,W4
integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, emethod
real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: winfo
integer, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: iinfo
integer, intent(in), dimension(nims:nime,njms:njme) :: imask
real, intent(in) :: cfield(cims:cime,ckms:ckme,cjms:cjme)
real, intent(out) :: nfield(nims:nime,nkms:nkme,njms:njme)
real, intent(In) :: evalue
integer :: j,i,a,nx,nz,ck,k
real :: weight
nx=min(nide-1,nite)-max(nids,nits)+1
nz=nkde-nkds+1
bigj: do j=max(njds,njts),min(njde-1,njte)
interploop: do i=max(nids,nits),min(nide-1,nite)
if(imask(i,j)/=0) cycle interploop
kinterploop: do k=nkds,nkde
a=1-mod(JJ(i,j),2)
weight=winfo(i,j,k)
ck=iinfo(i,j,k)
if(ck>1) then
nfield(i,k,j) = &
IKJCOPY(cfield,i,ck,j) * weight + &
IKJCOPY(cfield,i,ck-1,j) * (1.0-weight)
else
! Extrapolate
if(emethod==EConst) then
nfield(i,k,j)=evalue
elseif(emethod==ECopy) then
nfield(i,k,j)=IKJCOPY(cfield,i,1,j)
else ! Assume 2: linear extrapolation
nfield(i,k,j)=evalue * weight + &
IKJCOPY(cfield,i,1,j) * (1.0-weight)
endif
endif
enddo kinterploop
enddo interploop
end do bigj
end subroutine c2n_massikj
! ********************************************************************
! subs *_FULLDOM -- recalculates PD and PINT based on FIS if
! requested. Also, generates vertical interpolation arrays for use
! in later interpolations and interpolates T and Q while doing so.
! ********************************************************************
subroutine n2c_fulldom ( & 1,1
deta1,deta2, eta1,eta2, ptop,pdtop, &
cfis,cpint,ct,cpd,cq, &
cids, cide, cjds, cjde, ckds, ckde, &
cims, cime, cjms, cjme, ckms, ckme, &
cits, cite, cjts, cjte, ckts, ckte, &
nfis,npint,nt,npd,nq, &
ipos,jpos, &
! nri,nrj &
out_iinfo,out_winfo, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
implicit none
integer, intent(in) :: ipos,jpos
! integer, intent(in) :: nri,nrj
integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
integer, intent(in):: &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte
integer, intent(in):: &
cids, cide, cjds, cjde, ckds, ckde, &
cims, cime, cjms, cjme, ckms, ckme, &
cits, cite, cjts, cjte, ckts, ckte
real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
real, intent(inout), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_winfo
integer, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_iinfo
real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
real, intent(in), dimension(nims:nime,njms:njme,1) :: nPD
real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
real, intent(in) :: ptop,pdtop
real, dimension(1,nite-nits+1) :: inFIS,inPD,icFIS,icPD
real, dimension(nkde-1,nite-nits+1) :: inT0,inQ,icT,icQ, qinfo,winfo
real, dimension(nkde,nite-nits+1) :: inPINT,icPINT
integer, dimension(nkde-1,nite-nits+1) :: iinfo
integer :: nx,nz,k,i,a,j, istart,iend,jstart,jend, ni,nj,jprint,itest,jtest
character*255 :: message
logical bad
nx=min(cide-2,cite)-max(cids+1,cits)+1
nz=ckde-ckds+1
jstart=MAX(jpos+1,cjts)
jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
bigj: do j=jstart,jend
nj = (j-jpos)*nrj + 1
a=mod(j,2)
istart=MAX(ipos+a,cits)
iend=MIN(ipos+(nide-nids)/nri-1,cite)
nx=iend-istart+1
! STEP 1: Copy coarse and fine nest data into
! temporary arrays, reordering dimensions:
qtloop: do k=ckts,ckte-1
do i=istart,iend
ni = (i-ipos)*nri + 2 - a
UPCOPY(inT0,nT,i,j,k,ni,nj)
UPCOPY(inQ,nQ,i,j,k,ni,nj)
UPCOPY(inPINT,nPINT,i,j,k,ni,nj)
enddo
enddo qtloop
k=nkte
loop2d: do i=istart,iend
ni = (i-ipos)*nri + 2 - a
UPCOPY(inPINT,nPINT,i,j,k,ni,nj)
UPCOPY(inPD,nPD,i,j,1,ni,nj)
UPCOPY(inFIS,nFIS,i,j,1,ni,nj)
icPD(1,i-istart+1)=cPD(i,j,1)
! icPD(1,i-istart+1)=use_this_pd(i,j,1)
icFIS(1,i-istart+1)=cFIS(i,j,1)
enddo loop2d
! Step 2: Interpolate coarse grid to fine grid in reordered
! arrays:
call interp_T_PD_Q(nmm_method_linear, nmm_keep_pd, nx,nz, &
deta1,deta2,eta1,eta2,ptop,pdtop, &
inFIS,icFIS, inPINT,icPINT, inT0, icT, inPD,icPD, inQ,icQ, &
iinfo, winfo)
! Step 3: Copy back from reordered arrays to final nest arrays:
qtloop2: do k=ckts,max(ckte-1,ckde-1)
do i=istart,iend
cT(i,j,k)=icT(k,i-istart+1)
cQ(i,j,k)=icQ(k,i-istart+1)
enddo
enddo qtloop2
izloop: do k=ckts,max(ckte-1,ckde-1)
ixloop: do i=istart,iend
out_iinfo(i,j,k)=iinfo(k,i-istart+1)
out_winfo(i,j,k)=winfo(k,i-istart+1)
enddo ixloop
enddo izloop
iendloop: do i=istart,iend
out_iinfo(i,j,nkde)=nkde
out_winfo(i,j,nkde)=1.0
end do iendloop
k=nkte+1
loop2d2: do i=istart,iend
cPD(i,j,1)=icPD(1,i-istart+1)
enddo loop2d2
end do bigj
end subroutine n2c_fulldom
subroutine c2b_fulldom (II,JJ,W1,W2,W3,W4,& 1,2
deta1,deta2, eta1,eta2, ptop,pdtop, &
cfis,cpint,ct,cpd,cq, nfis, &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, &
ibxs, ibxe, ibys, ibye, &
wbxs, wbxe, wbys, wbye, &
pdbxs, pdbxe, pdbys, pdbye, &
tbxs, tbxe, tbys, tbye, &
qbxs, qbxe, qbys, qbye)
implicit none
integer, intent(in):: &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte
integer, parameter :: bdyw=1
! Domain info:
real, intent(in), dimension(nims:nime,njms:njme) :: W1,W2,W3,W4
integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
real, intent(in) :: ptop,pdtop
! Parent fields:
real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
! Nest terrain info:
real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
! T, Q, PINT, PD boundary info:
real,dimension(nims:nime,nkms:nkme,bdyw) :: tbys,tbye,qbys,qbye
real,dimension(njms:njme,nkms:nkme,bdyw) :: tbxs,tbxe,qbxs,qbxe
real,dimension(nims:nime,1,bdyw) :: pdbys,pdbye
real,dimension(njms:njme,1,bdyw) :: pdbxs,pdbxe
! Weights and indices:
real,dimension(nims:nime,nkms:nkme,bdyw) :: wbys,wbye
real,dimension(njms:njme,nkms:nkme,bdyw) :: wbxs,wbxe
integer,dimension(nims:nime,nkms:nkme,bdyw) :: ibys,ibye
integer,dimension(njms:njme,nkms:nkme,bdyw) :: ibxs,ibxe
integer :: i,j,k,a,b,nx,nz,used,j1,j2,used1
real, dimension(1,2*(nite-nits+5)+2*(njte-njts+5)) :: inFIS,inPD,icFIS,icPD
real, dimension(nkde-1,2*(nite-nits+5)+2*(njte-njts+5)) :: inT,inQ,icT,icQ,winfo
integer, dimension(nkde-1,2*(nite-nits+5)+2*(njte-njts+5)) :: iinfo
real, dimension(nkde,2*(nite-nits+5)+2*(njte-njts+5)) :: inPINT,icPINT
nx=min(nide-1,nite)-max(nids,nits)+1
nz=nkde-nkds+1
j1=max(njts-1,njds)
if(mod(j1,2)/=1) j1=j1+1
j2=min(njte+1,njde-1)
if(mod(j2,2)/=1) j2=j2-1
used=0
bdyloop: do b=1,4
if_xbdy: if(b==1 .or. b==2) then
if(b==1) then
if(nits/=1) cycle bdyloop
i=1
endif
if(b==2) then
if(nite<nide-1) cycle bdyloop
i=nide-1
endif
do j=j1,j2,2
used=used+1
do k=nkts,nkte-1
a=1-mod(JJ(i,j),2)
uFCOPY(icT,cT,i,j,k)
uFCOPY(icQ,cQ,i,j,k)
uFCOPY(icPINT,cPINT,i,j,k)
enddo
k=nkte
a=1-mod(JJ(i,j),2)
uFCOPY(icPINT,cPINT,i,j,k)
uFCOPY(icPD,cPD,i,j,1)
uFCOPY(icFIS,cFIS,i,j,1)
inFIS(1,used)=nFIS(i,j,1)
enddo
endif if_xbdy
if_ybdy: if(b==3 .or. b==4) then
if(b==3) then
if(njts/=1) cycle bdyloop
j=1
endif
if(b==4) then
if(njte<njde-1) cycle bdyloop
j=njde-1
endif
do i=max(nits-1,nids),min(nite+1,nide-1)
used=used+1
do k=nkts,nkte-1
a=1-mod(JJ(i,j),2)
uFCOPY(icT,cT,i,j,k)
uFCOPY(icQ,cQ,i,j,k)
uFCOPY(icPINT,cPINT,i,j,k)
enddo
k=nkte
a=1-mod(JJ(i,j),2)
uFCOPY(icPINT,cPINT,i,j,k)
uFCOPY(icPD,cPD,i,j,1)
uFCOPY(icFIS,cFIS,i,j,1)
inFIS(1,used)=nFIS(i,j,1)
enddo
endif if_ybdy
enddo bdyloop
if(used==0) then
! No boundary points on this tile so we are done
return
endif
call interp_T_PD_Q(nmm_method_linear, nmm_interp_pd, used,nz, &
deta1,deta2,eta1,eta2,ptop,pdtop, &
icFIS,inFIS, icPINT,inPINT, icT, inT, icPD,inPD, icQ,inQ, &
iinfo, winfo)
used1=used
used=0
if_bxs: if(nits==1) then
i=1
do j=j1,j2,2
used=used+1
do k=nkts,nkte-1
tbxs(j,k,1)=inT(k,used)
qbxs(j,k,1)=inQ(k,used)
ibxs(j,k,1)=iinfo(k,used)
wbxs(j,k,1)=winfo(k,used)
enddo
ibxs(j,nkde,1)=nkde
wbxs(j,nkde,1)=1.0
pdbxs(j,1,1)=inPD(1,used)
enddo
endif if_bxs
if_bxe: if(nite>=nide-1) then
i=nide-1
do j=j1,j2,2
used=used+1
do k=nkts,nkte-1
tbxe(j,k,1)=inT(k,used)
qbxe(j,k,1)=inQ(k,used)
ibxe(j,k,1)=iinfo(k,used)
wbxe(j,k,1)=winfo(k,used)
enddo
ibxe(j,nkde,1)=nkde
wbxe(j,nkde,1)=1.0
pdbxe(j,1,1)=inPD(1,used)
enddo
endif if_bxe
if_bys: if(njts==1) then
j=1
do i=max(nits-1,nids),min(nite+1,nide-1)
used=used+1
do k=nkts,nkte-1
tbys(i,k,1)=inT(k,used)
qbys(i,k,1)=inQ(k,used)
ibys(i,k,1)=iinfo(k,used)
wbys(i,k,1)=winfo(k,used)
enddo
ibys(i,nkde,1)=nkde
wbys(i,nkde,1)=1.0
pdbys(i,1,1)=inPD(1,used)
enddo
endif if_bys
if_bye: if(njte>=njde-1) then
j=njde-1
do i=max(nits-1,nids),min(nite+1,nide-1)
used=used+1
do k=nkts,nkte-1
tbye(i,k,1)=inT(k,used)
qbye(i,k,1)=inQ(k,used)
ibye(i,k,1)=iinfo(k,used)
wbye(i,k,1)=winfo(k,used)
enddo
ibye(i,nkde,1)=nkde
wbye(i,nkde,1)=1.0
pdbye(i,1,1)=inPD(1,used)
enddo
endif if_bye
if(used/=used1) then
call wrf_error_fatal('Number of input and output points does not match.')
endif
end subroutine c2b_fulldom
subroutine c2n_fulldom (II,JJ,W1,W2,W3,W4,& 1,2
deta1,deta2, eta1,eta2, ptop,pdtop, &
cfis,cpint,ct,cpd,cq, &
cims, cime, cjms, cjme, ckms, ckme, &
nfis,npint,nt,npd,nq, &
out_iinfo,out_winfo,imask,&
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
implicit none
integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte
real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_winfo
integer, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_iinfo
real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
real, intent(inout), dimension(nims:nime,njms:njme,1) :: nPD
integer, intent(in), dimension(nims:nime,njms:njme) :: imask
real, intent(in), dimension(nims:nime,njms:njme) :: &
W1,W2,W3,W4
integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
real, intent(in) :: ptop,pdtop
real, dimension(1,nite-nits+1) :: inFIS,inPD,icFIS,icPD
real, dimension(nkde-1,nite-nits+1) :: inT,inQ,icT,icQ,winfo
integer, dimension(nkde-1,nite-nits+1) :: iinfo
real, dimension(nkde,nite-nits+1) :: inPINT,icPINT
real :: pdcheck
integer :: i,j,k,a, nx,nz,used
logical :: badbad
nx=min(nide-1,nite)-max(nids,nits)+1
nz=nkde-nkds+1
bigj: do j=max(njds,njts),min(njde-1,njte)
! STEP 1: Copy coarse and fine nest data into
! temporary arrays, reordering dimensions:
used=0
iloop1: do i=max(nids,nits),min(nide-1,nite)
if(imask(i,j)/=0) cycle iloop1
used=used+1
qtloop: do k=nkts,nkte-1
a=1-mod(JJ(i,j),2)
uFCOPY(icT,cT,i,j,k)
uFCOPY(icQ,cQ,i,j,k)
uFCOPY(icPINT,cPINT,i,j,k)
enddo qtloop
k=nkte
a=1-mod(JJ(i,j),2)
uFCOPY(icPINT,cPINT,i,j,k)
uFCOPY(icPD,cPD,i,j,1)
uFCOPY(icFIS,cFIS,i,j,1)
! nPD(i,j,1)=use_this_pd(i,j,1) ! FIXME: remove this line
! if(nPD(i,j,1)> 1.3e5 .or. nPD(i,j,1)<2e4) then
! 3131 format('Invalid input nest PD in C2N interpolation: PD(',I0,',',I0,') = ',F0.7)
! write(0,3131) i,j,nPD(i,j,1)
! call wrf_error_fatal('Invalid input nest PD')
! endif
inPD(1,used)=nPD(i,j,1)
inFIS(1,used)=nFIS(i,j,1)
enddo iloop1
! Step 2: Interpolate coarse grid to fine grid in reordered
! arrays:
if(used==0) then
! No points in this row require interpolation, so we're done
! with this row:
cycle bigj
else
call interp_T_PD_Q(nmm_method_linear, nmm_interp_pd, used,nz, &
deta1,deta2,eta1,eta2,ptop,pdtop, &
icFIS,inFIS, icPINT,inPINT, icT, inT, icPD,inPD, icQ,inQ, &
iinfo, winfo)
endif
! Step 3: Copy back from reordered arrays to final nest arrays:
used=0
iloop2: do i=max(nids,nits),min(nide-1,nite)
if(imask(i,j)/=0) cycle iloop2
used=used+1
qtloop2: do k=nkts,nkte-1
nT(i,j,k)=inT(k,used)
nQ(i,j,k)=inQ(k,used)
nQ(i,j,k)=inQ(k,used)
enddo qtloop2
izloop: do k=nkts,nkte-1
out_iinfo(i,j,k)=iinfo(k,used)
out_winfo(i,j,k)=winfo(k,used)
enddo izloop
out_iinfo(i,j,nkde)=nkde
out_winfo(i,j,nkde)=1.0
k=nkte+1
a=1-mod(JJ(i,j),2)
nPD(i,j,1)=inPD(1,used)
pdcheck=npd(i,j,1)+ptop+pdtop
if(pdcheck<50000.0 .or. pdcheck>105000.0) then
3131 format('Invalid output nest PD in C2N interpolation: PD(',I0,',',I0,') = ',F0.7,' which is ',F0.7,' Pa')
write(0,3131) i,j,nPD(i,j,1),pdcheck
call wrf_error_fatal('Invalid output nest PD')
endif
enddo iloop2
enddo bigj
end subroutine c2n_fulldom
! ********************************************************************
! INTERP_T_PD_Q -- this is the actual vertical interpolation
! subroutine, the implementation of the *_fulldom subroutines. It
! takes arrays of atmospheric columns, with no knowledge of the
! horizontal layout of those columns. This routine will recalculate
! PD if requested, then generate interpolation information and
! interpolate T and Q, within each atmospheric column.
!
! Other variables are handled by the various other c2n, n2c and c2b
! routines in this module. PD, T and Q are handled specially since
! they are critical to the stability of the hydrostatic solver, and
! hence have more expensive extrapolation and mass rebalancing
! methods.
!
! Since this routine has no knowledge of the horizontal layout of
! the atmospheric columns it is given, the columns do not have to
! lie on an I- or J-row. The boundary interpolation routines take
! advantage of that and provide all boundary information to this
! routine in a single call. Of course, the caller must handle all
! horizontal interpolation.
! ********************************************************************
subroutine interp_T_PD_Q(method, pd_interp, nx, nz, & 3,4
deta1,deta2, eta1,eta2, ptop,pdtop, &
fisA,fisB, pintA,pintB, tA,tB, pdA,pdB, qA,qB, &
iinfo, winfo)
implicit none
integer, intent(in) :: pd_interp,method
! real, intent(in) :: dtA,dtB
! Coordinate system definitions must be the same for all domains:
real, intent(in) :: deta1(nz),deta2(nz),eta1(nz),eta2(nz),ptop,pdtop
integer, intent(in) :: nx,nz
! Surface height and mass field information for source (A):
real, intent(in) :: fisA(nx),tA(nz-1,nx),pdA(nx),qA(nz-1,nx),pintA(nz,nx)
! Surface height and mass field information for target (B):
real, intent(inout) :: fisB(nx),tB(nz-1,nx),pdB(nx),qB(nz-1,nx),pintB(nz,nx)
! Interpolation or extrapolation information for use in other
! calls later on:
real, intent(out) :: winfo(nz-1,nx)
integer, intent(out) :: iinfo(nz-1,nx)
! ==================== Local variables ====================
character*255 :: message
integer :: ix,iz,izA,izB,xpr
real :: zA,zB,znext,apelp,rtopp,dz,weight,A,B,pstd1,z,pstd2,pstd12
real :: tsfc(nx), slp(nx), zmslp(nx), z0mid(nx), zbelow(nx), &
tbelow(nx), RHbelow(nx)
real :: pb,pb1,pb2,pa,pa1,pa2,pnext,pa3, wnum,wdenom, QC, P0
! Constants from UPP params.f for RH calculation (all mks):
real, parameter :: PQ0=379.90516
real, parameter :: A2=17.2693882
real, parameter :: A3=273.16
real, parameter :: A4=35.86
real, parameter :: RHmin=1.0E-6 ! minimal RH bound
if(method/=nmm_method_linear) then
call wrf_error_fatal('only linear interpolation is supported')
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Step 1: calculate near-surface values !!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
pstd1=p_ref ! pstd(1) from base_state_parent
! pstd(2) from base_state_parent:
pstd2=eta1(2)*pdtop + eta2(2)*(p_ref-pdtop-ptop) + ptop
pstd12=exp((alog(pstd1)+alog(pstd2))*0.5)
do ix=1,nx
! These calculations are from base_state_parent:
APELP = (pintA(2,ix)+pintA(1,ix))
RTOPP = TRG*tA(1,ix)*(1.0+qA(1,ix)*P608)/APELP
DZ = RTOPP*(DETA1(1)*PDTOP+DETA2(1)*pdA(ix))
Z0MID(ix) = fisA(ix)/g + dz/2
zA=fisA(ix)/g
TSFC(ix) = TA(1,ix)*(1.+D608*QA(1,ix)) &
+ LAPSR*(zA + zA+dz)*0.5
A = LAPSR*zA/TSFC(ix)
SLP(ix) = PINTA(1,ix)*(1-A)**COEF2 ! sea level pressure
B = (pstd1/SLP(ix))**COEF3
ZMSLP(ix)= TSFC(ix)*LAPSI*(1.0 - B) ! Height at 1030. mb level
TBELOW(ix) = TA(1,ix) + LAPSR*(Z0MID(ix)-ZMSLP(ix))
! Calculate lowest level RH. This calculation is from UPP
! CALRH.f:
P0=pdA(ix)+ptop+pdtop ! Use hydrostatic lowest model level P
QC=PQ0/P0*EXP(A2*(TA(1,ix)-A3)/(TA(1,ix)-A4))
RHbelow(ix)=max(RHmin,min(1.0,QA(1,ix)/QC))
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Step 2: figure out the new surface pressures !!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if_pd_interp: if(pd_interp==nmm_keep_pd) then
! PD interpolation is turned off, so we use the original PD.
elseif(pd_interp==nmm_interp_pd) then
! PD interpolation is requested. As with the old base_state_parent,
! determine PD by interpolating or extrapolating the non-hydrostatic
! pressure field (PINT) in source grid to the surface height of the
! target grid.
xloop: do ix=1,nx
if(pintA(1,ix)>p_ref) then
! Cannot extrapolate PD below ground because pressure is
! unrealistically high. Follow base_state_parent method:
! when this happens, assign input pressure to output
! pressure.
WRITE(0,*)'WARNING: NESTED DOMAIN PRESSURE AT LOWEST LEVEL HIGHER THAN PSTD'
WRITE(0,*)'PINT(1),PD(1),PSTD(1)',pintA(1,ix),pdA(ix),p_ref
pdB(ix)=pdA(ix)
cycle xloop
endif
zA=fisA(ix)/g
zB=fisB(ix)/g
if(zB<zA) then
! Need to extrapolate PD below ground.
! Follow base_state_parent method:
! For pressures between pint(1) and 1030 mbar, interpolate
iz=1
weight=abs((zB-zmslp(ix))/(zA-zmslp(ix)))
pdB(ix) = exp(weight*log(pdA(ix)+pdtop+ptop) + &
(1.0-weight)*log(P_REF)) &
- pdtop - ptop
! write(0,*) 'WARNING: extrapolating below ground at ix=',ix
!1302 format('Extrap at ix=',I0,' zA=',F0.3,' zB=',F0.3,' zmslp=',F0.3,' weight=',F0.3,' pdB+pdtop+ptop=',F0.3,' pintA=',F0.3,' pdA(ix)+pdtop+ptop=',F0.3,' P_REF=',F0.3)
! write(0,1302) ix,zA,zB,zmslp(ix),weight,pdB(ix)+pdtop+ptop,pintA(1,ix),pdA(ix)+pdtop+ptop,P_REF
cycle xloop
endif
! PD in grid B lies above ground in grid A. Find out where by
! walking up from the ground until we find the two levels the
! pressure lies between.
znext=-9e9
z=zA
do iz=1,nz-1
! Calculations from base_state_parent:
APELP = (pintA(iz+1,ix)+pintA(iz,ix))
RTOPP = TRG*tA(iz,ix)*(1.0+qA(iz,ix)*P608)/APELP
DZ = RTOPP*(DETA1(iz)*PDTOP+DETA2(iz)*pdA(ix))
znext=z+dz
if(znext>zB) then
! Interpolate between these two levels
weight=(zB-z)/dz
pdB(ix) = exp(weight*log(pintA(iz+1,ix)) + &
(1.0-weight)*log(pintA(iz,ix))) &
- pdtop - ptop
! if(iz>1) then
! write(0,*) 'WARNING: interpolate above ground at ix=',ix
!1303 format('Interp at ix=',I0,' with iz=',I0,' zA=',F0.3,' zB=',F0.3,' zmslp=',F0.3,' weight=',F0.3,' pdB+pdtop+ptop=',F0.3,' pintA=',F0.3,' pdA(ix)+pdtop+ptop=',F0.3,' P_REF=',F0.3)
! write(0,1303) ix,iz,zA,zB,zmslp(ix),weight,pdB(ix)+pdtop+ptop,pintA(1,ix),pdA(ix)+pdtop+ptop,P_REF
! endif
cycle xloop
else
z=znext
endif
enddo
201 format('Target domain surface height ',F0.4,'m is higher than the model top ',F0.4,'m of the source domain. ABORTING')
write(message,201) zB,znext
call wrf_error_fatal(message)
enddo xloop
else
202 format('Invalid value ',I0,' for pd_interp in interp_T_PD_Q')
write(message,202) pd_interp
call wrf_error_fatal(message)
endif if_pd_interp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Step 3: interpolate T and Q in pressure space !!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
outer: do ix=1,nx
! For constant pressure levels, do a straight level-by-level copy:
iz=nz-1
copyloop: do while(eta2(iz) < 1e-5 .and. eta2(iz-1) < 1e-5)
tB(iz,ix)=tA(iz,ix)
qB(iz,ix)=qA(iz,ix)
! Save interpolation information
iinfo(iz,ix)=iz
winfo(iz,ix)=1.0
iz=iz-1
enddo copyloop
! For sigma levels where target domain lies within source,
! interpolate from top down, stopping when we go below ground
! in the source domain.
izA=iz
izB=iz
! FIXME: REMOVE THIS CHECK
if(iz>nz) then
! Make sure the entire vertical extent isn't sigma levels:
call wrf_error_fatal('ERROR: WRF-NMM does not support pure sigma levels (only sigma-pressure hybrid)')
endif
pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
pB1=log(eta1(izB+1)*pdtop + eta2(izB+1)*pdB(ix) + ptop)
pB=(pB2+pB1)*0.5
pA2=log(eta1(izA)*pdtop + eta2(izA)*pdA(ix) + ptop)
pA1=log(eta1(izA+1)*pdtop + eta2(izA+1)*pdA(ix) + ptop)
pA=(pA2+pA1)*0.5
! Find the lowest mass level izA that is above pB
interpinit: do while(izA>1)
pA3=log(eta1(izA-1)*pdtop + eta2(izA-1)*pdA(ix) + ptop)
pnext=(pA2+pA3)*0.5
wdenom=pnext-pA
wnum=pnext-pb
if(pA<=pB .and. wnum>1e-5) then
exit interpinit
else
pA1=pA2
pA2=pa3
izA=izA-1
pA=pnext
endif
enddo interpinit
! Loop over all remaining B points, interpolating to them.
interploop: do while(izB>0 .and. izA>1)
! Find the source domain levels that this target domain level lies between:
zinterp: do while(izA>1)
if(pnext>=pB) then
! We found the two levels, so interpolate and move on to
! the next target level:
weight=max(0.,wnum/wdenom)
tB(izB,ix)=weight*tA(izA,ix) + (1.0-weight)*tA(izA-1,ix)
qB(izB,ix)=weight*qA(izA,ix) + (1.0-weight)*qA(izA-1,ix)
! Save interpolation info
iinfo(izB,ix)=izA ! upper level
winfo(izB,ix)=weight ! linear interpolation weight
! We interpolated to a B point, so go to the next one:
pB1=pB2
izB=izB-1
if(izB>0) then
pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
pB=(pb1+pb2)*0.5
else
exit interploop
endif
wnum=pnext-pb
exit zinterp
else
izA=izA-1
pA1=pA2
pA2=pa3
pA=pnext
if(izA>1) then
pA3=log(eta1(izA-1)*pdtop + eta2(izA-1)*pdA(ix) + ptop)
pnext=(pA2+pA3)*0.5
wdenom=pnext-pa
wnum=pnext-pb
else
exit interploop
endif
endif
enddo zinterp
enddo interploop
! Follow what base_state_parent did for temperature:
! interpolate between the second to last level and P_REF using
! a lapse rate atmosphere. For Q, we use constant RH below
! ground, as suggested by Greg Thompson.
extraploop: do while(izB>=1)
! Decide extrapolation weight:
weight=(pB-pA)/(pstd12-pA)
! Extrapolate Q by copying lowest sigma layer value:
!qB(izB,ix)=qA(1,ix)
! Extrapolate Q by linearly interpolating between 0 and surface value:
!qB(izB,ix)=(1.0-weight)*qA(1,ix)
! Extrapolate T using a lapse rate atmosphere
tB(izB,ix)=weight*tbelow(ix) + (1.0-weight)*tA(1,ix)
! Extrapolate Q using constant RH below ground, as suggested
! by Greg Thompson. This is the inversion of the RH
! calculation used earlier in this function:
P0=eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop
QC=PQ0/P0*EXP(A2*(TB(izB,ix)-A3)/(TB(izB,ix)-A4))
QB(izB,ix)=QC*RHbelow(ix)
! Save extrapolation information
iinfo(izB,ix)=0
winfo(izB,ix)=weight
! Move to the next B level
izB=izB-1
if(izB>0) then
pB1=pB2
pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
pB=(pb1+pb2)*0.5
endif
enddo extraploop
enddo outer
end subroutine interp_T_PD_Q
end module module_interp_nmm
! ********************************************************************
! subs EXT_*_FULLDOM -- simple wrappers around module_interp_nmm's
! *_FULLDOM routines. These exist solely to allow module_dm to
! interface with this module -- this is needed because the WRF build
! scripts compile module_dm before module_interp_nmm.
! ********************************************************************
subroutine ext_c2n_fulldom (II,JJ,W1,W2,W3,W4,&,3
deta1,deta2, eta1,eta2, ptop,pdtop, &
cpint,ct,cpd,cq, &
cims, cime, cjms, cjme, ckms, ckme, &
npint,nt,npd,nq, &
out_iinfo,out_winfo,imask,&
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
! This subroutine is an alias for c2n_fulldom. It exists to allow
! the routine to be called from module_dm
use module_interp_store, only: parent_fis, nest_fis
use module_interp_nmm, only: c2n_fulldom
implicit none
integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte
real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD
real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_winfo
integer, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_iinfo
real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
real, intent(inout), dimension(nims:nime,njms:njme,1) :: nPD
integer, intent(in), dimension(nims:nime,njms:njme) :: imask
real, intent(in), dimension(nims:nime,njms:njme) :: &
W1,W2,W3,W4
integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
real, intent(in) :: ptop,pdtop
integer :: i,j,k
logical :: badbad
call c2n_fulldom(II,JJ,W1,W2,W3,W4,&
deta1,deta2, eta1,eta2, ptop,pdtop, &
parent_fis,cpint,ct,cpd,cq, &
cims, cime, cjms, cjme, ckms, ckme, &
nest_fis,npint,nt,npd,nq, &
out_iinfo,out_winfo,imask, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
end subroutine ext_c2n_fulldom
subroutine ext_n2c_fulldom ( &,3
deta1,deta2, eta1,eta2, ptop,pdtop, &
cpint,ct,cpd,cq, &
cids, cide, cjds, cjde, ckds, ckde, &
cims, cime, cjms, cjme, ckms, ckme, &
cits, cite, cjts, cjte, ckts, ckte, &
npint,nt,npd,nq, &
ipos,jpos, &
out_iinfo,out_winfo, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
! This subroutine is an alias for n2c_fulldom. It exists to allow
! the routine to be called from module_dm
use module_interp_store, only: parent_fis, nest_fis
use module_interp_nmm, only: n2c_fulldom
implicit none
integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme, &
cids, cide, cjds, cjde, ckds, ckde, &
cits, cite, cjts, cjte, ckts, ckte, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, &
ipos, jpos ! nest start loctions within parent
real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
real, intent(inout), dimension(cims:cime,cjms:cjme,1) :: cPD
real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_winfo
integer, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_iinfo
real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
real, intent(in), dimension(nims:nime,njms:njme,1) :: nPD
real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
real, intent(in) :: ptop,pdtop
call n2c_fulldom( &
deta1,deta2, eta1,eta2, ptop,pdtop, &
parent_fis,cpint,ct,cpd,cq, &
cids, cide, cjds, cjde, ckds, ckde, &
cims, cime, cjms, cjme, ckms, ckme, &
cits, cite, cjts, cjte, ckts, ckte, &
nest_fis,npint,nt,npd,nq, &
ipos,jpos, &
out_iinfo,out_winfo, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte)
end subroutine ext_n2c_fulldom
subroutine ext_c2b_fulldom (II,JJ,W1,W2,W3,W4,&,3
deta1,deta2, eta1,eta2, ptop,pdtop, &
cpint,ct,cpd,cq, &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, &
ibxs, ibxe, ibys, ibye, &
wbxs, wbxe, wbys, wbye, &
pdbxs, pdbxe, pdbys, pdbye, &
tbxs, tbxe, tbys, tbye, &
qbxs, qbxe, qbys, qbye)
use module_interp_nmm, only: c2b_fulldom
use module_interp_store, only: parent_fis, nest_fis
implicit none
integer, intent(in):: &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte
integer, parameter :: bdyw = 1
! Domain info:
real, intent(in), dimension(nims:nime,njms:njme) :: W1,W2,W3,W4
integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
real, intent(in) :: ptop,pdtop
! Parent fields:
real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD
! T, Q, PINT, PD boundary info:
real,dimension(nims:nime,nkms:nkme,bdyw) :: tbys,tbye,qbys,qbye
real,dimension(njms:njme,nkms:nkme,bdyw) :: tbxs,tbxe,qbxs,qbxe
real,dimension(nims:nime,1,bdyw) :: pdbys,pdbye
real,dimension(njms:njme,1,bdyw) :: pdbxs,pdbxe
! Weights and indices:
real,dimension(nims:nime,nkms:nkme,bdyw) :: wbys,wbye
real,dimension(njms:njme,nkms:nkme,bdyw) :: wbxs,wbxe
integer,dimension(nims:nime,nkms:nkme,bdyw) :: ibys,ibye
integer,dimension(njms:njme,nkms:nkme,bdyw) :: ibxs,ibxe
call c2b_fulldom (II,JJ,W1,W2,W3,W4,&
deta1,deta2, eta1,eta2, ptop,pdtop, &
parent_fis,cpint,ct,cpd,cq,nest_fis, &
cims, cime, cjms, cjme, ckms, ckme, &
nids, nide, njds, njde, nkds, nkde, &
nims, nime, njms, njme, nkms, nkme, &
nits, nite, njts, njte, nkts, nkte, &
ibxs, ibxe, ibys, ibye, &
wbxs, wbxe, wbys, wbye, &
pdbxs, pdbxe, pdbys, pdbye, &
tbxs, tbxe, tbys, tbye, &
qbxs, qbxe, qbys, qbye)
end subroutine ext_c2b_fulldom