share/module_interp_nmm.F

   1 ! FCOPY, ICOPY, etc.: C pre-processor macros for C->N and N->C
   2 ! horizontal interpolation.  These exist solely to reduce code
   3 ! complexity.
   4
   5 ! Copy from C array to N array:
   6 #define FCOPY(N,C,i,j,k) \
   7   N(k,i-nits+1)=W1(i,j)*C(II(i,j)  ,JJ(i,j)  ,k)  \
   8   +W2(i,j)*C(II(i,j)+1,JJ(i,j)  ,k)  \
   9   +W3(i,j)*C(II(i,j)+a,JJ(i,j)-1,k)  \
  10   +W4(i,j)*C(II(i,j)+a,JJ(i,j)+1,k)
  11
  12 ! Copy from C array to (k,used) indexed N array:
  13 #define uFCOPY(N,C,i,j,k) \
  14   N(k,used)=W1(i,j)*C(II(i,j)  ,JJ(i,j)  ,k)  \
  15   +W2(i,j)*C(II(i,j)+1,JJ(i,j)  ,k)  \
  16   +W3(i,j)*C(II(i,j)+a,JJ(i,j)-1,k)  \
  17   +W4(i,j)*C(II(i,j)+a,JJ(i,j)+1,k)
  18
  19 ! ICOPY: Grab from C array, but do not assign:
  20 #define ICOPY(C,i,j,k) \
  21    (W1(i,j)*C(II(i,j)  ,JJ(i,j)  ,k)  \
  22   + W2(i,j)*C(II(i,j)+1,JJ(i,j)  ,k)  \
  23   + W3(i,j)*C(II(i,j)+a,JJ(i,j)-1,k)  \
  24   + W4(i,j)*C(II(i,j)+a,JJ(i,j)+1,k))
  25
  26 ! IKJCOPY: grab from IKJ C array, do not assign:
  27 #define IKJCOPY(C,i,k,j) \
  28    (W1(i,j)*C(II(i,j)  ,k,JJ(i,j)  )  \
  29   + W2(i,j)*C(II(i,j)+1,k,JJ(i,j)  )  \
  30   + W3(i,j)*C(II(i,j)+a,k,JJ(i,j)-1)  \
  31   + W4(i,j)*C(II(i,j)+a,k,JJ(i,j)+1))
  32
  33 ! ICOPY2D: grab from IJ C array, do not assign:
  34 #define ICOPY2D(C,i,j) \
  35    (W1(i,j)*C(II(i,j)  ,JJ(i,j)  )  \
  36   + W2(i,j)*C(II(i,j)+1,JJ(i,j)  )  \
  37   + W3(i,j)*C(II(i,j)+a,JJ(i,j)-1)  \
  38   + W4(i,j)*C(II(i,j)+a,JJ(i,j)+1))
  39
  40 ! Copying from N array to C array:
  41 #define UPCOPY(C,N,i,j,k,ni,nj) \
  42   C(k,i-istart+1)=(N(ni,nj+2,k) \
  43  + N(ni-a  ,nj+1,k)+ N(ni+1-a,nj+1,k) \
  44  + N(ni-1,nj  ,k)+ N(ni,nj  ,k)  + N(ni+1,nj  ,k) \
  45  + N(ni-a  ,nj-1,k)+ N(ni+1-a,nj-1,k) \
  46  +  N(ni,nj-2,k) \
  47   ) / 9
  48
  49 ! Average to C points from N points without assignment:
  50 #define NGRAB(N,ni,nj,nk) \
  51 (N(ni,nj+2,nk) \
  52  + N(ni-a  ,nj+1,nk)+ N(ni+1-a,nj+1,nk)\
  53  + N(ni-1,nj  ,nk)  + N(ni,nj  ,nk)  + N(ni+1,nj  ,nk)\
  54  + N(ni-a  ,nj-1,nk)+ N(ni+1-a,nj-1,nk)\
  55  +  N(ni,nj-2,nk) \
  56   ) / 9
  57
  58 ! Average to C points from N points without assignment on an IKJ grid:
  59 #define NGRABIKJ(N,ni,nk,nj) \
  60 (N(ni,nk,nj+2) \
  61  + N(ni-a  ,nk,nj+1)+ N(ni+1-a,nk,nj+1)   \
  62  + N(ni-1,nk,nj  )  + N(ni,nk,nj  )  + N(ni+1,nk,nj  )\
  63  + N(ni-a  ,nk,nj-1)+ N(ni+1-a,nk,nj-1)\
  64  +  N(ni,nk,nj-2) \
  65   ) / 9
  66
  67 ! Average to C points from N points without assignment, no vertical levels:
  68 #define NGRAB2D(N,ni,nj) \
  69 (N(ni,nj+2) \
  70  + N(ni-a  ,nj+1)+ N(ni+1-a,nj+1)\
  71  + N(ni-1,nj  )  + N(ni,nj  )  + N(ni+1,nj  )\
  72  + N(ni-a  ,nj-1)+ N(ni+1-a,nj-1)\
  73  +  N(ni,nj-2) \
  74   ) / 9
  75
  76 ! Copying from N array to I array:
  77 #define N2ICOPY(C,N,i,j,k,ni,nj) \
  78   C(i,j,k)=( N(ni,nj+2,k) \
  79  + N(ni-a  ,nj+1,k)+ N(ni+1-a,nj+1,k) \
  80  + N(ni-1,nj  ,k)+ N(ni,nj  ,k)  + N(ni+1,nj  ,k) \
  81  + N(ni-a  ,nj-1,k)+ N(ni+1-a,nj-1,k) \
  82  +  N(ni,nj-2,k) \
  83   ) / 9
  84
  85 ! Maximum value from N array to C array:
  86 #define N2I_SET_MAX(C,N,i,j,k,ni,nj) \
  87   C(i,j,k)=max(C(i,j,k),max(N(ni,nj+2,k),   \
  88   max(N(ni-1,nj  ,k),max(N(ni,nj  ,k),  max(N(ni+1,nj  ,k),\
  89   max(N(ni-a  ,nj-1,k),max(N(ni+1-a,nj-1,k),\
  90  N(ni,nj-2,k) )))))))
  91
  92 ! Maximum value from N array to C array:
  93 #define N2I_SET_MAX_IJ(C,N,i,j,ni,nj) \
  94   C(i,j)=max(C(i,j), max(N(ni,nj+2),\
  95  max(N(ni-1,nj  ),max(N(ni,nj  ),  max(N(ni+1,nj  ),  \
  96    max(N(ni-a  ,nj-1),max(N(ni+1-a,nj-1), \
  97     N(ni,nj-2) )))))))
  98
  99
 100 module module_interp_nmm
 101   use module_model_constants, only: g, R_D, p608
 102   implicit none
 103
 104   private
 105   public :: find_kpres
 106
 107   public :: nmm_interp_pd, nmm_keep_pd, nmm_method_linear
 108
 109   public :: c2b_fulldom, c2n_fulldom, n2c_fulldom
 110   public :: n2c_max2d, n2c_max3d
 111   public :: c2b_fulldom_new, n2c_fulldom_new
 112   public :: c2b_mass, c2n_mass, n2c_mass
 113   public :: c2n_massikj, n2c_massikj
 114   public :: c2b_copy3d, c2n_copy3d, n2c_copy3d
 115
 116   public :: c2b_copy2d, c2n_copy2d, n2c_copy2d, c2n_copy2d_nomask
 117   public :: c2b_near2d, c2n_near2d, n2c_near2d
 118   public :: c2b_inear2d, c2n_inear2d, n2c_inear2d
 119
 120   public :: c2n_near3dikj, c2n_sst, c2n_near3d
 121
 122 ! unimplemented:  public :: nmm_method_quadratic, nmm_method_spline_log, nmm_method_copy, nmm_method_spline
 123
 124 !  integer, parameter :: nmm_method_copy = 5
 125 !  integer, parameter :: nmm_method_spline_log = 4
 126 !  integer, parameter :: nmm_method_spline = 3
 127 !  integer, parameter :: nmm_method_quadratic = 2
 128
 129   integer, parameter :: nmm_method_linear = 1
 130
 131   integer, parameter :: nmm_interp_pd = 1
 132   integer, parameter :: nmm_keep_pd = 0
 133
 134   ! Constants from the original base_state_parent:
 135   REAL, PARAMETER :: LAPSR=6.5E-3, GI=1./G,D608=0.608
 136   REAL, PARAMETER :: COEF3=287.05*GI*LAPSR, COEF2=-1./COEF3
 137   REAL, PARAMETER :: TRG=2.0*R_D*GI,LAPSI=1.0/LAPSR
 138   REAL, PARAMETER :: P_REF=103000.
 139
 140   integer, parameter :: EConst=0, ECopy=1, EExtrap=2 ! MUST match module_dm
 141
 142 contains
 143   ! ********************************************************************
 144   ! subs *_NEAR2D -- horizontally interpolate via nearest neighbor
 145   ! method, but do not vertically interpolate.  Only handles 2D H grid
 146   ! arrays.
 147   ! ********************************************************************
 148
 149   subroutine c2b_near2d (inear,jnear,cfield,         &
 150        fbxs,fbxe,fbys,fbye,                  &
 151        cims, cime, cjms, cjme,   &
 152        nids, nide, njds, njde,   &
 153        nims, nime, njms, njme,   &
 154        nits, nite, njts, njte)
 155     implicit none
 156     integer, parameter :: bdyw = 1
 157
 158     integer, intent(in):: &
 159          cims, cime, cjms, cjme,   &
 160          nids, nide, njds, njde,   &
 161          nims, nime, njms, njme,   &
 162          nits, nite, njts, njte
 163
 164     integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
 165     real, intent(in) :: cfield(cims:cime,cjms:cjme)
 166
 167     ! Output field boundary info:
 168     real,dimension(nims:nime,bdyw) :: fbys,fbye
 169     real,dimension(njms:njme,bdyw) :: fbxs,fbxe
 170
 171     integer :: j,i,a,nx,nz,ck,k,j1,j2,add
 172     real :: weight
 173
 174     nx=min(nide-1,nite)-max(nids,nits)+1
 175
 176     j1=max(njts-1,njds)
 177     if(mod(j1,2)/=1)   j1=j1+1
 178
 179     j2=min(njte+1,njde-1)
 180     if(mod(j2,2)/=1)   j2=j2-1
 181
 182     if(nits==1) then
 183        i=1
 184        do j=j1,j2,2
 185           fbxs(j,1)=cfield(inear(i,j),jnear(i,j))
 186        enddo
 187     endif
 188     if(nite>=nide-1) then
 189        i=nide-1
 190        do j=j1,j2,2
 191           fbxe(j,1)=cfield(inear(i,j),jnear(i,j))
 192        enddo
 193     endif
 194     if(njts==1) then
 195        j=1
 196        do i=max(nits-1,nids),min(nite+1,nide-1)
 197           fbys(i,1)=cfield(inear(i,j),jnear(i,j))
 198        enddo
 199     endif
 200     if(njte>=njde-1) then
 201        j=njde-1
 202        do i=max(nits-1,nids),min(nite+1,nide-1)
 203           fbye(i,1)=cfield(inear(i,j),jnear(i,j))
 204        enddo
 205     endif
 206   end subroutine c2b_near2d
 207
 208   subroutine n2c_near2d (&
 209        cfield,nfield,ipos,jpos,  &
 210        cids, cide, cjds, cjde,   &
 211        cims, cime, cjms, cjme,   &
 212        cits, cite, cjts, cjte,   &
 213        nids, nide, njds, njde,   &
 214        nims, nime, njms, njme,   &
 215        nits, nite, njts, njte)
 216     implicit none
 217     integer, intent(in) :: &
 218          cids, cide, cjds, cjde,   &
 219          cims, cime, cjms, cjme,   &
 220          cits, cite, cjts, cjte,   &
 221          nids, nide, njds, njde,   &
 222          nims, nime, njms, njme,   &
 223          nits, nite, njts, njte,   &
 224          ipos,jpos
 225     real, intent(out) :: cfield(cims:cime,cjms:cjme)
 226     real, intent(in) :: nfield(nims:nime,njms:njme)
 227
 228     integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
 229     integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j
 230
 231     jstart=MAX(jpos+1,cjts)
 232     jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
 233
 234     bigj: do j=jstart,jend
 235        a=mod(j,2)
 236        istart=MAX(ipos+a,cits)
 237        iend=MIN(ipos+(nide-nids)/nri-1,cite)
 238        nj = (j-jpos)*nrj + 1
 239
 240        iloop: do i=istart,iend
 241           ni = (i-ipos)*nri + 2 - a
 242           cfield(i,j)=nfield(ni,nj)
 243        enddo iloop
 244     enddo bigj
 245   end subroutine n2c_near2d
 246
 247     subroutine c2n_near2d (inear,jnear,    &
 248        cfield,nfield,imask,     &
 249        cims, cime, cjms, cjme,   &
 250        nids, nide, njds, njde,   &
 251        nims, nime, njms, njme,   &
 252        nits, nite, njts, njte)
 253     implicit none
 254
 255     integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
 256
 257     integer, intent(in):: cims, cime, cjms, cjme,   &
 258          nids, nide, njds, njde,   &
 259          nims, nime, njms, njme,   &
 260          nits, nite, njts, njte
 261     integer, intent(in), dimension(nims:nime,njms:njme) :: imask
 262     real, intent(in) :: cfield(cims:cime,cjms:cjme)
 263     real, intent(out) :: nfield(nims:nime,njms:njme)
 264
 265     integer :: j,i,a,nx
 266
 267     nx=min(nide-1,nite)-max(nids,nits)+1
 268
 269     bigj: do j=max(njds,njts),min(njde-1,njte)
 270        interploop: do i=max(nids,nits),min(nide-1,nite)
 271           if(imask(i,j)/=0) cycle interploop
 272           nfield(i,j)=cfield(inear(i,j),jnear(i,j))
 273        enddo interploop
 274     end do bigj
 275   end subroutine c2n_near2d
 276
 277     subroutine c2n_near3d (inear,jnear,    &
 278        cfield,nfield,imask,     &
 279        cims, cime, cjms, cjme, ckms, ckme,   &
 280        nids, nide, njds, njde, nkds, nkde,   &
 281        nims, nime, njms, njme, nkms, nkme,   &
 282        nits, nite, njts, njte, nkts, nkte)
 283     implicit none
 284
 285     integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
 286
 287     integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme,   &
 288          nids, nide, njds, njde, nkds, nkde,   &
 289          nims, nime, njms, njme, nkms, nkme,   &
 290          nits, nite, njts, njte, nkts, nkte
 291     integer, intent(in), dimension(nims:nime,njms:njme) :: imask
 292     real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
 293     real, intent(out) :: nfield(nims:nime,njms:njme,nkms:nkme)
 294
 295     integer :: j,i,a,nx,k
 296
 297     nx=min(nide-1,nite)-max(nids,nits)+1
 298
 299     bigk: do k=nkds,nkde
 300        bigj: do j=max(njds,njts),min(njde-1,njte)
 301           interploop: do i=max(nids,nits),min(nide-1,nite)
 302              if(imask(i,j)/=0) cycle interploop
 303              nfield(i,j,k)=cfield(inear(i,j),jnear(i,j),k)
 304           enddo interploop
 305        end do bigj
 306     enddo bigk
 307   end subroutine c2n_near3d
 308
 309     subroutine c2n_near3dikj (inear,jnear,    &
 310        cfield,nfield,imask,     &
 311        cims, cime, cjms, cjme, ckms, ckme,   &
 312        nids, nide, njds, njde, nkds, nkde,   &
 313        nims, nime, njms, njme, nkms, nkme,   &
 314        nits, nite, njts, njte, nkts, nkte)
 315     implicit none
 316
 317     integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
 318
 319     integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme,   &
 320          nids, nide, njds, njde, nkds, nkde,   &
 321          nims, nime, njms, njme, nkms, nkme,   &
 322          nits, nite, njts, njte, nkts, nkte
 323     integer, intent(in), dimension(nims:nime,njms:njme) :: imask
 324     real, intent(in) :: cfield(cims:cime,ckms:ckme,cjms:cjme)
 325     real, intent(out) :: nfield(nims:nime,nkms:nkme,njms:njme)
 326
 327     integer :: j,i,a,nx,k
 328
 329     nx=min(nide-1,nite)-max(nids,nits)+1
 330
 331     bigj: do j=max(njds,njts),min(njde-1,njte)
 332        bigk: do k=nkds,nkde
 333           interploop: do i=max(nids,nits),min(nide-1,nite)
 334              if(imask(i,j)/=0) cycle interploop
 335              nfield(i,k,j)=cfield(inear(i,j),k,jnear(i,j))
 336           enddo interploop
 337        enddo bigk
 338     end do bigj
 339   end subroutine c2n_near3dikj
 340
 341     subroutine c2n_sst (inear,jnear,    &
 342        cfield,nfield,     &
 343        cims, cime, cjms, cjme,   &
 344        nids, nide, njds, njde,   &
 345        nims, nime, njms, njme,   &
 346        nits, nite, njts, njte)
 347     implicit none
 348
 349     integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
 350
 351     integer, intent(in):: cims, cime, cjms, cjme,   &
 352          nids, nide, njds, njde,   &
 353          nims, nime, njms, njme,   &
 354          nits, nite, njts, njte
 355     real, intent(in) :: cfield(cims:cime,cjms:cjme)
 356     real, intent(out) :: nfield(nims:nime,njms:njme)
 357
 358     integer :: j,i,a,nx
 359
 360     nx=min(nide-1,nite)-max(nids,nits)+1
 361
 362     bigj: do j=max(njds,njts),min(njde-1,njte)
 363        interploop: do i=max(nids,nits),min(nide-1,nite)
 364           nfield(i,j)=cfield(inear(i,j),jnear(i,j))
 365        enddo interploop
 366     end do bigj
 367   end subroutine c2n_sst
 368
 369   ! ********************************************************************
 370   ! subs *_INEAR2D -- horizontally interpolate integers via nearest
 371   ! neighbor method, but do not vertically interpolate.  Only handles
 372   ! 2D H grid integer arrays.
 373   ! ********************************************************************
 374
 375   subroutine c2b_inear2d (inear,jnear,cfield,         &
 376        fbxs,fbxe,fbys,fbye,                  &
 377        cims, cime, cjms, cjme,   &
 378        nids, nide, njds, njde,   &
 379        nims, nime, njms, njme,   &
 380        nits, nite, njts, njte)
 381     implicit none
 382     integer, parameter :: bdyw = 1
 383
 384     integer, intent(in):: &
 385          cims, cime, cjms, cjme,   &
 386          nids, nide, njds, njde,   &
 387          nims, nime, njms, njme,   &
 388          nits, nite, njts, njte
 389
 390     integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
 391     integer, intent(in) :: cfield(cims:cime,cjms:cjme)
 392
 393     ! Output field boundary info:
 394     integer,dimension(nims:nime,bdyw) :: fbys,fbye
 395     integer,dimension(njms:njme,bdyw) :: fbxs,fbxe
 396
 397     integer :: j,i,a,nx,nz,ck,k,j1,j2,add
 398
 399     nx=min(nide-1,nite)-max(nids,nits)+1
 400
 401     j1=max(njts-1,njds)
 402     if(mod(j1,2)/=1)   j1=j1+1
 403
 404     j2=min(njte+1,njde-1)
 405     if(mod(j2,2)/=1)   j2=j2-1
 406
 407     if(nits==1) then
 408        i=1
 409        do j=j1,j2,2
 410           fbxs(j,1)=cfield(inear(i,j),jnear(i,j))
 411        enddo
 412     endif
 413     if(nite>=nide-1) then
 414        i=nide-1
 415        do j=j1,j2,2
 416           fbxe(j,1)=cfield(inear(i,j),jnear(i,j))
 417        enddo
 418     endif
 419     if(njts==1) then
 420        j=1
 421        do i=max(nits-1,nids),min(nite+1,nide-1)
 422           fbys(i,1)=cfield(inear(i,j),jnear(i,j))
 423        enddo
 424     endif
 425     if(njte>=njde-1) then
 426        j=njde-1
 427        do i=max(nits-1,nids),min(nite+1,nide-1)
 428           fbye(i,1)=cfield(inear(i,j),jnear(i,j))
 429        enddo
 430     endif
 431   end subroutine c2b_inear2d
 432
 433   subroutine n2c_inear2d (&
 434        cfield,nfield,ipos,jpos,  &
 435        cids, cide, cjds, cjde,   &
 436        cims, cime, cjms, cjme,   &
 437        cits, cite, cjts, cjte,   &
 438        nids, nide, njds, njde,   &
 439        nims, nime, njms, njme,   &
 440        nits, nite, njts, njte)
 441     implicit none
 442     integer, intent(in) :: &
 443          cids, cide, cjds, cjde,   &
 444          cims, cime, cjms, cjme,   &
 445          cits, cite, cjts, cjte,   &
 446          nids, nide, njds, njde,   &
 447          nims, nime, njms, njme,   &
 448          nits, nite, njts, njte,   &
 449          ipos,jpos
 450     integer, intent(out) :: cfield(cims:cime,cjms:cjme)
 451     integer, intent(in) :: nfield(nims:nime,njms:njme)
 452
 453     integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
 454     integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j
 455
 456     jstart=MAX(jpos+1,cjts)
 457     jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
 458
 459     bigj: do j=jstart,jend
 460        a=mod(j,2)
 461        istart=MAX(ipos+a,cits)
 462        iend=MIN(ipos+(nide-nids)/nri-1,cite)
 463        nj = (j-jpos)*nrj + 1
 464
 465        iloop: do i=istart,iend
 466           ni = (i-ipos)*nri + 2 - a
 467           cfield(i,j)=nfield(ni,nj)
 468        enddo iloop
 469     enddo bigj
 470   end subroutine n2c_inear2d
 471
 472     subroutine c2n_inear2d (inear,jnear,    &
 473        cfield,nfield,imask,     &
 474        cims, cime, cjms, cjme,   &
 475        nids, nide, njds, njde,   &
 476        nims, nime, njms, njme,   &
 477        nits, nite, njts, njte)
 478     implicit none
 479
 480     integer, intent(in), dimension(nims:nime,njms:njme) :: inear,jnear
 481
 482     integer, intent(in):: cims, cime, cjms, cjme,   &
 483          nids, nide, njds, njde,   &
 484          nims, nime, njms, njme,   &
 485          nits, nite, njts, njte
 486     integer, intent(in), dimension(nims:nime,njms:njme) :: imask
 487     integer, intent(in) :: cfield(cims:cime,cjms:cjme)
 488     integer, intent(out) :: nfield(nims:nime,njms:njme)
 489
 490     integer :: j,i,a,nx
 491
 492     nx=min(nide-1,nite)-max(nids,nits)+1
 493
 494     bigj: do j=max(njds,njts),min(njde-1,njte)
 495        interploop: do i=max(nids,nits),min(nide-1,nite)
 496           if(imask(i,j)/=0) cycle interploop
 497           nfield(i,j)=cfield(inear(i,j),jnear(i,j))
 498        enddo interploop
 499     end do bigj
 500   end subroutine c2n_inear2d
 501
 502   ! ********************************************************************
 503   ! subs *_COPY2D -- horizontally interpolate but do not vertically
 504   ! interpolate.  Only handles 2D arrays (H or V grid).
 505   ! ********************************************************************
 506
 507   subroutine c2b_copy2d (II,JJ,W1,W2,W3,W4,cfield,         &
 508        fbxs,fbxe,fbys,fbye,                  &
 509        cims, cime, cjms, cjme,   &
 510        nids, nide, njds, njde,   &
 511        nims, nime, njms, njme,   &
 512        nits, nite, njts, njte, hgrid)
 513     implicit none
 514     integer, parameter :: bdyw = 1
 515
 516     integer, intent(in):: &
 517          cims, cime, cjms, cjme,   &
 518          nids, nide, njds, njde,   &
 519          nims, nime, njms, njme,   &
 520          nits, nite, njts, njte
 521     logical, intent(in) :: hgrid
 522     real, intent(in), dimension(nims:nime,njms:njme) :: &
 523          W1,W2,W3,W4
 524     integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
 525     real, intent(in) :: cfield(cims:cime,cjms:cjme)
 526
 527     ! Output field boundary info:
 528     real,dimension(nims:nime,bdyw) :: fbys,fbye
 529     real,dimension(njms:njme,bdyw) :: fbxs,fbxe
 530
 531     integer :: j,i,a,nx,nz,ck,k,j1,j2,add
 532     real :: weight
 533
 534     nx=min(nide-1,nite)-max(nids,nits)+1
 535
 536     if(hgrid) then
 537        a=1
 538     else
 539        a=0
 540     endif
 541
 542     j1=max(njts-1,njds)
 543     if(mod(j1,2)/=a)   j1=j1+1
 544
 545     j2=min(njte+1,njde-1)
 546     if(mod(j2,2)/=a)   j2=j2-1
 547
 548     if(nits==1) then
 549        i=1
 550        do j=j1,j2,2
 551           if(hgrid) then
 552              a=1-mod(JJ(i,j),2)
 553           else
 554              a=mod(JJ(i,j),2)
 555           endif
 556           fbxs(j,1)=ICOPY2D(cfield,i,j)
 557        enddo
 558     endif
 559     if(nite>=nide-1) then
 560        i=nide-1
 561        do j=j1,j2,2
 562           if(hgrid) then
 563              a=1-mod(JJ(i,j),2)
 564           else
 565              a=mod(JJ(i,j),2)
 566           endif
 567           fbxe(j,1)=ICOPY2D(cfield,i,j)
 568        enddo
 569     endif
 570     if(njts==1) then
 571        j=1
 572        do i=max(nits-1,nids),min(nite+1,nide-1)
 573           if(hgrid) then
 574              a=1-mod(JJ(i,j),2)
 575           else
 576              a=mod(JJ(i,j),2)
 577           endif
 578           fbys(i,1)=ICOPY2D(cfield,i,j)
 579        enddo
 580     endif
 581     if(njte>=njde-1) then
 582        j=njde-1
 583        do i=max(nits-1,nids),min(nite+1,nide-1)
 584           if(hgrid) then
 585              a=1-mod(JJ(i,j),2)
 586           else
 587              a=mod(JJ(i,j),2)
 588           endif
 589           fbye(i,1)=ICOPY2D(cfield,i,j)
 590        enddo
 591     endif
 592   end subroutine c2b_copy2d
 593   subroutine n2c_copy2d (&
 594        cfield,nfield,ipos,jpos,  &
 595        cids, cide, cjds, cjde,   &
 596        cims, cime, cjms, cjme,   &
 597        cits, cite, cjts, cjte,   &
 598        nids, nide, njds, njde,   &
 599        nims, nime, njms, njme,   &
 600        nits, nite, njts, njte, &
 601        hgrid)
 602     implicit none
 603     integer, intent(in) :: &
 604          cids, cide, cjds, cjde,   &
 605          cims, cime, cjms, cjme,   &
 606          cits, cite, cjts, cjte,   &
 607          nids, nide, njds, njde,   &
 608          nims, nime, njms, njme,   &
 609          nits, nite, njts, njte,   &
 610          ipos,jpos
 611     real, intent(out) :: cfield(cims:cime,cjms:cjme)
 612     real, intent(in) :: nfield(nims:nime,njms:njme)
 613     logical, intent(in) :: hgrid
 614
 615     integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
 616     integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j
 617
 618     jstart=MAX(jpos+1,cjts)
 619     jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
 620
 621     bigj: do j=jstart,jend
 622           a=mod(j,2)
 623           if(.not.hgrid) then
 624              a=1-a
 625           endif
 626        istart=MAX(ipos+a,cits)
 627        iend=MIN(ipos+(nide-nids)/nri-1,cite)
 628        nj = (j-jpos)*nrj + 1
 629
 630        iloop: do i=istart,iend
 631           ni = (i-ipos)*nri + 2 - a
 632           cfield(i,j)=NGRAB2D(nfield,ni,nj)
 633        enddo iloop
 634     enddo bigj
 635  end subroutine n2c_copy2d
 636
 637   subroutine n2c_max2d (&
 638        cfield,nfield,ipos,jpos,  &
 639        cids, cide, cjds, cjde,   &
 640        cims, cime, cjms, cjme,   &
 641        cits, cite, cjts, cjte,   &
 642        nids, nide, njds, njde,   &
 643        nims, nime, njms, njme,   &
 644        nits, nite, njts, njte, &
 645        hgrid)
 646     implicit none
 647     integer, intent(in) :: &
 648          cids, cide, cjds, cjde,   &
 649          cims, cime, cjms, cjme,   &
 650          cits, cite, cjts, cjte,   &
 651          nids, nide, njds, njde,   &
 652          nims, nime, njms, njme,   &
 653          nits, nite, njts, njte,   &
 654          ipos,jpos
 655     real, intent(out) :: cfield(cims:cime,cjms:cjme)
 656     real, intent(in) :: nfield(nims:nime,njms:njme)
 657     logical, intent(in) :: hgrid
 658
 659     integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
 660     integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j
 661
 662     jstart=MAX(jpos+1,cjts)
 663     jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
 664
 665     bigj: do j=jstart,jend
 666           a=mod(j,2)
 667           if(.not.hgrid) then
 668              a=1-a
 669           endif
 670        istart=MAX(ipos+a,cits)
 671        iend=MIN(ipos+(nide-nids)/nri-1,cite)
 672        nj = (j-jpos)*nrj + 1
 673
 674        iloop: do i=istart,iend
 675           ni = (i-ipos)*nri + 2 - a
 676           N2I_SET_MAX_IJ(cfield,nfield,i,j,ni,nj)
 677        enddo iloop
 678     enddo bigj
 679   end subroutine n2c_max2d
 680
 681     subroutine c2n_copy2d (II,JJ,W1,W2,W3,W4,    &
 682        cfield,nfield,imask,     &
 683        cims, cime, cjms, cjme,   &
 684        nids, nide, njds, njde,   &
 685        nims, nime, njms, njme,   &
 686        nits, nite, njts, njte, hgrid)
 687     implicit none
 688     logical, intent(in) :: hgrid
 689     real, intent(in), dimension(nims:nime,njms:njme) :: &
 690          W1,W2,W3,W4
 691     integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
 692
 693     integer, intent(in):: cims, cime, cjms, cjme,   &
 694          nids, nide, njds, njde,   &
 695          nims, nime, njms, njme,   &
 696          nits, nite, njts, njte
 697     integer, intent(in), dimension(nims:nime,njms:njme) :: imask
 698     real, intent(in) :: cfield(cims:cime,cjms:cjme)
 699     real, intent(out) :: nfield(nims:nime,njms:njme)
 700
 701     integer :: j,i,a,nx
 702
 703     nx=min(nide-1,nite)-max(nids,nits)+1
 704
 705     bigj: do j=max(njds,njts),min(njde-1,njte)
 706        interploop: do i=max(nids,nits),min(nide-1,nite)
 707           if(imask(i,j)/=0) cycle interploop
 708           if(hgrid) then
 709              a=1-mod(JJ(i,j),2)
 710           else
 711              a=mod(JJ(i,j),2)
 712           endif
 713           nfield(i,j)=ICOPY2D(cfield,i,j)
 714        enddo interploop
 715     end do bigj
 716   end subroutine c2n_copy2d
 717
 718     subroutine c2n_copy2d_nomask (II,JJ,W1,W2,W3,W4,    &
 719        cfield,nfield,     &
 720        cims, cime, cjms, cjme,   &
 721        nids, nide, njds, njde,   &
 722        nims, nime, njms, njme,   &
 723        nits, nite, njts, njte, hgrid)
 724     implicit none
 725     logical, intent(in) :: hgrid
 726     real, intent(in), dimension(nims:nime,njms:njme) :: &
 727          W1,W2,W3,W4
 728     integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
 729
 730     integer, intent(in):: cims, cime, cjms, cjme,   &
 731          nids, nide, njds, njde,   &
 732          nims, nime, njms, njme,   &
 733          nits, nite, njts, njte
 734     real, intent(in) :: cfield(cims:cime,cjms:cjme)
 735     real, intent(out) :: nfield(nims:nime,njms:njme)
 736
 737     integer :: j,i,a,nx
 738
 739     nx=min(nide-1,nite)-max(nids,nits)+1
 740
 741     bigj: do j=max(njds,njts),min(njde-1,njte)
 742        interploop: do i=max(nids,nits),min(nide-1,nite)
 743           if(hgrid) then
 744              a=1-mod(JJ(i,j),2)
 745           else
 746              a=mod(JJ(i,j),2)
 747           endif
 748           nfield(i,j)=ICOPY2D(cfield,i,j)
 749        enddo interploop
 750     end do bigj
 751   end subroutine c2n_copy2d_nomask
 752
 753
 754   ! ********************************************************************
 755   ! subs *_COPY3D -- horizontally interpolate but do not vertically
 756   ! interpolate
 757   ! ********************************************************************
 758
 759   subroutine c2b_copy3d  (II,JJ,W1,W2,W3,W4,cfield,         &
 760        fbxs,fbxe,fbys,fbye,                  &
 761        cims, cime, cjms, cjme, ckms, ckme,   &
 762        nids, nide, njds, njde, nkds, nkde,   &
 763        nims, nime, njms, njme, nkms, nkme,   &
 764        nits, nite, njts, njte, nkts, nkte, hgrid)
 765     implicit none
 766     integer, parameter :: bdyw = 1
 767     logical, intent(in) :: hgrid
 768     integer, intent(in):: &
 769          cims, cime, cjms, cjme, ckms, ckme,   &
 770          nids, nide, njds, njde, nkds, nkde,   &
 771          nims, nime, njms, njme, nkms, nkme,   &
 772          nits, nite, njts, njte, nkts, nkte
 773
 774     real, intent(in), dimension(nims:nime,njms:njme) :: &
 775          W1,W2,W3,W4
 776     integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
 777     real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
 778
 779     ! Output field boundary info:
 780     real,dimension(nims:nime,nkms:nkme,bdyw) :: fbys,fbye
 781     real,dimension(njms:njme,nkms:nkme,bdyw) :: fbxs,fbxe
 782
 783     integer :: j,i,a,nx,nz,ck,k,j1,j2
 784     real :: weight
 785
 786     nx=min(nide-1,nite)-max(nids,nits)+1
 787     nz=nkde-nkds+1
 788
 789     j1=max(njts-1,njds)
 790     if(mod(j1,2)/=0)   j1=j1+1
 791
 792     j2=min(njte+1,njde-1)
 793     if(mod(j2,2)/=0)   j2=j2-1
 794
 795     if(nits==1) then
 796        i=1
 797        do j=j1,j2,2
 798           if(hgrid) then
 799              a=1-mod(JJ(i,j),2)
 800           else
 801              a=mod(JJ(i,j),2)
 802           endif
 803           kloop1: do k=nkds,min(nkde,nkte)
 804              fbxs(j,k,1)=ICOPY(cfield,i,j,k)
 805           enddo kloop1
 806        enddo
 807     endif
 808     if(nite>=nide-1) then
 809        i=nide-1
 810        do j=j1,j2,2
 811           if(hgrid) then
 812              a=1-mod(JJ(i,j),2)
 813           else
 814              a=mod(JJ(i,j),2)
 815           endif
 816           kloop2: do k=nkds,min(nkde,nkte)
 817              fbxe(j,k,1)=ICOPY(cfield,i,j,k)
 818           enddo kloop2
 819        enddo
 820     endif
 821     if(njts==1) then
 822        j=1
 823        do i=max(nits-1,nids),min(nite+1,nide-1)
 824           if(hgrid) then
 825              a=1-mod(JJ(i,j),2)
 826           else
 827              a=mod(JJ(i,j),2)
 828           endif
 829           kloop3: do k=nkts,min(nkde,nkte)
 830              fbys(i,k,1)=ICOPY(cfield,i,j,k)
 831           enddo kloop3
 832        enddo
 833     endif
 834     if(njte>=njde-1) then
 835        j=njde-1
 836        do i=max(nits-1,nids),min(nite+1,nide-1)
 837           if(hgrid) then
 838              a=1-mod(JJ(i,j),2)
 839           else
 840              a=mod(JJ(i,j),2)
 841           endif
 842           kloop4: do k=nkts,min(nkde,nkte)
 843              fbye(i,k,1)=ICOPY(cfield,i,j,k)
 844           enddo kloop4
 845        enddo
 846     endif
 847   end subroutine c2b_copy3d
 848
 849   subroutine n2c_copy3d (&
 850        cfield,nfield,ipos,jpos,  &
 851        cids, cide, cjds, cjde, ckds, ckde,   &
 852        cims, cime, cjms, cjme, ckms, ckme,   &
 853        cits, cite, cjts, cjte, ckts, ckte,   &
 854        nids, nide, njds, njde, nkds, nkde,   &
 855        nims, nime, njms, njme, nkms, nkme,   &
 856        nits, nite, njts, njte, nkts, nkte,   &
 857        hgrid)
 858     implicit none
 859     logical, intent(in) :: hgrid
 860     integer, intent(in) :: &
 861          cids, cide, cjds, cjde, ckds, ckde,   &
 862          cims, cime, cjms, cjme, ckms, ckme,   &
 863          cits, cite, cjts, cjte, ckts, ckte,   &
 864          nids, nide, njds, njde, nkds, nkde,   &
 865          nims, nime, njms, njme, nkms, nkme,   &
 866          nits, nite, njts, njte, nkts, nkte,   &
 867          ipos,jpos
 868     real, intent(out) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
 869     real, intent(in) :: nfield(nims:nime,njms:njme,nkms:nkme)
 870
 871     integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
 872     integer :: nx,nz,jstart,jend,istart,iend,ni,nj,a,i,j,k,nk
 873     real :: weight
 874
 875     nx=min(nide-2,nite)-max(nids+1,nits)+1
 876     nz=nkde-nkds+1
 877
 878     jstart=MAX(jpos+1,cjts)
 879     jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
 880
 881     bigj: do j=jstart,jend
 882        if(hgrid) then
 883           a=mod(j,2)
 884        else
 885           a=1-mod(j,2)
 886        endif
 887        istart=MAX(ipos+a,cits)
 888        iend=MIN(ipos+(nide-nids)/nri-1,cite)
 889        nj = (j-jpos)*nrj + 1
 890
 891        iloop: do i=istart,iend
 892           do k=nkds,nkde
 893              ni = (i-ipos)*nri + 2 - a
 894              cfield(i,j,k)=NGRAB(nfield,ni,nj,k)
 895           enddo
 896        enddo iloop
 897     enddo bigj
 898   end subroutine n2c_copy3d
 899
 900   subroutine n2c_max3d (&
 901        cfield,nfield,ipos,jpos,  &
 902        cids, cide, cjds, cjde, ckds, ckde,   &
 903        cims, cime, cjms, cjme, ckms, ckme,   &
 904        cits, cite, cjts, cjte, ckts, ckte,   &
 905        nids, nide, njds, njde, nkds, nkde,   &
 906        nims, nime, njms, njme, nkms, nkme,   &
 907        nits, nite, njts, njte, nkts, nkte,   &
 908        hgrid)
 909     implicit none
 910     logical, intent(in) :: hgrid
 911     integer, intent(in) :: &
 912          cids, cide, cjds, cjde, ckds, ckde,   &
 913          cims, cime, cjms, cjme, ckms, ckme,   &
 914          cits, cite, cjts, cjte, ckts, ckte,   &
 915          nids, nide, njds, njde, nkds, nkde,   &
 916          nims, nime, njms, njme, nkms, nkme,   &
 917          nits, nite, njts, njte, nkts, nkte,   &
 918          ipos,jpos
 919     real, intent(out) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
 920     real, intent(in) :: nfield(nims:nime,njms:njme,nkms:nkme)
 921
 922     integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
 923     integer :: nx,nz,jstart,jend,istart,iend,ni,nj,a,i,j,k,nk
 924     real :: weight
 925
 926     nx=min(nide-2,nite)-max(nids+1,nits)+1
 927     nz=nkde-nkds+1
 928
 929     jstart=MAX(jpos+1,cjts)
 930     jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
 931
 932     bigj: do j=jstart,jend
 933        if(hgrid) then
 934           a=mod(j,2)
 935        else
 936           a=1-mod(j,2)
 937        endif
 938        istart=MAX(ipos+a,cits)
 939        iend=MIN(ipos+(nide-nids)/nri-1,cite)
 940        nj = (j-jpos)*nrj + 1
 941
 942        iloop: do i=istart,iend
 943           do k=nkds,nkde
 944              ni = (i-ipos)*nri + 2 - a
 945              N2I_SET_MAX(cfield,nfield,i,j,k,ni,nj)
 946           enddo
 947        enddo iloop
 948     enddo bigj
 949   end subroutine n2c_max3d
 950
 951   subroutine c2n_copy3d (II,JJ,W1,W2,W3,W4,    &
 952        cfield,nfield,imask,      &
 953        cims, cime, cjms, cjme, ckms, ckme,   &
 954        nids, nide, njds, njde, nkds, nkde,   &
 955        nims, nime, njms, njme, nkms, nkme,   &
 956        nits, nite, njts, njte, nkts, nkte,   &
 957        hgrid)
 958     implicit none
 959     real, intent(in), dimension(nims:nime,njms:njme) :: &
 960          W1,W2,W3,W4
 961     integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
 962     logical, intent(in) :: hgrid
 963     integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme,   &
 964          nids, nide, njds, njde, nkds, nkde,   &
 965          nims, nime, njms, njme, nkms, nkme,   &
 966          nits, nite, njts, njte, nkts, nkte
 967     integer, intent(in), dimension(nims:nime,njms:njme) :: imask
 968     real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
 969     real, intent(out) :: nfield(nims:nime,njms:njme,nkms:nkme)
 970
 971     integer :: j,i,a,nx,nz,k
 972
 973     nx=min(nide-1,nite)-max(nids,nits)+1
 974     nz=nkde-nkds+1
 975
 976     bigj: do j=max(njds,njts),min(njde-1,njte)
 977        interploop: do i=max(nids,nits),min(nide-1,nite)
 978           if(imask(i,j)/=0) cycle interploop
 979           kinterploop: do k=nkds,nkde
 980              if(hgrid) then
 981                 a=1-mod(JJ(i,j),2)
 982              else
 983                 a=mod(JJ(i,j),2)
 984              endif
 985              nfield(i,j,k) = ICOPY(cfield,i,j,k)
 986           enddo kinterploop
 987        enddo interploop
 988     end do bigj
 989   end subroutine c2n_copy3d
 990
 991   ! ********************************************************************
 992   ! subs *_MASS -- horizontally and vertically interpolate.  Vertical
 993   ! interpolation uses the iinfo and winfo arrays (and equivalent
 994   ! boundary arrays) generated in the *_fulldom subroutines.
 995   ! ********************************************************************
 996
 997   subroutine c2b_mass (II,JJ,W1,W2,W3,W4,cfield,         &
 998        ibxs,ibxe,ibys,ibye,                  &
 999        wbxs,wbxe,wbys,wbye,                  &
1000        fbxs,fbxe,fbys,fbye,                  &
1001        emethod,evalue,                       &
1002        cims, cime, cjms, cjme, ckms, ckme,   &
1003        nids, nide, njds, njde, nkds, nkde,   &
1004        nims, nime, njms, njme, nkms, nkme,   &
1005        nits, nite, njts, njte, nkts, nkte)
1006     use module_interp_store, only: kpres
1007     implicit none
1008     integer, parameter :: bdyw = 1
1009
1010     integer, intent(in):: &
1011          cims, cime, cjms, cjme, ckms, ckme,   &
1012          nids, nide, njds, njde, nkds, nkde,   &
1013          nims, nime, njms, njme, nkms, nkme,   &
1014          nits, nite, njts, njte, nkts, nkte
1015
1016     integer, intent(in) :: emethod ! extrapolation method
1017     real, intent(in) :: evalue
1018     real, intent(in), dimension(nims:nime,njms:njme) :: &
1019          W1,W2,W3,W4
1020     integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
1021     real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
1022
1023     ! Output field boundary info:
1024     real,dimension(nims:nime,nkms:nkme,bdyw) :: fbys,fbye
1025     real,dimension(njms:njme,nkms:nkme,bdyw) :: fbxs,fbxe
1026
1027     ! Weights and indices:
1028     real,dimension(nims:nime,nkms:nkme,bdyw) :: wbys,wbye
1029     real,dimension(njms:njme,nkms:nkme,bdyw) :: wbxs,wbxe
1030     integer,dimension(nims:nime,nkms:nkme,bdyw) :: ibys,ibye
1031     integer,dimension(njms:njme,nkms:nkme,bdyw) :: ibxs,ibxe
1032
1033     integer :: j,i,a,nx,nz,ck,k,j1,j2
1034     real :: weight
1035
1036     nx=min(nide-1,nite)-max(nids,nits)+1
1037     nz=nkde-nkds+1
1038
1039     j1=max(njts-1,njds)
1040     if(mod(j1,2)/=1)   j1=j1+1
1041
1042     j2=min(njte+1,njde-1)
1043     if(mod(j2,2)/=1)   j2=j2-1
1044
1045     if(nits==1) then
1046        i=1
1047        do j=j1,j2,2
1048           a=1-mod(JJ(i,j),2)
1049           kcopy1: do k=min(nkde,nkte),kpres+1,-1
1050              fbxs(j,k,1) = ICOPY(cfield,i,j,k)
1051           enddo kcopy1
1052           kloop1: do k=kpres,nkds,-1
1053              weight=wbxs(j,k,1)
1054              ck=ibxs(j,k,1)
1055
1056              if(ck>1) then
1057                 fbxs(j,k,1) = &
1058                      ICOPY(cfield,i,j,ck)   * weight + &
1059                      ICOPY(cfield,i,j,ck-1) * (1.0-weight)
1060              else
1061                 ! Extrapolate
1062                 if(emethod==EConst) then
1063                    fbxs(j,k,1)=evalue
1064                 else if(emethod==ECopy) then
1065                    fbxs(j,k,1)=ICOPY(cfield,i,j,1)
1066                 else ! assume 2: linear extrap
1067                    fbxs(j,k,1)=evalue * weight + &
1068                                ICOPY(cfield,i,j,1) * (1.0-weight)
1069                 endif
1070              endif
1071           enddo kloop1
1072        enddo
1073     endif
1074     if(nite>=nide-1) then
1075        i=nide-1
1076        do j=j1,j2,2
1077           a=1-mod(JJ(i,j),2)
1078           kcopy2: do k=min(nkde,nkte),kpres+1,-1
1079              fbxe(j,k,1)=ICOPY(cfield,i,j,k)
1080           enddo kcopy2
1081           kloop2: do k=kpres,nkds,-1
1082              weight=wbxe(j,k,1)
1083              ck=ibxe(j,k,1)
1084
1085              if(ck>1) then
1086                 fbxe(j,k,1) = &
1087                      ICOPY(cfield,i,j,ck)   * weight + &
1088                      ICOPY(cfield,i,j,ck-1) * (1.0-weight)
1089              else
1090                 ! Extrapolate
1091                 if(emethod==EConst) then
1092                    fbxe(j,k,1)=evalue
1093                 else if(emethod==ECopy) then
1094                    fbxe(j,k,1)=ICOPY(cfield,i,j,1)
1095                 else ! assume 2: linear extrap
1096                    fbxe(j,k,1)=evalue * weight + &
1097                         ICOPY(cfield,i,j,1) * (1.0-weight)
1098                 endif
1099              endif
1100           enddo kloop2
1101        enddo
1102     endif
1103     if(njts==1) then
1104        j=1
1105        do i=max(nits-1,nids),min(nite+1,nide-1)
1106           a=1-mod(JJ(i,j),2)
1107           kcopy3: do k=min(nkde,nkte),kpres+1,-1
1108              fbys(i,k,1) = ICOPY(cfield,i,j,k)
1109           enddo kcopy3
1110           kloop3: do k=kpres,nkds,-1
1111              weight=wbys(i,k,1)
1112              ck=ibys(i,k,1)
1113
1114              if(ck>1) then
1115                 fbys(i,k,1) = &
1116                      ICOPY(cfield,i,j,ck)   * weight + &
1117                      ICOPY(cfield,i,j,ck-1) * (1.0-weight)
1118              else
1119                 ! Extrapolate
1120                 if(emethod==EConst) then
1121                    fbys(i,k,1)=evalue
1122                 else if(emethod==ECopy) then
1123                    fbys(i,k,1)=ICOPY(cfield,i,j,1)
1124                 else ! assume 2: linear extrap
1125                    fbys(i,k,1)=evalue*weight + &
1126                         ICOPY(cfield,i,j,1)*(1.0-weight)
1127                 endif
1128              endif
1129           enddo kloop3
1130        enddo
1131     endif
1132     if(njte>=njde-1) then
1133        j=njde-1
1134        do i=max(nits-1,nids),min(nite+1,nide-1)
1135           a=1-mod(JJ(i,j),2)
1136           kcopy4: do k=min(nkde,nkte),kpres+1,-1
1137              fbye(i,k,1) = ICOPY(cfield,i,j,k)
1138           enddo kcopy4
1139           kloop4: do k=kpres,nkds,-1
1140              weight=wbye(i,k,1)
1141              ck=ibye(i,k,1)
1142
1143              if(ck>1) then
1144                 fbye(i,k,1) = &
1145                      ICOPY(cfield,i,j,ck)   * weight + &
1146                      ICOPY(cfield,i,j,ck-1) * (1.0-weight)
1147              else
1148                 if(emethod==EConst) then
1149                    fbye(i,k,1)=evalue
1150                 else if(emethod==ECopy) then
1151                    fbye(i,k,1)=ICOPY(cfield,i,j,1)
1152                 else ! assume 2: linear extrap
1153                    fbye(i,k,1)=evalue*weight + &
1154                         ICOPY(cfield,i,j,1)*(1.0-weight)
1155                 endif
1156              endif
1157           enddo kloop4
1158        enddo
1159     endif
1160
1161   end subroutine c2b_mass
1162
1163   subroutine n2c_mass (&
1164        cfield,nfield,iinfo,winfo,ipos,jpos,  &
1165        emethod, evalue,                      &
1166        cids, cide, cjds, cjde, ckds, ckde,   &
1167        cims, cime, cjms, cjme, ckms, ckme,   &
1168        cits, cite, cjts, cjte, ckts, ckte,   &
1169        nids, nide, njds, njde, nkds, nkde,   &
1170        nims, nime, njms, njme, nkms, nkme,   &
1171        nits, nite, njts, njte, nkts, nkte)
1172     use module_interp_store, only: kpres
1173     implicit none
1174     integer, intent(in) :: &
1175          cids, cide, cjds, cjde, ckds, ckde,   &
1176          cims, cime, cjms, cjme, ckms, ckme,   &
1177          cits, cite, cjts, cjte, ckts, ckte,   &
1178          nids, nide, njds, njde, nkds, nkde,   &
1179          nims, nime, njms, njme, nkms, nkme,   &
1180          nits, nite, njts, njte, nkts, nkte,   &
1181          ipos,jpos, emethod
1182     real, intent(in) :: evalue
1183     real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: winfo
1184     integer, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: iinfo
1185     real, intent(out) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
1186     real, intent(in) :: nfield(nims:nime,njms:njme,nkms:nkme)
1187
1188     integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
1189     integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j,k,nk
1190     real :: weight
1191
1192     nz=nkde-nkds+1
1193
1194     jstart=MAX(jpos+1,cjts)
1195     jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
1196
1197     kcopy: do k=ckde,kpres+1,-1
1198        jcopy: do j=jstart,jend
1199           a=mod(j,2)
1200           istart=MAX(ipos+a,cits)
1201           iend=MIN(ipos+(nide-nids)/nri-1,cite)
1202           nj = (j-jpos)*nrj + 1
1203           icopy: do i=istart,iend
1204              ni = (i-ipos)*nri + 2 - a
1205              cfield(i,j,k) = NGRAB(nfield,ni,nj,k)
1206           enddo icopy
1207        enddo jcopy
1208     enddo kcopy
1209
1210     bigj: do j=jstart,jend
1211           a=mod(j,2)
1212        istart=MAX(ipos+a,cits)
1213        iend=MIN(ipos+(nide-nids)/nri-1,cite)
1214        nj = (j-jpos)*nrj + 1
1215
1216        iloop: do i=istart,iend
1217           ni = (i-ipos)*nri + 2 - a
1218           kinterploop: do k=kpres,ckds,-1
1219              weight=winfo(i,j,k)
1220              nk=iinfo(i,j,k)
1221
1222              if(nk>1) then
1223                 cfield(i,j,k) = &
1224                      NGRAB(nfield,ni,nj,nk)   * weight + &
1225 ! pjj/cray - source line limit in Cray compiler
1226              NGRAB(nfield,ni,nj,nk-1) &
1227              * (1.0-weight)
1228              else
1229                 ! Extrapolate
1230                 if(emethod==EConst) then
1231                    cfield(i,j,k)=evalue
1232                 elseif(emethod==ECopy) then
1233                    cfield(i,j,k)=NGRAB(nfield,ni,nj,1)
1234                 else ! Assume 2: linear extrap
1235                    cfield(i,j,k)=evalue*weight + &
1236                         NGRAB(nfield,ni,nj,1)*(1.0-weight)
1237                 endif
1238              endif
1239           end do kinterploop
1240        enddo iloop
1241     enddo bigj
1242   end subroutine n2c_mass
1243
1244   subroutine n2c_massikj (&
1245        cfield,nfield,iinfo,winfo,ipos,jpos,  &
1246        emethod, evalue,                      &
1247        cids, cide, cjds, cjde, ckds, ckde,   &
1248        cims, cime, cjms, cjme, ckms, ckme,   &
1249        cits, cite, cjts, cjte, ckts, ckte,   &
1250        nids, nide, njds, njde, nkds, nkde,   &
1251        nims, nime, njms, njme, nkms, nkme,   &
1252        nits, nite, njts, njte, nkts, nkte)
1253     use module_interp_store, only: kpres
1254     implicit none
1255     integer, intent(in) :: &
1256          cids, cide, cjds, cjde, ckds, ckde,   &
1257          cims, cime, cjms, cjme, ckms, ckme,   &
1258          cits, cite, cjts, cjte, ckts, ckte,   &
1259          nids, nide, njds, njde, nkds, nkde,   &
1260          nims, nime, njms, njme, nkms, nkme,   &
1261          nits, nite, njts, njte, nkts, nkte,   &
1262          ipos,jpos, emethod
1263     real, intent(in) :: evalue
1264     real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: winfo
1265     integer, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: iinfo
1266     real, intent(out) :: cfield(cims:cime,ckms:ckme,cjms:cjme)
1267     real, intent(in) :: nfield(nims:nime,nkms:nkme,njms:njme)
1268
1269     integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
1270     integer :: nz,jstart,jend,istart,iend,ni,nj,a,i,j,k,nk
1271     real :: weight
1272
1273     nz=nkde-nkds+1
1274
1275     jstart=MAX(jpos+1,cjts)
1276     jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
1277
1278     bigj: do j=jstart,jend
1279           a=mod(j,2)
1280        istart=MAX(ipos+a,cits)
1281        iend=MIN(ipos+(nide-nids)/nri-1,cite)
1282        nj = (j-jpos)*nrj + 1
1283
1284        kcopyloop: do k=nkde,kpres+1,-1
1285           icopyloop: do i=istart,iend
1286              ni = (i-ipos)*nri + 2 - a
1287              cfield(i,k,j) = NGRABIKJ(nfield,ni,k,nj)
1288           enddo icopyloop
1289        enddo kcopyloop
1290
1291        kinterploop: do k=kpres+1,nkds,-1
1292           iloop: do i=istart,iend
1293              ni = (i-ipos)*nri + 2 - a
1294              weight=winfo(i,j,k)
1295              nk=iinfo(i,j,k)
1296
1297              if(nk>1) then
1298                 cfield(i,k,j) = &
1299                      NGRABIKJ(nfield,ni,nk,nj) * weight + &
1300 ! pjj/cray - source line limit in Cray compiler
1301              NGRABIKJ(nfield,ni,nk-1,nj) &
1302              * (1.0-weight)
1303              else
1304                 ! Extrapolate
1305                 if(emethod==EConst) then
1306                    cfield(i,k,j)=evalue
1307                 elseif(emethod==ECopy) then
1308                    cfield(i,k,j)=NGRABIKJ(nfield,ni,1,nj)
1309                 else ! Assume 2: linear extrap
1310                    cfield(i,k,j)=evalue*weight + &
1311                         NGRABIKJ(nfield,ni,1,nj)*(1.0-weight)
1312                 endif
1313              endif
1314           enddo iloop
1315        end do kinterploop
1316     enddo bigj
1317   end subroutine n2c_massikj
1318
1319   subroutine c2n_mass (II,JJ,W1,W2,W3,W4,    &
1320        cfield,nfield,iinfo,winfo,imask,      &
1321        emethod, evalue,                      &
1322        cims, cime, cjms, cjme, ckms, ckme,   &
1323        nids, nide, njds, njde, nkds, nkde,   &
1324        nims, nime, njms, njme, nkms, nkme,   &
1325        nits, nite, njts, njte, nkts, nkte)
1326     use module_interp_store, only: kpres
1327     implicit none
1328     real, intent(in), dimension(nims:nime,njms:njme) :: &
1329          W1,W2,W3,W4
1330     integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
1331
1332     integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme,   &
1333          nids, nide, njds, njde, nkds, nkde,   &
1334          nims, nime, njms, njme, nkms, nkme,   &
1335          nits, nite, njts, njte, nkts, nkte, emethod
1336     real, intent(in) :: evalue
1337     real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: winfo
1338     integer, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: iinfo
1339     integer, intent(in), dimension(nims:nime,njms:njme) :: imask
1340     real, intent(in) :: cfield(cims:cime,cjms:cjme,ckms:ckme)
1341     real, intent(out) :: nfield(nims:nime,njms:njme,nkms:nkme)
1342     character*255 :: message
1343     integer :: j,i,a,nx,nz,ck,k
1344     real :: weight
1345
1346     nx=min(nide-1,nite)-max(nids,nits)+1
1347     nz=nkde-nkds+1
1348     if(kpres<=nkds .or. kpres>=nkde) then
1349        call wrf_error_fatal('invalid kpres: outside domain bounds')
1350     end if
1351     bigj: do j=max(njds,njts),min(njde-1,njte)
1352        interploop: do i=max(nids,nits),min(nide-1,nite)
1353           if(imask(i,j)/=0) cycle interploop
1354           a=1-mod(JJ(i,j),2)
1355           kcopyloop: do k=nkde,kpres+1,-1
1356              nfield(i,j,k)=ICOPY(cfield,i,j,k)
1357              !weight=winfo(i,j,k)
1358              ck=iinfo(i,j,k)
1359 !              if(1.-weight>1e-5) then
1360 ! 2000            format(I0,", ",I0,", ",I0,": invalid weight at constant pressure level: w=",F0.5," i=",I0)
1361 !                 write(message,2000) i,j,k, weight,ck
1362 !                 call wrf_error_fatal(message)
1363 !              endif
1364 !              if(ck/=k) then
1365 ! 2001            format(I0,", ",I0,", ",I0,": invalid iinfo at constant pressure level: w=",F0.5," i=",I0)
1366 !                 write(message,2000) i,j,k, weight,ck
1367 !                 call wrf_error_fatal(message)
1368 !              endif
1369           enddo kcopyloop
1370           kinterploop: do k=kpres,nkds,-1
1371              weight=winfo(i,j,k)
1372              ck=iinfo(i,j,k)
1373
1374              if(ck>1) then
1375                 nfield(i,j,k) = &
1376                      ICOPY(cfield,i,j,ck)   * weight + &
1377                      ICOPY(cfield,i,j,ck-1) * (1.0-weight)
1378              else
1379                 ! Extrapolate
1380                 if(emethod==EConst) then
1381                    nfield(i,j,k)=evalue
1382                 elseif(emethod==ECopy) then
1383                    nfield(i,j,k)=ICOPY(cfield,i,j,1)
1384                 else ! Assume 2: linear extrap
1385                    nfield(i,j,k)=evalue*weight + &
1386                         ICOPY(cfield,i,j,1)*(1.0-weight)
1387                 endif
1388              endif
1389           enddo kinterploop
1390        enddo interploop
1391     end do bigj
1392
1393   end subroutine c2n_mass
1394
1395   subroutine c2n_massikj (II,JJ,W1,W2,W3,W4, &
1396        cfield,nfield,iinfo,winfo,imask,      &
1397        emethod, evalue,                      &
1398        cims, cime, cjms, cjme, ckms, ckme,   &
1399        nids, nide, njds, njde, nkds, nkde,   &
1400        nims, nime, njms, njme, nkms, nkme,   &
1401        nits, nite, njts, njte, nkts, nkte)
1402     use module_interp_store, only: kpres
1403     implicit none
1404     real, intent(in), dimension(nims:nime,njms:njme) :: &
1405          W1,W2,W3,W4
1406     integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
1407
1408     integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme,   &
1409          nids, nide, njds, njde, nkds, nkde,   &
1410          nims, nime, njms, njme, nkms, nkme,   &
1411          nits, nite, njts, njte, nkts, nkte, emethod
1412     real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: winfo
1413     integer, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: iinfo
1414     integer, intent(in), dimension(nims:nime,njms:njme) :: imask
1415     real, intent(in) :: cfield(cims:cime,ckms:ckme,cjms:cjme)
1416     real, intent(out) :: nfield(nims:nime,nkms:nkme,njms:njme)
1417     real, intent(In) :: evalue
1418
1419     integer :: j,i,a,nx,nz,ck,k
1420     real :: weight
1421
1422     nx=min(nide-1,nite)-max(nids,nits)+1
1423     nz=nkde-nkds+1
1424
1425     bigj: do j=max(njds,njts),min(njde-1,njte)
1426        interploop: do i=max(nids,nits),min(nide-1,nite)
1427           if(imask(i,j)/=0) cycle interploop
1428           a=1-mod(JJ(i,j),2)
1429           kcopyloop: do k=nkde,kpres+1,-1
1430              nfield(i,k,j)=IKJCOPY(cfield,i,k,j)
1431           end do kcopyloop
1432           kinterploop: do k=kpres,nkds,-1
1433              weight=winfo(i,j,k)
1434              ck=iinfo(i,j,k)
1435
1436              if(ck>1) then
1437                 nfield(i,k,j) = &
1438                      IKJCOPY(cfield,i,ck,j)   * weight + &
1439                      IKJCOPY(cfield,i,ck-1,j) * (1.0-weight)
1440              else
1441                 ! Extrapolate
1442                 if(emethod==EConst) then
1443                    nfield(i,k,j)=evalue
1444                 elseif(emethod==ECopy) then
1445                    nfield(i,k,j)=IKJCOPY(cfield,i,1,j)
1446                 else ! Assume 2: linear extrapolation
1447                    nfield(i,k,j)=evalue * weight + &
1448                                IKJCOPY(cfield,i,1,j) * (1.0-weight)
1449                 endif
1450              endif
1451           enddo kinterploop
1452        enddo interploop
1453     end do bigj
1454
1455   end subroutine c2n_massikj
1456
1457   subroutine find_kpres(kpres, eta2, kds,kde, kms,kme)
1458     ! Find the level at which the pressure/sigma transition occurs.
1459     ! All levels above this are pure pressure levels, whereas kpres
1460     ! and below are sigma.
1461     integer, intent(out) :: kpres
1462     real, intent(in) :: eta2(kms:kme)
1463     integer, intent(in) :: kds,kde, kms,kme
1464
1465     integer :: k
1466
1467     k=kde-1
1468     do while(eta2(k) < 1e-5 .and. eta2(k-1) < 1e-5)
1469        k=k-1
1470        if(k<=kds) then
1471           call wrf_error_fatal('New NMM interpolation routines do not work in a constant pressure space.')
1472        endif
1473     enddo
1474     kpres=k
1475   end subroutine find_kpres
1476
1477   ! ********************************************************************
1478   ! subs *_FULLDOM -- recalculates PD and PINT based on FIS if
1479   ! requested.  Also, generates vertical interpolation arrays for use
1480   ! in later interpolations and interpolates T and Q while doing so.
1481   ! ********************************************************************
1482
1483   subroutine n2c_fulldom  (                  &
1484        deta1,deta2, eta1,eta2, ptop,pdtop,   &
1485        cfis,cpint,ct,cpd,cq,                 &
1486        cids, cide, cjds, cjde, ckds, ckde,   &
1487        cims, cime, cjms, cjme, ckms, ckme,   &
1488        cits, cite, cjts, cjte, ckts, ckte,   &
1489        nfis,npint,nt,npd,nq,                 &
1490        ipos,jpos,                            &
1491 !       nri,nrj                               &
1492        out_iinfo,out_winfo,                  &
1493        nids, nide, njds, njde, nkds, nkde,   &
1494        nims, nime, njms, njme, nkms, nkme,   &
1495        nits, nite, njts, njte, nkts, nkte, kpres)
1496
1497     implicit none
1498
1499     integer, intent(in) :: ipos,jpos
1500 !    integer, intent(in) :: nri,nrj
1501     integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
1502
1503     integer, intent(in):: &
1504          nids, nide, njds, njde, nkds, nkde,   &
1505          nims, nime, njms, njme, nkms, nkme,   &
1506          nits, nite, njts, njte, nkts, nkte,   kpres
1507     integer, intent(in):: &
1508          cids, cide, cjds, cjde, ckds, ckde,   &
1509          cims, cime, cjms, cjme, ckms, ckme,   &
1510          cits, cite, cjts, cjte, ckts, ckte
1511     real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
1512     real, intent(inout), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
1513
1514     real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_winfo
1515     integer, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_iinfo
1516
1517     real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
1518     real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
1519     real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
1520     real, intent(in), dimension(nims:nime,njms:njme,1) :: nPD
1521
1522     real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
1523     real, intent(in) :: ptop,pdtop
1524
1525     real, dimension(1,nite-nits+1) :: inFIS,inPD,icFIS,icPD
1526     real, dimension(nkde-1,nite-nits+1) :: inT0,inQ,icT,icQ, qinfo,winfo
1527     real, dimension(nkde,nite-nits+1) :: inPINT,icPINT
1528     integer, dimension(nkde-1,nite-nits+1) :: iinfo
1529     integer :: nx,nz,k,i,a,j, istart,iend,jstart,jend, ni,nj,jprint,itest,jtest
1530     character*255 :: message
1531     logical bad, warned
1532
1533     warned=.false. ! allow one P>PSTD message
1534
1535     nx=min(cide-2,cite)-max(cids+1,cits)+1
1536     nz=ckde-ckds+1
1537
1538     jstart=MAX(jpos+1,cjts)
1539     jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
1540
1541     bigj: do j=jstart,jend
1542        nj = (j-jpos)*nrj + 1
1543
1544        a=mod(j,2)
1545        istart=MAX(ipos+a,cits)
1546        iend=MIN(ipos+(nide-nids)/nri-1,cite)
1547        nx=iend-istart+1
1548
1549        ! STEP 1: Copy coarse and fine nest data into
1550        ! temporary arrays, reordering dimensions:
1551        qtloop: do k=ckts,ckte-1
1552           do i=istart,iend
1553              ni = (i-ipos)*nri + 2 - a
1554
1555              UPCOPY(inT0,nT,i,j,k,ni,nj)
1556              UPCOPY(inQ,nQ,i,j,k,ni,nj)
1557              UPCOPY(inPINT,nPINT,i,j,k,ni,nj)
1558           enddo
1559        enddo qtloop
1560
1561        k=nkte
1562        loop2d: do i=istart,iend
1563           ni = (i-ipos)*nri + 2 - a
1564
1565           UPCOPY(inPINT,nPINT,i,j,k,ni,nj)
1566           UPCOPY(inPD,nPD,i,j,1,ni,nj)
1567           UPCOPY(inFIS,nFIS,i,j,1,ni,nj)
1568
1569           icPD(1,i-istart+1)=cPD(i,j,1)
1570 !          icPD(1,i-istart+1)=use_this_pd(i,j,1)
1571           icFIS(1,i-istart+1)=cFIS(i,j,1)
1572        enddo loop2d
1573
1574        ! Step 2: Interpolate coarse grid to fine grid in reordered
1575        ! arrays:
1576        call interp_T_PD_Q(nmm_method_linear, nmm_keep_pd, nx,nz, &
1577             deta1,deta2,eta1,eta2,ptop,pdtop, kpres, &
1578             inFIS,icFIS, inPINT,icPINT, inT0, icT, inPD,icPD, inQ,icQ, &
1579             iinfo, winfo, warned)
1580
1581        ! Step 3: Copy back from reordered arrays to final nest arrays:
1582
1583        qtloop2: do k=ckts,max(ckte-1,ckde-1)
1584           do i=istart,iend
1585              cT(i,j,k)=icT(k,i-istart+1)
1586              cQ(i,j,k)=icQ(k,i-istart+1)
1587           enddo
1588        enddo qtloop2
1589
1590        izloop: do k=ckts,max(ckte-1,ckde-1)
1591           ixloop: do i=istart,iend
1592              out_iinfo(i,j,k)=iinfo(k,i-istart+1)
1593              out_winfo(i,j,k)=winfo(k,i-istart+1)
1594           enddo ixloop
1595        enddo izloop
1596
1597        iendloop: do i=istart,iend
1598           out_iinfo(i,j,nkde)=nkde
1599           out_winfo(i,j,nkde)=1.0
1600        end do iendloop
1601
1602        k=nkte+1
1603        loop2d2: do i=istart,iend
1604           cPD(i,j,1)=icPD(1,i-istart+1)
1605        enddo loop2d2
1606     end do bigj
1607   end subroutine n2c_fulldom
1608
1609   subroutine n2c_fulldom_new (               &
1610        deta1,deta2, eta1,eta2, ptop,pdtop,   &
1611        cfis,cpint,ct,cpd,cq,                 &
1612        cids, cide, cjds, cjde, ckds, ckde,   &
1613        cims, cime, cjms, cjme, ckms, ckme,   &
1614        cits, cite, cjts, cjte, ckts, ckte,   &
1615        nfis,npint,nt,npd,nq,                 &
1616        ipos,jpos,                            &
1617 !       nri,nrj                               &
1618        out_iinfo,out_winfo,                  &
1619        nids, nide, njds, njde, nkds, nkde,   &
1620        nims, nime, njms, njme, nkms, nkme,   &
1621        nits, nite, njts, njte, nkts, nkte, kpres)
1622
1623     implicit none
1624
1625     integer, intent(in) :: ipos,jpos
1626 !    integer, intent(in) :: nri,nrj
1627     integer, parameter :: nri=3, nrj=3 ! parent:nest ratio, must be 3 in NMM
1628
1629     integer, intent(in):: &
1630          nids, nide, njds, njde, nkds, nkde,   &
1631          nims, nime, njms, njme, nkms, nkme,   &
1632          nits, nite, njts, njte, nkts, nkte,   kpres
1633     integer, intent(in):: &
1634          cids, cide, cjds, cjde, ckds, ckde,   &
1635          cims, cime, cjms, cjme, ckms, ckme,   &
1636          cits, cite, cjts, cjte, ckts, ckte
1637     real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
1638     real, intent(inout), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
1639
1640     real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_winfo
1641     integer, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_iinfo
1642
1643     real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
1644     real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
1645     real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
1646     real, intent(in), dimension(nims:nime,njms:njme,1) :: nPD
1647
1648     real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
1649     real, intent(in) :: ptop,pdtop
1650
1651     real, dimension(1,nite-nits+1) :: inFIS,inPD,icFIS,icPD
1652     real, dimension(kpres+1,nite-nits+1) :: inT0,inQ,icT,icQ, qinfo,winfo
1653     real, dimension(kpres+2,nite-nits+1) :: inPINT,icPINT
1654     integer, dimension(kpres+1,nite-nits+1) :: iinfo
1655     integer :: nx,nz,k,i,a,j, istart,iend,jstart,jend, ni,nj,jprint,itest,jtest
1656     character*255 :: message
1657     logical bad, warned
1658
1659     warned=.false. ! Allow one P>PSTD message
1660
1661     nx=min(cide-2,cite)-max(cids+1,cits)+1
1662     nz=ckde-ckds+1
1663
1664     jstart=MAX(jpos+1,cjts)
1665     jend=MIN(jpos+(njde-njds)/nrj-1,cjte)
1666
1667     bigj: do j=jstart,jend
1668        nj = (j-jpos)*nrj + 1
1669
1670        a=mod(j,2)
1671        istart=MAX(ipos+a,cits)
1672        iend=MIN(ipos+(nide-nids)/nri-1,cite)
1673        nx=iend-istart+1
1674
1675        ! STEP 1: Copy coarse and fine nest data into
1676        ! temporary arrays, reordering dimensions:
1677        qtloop: do k=ckts,kpres+1
1678           do i=istart,iend
1679              ni = (i-ipos)*nri + 2 - a
1680
1681              UPCOPY(inT0,nT,i,j,k,ni,nj)
1682              UPCOPY(inQ,nQ,i,j,k,ni,nj)
1683              UPCOPY(inPINT,nPINT,i,j,k,ni,nj)
1684           enddo
1685        enddo qtloop
1686
1687        k=kpres+1
1688        loop2d: do i=istart,iend
1689           ni = (i-ipos)*nri + 2 - a
1690
1691           UPCOPY(inPINT,nPINT,i,j,k,ni,nj)
1692           UPCOPY(inPD,nPD,i,j,1,ni,nj)
1693           UPCOPY(inFIS,nFIS,i,j,1,ni,nj)
1694
1695           icPD(1,i-istart+1)=cPD(i,j,1)
1696 !          icPD(1,i-istart+1)=use_this_pd(i,j,1)
1697           icFIS(1,i-istart+1)=cFIS(i,j,1)
1698        enddo loop2d
1699
1700        ! Step 2: Interpolate coarse grid to fine grid in reordered
1701        ! arrays:
1702        call interp_T_PD_Q_kpres(nmm_method_linear, nmm_keep_pd, nx,nz, &
1703             deta1,deta2,eta1,eta2,ptop,pdtop, kpres, kpres+2, &
1704             inFIS,icFIS, inPINT,icPINT, inT0, icT, inPD,icPD, inQ,icQ, &
1705             iinfo, winfo, warned)
1706
1707        ! Step 3: Copy back from reordered arrays to final nest arrays:
1708
1709        qtloop2: do k=ckts,kpres+1
1710           do i=istart,iend
1711              cT(i,j,k)=icT(k,i-istart+1)
1712              cQ(i,j,k)=icQ(k,i-istart+1)
1713           enddo
1714        enddo qtloop2
1715
1716        izloop: do k=ckts,kpres+1
1717           ixloop: do i=istart,iend
1718              out_iinfo(i,j,k)=iinfo(k,i-istart+1)
1719              out_winfo(i,j,k)=winfo(k,i-istart+1)
1720           enddo ixloop
1721        enddo izloop
1722
1723        k=nkte+1
1724        loop2d2: do i=istart,iend
1725           cPD(i,j,1)=icPD(1,i-istart+1)
1726        enddo loop2d2
1727     end do bigj
1728
1729     kcopy: do k=kpres+1,ckde-1
1730        jcopy: do j=jstart,jend
1731           nj = (j-jpos)*nrj + 1
1732
1733           a=mod(j,2)
1734           istart=MAX(ipos+a,cits)
1735           iend=MIN(ipos+(nide-nids)/nri-1,cite)
1736           icopy: do i=istart,iend
1737              ni = (i-ipos)*nri + 2 - a
1738              cT(i,j,k)=NGRAB(nT,ni,nj,k)
1739              cQ(i,j,k)=NGRAB(nQ,ni,nj,k)
1740              out_iinfo(i,j,k)=k
1741              out_winfo(i,j,k)=1.0
1742           enddo icopy
1743           if(k==ckde-1) then
1744              icopy2: do i=istart,iend
1745                 out_iinfo(i,j,nkde)=nkde
1746                 out_winfo(i,j,nkde)=1.0
1747              enddo icopy2
1748           endif
1749        end do jcopy
1750     end do kcopy
1751   end subroutine n2c_fulldom_new
1752
1753   subroutine c2b_fulldom  (II,JJ,W1,W2,W3,W4,&
1754        deta1,deta2, eta1,eta2, ptop,pdtop,   &
1755        cfis,cpint,ct,cpd,cq,    nfis,        &
1756        cims, cime, cjms, cjme, ckms, ckme,   &
1757        nids, nide, njds, njde, nkds, nkde,   &
1758        nims, nime, njms, njme, nkms, nkme,   &
1759        nits, nite, njts, njte, nkts, nkte,   &
1760        kpres,                                &
1761        ibxs,    ibxe,    ibys,    ibye,      &
1762        wbxs,    wbxe,    wbys,    wbye,      &
1763        pdbxs,   pdbxe,   pdbys,   pdbye,     &
1764        tbxs,    tbxe,    tbys,    tbye,      &
1765        qbxs,    qbxe,    qbys,    qbye)
1766     implicit none
1767     integer, intent(in):: &
1768          cims, cime, cjms, cjme, ckms, ckme,   &
1769          nids, nide, njds, njde, nkds, nkde,   &
1770          nims, nime, njms, njme, nkms, nkme,   &
1771          nits, nite, njts, njte, nkts, nkte
1772
1773     integer, parameter :: bdyw=1
1774
1775     ! Domain info:
1776     real, intent(in), dimension(nims:nime,njms:njme) :: W1,W2,W3,W4
1777     integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
1778     real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
1779     real, intent(in) :: ptop,pdtop
1780     integer, intent(in) :: kpres
1781
1782     ! Parent fields:
1783     real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
1784     real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
1785
1786     ! Nest terrain info:
1787     real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
1788
1789     ! T, Q, PINT, PD boundary info:
1790     real,dimension(nims:nime,nkms:nkme,bdyw) :: tbys,tbye,qbys,qbye
1791     real,dimension(njms:njme,nkms:nkme,bdyw) :: tbxs,tbxe,qbxs,qbxe
1792     real,dimension(nims:nime,1,bdyw) :: pdbys,pdbye
1793     real,dimension(njms:njme,1,bdyw) :: pdbxs,pdbxe
1794
1795     ! Weights and indices:
1796     real,dimension(nims:nime,nkms:nkme,bdyw) :: wbys,wbye
1797     real,dimension(njms:njme,nkms:nkme,bdyw) :: wbxs,wbxe
1798     integer,dimension(nims:nime,nkms:nkme,bdyw) :: ibys,ibye
1799     integer,dimension(njms:njme,nkms:nkme,bdyw) :: ibxs,ibxe
1800
1801     integer :: i,j,k,a,b,nx,nz,used,j1,j2,used1
1802
1803     real, dimension(1,2*(nite-nits+5)+2*(njte-njts+5)) :: inFIS,inPD,icFIS,icPD
1804     real, dimension(nkde-1,2*(nite-nits+5)+2*(njte-njts+5)) :: inT,inQ,icT,icQ,winfo
1805     integer, dimension(nkde-1,2*(nite-nits+5)+2*(njte-njts+5)) :: iinfo
1806     real, dimension(nkde,2*(nite-nits+5)+2*(njte-njts+5)) :: inPINT,icPINT
1807     logical :: warned
1808
1809     warned=.false. ! Allow one P>PSTD message
1810
1811     nx=min(nide-1,nite)-max(nids,nits)+1
1812     nz=nkde-nkds+1
1813
1814     j1=max(njts-1,njds)
1815     if(mod(j1,2)/=1)   j1=j1+1
1816
1817     j2=min(njte+1,njde-1)
1818     if(mod(j2,2)/=1)   j2=j2-1
1819
1820     used=0
1821     bdyloop: do b=1,4
1822        if_xbdy: if(b==1 .or. b==2) then
1823           if(b==1) then
1824              if(nits/=1) cycle bdyloop
1825              i=1
1826           endif
1827           if(b==2) then
1828              if(nite<nide-1) cycle bdyloop
1829              i=nide-1
1830           endif
1831           do j=j1,j2,2
1832              a=1-mod(JJ(i,j),2)
1833              used=used+1
1834              do k=nkts,nkte-1
1835                 uFCOPY(icT,cT,i,j,k)
1836                 uFCOPY(icQ,cQ,i,j,k)
1837                 uFCOPY(icPINT,cPINT,i,j,k)
1838              enddo
1839
1840              k=nkte
1841              a=1-mod(JJ(i,j),2)
1842
1843              uFCOPY(icPINT,cPINT,i,j,k)
1844              uFCOPY(icPD,cPD,i,j,1)
1845              uFCOPY(icFIS,cFIS,i,j,1)
1846
1847              inFIS(1,used)=nFIS(i,j,1)
1848
1849           enddo
1850        endif if_xbdy
1851        if_ybdy: if(b==3 .or. b==4) then
1852           if(b==3) then
1853              if(njts/=1) cycle bdyloop
1854              j=1
1855           endif
1856           if(b==4) then
1857              if(njte<njde-1) cycle bdyloop
1858              j=njde-1
1859           endif
1860           do i=max(nits-1,nids),min(nite+1,nide-1)
1861              used=used+1
1862              do k=nkts,nkte-1
1863                 a=1-mod(JJ(i,j),2)
1864                 uFCOPY(icT,cT,i,j,k)
1865                 uFCOPY(icQ,cQ,i,j,k)
1866                 uFCOPY(icPINT,cPINT,i,j,k)
1867              enddo
1868
1869              k=nkte
1870              a=1-mod(JJ(i,j),2)
1871
1872              uFCOPY(icPINT,cPINT,i,j,k)
1873              uFCOPY(icPD,cPD,i,j,1)
1874              uFCOPY(icFIS,cFIS,i,j,1)
1875
1876              inFIS(1,used)=nFIS(i,j,1)
1877           enddo
1878        endif if_ybdy
1879     enddo bdyloop
1880
1881     if(used==0) then
1882        ! No boundary points on this tile so we are done
1883        return
1884     endif
1885
1886     call interp_T_PD_Q_kpres(nmm_method_linear, nmm_interp_pd, used,nz, &
1887          deta1,deta2,eta1,eta2,ptop,pdtop, kpres, nz, &
1888          icFIS,inFIS, icPINT,inPINT, icT, inT, icPD,inPD, icQ,inQ, &
1889          iinfo, winfo, warned)
1890
1891     used1=used
1892     used=0
1893
1894     if_bxs: if(nits==1) then
1895        i=1
1896        do j=j1,j2,2
1897           used=used+1
1898           do k=nkts,nkte-1
1899              tbxs(j,k,1)=inT(k,used)
1900              qbxs(j,k,1)=inQ(k,used)
1901              ibxs(j,k,1)=iinfo(k,used)
1902              wbxs(j,k,1)=winfo(k,used)
1903           enddo
1904           ibxs(j,nkde,1)=nkde
1905           wbxs(j,nkde,1)=1.0
1906           pdbxs(j,1,1)=inPD(1,used)
1907        enddo
1908     endif if_bxs
1909
1910     if_bxe: if(nite>=nide-1) then
1911        i=nide-1
1912        do j=j1,j2,2
1913           used=used+1
1914           do k=nkts,nkte-1
1915              tbxe(j,k,1)=inT(k,used)
1916              qbxe(j,k,1)=inQ(k,used)
1917              ibxe(j,k,1)=iinfo(k,used)
1918              wbxe(j,k,1)=winfo(k,used)
1919           enddo
1920           ibxe(j,nkde,1)=nkde
1921           wbxe(j,nkde,1)=1.0
1922           pdbxe(j,1,1)=inPD(1,used)
1923        enddo
1924     endif if_bxe
1925
1926     if_bys: if(njts==1) then
1927        j=1
1928        do i=max(nits-1,nids),min(nite+1,nide-1)
1929           used=used+1
1930           do k=nkts,nkte-1
1931              tbys(i,k,1)=inT(k,used)
1932              qbys(i,k,1)=inQ(k,used)
1933              ibys(i,k,1)=iinfo(k,used)
1934              wbys(i,k,1)=winfo(k,used)
1935           enddo
1936           ibys(i,nkde,1)=nkde
1937           wbys(i,nkde,1)=1.0
1938           pdbys(i,1,1)=inPD(1,used)
1939        enddo
1940     endif if_bys
1941
1942     if_bye: if(njte>=njde-1) then
1943        j=njde-1
1944        do i=max(nits-1,nids),min(nite+1,nide-1)
1945           used=used+1
1946           do k=nkts,nkte-1
1947              tbye(i,k,1)=inT(k,used)
1948              qbye(i,k,1)=inQ(k,used)
1949              ibye(i,k,1)=iinfo(k,used)
1950              wbye(i,k,1)=winfo(k,used)
1951           enddo
1952           ibye(i,nkde,1)=nkde
1953           wbye(i,nkde,1)=1.0
1954           pdbye(i,1,1)=inPD(1,used)
1955        enddo
1956     endif if_bye
1957
1958     if(used/=used1) then
1959        call wrf_error_fatal('Number of input and output points does not match.')
1960     endif
1961
1962   end subroutine c2b_fulldom
1963
1964   ! ####################################################################
1965
1966   subroutine c2b_fulldom_new  (II,JJ,W1,W2,W3,W4,&
1967        deta1,deta2, eta1,eta2, ptop,pdtop,   &
1968        cfis,cpint,ct,cpd,cq,    nfis,        &
1969        cims, cime, cjms, cjme, ckms, ckme,   &
1970        nids, nide, njds, njde, nkds, nkde,   &
1971        nims, nime, njms, njme, nkms, nkme,   &
1972        nits, nite, njts, njte, nkts, nkte,   &
1973        kpres,                                &
1974        ibxs,    ibxe,    ibys,    ibye,      &
1975        wbxs,    wbxe,    wbys,    wbye,      &
1976        pdbxs,   pdbxe,   pdbys,   pdbye,     &
1977        tbxs,    tbxe,    tbys,    tbye,      &
1978        qbxs,    qbxe,    qbys,    qbye)
1979     implicit none
1980     integer, intent(in):: &
1981          cims, cime, cjms, cjme, ckms, ckme,   &
1982          nids, nide, njds, njde, nkds, nkde,   &
1983          nims, nime, njms, njme, nkms, nkme,   &
1984          nits, nite, njts, njte, nkts, nkte
1985
1986     integer, parameter :: bdyw=1
1987
1988     ! Domain info:
1989     real, intent(in), dimension(nims:nime,njms:njme) :: W1,W2,W3,W4
1990     integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
1991     real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
1992     real, intent(in) :: ptop,pdtop
1993     integer, intent(in) :: kpres
1994
1995     ! Parent fields:
1996     real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
1997     real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
1998
1999     ! Nest terrain info:
2000     real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
2001
2002     ! T, Q, PINT, PD boundary info:
2003     real,dimension(nims:nime,nkms:nkme,bdyw) :: tbys,tbye,qbys,qbye
2004     real,dimension(njms:njme,nkms:nkme,bdyw) :: tbxs,tbxe,qbxs,qbxe
2005     real,dimension(nims:nime,1,bdyw) :: pdbys,pdbye
2006     real,dimension(njms:njme,1,bdyw) :: pdbxs,pdbxe
2007
2008     ! Weights and indices:
2009     real,dimension(nims:nime,nkms:nkme,bdyw) :: wbys,wbye
2010     real,dimension(njms:njme,nkms:nkme,bdyw) :: wbxs,wbxe
2011     integer,dimension(nims:nime,nkms:nkme,bdyw) :: ibys,ibye
2012     integer,dimension(njms:njme,nkms:nkme,bdyw) :: ibxs,ibxe
2013
2014     integer :: i,j,k,a,b,nx,nz,used,j1,j2,used1
2015
2016     real, dimension(1,2*(nite-nits+5)+2*(njte-njts+5)) :: inFIS,inPD,icFIS,icPD
2017     real, dimension(kpres+1,2*(nite-nits+5)+2*(njte-njts+5)) :: inT,inQ,icT,icQ,winfo
2018     integer, dimension(kpres+1,2*(nite-nits+5)+2*(njte-njts+5)) :: iinfo
2019     real, dimension(kpres+2,2*(nite-nits+5)+2*(njte-njts+5)) :: inPINT,icPINT
2020     logical :: warned
2021
2022     warned=.false. ! Allow one P>PSTD message
2023
2024     nx=min(nide-1,nite)-max(nids,nits)+1
2025     nz=nkde-nkds+1
2026
2027     j1=max(njts-1,njds)
2028     if(mod(j1,2)/=1)   j1=j1+1
2029
2030     j2=min(njte+1,njde-1)
2031     if(mod(j2,2)/=1)   j2=j2-1
2032
2033     used=0
2034     bdyloop: do b=1,4
2035        if_xbdy: if(b==1 .or. b==2) then
2036           if(b==1) then
2037              if(nits/=1) cycle bdyloop
2038              i=1
2039           endif
2040           if(b==2) then
2041              if(nite<nide-1) cycle bdyloop
2042              i=nide-1
2043           endif
2044           do j=j1,j2,2
2045              used=used+1
2046              do k=nkts,kpres+1
2047                 a=1-mod(JJ(i,j),2)
2048                 uFCOPY(icT,cT,i,j,k)
2049                 uFCOPY(icQ,cQ,i,j,k)
2050                 uFCOPY(icPINT,cPINT,i,j,k)
2051              enddo
2052
2053              k=kpres+2
2054              a=1-mod(JJ(i,j),2)
2055
2056              uFCOPY(icPINT,cPINT,i,j,k)
2057              uFCOPY(icPD,cPD,i,j,1)
2058              uFCOPY(icFIS,cFIS,i,j,1)
2059
2060              inFIS(1,used)=nFIS(i,j,1)
2061
2062           enddo
2063        endif if_xbdy
2064        if_ybdy: if(b==3 .or. b==4) then
2065           if(b==3) then
2066              if(njts/=1) cycle bdyloop
2067              j=1
2068           endif
2069           if(b==4) then
2070              if(njte<njde-1) cycle bdyloop
2071              j=njde-1
2072           endif
2073           do i=max(nits-1,nids),min(nite+1,nide-1)
2074              used=used+1
2075              do k=nkts,kpres+1
2076                 a=1-mod(JJ(i,j),2)
2077                 uFCOPY(icT,cT,i,j,k)
2078                 uFCOPY(icQ,cQ,i,j,k)
2079                 uFCOPY(icPINT,cPINT,i,j,k)
2080              enddo
2081
2082              k=kpres+2
2083              a=1-mod(JJ(i,j),2)
2084
2085              uFCOPY(icPINT,cPINT,i,j,k)
2086              uFCOPY(icPD,cPD,i,j,1)
2087              uFCOPY(icFIS,cFIS,i,j,1)
2088
2089              inFIS(1,used)=nFIS(i,j,1)
2090           enddo
2091        endif if_ybdy
2092     enddo bdyloop
2093
2094     if(used==0) then
2095        ! No boundary points on this tile so we are done
2096        return
2097     endif
2098
2099     call interp_T_PD_Q_kpres(nmm_method_linear, nmm_interp_pd, used,nz, &
2100          deta1,deta2,eta1,eta2,ptop,pdtop, kpres, kpres+2, &
2101          icFIS,inFIS, icPINT,inPINT, icT, inT, icPD,inPD, icQ,inQ, &
2102          iinfo, winfo, warned)
2103
2104     used1=used
2105     used=0
2106
2107     if_bxs: if(nits==1) then
2108        i=1
2109        do j=j1,j2,2
2110           used=used+1
2111           do k=nkts,kpres
2112              tbxs(j,k,1)=inT(k,used)
2113              qbxs(j,k,1)=inQ(k,used)
2114              ibxs(j,k,1)=iinfo(k,used)
2115              wbxs(j,k,1)=winfo(k,used)
2116           enddo
2117           ibxs(j,nkde,1)=nkde
2118           wbxs(j,nkde,1)=1.0
2119           pdbxs(j,1,1)=inPD(1,used)
2120        enddo
2121        do k=kpres+1,nkde-1
2122           do j=j1,j2,2
2123              a=1-mod(JJ(i,j),2)
2124              tbxs(j,k,1)=ICOPY(cT,i,j,k)
2125              qbxs(j,k,1)=ICOPY(cQ,i,j,k)
2126              ibxs(j,k,1)=k
2127              wbxs(j,k,1)=1.0
2128           enddo
2129        enddo
2130     endif if_bxs
2131
2132     if_bxe: if(nite>=nide-1) then
2133        i=nide-1
2134        do j=j1,j2,2
2135           used=used+1
2136           do k=nkts,kpres
2137              tbxe(j,k,1)=inT(k,used)
2138              qbxe(j,k,1)=inQ(k,used)
2139              ibxe(j,k,1)=iinfo(k,used)
2140              wbxe(j,k,1)=winfo(k,used)
2141           enddo
2142           ibxe(j,nkde,1)=nkde
2143           wbxe(j,nkde,1)=1.0
2144           pdbxe(j,1,1)=inPD(1,used)
2145        enddo
2146        do k=kpres+1,nkde-1
2147           do j=j1,j2,2
2148              a=1-mod(JJ(i,j),2)
2149              tbxe(j,k,1)=ICOPY(cT,i,j,k)
2150              qbxe(j,k,1)=ICOPY(cQ,i,j,k)
2151              ibxe(j,k,1)=k
2152              wbxe(j,k,1)=1.0
2153           enddo
2154        enddo
2155     endif if_bxe
2156
2157     if_bys: if(njts==1) then
2158        j=1
2159        do i=max(nits-1,nids),min(nite+1,nide-1)
2160           used=used+1
2161           do k=nkts,kpres
2162              tbys(i,k,1)=inT(k,used)
2163              qbys(i,k,1)=inQ(k,used)
2164              ibys(i,k,1)=iinfo(k,used)
2165              wbys(i,k,1)=winfo(k,used)
2166           enddo
2167           ibys(i,nkde,1)=nkde
2168           wbys(i,nkde,1)=1.0
2169           pdbys(i,1,1)=inPD(1,used)
2170        enddo
2171        do k=kpres+1,nkde-1
2172           do i=max(nits-1,nids),min(nite+1,nide-1)
2173              a=1-mod(JJ(i,j),2)
2174              tbys(i,k,1)=ICOPY(cT,i,j,k)
2175              qbys(i,k,1)=ICOPY(cQ,i,j,k)
2176              ibys(i,k,1)=k
2177              wbys(i,k,1)=1.0
2178           enddo
2179        enddo
2180     endif if_bys
2181
2182     if_bye: if(njte>=njde-1) then
2183        j=njde-1
2184        do i=max(nits-1,nids),min(nite+1,nide-1)
2185           used=used+1
2186           do k=nkts,kpres
2187              tbye(i,k,1)=inT(k,used)
2188              qbye(i,k,1)=inQ(k,used)
2189              ibye(i,k,1)=iinfo(k,used)
2190              wbye(i,k,1)=winfo(k,used)
2191           enddo
2192           ibye(i,nkde,1)=nkde
2193           wbye(i,nkde,1)=1.0
2194           pdbye(i,1,1)=inPD(1,used)
2195        enddo
2196        do k=kpres+1,nkde-1
2197           do i=max(nits-1,nids),min(nite+1,nide-1)
2198              a=1-mod(JJ(i,j),2)
2199              tbye(i,k,1)=ICOPY(cT,i,j,k)
2200              qbye(i,k,1)=ICOPY(cQ,i,j,k)
2201              ibye(i,k,1)=k
2202              wbye(i,k,1)=1.0
2203           enddo
2204        enddo
2205     endif if_bye
2206
2207     if(used/=used1) then
2208        call wrf_error_fatal('Number of input and output points does not match.')
2209     endif
2210
2211   end subroutine c2b_fulldom_new
2212
2213   ! ####################################################################
2214
2215   subroutine c2n_fulldom  (II,JJ,W1,W2,W3,W4,&
2216        deta1,deta2, eta1,eta2, ptop,pdtop,   &
2217        cfis,cpint,ct,cpd,cq,                 &
2218        cims, cime, cjms, cjme, ckms, ckme,   &
2219        nfis,npint,nt,npd,nq,                 &
2220        out_iinfo,out_winfo,imask,&
2221        nids, nide, njds, njde, nkds, nkde,   &
2222        nims, nime, njms, njme, nkms, nkme,   &
2223        nits, nite, njts, njte, nkts, nkte, kpres)
2224     implicit none
2225     integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme,   &
2226          nids, nide, njds, njde, nkds, nkde,   &
2227          nims, nime, njms, njme, nkms, nkme,   &
2228          nits, nite, njts, njte, nkts, nkte, kpres
2229     real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
2230     real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD,cFIS
2231
2232     real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_winfo
2233     integer, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_iinfo
2234
2235     real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
2236     real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
2237     real, intent(in), dimension(nims:nime,njms:njme,1) :: nFIS
2238     real, intent(inout), dimension(nims:nime,njms:njme,1) :: nPD
2239     integer, intent(in), dimension(nims:nime,njms:njme) :: imask
2240
2241     real, intent(in), dimension(nims:nime,njms:njme) :: &
2242          W1,W2,W3,W4
2243     integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
2244
2245     real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
2246     real, intent(in) :: ptop,pdtop
2247
2248     real, dimension(1,nite-nits+1) :: inFIS,inPD,icFIS,icPD
2249     real, dimension(nkde-1,nite-nits+1) :: inT,inQ,icT,icQ,winfo
2250     integer, dimension(nkde-1,nite-nits+1) :: iinfo
2251     real, dimension(nkde,nite-nits+1) :: inPINT,icPINT
2252
2253     real :: pdcheck
2254     integer :: i,j,k,a, nx,nz,used
2255     logical :: badbad, warned
2256
2257     warned=.false. ! Allow one P>PSTD message
2258
2259     nx=min(nide-1,nite)-max(nids,nits)+1
2260     nz=nkde-nkds+1
2261
2262     bigj: do j=max(njds,njts),min(njde-1,njte)
2263        ! STEP 1: Copy coarse and fine nest data into
2264        ! temporary arrays, reordering dimensions:
2265        used=0
2266        iloop1: do i=max(nids,nits),min(nide-1,nite)
2267           if(imask(i,j)/=0) cycle iloop1
2268           used=used+1
2269
2270           qtloop: do k=nkts,nkte-1
2271              a=1-mod(JJ(i,j),2)
2272              uFCOPY(icT,cT,i,j,k)
2273              uFCOPY(icQ,cQ,i,j,k)
2274              uFCOPY(icPINT,cPINT,i,j,k)
2275           enddo qtloop
2276
2277           k=nkte
2278           a=1-mod(JJ(i,j),2)
2279
2280           uFCOPY(icPINT,cPINT,i,j,k)
2281           uFCOPY(icPD,cPD,i,j,1)
2282           uFCOPY(icFIS,cFIS,i,j,1)
2283
2284 !           nPD(i,j,1)=use_this_pd(i,j,1) ! FIXME: remove this line
2285  !           if(nPD(i,j,1)> 1.3e5 .or. nPD(i,j,1)<2e4) then
2286  ! 3131         format('Invalid input nest PD in C2N interpolation: PD(',I0,',',I0,') = ',F0.7)
2287  !              write(0,3131) i,j,nPD(i,j,1)
2288  !              call wrf_error_fatal('Invalid input nest PD')
2289  !           endif
2290           inPD(1,used)=nPD(i,j,1)
2291           inFIS(1,used)=nFIS(i,j,1)
2292        enddo iloop1
2293
2294        ! Step 2: Interpolate coarse grid to fine grid in reordered
2295        ! arrays:
2296
2297        if(used==0) then
2298           ! No points in this row require interpolation, so we're done
2299           ! with this row:
2300           cycle bigj
2301        else
2302           call interp_T_PD_Q(nmm_method_linear, nmm_interp_pd, used,nz, &
2303                deta1,deta2,eta1,eta2,ptop,pdtop,kpres, &
2304                icFIS,inFIS, icPINT,inPINT, icT, inT, icPD,inPD, icQ,inQ, &
2305                iinfo, winfo, warned)
2306        endif
2307
2308        ! Step 3: Copy back from reordered arrays to final nest arrays:
2309        used=0
2310        iloop2: do i=max(nids,nits),min(nide-1,nite)
2311           if(imask(i,j)/=0) cycle iloop2
2312           used=used+1
2313
2314           qtloop2: do k=nkts,nkte-1
2315              nT(i,j,k)=inT(k,used)
2316              nQ(i,j,k)=inQ(k,used)
2317              nQ(i,j,k)=inQ(k,used)
2318           enddo qtloop2
2319
2320           izloop: do k=nkts,nkte-1
2321              out_iinfo(i,j,k)=iinfo(k,used)
2322              out_winfo(i,j,k)=winfo(k,used)
2323           enddo izloop
2324           out_iinfo(i,j,nkde)=nkde
2325           out_winfo(i,j,nkde)=1.0
2326
2327           k=nkte+1
2328           a=1-mod(JJ(i,j),2)
2329           nPD(i,j,1)=inPD(1,used)
2330           pdcheck=npd(i,j,1)+ptop+pdtop
2331  !          if(pdcheck<50000.0 .or. pdcheck>105000.0) then
2332  ! 3131         format('Invalid output nest PD in C2N interpolation: PD(',I0,',',I0,') = ',F0.7,' which is ',F0.7,' Pa')
2333  !              write(0,3131) i,j,nPD(i,j,1),pdcheck
2334  !              call wrf_error_fatal('Invalid output nest PD')
2335  !          endif
2336        enddo iloop2
2337     enddo bigj
2338
2339   end subroutine c2n_fulldom
2340
2341   ! ********************************************************************
2342   ! INTERP_T_PD_Q -- this is the actual vertical interpolation
2343   ! subroutine, the implementation of the *_fulldom subroutines.  It
2344   ! takes arrays of atmospheric columns, with no knowledge of the
2345   ! horizontal layout of those columns.  This routine will recalculate
2346   ! PD if requested, then generate interpolation information and
2347   ! interpolate T and Q, within each atmospheric column.
2348   !
2349   ! Other variables are handled by the various other c2n, n2c and c2b
2350   ! routines in this module.  PD, T and Q are handled specially since
2351   ! they are critical to the stability of the hydrostatic solver, and
2352   ! hence have more expensive extrapolation and mass rebalancing
2353   ! methods.
2354   !
2355   ! Since this routine has no knowledge of the horizontal layout of
2356   ! the atmospheric columns it is given, the columns do not have to
2357   ! lie on an I- or J-row.  The boundary interpolation routines take
2358   ! advantage of that and provide all boundary information to this
2359   ! routine in a single call.  Of course, the caller must handle all
2360   ! horizontal interpolation.
2361   ! ********************************************************************
2362
2363   subroutine interp_T_PD_Q(method, pd_interp, nx, nz, &
2364        deta1,deta2, eta1,eta2, ptop,pdtop, kpres,     &
2365        fisA,fisB, pintA,pintB, tA,tB, pdA,pdB, qA,qB, &
2366        iinfo, winfo, warned)
2367     implicit none
2368
2369     integer, intent(in) :: pd_interp,method, kpres
2370     !      real, intent(in) :: dtA,dtB
2371
2372     ! Coordinate system definitions must be the same for all domains:
2373     real, intent(in) :: deta1(nz),deta2(nz),eta1(nz),eta2(nz),ptop,pdtop
2374     integer, intent(in) :: nx,nz
2375
2376     ! Surface height and mass field information for source (A):
2377     real, intent(in) :: fisA(nx),tA(nz-1,nx),pdA(nx),qA(nz-1,nx),pintA(nz,nx)
2378
2379     ! Surface height and mass field information for target (B):
2380     real, intent(inout) :: fisB(nx),tB(nz-1,nx),pdB(nx),qB(nz-1,nx),pintB(nz,nx)
2381
2382     ! Interpolation or extrapolation information for use in other
2383     ! calls later on:
2384     real, intent(out) :: winfo(nz-1,nx)
2385     integer, intent(out) :: iinfo(nz-1,nx)
2386     logical, intent(inout) :: warned
2387
2388     ! ==================== Local variables ====================
2389
2390     character*255 :: message
2391     integer :: ix,iz,izA,izB,xpr
2392     real :: zA,zB,znext,apelp,rtopp,dz,weight,A,B,pstd1,z,pstd2,pstd12
2393     real :: tsfc(nx), slp(nx), zmslp(nx), z0mid(nx), zbelow(nx), &
2394             tbelow(nx), RHbelow(nx)
2395     real :: pb,pb1,pb2,pa,pa1,pa2,pnext,pa3, wnum,wdenom, QC, P0
2396
2397     ! Constants from UPP params.f for RH calculation (all mks):
2398     real, parameter :: PQ0=379.90516
2399     real, parameter :: A2=17.2693882
2400     real, parameter :: A3=273.16
2401     real, parameter :: A4=35.86
2402     real, parameter :: RHmin=1.0E-6     ! minimal RH bound
2403
2404     if(method/=nmm_method_linear) then
2405        call wrf_error_fatal('only linear interpolation is supported')
2406     endif
2407
2408     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2409     !! Step 1: calculate near-surface values !!!!!!!!!!!!!!!!!!!!!
2410     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2411
2412     pstd1=p_ref ! pstd(1) from base_state_parent
2413     ! pstd(2) from base_state_parent:
2414     pstd2=eta1(2)*pdtop + eta2(2)*(p_ref-pdtop-ptop) + ptop
2415     pstd12=exp((alog(pstd1)+alog(pstd2))*0.5)
2416
2417     do ix=1,nx
2418        ! These calculations are from base_state_parent:
2419        APELP    = (pintA(2,ix)+pintA(1,ix))
2420        RTOPP    = TRG*tA(1,ix)*(1.0+qA(1,ix)*P608)/APELP
2421        DZ       = RTOPP*(DETA1(1)*PDTOP+DETA2(1)*pdA(ix))
2422
2423        Z0MID(ix) = fisA(ix)/g + dz/2
2424
2425        zA=fisA(ix)/g
2426
2427        TSFC(ix) = TA(1,ix)*(1.+D608*QA(1,ix)) &
2428             + LAPSR*(zA + zA+dz)*0.5
2429
2430        A         = LAPSR*zA/TSFC(ix)
2431        SLP(ix)  = PINTA(1,ix)*(1-A)**COEF2    ! sea level pressure
2432        B         = (pstd1/SLP(ix))**COEF3
2433        ZMSLP(ix)= TSFC(ix)*LAPSI*(1.0 - B)   ! Height at 1030. mb level
2434
2435        TBELOW(ix) = TA(1,ix) + LAPSR*(Z0MID(ix)-ZMSLP(ix))
2436
2437        ! Calculate lowest level RH.  This calculation is from UPP
2438        ! CALRH.f:
2439        P0=pdA(ix)+ptop+pdtop ! Use hydrostatic lowest model level P
2440        QC=PQ0/P0*EXP(A2*(TA(1,ix)-A3)/(TA(1,ix)-A4))
2441        RHbelow(ix)=max(RHmin,min(1.0,QA(1,ix)/QC))
2442     enddo
2443
2444
2445
2446     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2447     !! Step 2: figure out the new surface pressures !!!!!!!!!!!!!!
2448     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2449
2450
2451
2452     if_pd_interp: if(pd_interp==nmm_keep_pd) then
2453        ! PD interpolation is turned off, so we use the original PD.
2454     elseif(pd_interp==nmm_interp_pd) then
2455        ! PD interpolation is requested.  As with the old base_state_parent,
2456        ! determine PD by interpolating or extrapolating the non-hydrostatic
2457        ! pressure field (PINT) in source grid to the surface height of the
2458        ! target grid.
2459        xloop: do ix=1,nx
2460           if(pintA(1,ix)>p_ref) then
2461              ! Cannot extrapolate PD below ground because pressure is
2462              ! unrealistically high.  Follow base_state_parent method:
2463              ! when this happens, assign input pressure to output
2464              ! pressure.
2465              if(.not.warned) then
2466                 call wrf_message2('WARNING: NESTED DOMAIN PRESSURE AT LOWEST LEVEL HIGHER THAN PSTD')
2467                 write(message,*) 'PINT(1),PD(1),PSTD(1)',pintA(1,ix),pdA(ix),p_ref
2468                 call wrf_message2(message)
2469                 warned=.true.
2470              endif
2471              pdB(ix)=pdA(ix)
2472              cycle xloop
2473           endif
2474
2475           zA=fisA(ix)/g
2476           zB=fisB(ix)/g
2477
2478           if(zB<zA) then
2479              ! Need to extrapolate PD below ground.
2480              ! Follow base_state_parent method:
2481              !   For pressures between pint(1) and 1030 mbar, interpolate
2482              iz=1
2483              weight=abs((zB-zmslp(ix))/(zA-zmslp(ix)))
2484              pdB(ix) = exp(weight*log(pdA(ix)+pdtop+ptop) + &
2485                         (1.0-weight)*log(P_REF)) &
2486                         - pdtop - ptop
2487 !             write(0,*) 'WARNING: extrapolating below ground at ix=',ix
2488 !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)
2489 !             write(0,1302) ix,zA,zB,zmslp(ix),weight,pdB(ix)+pdtop+ptop,pintA(1,ix),pdA(ix)+pdtop+ptop,P_REF
2490              cycle xloop
2491           endif
2492
2493           ! PD in grid B lies above ground in grid A.  Find out where by
2494           ! walking up from the ground until we find the two levels the
2495           ! pressure lies between.
2496           znext=-9e9
2497           z=zA
2498           do iz=1,nz-1
2499              ! Calculations from base_state_parent:
2500              APELP    = (pintA(iz+1,ix)+pintA(iz,ix))
2501              RTOPP    = TRG*tA(iz,ix)*(1.0+qA(iz,ix)*P608)/APELP
2502              DZ       = RTOPP*(DETA1(iz)*PDTOP+DETA2(iz)*pdA(ix))
2503              znext=z+dz
2504
2505              if(znext>zB) then
2506                 ! Interpolate between these two levels
2507                 weight=(zB-z)/dz
2508                 pdB(ix) = exp(weight*log(pintA(iz+1,ix)) + &
2509                               (1.0-weight)*log(pintA(iz,ix))) &
2510                       - pdtop - ptop
2511 !                if(iz>1) then
2512 !                   write(0,*) 'WARNING: interpolate above ground at ix=',ix
2513 !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)
2514 !                   write(0,1303) ix,iz,zA,zB,zmslp(ix),weight,pdB(ix)+pdtop+ptop,pintA(1,ix),pdA(ix)+pdtop+ptop,P_REF
2515 !                endif
2516                 cycle xloop
2517              else
2518                 z=znext
2519              endif
2520
2521           enddo
2522
2523 201       format('interp_T_PD_Q: Target domain surface height ',F0.4,'m is higher than the model top ',F0.4,'m of the source domain.  ABORTING')
2524           write(message,201) zB,znext
2525           call wrf_error_fatal(message)
2526        enddo xloop
2527     else
2528 202    format('Invalid value ',I0,' for pd_interp in interp_T_PD_Q')
2529        write(message,202) pd_interp
2530        call wrf_error_fatal(message)
2531     endif if_pd_interp
2532
2533
2534
2535     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2536     !! Step 3: interpolate T and Q in pressure space !!!!!!!!!!!!!
2537     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2538
2539
2540
2541     outer: do ix=1,nx
2542        ! For constant pressure levels, do a straight level-by-level copy:
2543        iz=nz-1
2544        copyloop: do while(iz>kpres)
2545           tB(iz,ix)=tA(iz,ix)
2546           qB(iz,ix)=qA(iz,ix)
2547
2548           ! Save interpolation information
2549           iinfo(iz,ix)=iz
2550           winfo(iz,ix)=1.0
2551
2552           iz=iz-1
2553        enddo copyloop
2554
2555        ! For sigma levels where target domain lies within source,
2556        ! interpolate from top down, stopping when we go below ground
2557        ! in the source domain.
2558        izA=iz
2559        izB=iz
2560
2561        ! FIXME: REMOVE THIS CHECK
2562        if(iz>nz) then
2563           ! Make sure the entire vertical extent isn't sigma levels:
2564           call wrf_error_fatal('ERROR: WRF-NMM does not support pure sigma levels (only sigma-pressure hybrid)')
2565        endif
2566
2567        pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2568        pB1=log(eta1(izB+1)*pdtop + eta2(izB+1)*pdB(ix) + ptop)
2569        pB=(pB2+pB1)*0.5
2570
2571        pA2=log(eta1(izA)*pdtop + eta2(izA)*pdA(ix) + ptop)
2572        pA1=log(eta1(izA+1)*pdtop + eta2(izA+1)*pdA(ix) + ptop)
2573        pA=(pA2+pA1)*0.5
2574
2575        ! Find the lowest mass level izA that is above pB
2576        interpinit: do while(izA>1)
2577           pA3=log(eta1(izA-1)*pdtop + eta2(izA-1)*pdA(ix) + ptop)
2578           pnext=(pA2+pA3)*0.5
2579           wdenom=pnext-pA
2580           wnum=pnext-pb
2581           if(pA<=pB .and. wnum>1e-5) then
2582              exit interpinit
2583           else
2584              pA1=pA2
2585              pA2=pa3
2586              izA=izA-1
2587              pA=pnext
2588           endif
2589        enddo interpinit
2590
2591        ! Loop over all remaining B points, interpolating to them.
2592        interploop: do while(izB>0 .and. izA>1)
2593           ! Find the source domain levels that this target domain level lies between:
2594           zinterp: do while(izA>1)
2595              if(pnext>=pB) then
2596                 ! We found the two levels, so interpolate and move on to
2597                 ! the next target level:
2598                 weight=max(0.,wnum/wdenom)
2599                 tB(izB,ix)=weight*tA(izA,ix) + (1.0-weight)*tA(izA-1,ix)
2600                 qB(izB,ix)=weight*qA(izA,ix) + (1.0-weight)*qA(izA-1,ix)
2601
2602                 ! Save interpolation info
2603                 iinfo(izB,ix)=izA    ! upper level
2604                 winfo(izB,ix)=weight ! linear interpolation weight
2605
2606                 ! We interpolated to a B point, so go to the next one:
2607                 pB1=pB2
2608                 izB=izB-1
2609                 if(izB>0) then
2610                    pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2611                    pB=(pb1+pb2)*0.5
2612                 else
2613                    exit interploop
2614                 endif
2615                 wnum=pnext-pb
2616
2617                 exit zinterp
2618              else
2619                 izA=izA-1
2620                 pA1=pA2
2621                 pA2=pa3
2622                 pA=pnext
2623                 if(izA>1) then
2624                    pA3=log(eta1(izA-1)*pdtop + eta2(izA-1)*pdA(ix) + ptop)
2625                    pnext=(pA2+pA3)*0.5
2626                    wdenom=pnext-pa
2627                    wnum=pnext-pb
2628                 else
2629                    exit interploop
2630                 endif
2631              endif
2632           enddo zinterp
2633        enddo interploop
2634
2635        ! Follow what base_state_parent did for temperature:
2636        ! interpolate between the second to last level and P_REF using
2637        ! a lapse rate atmosphere.  For Q, we use constant RH below
2638        ! ground, as suggested by Greg Thompson.
2639        extraploop: do while(izB>=1)
2640           ! Decide extrapolation weight:
2641           weight=(pB-pA)/(pstd12-pA)
2642
2643           ! Extrapolate Q by copying lowest sigma layer value:
2644           !qB(izB,ix)=qA(1,ix)
2645
2646           ! Extrapolate Q by linearly interpolating between 0 and surface value:
2647           !qB(izB,ix)=(1.0-weight)*qA(1,ix)
2648
2649           ! Extrapolate T using a lapse rate atmosphere
2650           tB(izB,ix)=weight*tbelow(ix) + (1.0-weight)*tA(1,ix)
2651
2652           ! Extrapolate Q using constant RH below ground, as suggested
2653           ! by Greg Thompson.  This is the inversion of the RH
2654           ! calculation used earlier in this function:
2655           P0=eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop
2656           QC=PQ0/P0*EXP(A2*(TB(izB,ix)-A3)/(TB(izB,ix)-A4))
2657           QB(izB,ix)=QC*RHbelow(ix)
2658
2659           ! Save extrapolation information
2660           iinfo(izB,ix)=0
2661           winfo(izB,ix)=weight
2662
2663           ! Move to the next B level
2664           izB=izB-1
2665           if(izB>0) then
2666              pB1=pB2
2667              pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2668              pB=(pb1+pb2)*0.5
2669           endif
2670        enddo extraploop
2671     enddo outer
2672   end subroutine interp_T_PD_Q
2673
2674   subroutine interp_T_PD_Q_kpres(method, pd_interp, nx, nz, &
2675        deta1,deta2, eta1,eta2, ptop,pdtop, kpres, nz2,     &
2676        fisA,fisB, pintA,pintB, tA,tB, pdA,pdB, qA,qB, &
2677        iinfo, winfo, warned)
2678     implicit none
2679
2680     integer, intent(in) :: pd_interp,method, kpres, nz2
2681     !      real, intent(in) :: dtA,dtB
2682
2683     ! Coordinate system definitions must be the same for all domains:
2684     real, intent(in) :: deta1(nz),deta2(nz),eta1(nz),eta2(nz),ptop,pdtop
2685     integer, intent(in) :: nx,nz
2686
2687     ! Surface height and mass field information for source (A):
2688     real, intent(in) :: fisA(nx),tA(nz2-1,nx),pdA(nx),qA(nz2-1,nx),pintA(nz2,nx)
2689
2690     ! Surface height and mass field information for target (B):
2691     real, intent(inout) :: fisB(nx),tB(nz2-1,nx),pdB(nx),qB(nz2-1,nx),pintB(nz2,nx)
2692
2693     ! Interpolation or extrapolation information for use in other
2694     ! calls later on:
2695     real, intent(out) :: winfo(nz2-1,nx)
2696     integer, intent(out) :: iinfo(nz2-1,nx)
2697
2698     logical, intent(inout) :: warned
2699
2700     ! ==================== Local variables ====================
2701
2702     character*255 :: message
2703     integer :: ix,iz,izA,izB,xpr
2704     real :: zA,zB,znext,apelp,rtopp,dz,weight,A,B,pstd1,z,pstd2,pstd12
2705     real :: tsfc(nx), slp(nx), zmslp(nx), z0mid(nx), zbelow(nx), &
2706             tbelow(nx), RHbelow(nx)
2707     real :: pb,pb1,pb2,pa,pa1,pa2,pnext,pa3, wnum,wdenom, QC, P0
2708
2709     ! Constants from UPP params.f for RH calculation (all mks):
2710     real, parameter :: PQ0=379.90516
2711     real, parameter :: A2=17.2693882
2712     real, parameter :: A3=273.16
2713     real, parameter :: A4=35.86
2714     real, parameter :: RHmin=1.0E-6     ! minimal RH bound
2715
2716     if(method/=nmm_method_linear) then
2717        call wrf_error_fatal('only linear interpolation is supported')
2718     endif
2719
2720     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2721     !! Step 1: calculate near-surface values !!!!!!!!!!!!!!!!!!!!!
2722     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2723
2724
2725
2726     pstd1=p_ref ! pstd(1) from base_state_parent
2727     ! pstd(2) from base_state_parent:
2728     pstd2=eta1(2)*pdtop + eta2(2)*(p_ref-pdtop-ptop) + ptop
2729     pstd12=exp((alog(pstd1)+alog(pstd2))*0.5)
2730
2731     do ix=1,nx
2732        ! These calculations are from base_state_parent:
2733        APELP    = (pintA(2,ix)+pintA(1,ix))
2734        RTOPP    = TRG*tA(1,ix)*(1.0+qA(1,ix)*P608)/APELP
2735        DZ       = RTOPP*(DETA1(1)*PDTOP+DETA2(1)*pdA(ix))
2736
2737        Z0MID(ix) = fisA(ix)/g + dz/2
2738
2739        zA=fisA(ix)/g
2740
2741        TSFC(ix) = TA(1,ix)*(1.+D608*QA(1,ix)) &
2742             + LAPSR*(zA + zA+dz)*0.5
2743
2744        A         = LAPSR*zA/TSFC(ix)
2745        SLP(ix)  = PINTA(1,ix)*(1-A)**COEF2    ! sea level pressure
2746        B         = (pstd1/SLP(ix))**COEF3
2747        ZMSLP(ix)= TSFC(ix)*LAPSI*(1.0 - B)   ! Height at 1030. mb level
2748
2749        TBELOW(ix) = TA(1,ix) + LAPSR*(Z0MID(ix)-ZMSLP(ix))
2750
2751        ! Calculate lowest level RH.  This calculation is from UPP
2752        ! CALRH.f:
2753        P0=pdA(ix)+ptop+pdtop ! Use hydrostatic lowest model level P
2754        QC=PQ0/P0*EXP(A2*(TA(1,ix)-A3)/(TA(1,ix)-A4))
2755        RHbelow(ix)=max(RHmin,min(1.0,QA(1,ix)/QC))
2756     enddo
2757
2758
2759
2760     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2761     !! Step 2: figure out the new surface pressures !!!!!!!!!!!!!!
2762     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2763
2764
2765
2766     if_pd_interp: if(pd_interp==nmm_keep_pd) then
2767        ! PD interpolation is turned off, so we use the original PD.
2768     elseif(pd_interp==nmm_interp_pd) then
2769        ! PD interpolation is requested.  As with the old base_state_parent,
2770        ! determine PD by interpolating or extrapolating the non-hydrostatic
2771        ! pressure field (PINT) in source grid to the surface height of the
2772        ! target grid.
2773        xloop: do ix=1,nx
2774           if(pintA(1,ix)>p_ref) then
2775              ! Cannot extrapolate PD below ground because pressure is
2776              ! unrealistically high.  Follow base_state_parent method:
2777              ! when this happens, assign input pressure to output
2778              ! pressure.
2779              if(.not.warned) then
2780                 call wrf_message2('WARNING: NESTED DOMAIN PRESSURE AT LOWEST LEVEL HIGHER THAN PSTD')
2781                 WRITE(message,*) 'PINT(1),PD(1),PSTD(1)',pintA(1,ix),pdA(ix),p_ref
2782                 call wrf_message2(message)
2783                 warned=.true.
2784              endif
2785              pdB(ix)=pdA(ix)
2786              cycle xloop
2787           endif
2788
2789           zA=fisA(ix)/g
2790           zB=fisB(ix)/g
2791
2792           if(zB<zA) then
2793              ! Need to extrapolate PD below ground.
2794              ! Follow base_state_parent method:
2795              !   For pressures between pint(1) and 1030 mbar, interpolate
2796              iz=1
2797              weight=abs((zB-zmslp(ix))/(zA-zmslp(ix)))
2798              pdB(ix) = exp(weight*log(pdA(ix)+pdtop+ptop) + &
2799                         (1.0-weight)*log(P_REF)) &
2800                         - pdtop - ptop
2801 !             write(0,*) 'WARNING: extrapolating below ground at ix=',ix
2802 !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)
2803 !             write(0,1302) ix,zA,zB,zmslp(ix),weight,pdB(ix)+pdtop+ptop,pintA(1,ix),pdA(ix)+pdtop+ptop,P_REF
2804              cycle xloop
2805           endif
2806
2807           ! PD in grid B lies above ground in grid A.  Find out where by
2808           ! walking up from the ground until we find the two levels the
2809           ! pressure lies between.
2810           znext=-9e9
2811           z=zA
2812           do iz=1,nz2-1
2813              ! Calculations from base_state_parent:
2814              APELP    = (pintA(iz+1,ix)+pintA(iz,ix))
2815              RTOPP    = TRG*tA(iz,ix)*(1.0+qA(iz,ix)*P608)/APELP
2816              DZ       = RTOPP*(DETA1(iz)*PDTOP+DETA2(iz)*pdA(ix))
2817              znext=z+dz
2818
2819              if(znext>zB) then
2820                 ! Interpolate between these two levels
2821                 weight=(zB-z)/dz
2822                 pdB(ix) = exp(weight*log(pintA(iz+1,ix)) + &
2823                               (1.0-weight)*log(pintA(iz,ix))) &
2824                       - pdtop - ptop
2825 !                if(iz>1) then
2826 !                   write(0,*) 'WARNING: interpolate above ground at ix=',ix
2827 !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)
2828 !                   write(0,1303) ix,iz,zA,zB,zmslp(ix),weight,pdB(ix)+pdtop+ptop,pintA(1,ix),pdA(ix)+pdtop+ptop,P_REF
2829 !                endif
2830                 cycle xloop
2831              else
2832                 z=znext
2833              endif
2834
2835           enddo
2836
2837 201       format('interp_T_PD_Q_kpres: Target domain surface height ',F0.4,'m is higher than the model top ',F0.4,'m of the source domain.  ABORTING')
2838           write(message,201) zB,znext
2839           call wrf_error_fatal(message)
2840        enddo xloop
2841     else
2842 202    format('Invalid value ',I0,' for pd_interp in interp_T_PD_Q')
2843        write(message,202) pd_interp
2844        call wrf_error_fatal(message)
2845     endif if_pd_interp
2846
2847
2848
2849     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2850     !! Step 3: interpolate T and Q in pressure space !!!!!!!!!!!!!
2851     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2852
2853
2854
2855     outer: do ix=1,nx
2856        ! For constant pressure levels, do a straight level-by-level copy:
2857        iz=nz2-1
2858        copyloop: do while(iz>kpres)
2859           tB(iz,ix)=tA(iz,ix)
2860           qB(iz,ix)=qA(iz,ix)
2861
2862           ! Save interpolation information
2863           iinfo(iz,ix)=iz
2864           winfo(iz,ix)=1.0
2865
2866           iz=iz-1
2867        enddo copyloop
2868
2869        ! For sigma levels where target domain lies within source,
2870        ! interpolate from top down, stopping when we go below ground
2871        ! in the source domain.
2872        izA=iz
2873        izB=iz
2874
2875        ! FIXME: REMOVE THIS CHECK
2876        if(iz>nz2) then
2877           ! Make sure the entire vertical extent isn't sigma levels:
2878           call wrf_error_fatal('ERROR: WRF-NMM does not support pure sigma levels (only sigma-pressure hybrid)')
2879        endif
2880
2881        pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2882        pB1=log(eta1(izB+1)*pdtop + eta2(izB+1)*pdB(ix) + ptop)
2883        pB=(pB2+pB1)*0.5
2884
2885        pA2=log(eta1(izA)*pdtop + eta2(izA)*pdA(ix) + ptop)
2886        pA1=log(eta1(izA+1)*pdtop + eta2(izA+1)*pdA(ix) + ptop)
2887        pA=(pA2+pA1)*0.5
2888
2889        ! Find the lowest mass level izA that is above pB
2890        interpinit: do while(izA>1)
2891           pA3=log(eta1(izA-1)*pdtop + eta2(izA-1)*pdA(ix) + ptop)
2892           pnext=(pA2+pA3)*0.5
2893           wdenom=pnext-pA
2894           wnum=pnext-pb
2895           if(pA<=pB .and. wnum>1e-5) then
2896              exit interpinit
2897           else
2898              pA1=pA2
2899              pA2=pa3
2900              izA=izA-1
2901              pA=pnext
2902           endif
2903        enddo interpinit
2904
2905        ! Loop over all remaining B points, interpolating to them.
2906        interploop: do while(izB>0 .and. izA>1)
2907           ! Find the source domain levels that this target domain level lies between:
2908           zinterp: do while(izA>1)
2909              if(pnext>=pB) then
2910                 ! We found the two levels, so interpolate and move on to
2911                 ! the next target level:
2912                 weight=max(0.,wnum/wdenom)
2913                 tB(izB,ix)=weight*tA(izA,ix) + (1.0-weight)*tA(izA-1,ix)
2914                 qB(izB,ix)=weight*qA(izA,ix) + (1.0-weight)*qA(izA-1,ix)
2915
2916                 ! Save interpolation info
2917                 iinfo(izB,ix)=izA    ! upper level
2918                 winfo(izB,ix)=weight ! linear interpolation weight
2919
2920                 ! We interpolated to a B point, so go to the next one:
2921                 pB1=pB2
2922                 izB=izB-1
2923                 if(izB>0) then
2924                    pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2925                    pB=(pb1+pb2)*0.5
2926                 else
2927                    exit interploop
2928                 endif
2929                 wnum=pnext-pb
2930
2931                 exit zinterp
2932              else
2933                 izA=izA-1
2934                 pA1=pA2
2935                 pA2=pa3
2936                 pA=pnext
2937                 if(izA>1) then
2938                    pA3=log(eta1(izA-1)*pdtop + eta2(izA-1)*pdA(ix) + ptop)
2939                    pnext=(pA2+pA3)*0.5
2940                    wdenom=pnext-pa
2941                    wnum=pnext-pb
2942                 else
2943                    exit interploop
2944                 endif
2945              endif
2946           enddo zinterp
2947        enddo interploop
2948
2949        ! Follow what base_state_parent did for temperature:
2950        ! interpolate between the second to last level and P_REF using
2951        ! a lapse rate atmosphere.  For Q, we use constant RH below
2952        ! ground, as suggested by Greg Thompson.
2953        extraploop: do while(izB>=1)
2954           ! Decide extrapolation weight:
2955           weight=(pB-pA)/(pstd12-pA)
2956
2957           ! Extrapolate Q by copying lowest sigma layer value:
2958           !qB(izB,ix)=qA(1,ix)
2959
2960           ! Extrapolate Q by linearly interpolating between 0 and surface value:
2961           !qB(izB,ix)=(1.0-weight)*qA(1,ix)
2962
2963           ! Extrapolate T using a lapse rate atmosphere
2964           tB(izB,ix)=weight*tbelow(ix) + (1.0-weight)*tA(1,ix)
2965
2966           ! Extrapolate Q using constant RH below ground, as suggested
2967           ! by Greg Thompson.  This is the inversion of the RH
2968           ! calculation used earlier in this function:
2969           P0=eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop
2970           QC=PQ0/P0*EXP(A2*(TB(izB,ix)-A3)/(TB(izB,ix)-A4))
2971           QB(izB,ix)=QC*RHbelow(ix)
2972
2973           ! Save extrapolation information
2974           iinfo(izB,ix)=0
2975           winfo(izB,ix)=weight
2976
2977           ! Move to the next B level
2978           izB=izB-1
2979           if(izB>0) then
2980              pB1=pB2
2981              pB2=log(eta1(izB)*pdtop + eta2(izB)*pdB(ix) + ptop)
2982              pB=(pb1+pb2)*0.5
2983           endif
2984        enddo extraploop
2985     enddo outer
2986   end subroutine interp_T_PD_Q_kpres
2987 end module module_interp_nmm
2988
2989 ! ********************************************************************
2990 ! subs EXT_*_FULLDOM -- simple wrappers around module_interp_nmm's
2991 ! *_FULLDOM routines.  These exist solely to allow module_dm to
2992 ! interface with this module -- this is needed because the WRF build
2993 ! scripts compile module_dm before module_interp_nmm.
2994 ! ********************************************************************
2995
2996 subroutine ext_c2n_fulldom  (II,JJ,W1,W2,W3,W4,&
2997        deta1,deta2, eta1,eta2, ptop,pdtop,   &
2998        cpint,ct,cpd,cq,                 &
2999        cims, cime, cjms, cjme, ckms, ckme,   &
3000        npint,nt,npd,nq,                 &
3001        out_iinfo,out_winfo,imask,&
3002        nids, nide, njds, njde, nkds, nkde,   &
3003        nims, nime, njms, njme, nkms, nkme,   &
3004        nits, nite, njts, njte, nkts, nkte)
3005   ! This subroutine is an alias for c2n_fulldom.  It exists to allow
3006   ! the routine to be called from module_dm
3007   use module_interp_store, only: parent_fis, nest_fis, kpres
3008   use module_interp_nmm, only: c2n_fulldom, find_kpres
3009   implicit none
3010   integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme,   &
3011        nids, nide, njds, njde, nkds, nkde,   &
3012        nims, nime, njms, njme, nkms, nkme,   &
3013        nits, nite, njts, njte, nkts, nkte
3014   real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
3015   real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD
3016
3017   real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_winfo
3018   integer, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: out_iinfo
3019
3020   real, intent(out), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
3021   real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
3022   real, intent(inout), dimension(nims:nime,njms:njme,1) :: nPD
3023   integer, intent(in), dimension(nims:nime,njms:njme) :: imask
3024
3025   real, intent(in), dimension(nims:nime,njms:njme) :: &
3026        W1,W2,W3,W4
3027   integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
3028
3029   real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
3030   real, intent(in) :: ptop,pdtop
3031   integer :: i,j,k
3032   logical :: badbad
3033
3034   call find_kpres(kpres,eta2,nkds,nkde,nkms,nkme)
3035   call c2n_fulldom(II,JJ,W1,W2,W3,W4,&
3036        deta1,deta2, eta1,eta2, ptop,pdtop,   &
3037        parent_fis,cpint,ct,cpd,cq,           &
3038        cims, cime, cjms, cjme, ckms, ckme,   &
3039        nest_fis,npint,nt,npd,nq,             &
3040        out_iinfo,out_winfo,imask,            &
3041        nids, nide, njds, njde, nkds, nkde,   &
3042        nims, nime, njms, njme, nkms, nkme,   &
3043        nits, nite, njts, njte, nkts, nkte,   &
3044        kpres)
3045
3046 end subroutine ext_c2n_fulldom
3047
3048 subroutine ext_n2c_fulldom  ( &
3049        deta1,deta2, eta1,eta2, ptop,pdtop,   &
3050        cpint,ct,cpd,cq,                 &
3051        cids, cide, cjds, cjde, ckds, ckde,   &
3052        cims, cime, cjms, cjme, ckms, ckme,   &
3053        cits, cite, cjts, cjte, ckts, ckte,   &
3054        npint,nt,npd,nq,                 &
3055        ipos,jpos,                            &
3056        out_iinfo,out_winfo,      &
3057        nids, nide, njds, njde, nkds, nkde,   &
3058        nims, nime, njms, njme, nkms, nkme,   &
3059        nits, nite, njts, njte, nkts, nkte)
3060   ! This subroutine is an alias for n2c_fulldom.  It exists to allow
3061   ! the routine to be called from module_dm
3062   use module_interp_store, only: parent_fis, nest_fis, kpres
3063   use module_interp_nmm, only: n2c_fulldom, find_kpres, n2c_fulldom_new
3064   implicit none
3065   integer, intent(in):: cims, cime, cjms, cjme, ckms, ckme,   &
3066        cids, cide, cjds, cjde, ckds, ckde,   &
3067        cits, cite, cjts, cjte, ckts, ckte,   &
3068        nids, nide, njds, njde, nkds, nkde,   &
3069        nims, nime, njms, njme, nkms, nkme,   &
3070        nits, nite, njts, njte, nkts, nkte,   &
3071        ipos, jpos ! nest start loctions within parent
3072   real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
3073   real, intent(inout), dimension(cims:cime,cjms:cjme,1) :: cPD
3074
3075   real, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_winfo
3076   integer, intent(out), dimension(cims:cime,cjms:cjme,ckms:ckme) :: out_iinfo
3077
3078   real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nT,nQ
3079   real, intent(in), dimension(nims:nime,njms:njme,nkms:nkme) :: nPINT
3080   real, intent(in), dimension(nims:nime,njms:njme,1) :: nPD
3081
3082   real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
3083   real, intent(in) :: ptop,pdtop
3084
3085   call find_kpres(kpres,eta2,nkds,nkde,nkms,nkme)
3086   call n2c_fulldom_new( &
3087        deta1,deta2, eta1,eta2, ptop,pdtop,   &
3088        parent_fis,cpint,ct,cpd,cq,           &
3089        cids, cide, cjds, cjde, ckds, ckde,   &
3090        cims, cime, cjms, cjme, ckms, ckme,   &
3091        cits, cite, cjts, cjte, ckts, ckte,   &
3092        nest_fis,npint,nt,npd,nq,             &
3093        ipos,jpos,                            &
3094        out_iinfo,out_winfo,                  &
3095        nids, nide, njds, njde, nkds, nkde,   &
3096        nims, nime, njms, njme, nkms, nkme,   &
3097        nits, nite, njts, njte, nkts, nkte, kpres)
3098 end subroutine ext_n2c_fulldom
3099
3100 subroutine ext_c2b_fulldom  (II,JJ,W1,W2,W3,W4,&
3101      deta1,deta2, eta1,eta2, ptop,pdtop,   &
3102      cpint,ct,cpd,cq,                      &
3103      cims, cime, cjms, cjme, ckms, ckme,   &
3104      nids, nide, njds, njde, nkds, nkde,   &
3105      nims, nime, njms, njme, nkms, nkme,   &
3106      nits, nite, njts, njte, nkts, nkte,   &
3107      ibxs,    ibxe,    ibys,    ibye,      &
3108      wbxs,    wbxe,    wbys,    wbye,      &
3109      pdbxs,   pdbxe,   pdbys,   pdbye,     &
3110      tbxs,    tbxe,    tbys,    tbye,      &
3111      qbxs,    qbxe,    qbys,    qbye)
3112   use module_interp_nmm, only: c2b_fulldom, find_kpres, c2b_fulldom_new
3113   use module_interp_store, only: parent_fis, nest_fis, kpres
3114   implicit none
3115   integer, intent(in):: &
3116        cims, cime, cjms, cjme, ckms, ckme,   &
3117        nids, nide, njds, njde, nkds, nkde,   &
3118        nims, nime, njms, njme, nkms, nkme,   &
3119        nits, nite, njts, njte, nkts, nkte
3120   integer, parameter :: bdyw = 1
3121   ! Domain info:
3122   real, intent(in), dimension(nims:nime,njms:njme) :: W1,W2,W3,W4
3123   integer, intent(in), dimension(nims:nime,njms:njme) :: II,JJ
3124   real, intent(in), dimension(nkms:nkme) :: eta1,eta2,deta1,deta2
3125   real, intent(in) :: ptop,pdtop
3126
3127   ! Parent fields:
3128   real, intent(in), dimension(cims:cime,cjms:cjme,ckms:ckme) :: cT,cQ,cPINT
3129   real, intent(in), dimension(cims:cime,cjms:cjme,1) :: cPD
3130
3131   ! T, Q, PINT, PD boundary info:
3132   real,dimension(nims:nime,nkms:nkme,bdyw) :: tbys,tbye,qbys,qbye
3133   real,dimension(njms:njme,nkms:nkme,bdyw) :: tbxs,tbxe,qbxs,qbxe
3134   real,dimension(nims:nime,1,bdyw) :: pdbys,pdbye
3135   real,dimension(njms:njme,1,bdyw) :: pdbxs,pdbxe
3136
3137   ! Weights and indices:
3138   real,dimension(nims:nime,nkms:nkme,bdyw) :: wbys,wbye
3139   real,dimension(njms:njme,nkms:nkme,bdyw) :: wbxs,wbxe
3140   integer,dimension(nims:nime,nkms:nkme,bdyw) :: ibys,ibye
3141   integer,dimension(njms:njme,nkms:nkme,bdyw) :: ibxs,ibxe
3142
3143   call find_kpres(kpres,eta2,nkds,nkde,nkms,nkme)
3144   call c2b_fulldom_new    (II,JJ,W1,W2,W3,W4,&
3145        deta1,deta2, eta1,eta2, ptop,pdtop,   &
3146        parent_fis,cpint,ct,cpd,cq,nest_fis,  &
3147        cims, cime, cjms, cjme, ckms, ckme,   &
3148        nids, nide, njds, njde, nkds, nkde,   &
3149        nims, nime, njms, njme, nkms, nkme,   &
3150        nits, nite, njts, njte, nkts, nkte,   &
3151        kpres,                                &
3152        ibxs,    ibxe,    ibys,    ibye,      &
3153        wbxs,    wbxe,    wbys,    wbye,      &
3154        pdbxs,   pdbxe,   pdbys,   pdbye,     &
3155        tbxs,    tbxe,    tbys,    tbye,      &
3156        qbxs,    qbxe,    qbys,    qbye)
3157 end subroutine ext_c2b_fulldom