chem/module_interpolate.F

   1 module module_interpolate
   2 ! Description:
   3 !   Routines for interpolation of data in latitude, longitude and time.
   4 !
   5 ! Author: Gathered from a number of places and put into the current format by Jim Edwards
   6 !
   7 ! Modules Used:
   8 !
   9 ! use shr_kind_mod,   only: r8 => shr_kind_r8
  10 ! use abortutils,     only: endrun
  11 ! use scamMod,        only: single_column
  12 ! use cam_logfile,    only: iulog
  13
  14   implicit none
  15
  16   private
  17
  18   integer,parameter :: r8 = selected_real_kind(12) ! 8 byte real
  19 ! integer,parameter :: r8 = selected_real_kind(8) ! 8 byte real
  20 !
  21 ! Public Methods:
  22 !
  23   public :: interp_type, lininterp, vertinterp, bilin, get_timeinterp_factors
  24   public :: lininterp_init, lininterp_finish
  25   type interp_type
  26      real(r8), pointer :: wgts(:)
  27      real(r8), pointer :: wgtn(:)
  28      integer, pointer  :: jjm(:)
  29      integer, pointer  :: jjp(:)
  30   end type interp_type
  31   interface lininterp
  32      module procedure lininterp_original
  33      module procedure lininterp_full1d
  34      module procedure lininterp1d
  35      module procedure lininterp2d2d
  36      module procedure lininterp2d1d
  37      module procedure lininterp3d2d
  38   end interface
  39
  40
  41 contains
  42   subroutine lininterp_full1d (arrin, yin, nin, arrout, yout, nout)
  43     integer, intent(in) :: nin, nout
  44     real(r8), intent(in) :: arrin(nin), yin(nin), yout(nout)
  45     real(r8), intent(out) :: arrout(nout)
  46     type (interp_type) :: interp_wgts
  47
  48     call lininterp_init(yin, nin, yout, nout, 1, interp_wgts)
  49     call lininterp1d(arrin, nin, arrout, nout, interp_wgts)
  50     call lininterp_finish(interp_wgts)
  51
  52   end subroutine lininterp_full1d
  53
  54   subroutine lininterp_init(yin, nin, yout, nout, extrap_method, interp_wgts, &
  55        cyclicmin, cyclicmax)
  56 !
  57 ! Description:
  58 !   Initialize a variable of type(interp_type) with weights for linear interpolation.
  59 !       this variable can then be used in calls to lininterp1d and lininterp2d.
  60 !   yin is a 1d array of length nin of locations to interpolate from - this array must
  61 !       be monotonic but can be increasing or decreasing
  62 !   yout is a 1d array of length nout of locations to interpolate to, this array need
  63 !       not be ordered
  64 !   extrap_method determines how to handle yout points beyond the bounds of yin
  65 !       if 0 set values outside output grid to 0
  66 !       if 1 set to boundary value
  67 !       if 2 set to cyclic boundaries
  68 !         optional values cyclicmin and cyclicmax can be used to set the bounds of the
  69 !         cyclic mapping - these default to 0 and 360.
  70 !
  71
  72     integer, intent(in) :: nin
  73     integer, intent(in) :: nout
  74     real(r8), intent(in) :: yin(:)           ! input mesh
  75     real(r8), intent(in) :: yout(:)         ! output mesh
  76     integer, intent(in) :: extrap_method       ! if 0 set values outside output grid to 0
  77                                                ! if 1 set to boundary value
  78                                                ! if 2 set to cyclic boundaries
  79     real(r8), intent(in), optional :: cyclicmin, cyclicmax
  80
  81     type (interp_type), intent(out) :: interp_wgts
  82
  83     real(r8) :: cmin, cmax
  84     real(r8) :: extrap
  85     real(r8) :: dyinwrap
  86     real(r8) :: ratio
  87     real(r8) :: avgdyin
  88     integer :: i, j, icount
  89     integer :: jj
  90     real(r8), pointer :: wgts(:)
  91     real(r8), pointer :: wgtn(:)
  92     integer, pointer :: jjm(:)
  93     integer, pointer :: jjp(:)
  94     logical :: increasing
  95     character(len=132) :: message
  96     !
  97     ! Check validity of input coordinate arrays: must be monotonically increasing,
  98     ! and have a total of at least 2 elements
  99     !
 100     if (nin.lt.2) then
 101        call wrf_abort  !('LININTERP: Must have at least 2 input points for interpolation')
 102     end if
 103     if(present(cyclicmin)) then
 104        cmin=cyclicmin
 105     else
 106        cmin=0_r8
 107     end if
 108     if(present(cyclicmax)) then
 109        cmax=cyclicmax
 110     else
 111        cmax=360_r8
 112     end if
 113     if(cmax<=cmin) then
 114        call wrf_abort  !('LININTERP: cyclic min value must be < max value')
 115     end if
 116     increasing=.true.
 117     icount = 0
 118     do j=1,nin-1
 119        if (yin(j).gt.yin(j+1)) icount = icount + 1
 120     end do
 121     if(icount.eq.nin-1) then
 122        increasing = .false.
 123        icount=0
 124     endif
 125     if (icount.gt.0) then
 126        call wrf_abort  !('LININTERP: Non-monotonic input coordinate array found')
 127     end if
 128     allocate(interp_wgts%jjm(nout), &
 129          interp_wgts%jjp(nout), &
 130          interp_wgts%wgts(nout), &
 131          interp_wgts%wgtn(nout))
 132
 133     jjm => interp_wgts%jjm
 134     jjp => interp_wgts%jjp
 135     wgts =>  interp_wgts%wgts
 136     wgtn =>  interp_wgts%wgtn
 137
 138     !
 139     ! Initialize index arrays for later checking
 140     !
 141     jjm = 0
 142     jjp = 0
 143
 144     extrap = 0.
 145     if(extrap_method.eq.0) then
 146        !
 147        ! For values which extend beyond boundaries, set weights
 148        ! such that values will be 0.
 149        !
 150        do j=1,nout
 151           if(increasing) then
 152              if (yout(j).lt.yin(1)) then
 153                 jjm(j) = 1
 154                 jjp(j) = 1
 155                 wgts(j) = 0.
 156                 wgtn(j) = 0.
 157                 extrap = extrap + 1.
 158              else if (yout(j).gt.yin(nin)) then
 159                 jjm(j) = nin
 160                 jjp(j) = nin
 161                 wgts(j) = 0.
 162                 wgtn(j) = 0.
 163                 extrap = extrap + 1.
 164              end if
 165           else
 166              if (yout(j).gt.yin(1)) then
 167                 jjm(j) = 1
 168                 jjp(j) = 1
 169                 wgts(j) = 0.
 170                 wgtn(j) = 0.
 171                 extrap = extrap + 1.
 172              else if (yout(j).lt.yin(nin)) then
 173                 jjm(j) = nin
 174                 jjp(j) = nin
 175                 wgts(j) = 0.
 176                 wgtn(j) = 0.
 177                 extrap = extrap + 1.
 178              end if
 179           end if
 180        end do
 181     else if(extrap_method.eq.1) then
 182        !
 183        ! For values which extend beyond boundaries, set weights
 184        ! such that values will just be copied.
 185        !
 186        do j=1,nout
 187           if(increasing) then
 188              if (yout(j).le.yin(1)) then
 189                 jjm(j) = 1
 190                 jjp(j) = 1
 191                 wgts(j) = 1.
 192                 wgtn(j) = 0.
 193                 extrap = extrap + 1.
 194              else if (yout(j).gt.yin(nin)) then
 195                 jjm(j) = nin
 196                 jjp(j) = nin
 197                 wgts(j) = 1.
 198                 wgtn(j) = 0.
 199                 extrap = extrap + 1.
 200              end if
 201           else
 202              if (yout(j).gt.yin(1)) then
 203                 jjm(j) = 1
 204                 jjp(j) = 1
 205                 wgts(j) = 1.
 206                 wgtn(j) = 0.
 207                 extrap = extrap + 1.
 208              else if (yout(j).le.yin(nin)) then
 209                 jjm(j) = nin
 210                 jjp(j) = nin
 211                 wgts(j) = 1.
 212                 wgtn(j) = 0.
 213                 extrap = extrap + 1.
 214              end if
 215           end if
 216        end do
 217     else if(extrap_method.eq.2) then
 218        !
 219        ! For values which extend beyond boundaries, set weights
 220        ! for circular boundaries
 221        !
 222        dyinwrap = yin(1) + (cmax-cmin) - yin(nin)
 223        avgdyin = abs(yin(nin)-yin(1))/(nin-1.)
 224        ratio = dyinwrap/avgdyin
 225        if (ratio < 0.9 .or. ratio > 1.1) then
 226           write(message,*) 'Lininterp: Bad dyinwrap value =',dyinwrap,&
 227                ' avg=', avgdyin, yin(1),yin(nin)
 228           call wrf_message( trim(message) )
 229           call wrf_abort  !('interpolate_data')
 230        end if
 231
 232        do j=1,nout
 233           if(increasing) then
 234              if (yout(j) <= yin(1)) then
 235                 jjm(j) = nin
 236                 jjp(j) = 1
 237                 wgts(j) = (yin(1)-yout(j))/dyinwrap
 238                 wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap
 239              else if (yout(j) > yin(nin)) then
 240                 jjm(j) = nin
 241                 jjp(j) = 1
 242                 wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap
 243                 wgtn(j) = (yout(j)-yin(nin))/dyinwrap
 244              end if
 245           else
 246              if (yout(j) > yin(1)) then
 247                 jjm(j) = nin
 248                 jjp(j) = 1
 249                 wgts(j) = (yin(1)-yout(j))/dyinwrap
 250                 wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap
 251              else if (yout(j) <= yin(nin)) then
 252                 jjm(j) = nin
 253                 jjp(j) = 1
 254                 wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap
 255                 wgtn(j) = (yout(j)+(cmax-cmin)-yin(nin))/dyinwrap
 256              end if
 257
 258           endif
 259        end do
 260     end if
 261
 262     !
 263     ! Loop though output indices finding input indices and weights
 264     !
 265     if(increasing) then
 266        do j=1,nout
 267           do jj=1,nin-1
 268              if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then
 269                 jjm(j) = jj
 270                 jjp(j) = jj + 1
 271                 wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj))
 272                 wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj))
 273                 exit
 274              end if
 275           end do
 276        end do
 277     else
 278        do j=1,nout
 279           do jj=1,nin-1
 280              if (yout(j).le.yin(jj) .and. yout(j).gt.yin(jj+1)) then
 281                 jjm(j) = jj
 282                 jjp(j) = jj + 1
 283                 wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj))
 284                 wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj))
 285                 exit
 286              end if
 287           end do
 288        end do
 289     end if
 290
 291 #ifndef SPMD
 292     !
 293     ! Check grid overlap
 294     !
 295     extrap = 100.*extrap/real(nout,r8)
 296     if (extrap.gt.50. ) then
 297        write(message,*) 'interpolate_data:','yout=',minval(yout),maxval(yout),increasing,nout
 298        call wrf_message( trim(message) )
 299        write(message,*) 'interpolate_data:','yin=',yin(1),yin(nin)
 300        call wrf_message( trim(message) )
 301        write(message,*) 'interpolate_data:',extrap,' % of output grid will have to be extrapolated'
 302        call wrf_message( trim(message) )
 303        call wrf_abort  !('interpolate_data: ')
 304     end if
 305 #endif
 306
 307     !
 308     ! Check that interp/extrap points have been found for all outputs
 309     !
 310     icount = 0
 311     do j=1,nout
 312        if (jjm(j).eq.0 .or. jjp(j).eq.0) icount = icount + 1
 313        ratio=wgts(j)+wgtn(j)
 314        if((ratio<0.9.or.ratio>1.1).and.extrap_method.ne.0) then
 315           write(message,*) j, wgts(j),wgtn(j),jjm(j),jjp(j), increasing,extrap_method
 316           call wrf_message( trim(message) )
 317           call wrf_abort  !('Bad weight computed in LININTERP_init')
 318        end if
 319     end do
 320     if (icount.gt.0) then
 321        call wrf_abort  !('LININTERP: Point found without interp indices')
 322     end if
 323
 324   end subroutine lininterp_init
 325
 326   subroutine lininterp1d (arrin, n1, arrout, m1, interp_wgts)
 327     !-----------------------------------------------------------------------
 328     !
 329     ! Purpose: Do a linear interpolation from input mesh to output
 330     !          mesh with weights as set in lininterp_init.
 331     !
 332     !
 333     ! Author: Jim Edwards
 334     !
 335     !-----------------------------------------------------------------------
 336     !-----------------------------------------------------------------------
 337     implicit none
 338     !-----------------------------------------------------------------------
 339     !
 340     ! Arguments
 341     !
 342     integer, intent(in) :: n1                 ! number of input latitudes
 343     integer, intent(in) :: m1                ! number of output latitudes
 344
 345     real(r8), intent(in) :: arrin(n1)    ! input array of values to interpolate
 346     type(interp_type), intent(in) :: interp_wgts
 347     real(r8), intent(out) :: arrout(m1) ! interpolated array
 348
 349     !
 350     ! Local workspace
 351     !
 352     integer j                ! latitude indices
 353     integer, pointer :: jjm(:)
 354     integer, pointer :: jjp(:)
 355
 356     real(r8), pointer :: wgts(:)
 357     real(r8), pointer :: wgtn(:)
 358
 359
 360     jjm => interp_wgts%jjm
 361     jjp => interp_wgts%jjp
 362     wgts =>  interp_wgts%wgts
 363     wgtn =>  interp_wgts%wgtn
 364
 365     !
 366     ! Do the interpolation
 367     !
 368     do j=1,m1
 369       arrout(j) = arrin(jjm(j))*wgts(j) + arrin(jjp(j))*wgtn(j)
 370     end do
 371
 372     return
 373   end subroutine lininterp1d
 374
 375   subroutine lininterp2d2d(arrin, n1, n2, arrout, m1, m2, wgt1, wgt2)
 376     implicit none
 377     !-----------------------------------------------------------------------
 378     !
 379     ! Arguments
 380     !
 381     integer, intent(in) :: n1, n2, m1, m2
 382     real(r8), intent(in) :: arrin(n1,n2)    ! input array of values to interpolate
 383     type(interp_type), intent(in) :: wgt1, wgt2
 384     real(r8), intent(out) :: arrout(m1,m2) ! interpolated array
 385     !
 386     ! locals
 387     !
 388     integer i,j                ! indices
 389     integer, pointer :: iim(:), jjm(:)
 390     integer, pointer :: iip(:), jjp(:)
 391
 392     real(r8), pointer :: wgts1(:), wgts2(:)
 393     real(r8), pointer :: wgtn1(:), wgtn2(:)
 394
 395     real(r8) :: arrtmp(n1,m2)
 396
 397
 398     jjm => wgt2%jjm
 399     jjp => wgt2%jjp
 400     wgts2 => wgt2%wgts
 401     wgtn2 => wgt2%wgtn
 402
 403     iim => wgt1%jjm
 404     iip => wgt1%jjp
 405     wgts1 => wgt1%wgts
 406     wgtn1 => wgt1%wgtn
 407
 408     do j=1,m2
 409       do i=1,n1
 410         arrtmp(i,j) = arrin(i,jjm(j))*wgts2(j) + arrin(i,jjp(j))*wgtn2(j)
 411       end do
 412     end do
 413
 414     do j=1,m2
 415       do i=1,m1
 416         arrout(i,j) = arrtmp(iim(i),j)*wgts1(i) + arrtmp(iip(i),j)*wgtn1(i)
 417       end do
 418     end do
 419
 420   end subroutine lininterp2d2d
 421   subroutine lininterp2d1d(arrin, n1, n2, arrout, m1, wgt1, wgt2, fldname)
 422     implicit none
 423     !-----------------------------------------------------------------------
 424     !
 425     ! Arguments
 426     !
 427     integer, intent(in) :: n1, n2, m1
 428     real(r8), intent(in) :: arrin(n1,n2)    ! input array of values to interpolate
 429     type(interp_type), intent(in) :: wgt1, wgt2
 430     real(r8), intent(out) :: arrout(m1) ! interpolated array
 431     character(len=*), intent(in), optional :: fldname(:)
 432     !
 433     ! locals
 434     !
 435     integer i                ! indices
 436     integer, pointer :: iim(:), jjm(:)
 437     integer, pointer :: iip(:), jjp(:)
 438
 439     real(r8), pointer :: wgts(:), wgte(:)
 440     real(r8), pointer :: wgtn(:), wgtw(:)
 441
 442     jjm => wgt2%jjm
 443     jjp => wgt2%jjp
 444     wgts => wgt2%wgts
 445     wgtn => wgt2%wgtn
 446
 447     iim => wgt1%jjm
 448     iip => wgt1%jjp
 449     wgtw => wgt1%wgts
 450     wgte => wgt1%wgtn
 451
 452     do i=1,m1
 453        arrout(i) = arrin(iim(i),jjm(i))*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i))*wgte(i)*wgts(i) + &
 454                    arrin(iim(i),jjp(i))*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i))*wgte(i)*wgtn(i)
 455     end do
 456
 457
 458   end subroutine lininterp2d1d
 459   subroutine lininterp3d2d(arrin, n1, n2, n3, arrout, m1, len1, wgt1, wgt2)
 460     implicit none
 461     !-----------------------------------------------------------------------
 462     !
 463     ! Arguments
 464     !
 465     integer, intent(in) :: n1, n2, n3, m1, len1   ! m1 is to len1 as ncols is to pcols
 466     real(r8), intent(in) :: arrin(n1,n2,n3)    ! input array of values to interpolate
 467     type(interp_type), intent(in) :: wgt1, wgt2
 468     real(r8), intent(out) :: arrout(len1, n3) ! interpolated array
 469
 470     !
 471     ! locals
 472     !
 473     integer i, k               ! indices
 474     integer, pointer :: iim(:), jjm(:)
 475     integer, pointer :: iip(:), jjp(:)
 476
 477     real(r8), pointer :: wgts(:), wgte(:)
 478     real(r8), pointer :: wgtn(:), wgtw(:)
 479
 480     jjm => wgt2%jjm
 481     jjp => wgt2%jjp
 482     wgts => wgt2%wgts
 483     wgtn => wgt2%wgtn
 484
 485     iim => wgt1%jjm
 486     iip => wgt1%jjp
 487     wgtw => wgt1%wgts
 488     wgte => wgt1%wgtn
 489
 490     do k=1,n3
 491        do i=1,m1
 492           arrout(i,k) = arrin(iim(i),jjm(i),k)*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i),k)*wgte(i)*wgts(i) + &
 493                arrin(iim(i),jjp(i),k)*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i),k)*wgte(i)*wgtn(i)
 494        end do
 495     end do
 496
 497   end subroutine lininterp3d2d
 498
 499
 500
 501
 502   subroutine lininterp_finish(interp_wgts)
 503     type(interp_type) :: interp_wgts
 504
 505     deallocate(interp_wgts%jjm, &
 506          interp_wgts%jjp, &
 507          interp_wgts%wgts, &
 508          interp_wgts%wgtn)
 509
 510     nullify(interp_wgts%jjm, &
 511          interp_wgts%jjp, &
 512          interp_wgts%wgts, &
 513          interp_wgts%wgtn)
 514   end subroutine lininterp_finish
 515
 516   subroutine lininterp_original (arrin, yin, nlev, nlatin, arrout, &
 517        yout, nlatout)
 518     !-----------------------------------------------------------------------
 519     !
 520     ! Purpose: Do a linear interpolation from input mesh defined by yin to output
 521     !          mesh defined by yout.  Where extrapolation is necessary, values will
 522     !          be copied from the extreme edge of the input grid.  Vectorization is over
 523     !          the vertical (nlev) dimension.
 524     !
 525     ! Method: Check validity of input, then determine weights, then do the N-S interpolation.
 526     !
 527     ! Author: Jim Rosinski
 528     ! Modified: Jim Edwards so that there is no requirement of monotonacity on the yout array
 529     !
 530     !-----------------------------------------------------------------------
 531     implicit none
 532     !-----------------------------------------------------------------------
 533     !
 534     ! Arguments
 535     !
 536     integer, intent(in) :: nlev                   ! number of vertical levels
 537     integer, intent(in) :: nlatin                 ! number of input latitudes
 538     integer, intent(in) :: nlatout                ! number of output latitudes
 539
 540     real(r8), intent(in) :: arrin(nlev,nlatin)    ! input array of values to interpolate
 541     real(r8), intent(in) :: yin(nlatin)           ! input mesh
 542     real(r8), intent(in) :: yout(nlatout)         ! output mesh
 543
 544     real(r8), intent(out) :: arrout(nlev,nlatout) ! interpolated array
 545     !
 546     ! Local workspace
 547     !
 548     integer j, jj              ! latitude indices
 549     integer jjprev             ! latitude indices
 550     integer k                  ! level index
 551     integer icount             ! number of values
 552
 553     real(r8) extrap            ! percent grid non-overlap
 554     !
 555     ! Dynamic
 556     !
 557     integer :: jjm(nlatout)
 558     integer :: jjp(nlatout)
 559
 560     real(r8) :: wgts(nlatout)
 561     real(r8) :: wgtn(nlatout)
 562     !
 563     ! Check validity of input coordinate arrays: must be monotonically increasing,
 564     ! and have a total of at least 2 elements
 565     !
 566     if (nlatin.lt.2) then
 567        call wrf_abort  !('LININTERP: Must have at least 2 input points for interpolation')
 568     end if
 569
 570     icount = 0
 571     do j=1,nlatin-1
 572        if (yin(j).gt.yin(j+1)) icount = icount + 1
 573     end do
 574
 575
 576     if (icount.gt.0) then
 577        call wrf_abort  !('LININTERP: Non-monotonic coordinate array(s) found')
 578     end if
 579     !
 580     ! Initialize index arrays for later checking
 581     !
 582     do j=1,nlatout
 583        jjm(j) = 0
 584        jjp(j) = 0
 585     end do
 586     !
 587     ! For values which extend beyond N and S boundaries, set weights
 588     ! such that values will just be copied.
 589     !
 590     extrap = 0.
 591
 592     do j=1,nlatout
 593        if (yout(j).le.yin(1)) then
 594           jjm(j) = 1
 595           jjp(j) = 1
 596           wgts(j) = 1.
 597           wgtn(j) = 0.
 598           extrap=extrap+1.
 599        else if (yout(j).gt.yin(nlatin)) then
 600           jjm(j) = nlatin
 601           jjp(j) = nlatin
 602           wgts(j) = 1.
 603           wgtn(j) = 0.
 604           extrap=extrap+1.
 605        endif
 606     end do
 607
 608     !
 609     ! Loop though output indices finding input indices and weights
 610     !
 611     do j=1,nlatout
 612        do jj=1,nlatin-1
 613           if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then
 614              jjm(j) = jj
 615              jjp(j) = jj + 1
 616              wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj))
 617              wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj))
 618              exit
 619           end if
 620        end do
 621     end do
 622     !
 623     ! Check that interp/extrap points have been found for all outputs
 624     !
 625     icount = 0
 626     do j=1,nlatout
 627        if (jjm(j).eq.0 .or. jjp(j).eq.0) then
 628           icount = icount + 1
 629        end if
 630     end do
 631     if (icount.gt.0) then
 632        call wrf_abort  !('LININTERP: Point found without interp indices')
 633     end if
 634     !
 635     ! Do the interpolation
 636     !
 637     do j=1,nlatout
 638        do k=1,nlev
 639           arrout(k,j) = arrin(k,jjm(j))*wgts(j) + arrin(k,jjp(j))*wgtn(j)
 640        end do
 641     end do
 642
 643     return
 644   end subroutine lininterp_original
 645
 646
 647   subroutine bilin (arrin, xin, yin, nlondin, nlonin, &
 648        nlevdin, nlev, nlatin, arrout, xout, &
 649        yout, nlondout, nlonout, nlevdout, nlatout)
 650     !-----------------------------------------------------------------------
 651     !
 652     ! Purpose:
 653     !
 654     ! Do a bilinear interpolation from input mesh defined by xin, yin to output
 655     ! mesh defined by xout, yout.  Circularity is assumed in the x-direction so
 656     ! input x-grid must be in degrees east and must start from Greenwich.  When
 657     ! extrapolation is necessary in the N-S direction, values will be copied
 658     ! from the extreme edge of the input grid.  Vectorization is over the
 659     ! longitude dimension.  Input grid is assumed rectangular. Output grid
 660     ! is assumed ragged, i.e. xout is a 2-d mesh.
 661     !
 662     ! Method: Interpolate first in longitude, then in latitude.
 663     !
 664     ! Author: Jim Rosinski
 665     !
 666     !-----------------------------------------------------------------------
 667 !   use shr_kind_mod, only: r8 => shr_kind_r8
 668 !   use abortutils,   only: endrun
 669     !-----------------------------------------------------------------------
 670     implicit none
 671     !-----------------------------------------------------------------------
 672     integer,parameter :: r8 = selected_real_kind(12) ! 8 byte real
 673     !
 674     ! Input arguments
 675     !
 676     integer, intent(in) :: nlondin                        ! longitude dimension of input grid
 677     integer, intent(in) :: nlonin                         ! number of real longitudes (input)
 678     integer, intent(in) :: nlevdin                        ! vertical dimension of input grid
 679     integer, intent(in) :: nlev                           ! number of vertical levels
 680     integer, intent(in) :: nlatin                         ! number of input latitudes
 681     integer, intent(in) :: nlatout                        ! number of output latitudes
 682     integer, intent(in) :: nlondout                       ! longitude dimension of output grid
 683     integer, intent(in) :: nlonout(nlatout)               ! number of output longitudes per lat
 684     integer, intent(in) :: nlevdout                       ! vertical dimension of output grid
 685
 686     real(r8), intent(in) :: arrin(nlondin,nlevdin,nlatin) ! input array of values to interpolate
 687     real(r8), intent(in) :: xin(nlondin)                  ! input x mesh
 688     real(r8), intent(in) :: yin(nlatin)                   ! input y mesh
 689     real(r8), intent(in) :: xout(nlondout,nlatout)        ! output x mesh
 690     real(r8), intent(in) :: yout(nlatout)                 ! output y mesh
 691     !
 692     ! Output arguments
 693     !
 694     real(r8), intent(out) :: arrout(nlondout,nlevdout,nlatout) ! interpolated array
 695     !
 696     ! Local workspace
 697     !
 698     integer :: i, ii, iw, ie, iiprev ! longitude indices
 699     integer :: j, jj, js, jn, jjprev ! latitude indices
 700     integer :: k                     ! level index
 701     integer :: icount                ! number of bad values
 702
 703     real(r8) :: extrap               ! percent grid non-overlap
 704     real(r8) :: dxinwrap             ! delta-x on input grid for 2-pi
 705     real(r8) :: avgdxin              ! avg input delta-x
 706     real(r8) :: ratio                ! compare dxinwrap to avgdxin
 707     real(r8) :: sum                  ! sum of weights (used for testing)
 708     !
 709     ! Dynamic
 710     !
 711     integer :: iim(nlondout)         ! interpolation index to the left
 712     integer :: iip(nlondout)         ! interpolation index to the right
 713     integer :: jjm(nlatout)          ! interpolation index to the south
 714     integer :: jjp(nlatout)          ! interpolation index to the north
 715
 716     real(r8) :: wgts(nlatout)        ! interpolation weight to the north
 717     real(r8) :: wgtn(nlatout)        ! interpolation weight to the north
 718     real(r8) :: wgte(nlondout)       ! interpolation weight to the north
 719     real(r8) :: wgtw(nlondout)       ! interpolation weight to the north
 720     real(r8) :: igrid(nlonin)        ! interpolation weight to the north
 721
 722     character(len=132) :: message
 723     !
 724     ! Check validity of input coordinate arrays: must be monotonically increasing,
 725     ! and have a total of at least 2 elements
 726     !
 727     if (nlonin < 2 .or. nlatin < 2) then
 728        call wrf_abort  ! ('BILIN: Must have at least 2 input points for interpolation')
 729     end if
 730
 731     if (xin(1) < 0._r8 .or. xin(nlonin) > 360._r8) then
 732        call wrf_abort  ! ('BILIN: Input x-grid must be between 0 and 360')
 733     end if
 734
 735     icount = 0
 736     do i=1,nlonin-1
 737        if (xin(i) >= xin(i+1)) icount = icount + 1
 738     end do
 739
 740     do j=1,nlatin-1
 741        if (yin(j) >= yin(j+1)) icount = icount + 1
 742     end do
 743
 744     do j=1,nlatout-1
 745        if (yout(j) >= yout(j+1)) icount = icount + 1
 746     end do
 747
 748     do j=1,nlatout
 749        do i=1,nlonout(j)-1
 750           if (xout(i,j) >= xout(i+1,j)) icount = icount + 1
 751        end do
 752     end do
 753
 754     if (icount > 0) then
 755        call wrf_abort  ! ('BILIN: Non-monotonic coordinate array(s) found')
 756     end if
 757
 758     if (yout(nlatout) <= yin(1) .or. yout(1) >= yin(nlatin)) then
 759        call wrf_abort  ! ('BILIN: No overlap in y-coordinates')
 760     end if
 761
 762     do j=1,nlatout
 763        if (xout(1,j) < 0._r8 .or. xout(nlonout(j),j) > 360._r8) then
 764           call wrf_abort  ! ('BILIN: Output x-grid must be between 0 and 360')
 765        end if
 766
 767        if (xout(nlonout(j),j) <= xin(1) .or.  &
 768             xout(1,j)          >= xin(nlonin)) then
 769           call wrf_abort  ! ('BILIN: No overlap in x-coordinates')
 770        end if
 771     end do
 772     !
 773     ! Initialize index arrays for later checking
 774     !
 775     do j=1,nlatout
 776        jjm(j) = 0
 777        jjp(j) = 0
 778     end do
 779     !
 780     ! For values which extend beyond N and S boundaries, set weights
 781     ! such that values will just be copied.
 782     !
 783     do js=1,nlatout
 784        if (yout(js) > yin(1)) exit
 785        jjm(js) = 1
 786        jjp(js) = 1
 787        wgts(js) = 1._r8
 788        wgtn(js) = 0._r8
 789     end do
 790
 791     do jn=nlatout,1,-1
 792        if (yout(jn) <= yin(nlatin)) exit
 793        jjm(jn) = nlatin
 794        jjp(jn) = nlatin
 795        wgts(jn) = 1._r8
 796        wgtn(jn) = 0._r8
 797     end do
 798     !
 799     ! Loop though output indices finding input indices and weights
 800     !
 801     jjprev = 1
 802     do j=js,jn
 803        do jj=jjprev,nlatin-1
 804           if (yout(j) > yin(jj) .and. yout(j) <= yin(jj+1)) then
 805              jjm(j) = jj
 806              jjp(j) = jj + 1
 807              wgts(j) = (yin(jj+1) - yout(j)) / (yin(jj+1) - yin(jj))
 808              wgtn(j) = (yout(j)   - yin(jj)) / (yin(jj+1) - yin(jj))
 809              goto 30
 810           end if
 811        end do
 812        call wrf_abort  ! ('BILIN: Failed to find interp values')
 813 30     jjprev = jj
 814     end do
 815
 816     dxinwrap = xin(1) + 360._r8 - xin(nlonin)
 817     !
 818     ! Check for sane dxinwrap values.  Allow to differ no more than 10% from avg
 819     !
 820     avgdxin = (xin(nlonin)-xin(1))/(nlonin-1._r8)
 821     ratio = dxinwrap/avgdxin
 822     if (ratio < 0.9_r8 .or. ratio > 1.1_r8) then
 823        write(message,*)'BILIN: Insane dxinwrap value =',dxinwrap,' avg=', avgdxin
 824        call wrf_message( trim(message) )
 825        call wrf_abort  !
 826     end if
 827     !
 828     ! Check grid overlap
 829     ! Do not do on spmd since distributed output grid may be expected to fail this test
 830 #ifndef SPMD
 831     extrap = 100._r8*((js - 1._r8) + real(nlatout - jn,r8))/nlatout
 832     if (extrap > 20._r8) then
 833        write(message,*)'BILIN:',extrap,' % of N/S output grid will have to be extrapolated'
 834        call wrf_message( trim(message) )
 835     end if
 836 #endif
 837     !
 838     ! Check that interp/extrap points have been found for all outputs, and that
 839     ! they are reasonable (i.e. within 32-bit roundoff)
 840     !
 841     icount = 0
 842     do j=1,nlatout
 843        if (jjm(j) == 0 .or. jjp(j) == 0) icount = icount + 1
 844        sum = wgts(j) + wgtn(j)
 845        if (sum < 0.99999_r8 .or. sum > 1.00001_r8) icount = icount + 1
 846        if (wgts(j) < 0._r8 .or. wgts(j) > 1._r8) icount = icount + 1
 847        if (wgtn(j) < 0._r8 .or. wgtn(j) > 1._r8) icount = icount + 1
 848     end do
 849
 850     if (icount > 0) then
 851        call wrf_abort  ! ('BILIN: something bad in latitude indices or weights')
 852     end if
 853     !
 854     ! Do the bilinear interpolation
 855     !
 856     do j=1,nlatout
 857        !
 858        ! Initialize index arrays for later checking
 859        !
 860        do i=1,nlondout
 861           iim(i) = 0
 862           iip(i) = 0
 863        end do
 864        !
 865        ! For values which extend beyond E and W boundaries, set weights such that
 866        ! values will be interpolated between E and W edges of input grid.
 867        !
 868        do iw=1,nlonout(j)
 869           if (xout(iw,j) > xin(1)) exit
 870           iim(iw) = nlonin
 871           iip(iw) = 1
 872           wgtw(iw) = (xin(1)        - xout(iw,j))   /dxinwrap
 873           wgte(iw) = (xout(iw,j)+360._r8 - xin(nlonin))/dxinwrap
 874        end do
 875
 876        do ie=nlonout(j),1,-1
 877           if (xout(ie,j) <= xin(nlonin)) exit
 878           iim(ie) = nlonin
 879           iip(ie) = 1
 880           wgtw(ie) = (xin(1)+360._r8 - xout(ie,j))   /dxinwrap
 881           wgte(ie) = (xout(ie,j)    - xin(nlonin))/dxinwrap
 882        end do
 883        !
 884        ! Loop though output indices finding input indices and weights
 885        !
 886        iiprev = 1
 887        do i=iw,ie
 888           do ii=iiprev,nlonin-1
 889              if (xout(i,j) > xin(ii) .and. xout(i,j) <= xin(ii+1)) then
 890                 iim(i) = ii
 891                 iip(i) = ii + 1
 892                 wgtw(i) = (xin(ii+1) - xout(i,j)) / (xin(ii+1) - xin(ii))
 893                 wgte(i) = (xout(i,j) - xin(ii))   / (xin(ii+1) - xin(ii))
 894                 goto 60
 895              end if
 896           end do
 897           call wrf_abort  ! ('BILIN: Failed to find interp values')
 898 60        iiprev = ii
 899        end do
 900
 901        icount = 0
 902        do i=1,nlonout(j)
 903           if (iim(i) == 0 .or. iip(i) == 0) icount = icount + 1
 904           sum = wgtw(i) + wgte(i)
 905           if (sum < 0.99999_r8 .or. sum > 1.00001_r8) icount = icount + 1
 906           if (wgtw(i) < 0._r8 .or. wgtw(i) > 1._r8) icount = icount + 1
 907           if (wgte(i) < 0._r8 .or. wgte(i) > 1._r8) icount = icount + 1
 908        end do
 909
 910        if (icount > 0) then
 911           write(message,*)'BILIN: j=',j,' Something bad in longitude indices or weights'
 912           call wrf_message( trim(message) )
 913           call wrf_abort  !
 914        end if
 915        !
 916        ! Do the interpolation, 1st in longitude then latitude
 917        !
 918        do k=1,nlev
 919           do i=1,nlonin
 920              igrid(i) = arrin(i,k,jjm(j))*wgts(j) + arrin(i,k,jjp(j))*wgtn(j)
 921           end do
 922
 923           do i=1,nlonout(j)
 924              arrout(i,k,j) = igrid(iim(i))*wgtw(i) + igrid(iip(i))*wgte(i)
 925           end do
 926        end do
 927     end do
 928
 929
 930     return
 931   end subroutine bilin
 932
 933   subroutine vertinterp(ncol, ncold, nlev, pmid, pout, arrin, arrout)
 934
 935     !-----------------------------------------------------------------------
 936     !
 937     ! Purpose:
 938     ! Vertically interpolate input array to output pressure level
 939     ! Copy values at boundaries.
 940     !
 941     ! Method:
 942     !
 943     ! Author:
 944     !
 945     !-----------------------------------------------------------------------
 946
 947     implicit none
 948
 949     !------------------------------Arguments--------------------------------
 950     integer , intent(in)  :: ncol              ! column dimension
 951     integer , intent(in)  :: ncold             ! declared column dimension
 952     integer , intent(in)  :: nlev              ! vertical dimension
 953     real(r8), intent(in)  :: pmid(ncold,nlev)  ! input level pressure levels
 954     real(r8), intent(in)  :: pout              ! output pressure level
 955     real(r8), intent(in)  :: arrin(ncold,nlev) ! input  array
 956     real(r8), intent(out) :: arrout(ncold)     ! output array (interpolated)
 957     !--------------------------------------------------------------------------
 958
 959     !---------------------------Local variables-----------------------------
 960     integer i,k               ! indices
 961     integer kupper(ncold)     ! Level indices for interpolation
 962     real(r8) dpu              ! upper level pressure difference
 963     real(r8) dpl              ! lower level pressure difference
 964     logical found(ncold)      ! true if input levels found
 965     logical error             ! error flag
 966     !-----------------------------------------------------------------
 967     !
 968     ! Initialize index array and logical flags
 969     !
 970     do i=1,ncol
 971        found(i)  = .false.
 972        kupper(i) = 1
 973     end do
 974     error = .false.
 975     !
 976     ! Store level indices for interpolation.
 977     ! If all indices for this level have been found,
 978     ! do the interpolation
 979     !
 980     do k=1,nlev-1
 981        do i=1,ncol
 982           if ((.not. found(i)) .and. pmid(i,k)<pout .and. pout<=pmid(i,k+1)) then
 983              found(i) = .true.
 984              kupper(i) = k
 985           end if
 986        end do
 987     end do
 988     !
 989     ! If we've fallen through the k=1,nlev-1 loop, we cannot interpolate and
 990     ! must extrapolate from the bottom or top data level for at least some
 991     ! of the longitude points.
 992     !
 993     do i=1,ncol
 994        if (pout <= pmid(i,1)) then
 995           arrout(i) = arrin(i,1)
 996        else if (pout >= pmid(i,nlev)) then
 997           arrout(i) = arrin(i,nlev)
 998        else if (found(i)) then
 999           dpu = pout - pmid(i,kupper(i))
1000           dpl = pmid(i,kupper(i)+1) - pout
1001           arrout(i) = (arrin(i,kupper(i)  )*dpl + arrin(i,kupper(i)+1)*dpu)/(dpl + dpu)
1002        else
1003           error = .true.
1004        end if
1005     end do
1006     !
1007     ! Error check
1008     !
1009     if (error) then
1010        call wrf_abort  ! ('VERTINTERP: ERROR FLAG')
1011     end if
1012
1013     return
1014   end subroutine vertinterp
1015
1016   subroutine get_timeinterp_factors (cycflag, np1, cdayminus, cdayplus, cday, &
1017        fact1, fact2, str)
1018     !---------------------------------------------------------------------------
1019     !
1020     ! Purpose: Determine time interpolation factors (normally for a boundary dataset)
1021     !          for linear interpolation.
1022     !
1023     ! Method:  Assume 365 days per year.  Output variable fact1 will be the weight to
1024     !          apply to data at calendar time "cdayminus", and fact2 the weight to apply
1025     !          to data at time "cdayplus".  Combining these values will produce a result
1026     !          valid at time "cday".  Output arguments fact1 and fact2 will be between
1027     !          0 and 1, and fact1 + fact2 = 1 to roundoff.
1028     !
1029     ! Author:  Jim Rosinski
1030     !
1031     !---------------------------------------------------------------------------
1032     implicit none
1033     !
1034     ! Arguments
1035     !
1036     logical, intent(in) :: cycflag             ! flag indicates whether dataset is being cycled yearly
1037
1038     integer, intent(in) :: np1                 ! index points to forward time slice matching cdayplus
1039
1040     real(r8), intent(in) :: cdayminus          ! calendar day of rearward time slice
1041     real(r8), intent(in) :: cdayplus           ! calendar day of forward time slice
1042     real(r8), intent(in) :: cday               ! calenar day to be interpolated to
1043     real(r8), intent(out) :: fact1             ! time interpolation factor to apply to rearward time slice
1044     real(r8), intent(out) :: fact2             ! time interpolation factor to apply to forward time slice
1045
1046     character(len=*), intent(in) :: str        ! string to be added to print in case of error (normally the callers name)
1047     !
1048     ! Local workspace
1049     !
1050     real(r8) :: deltat                         ! time difference (days) between cdayminus and cdayplus
1051     real(r8), parameter :: daysperyear = 365.  ! number of days in a year
1052
1053     character(len=132) :: message
1054     !
1055     ! Initial sanity checks
1056     !
1057     if (np1 == 1 .and. .not. cycflag) then
1058        call wrf_abort  ! ('GETFACTORS:'//str//' cycflag false and forward month index = Jan. not allowed')
1059     end if
1060
1061     if (np1 < 1) then
1062        call wrf_abort  ! ('GETFACTORS:'//str//' input arg np1 must be > 0')
1063     end if
1064
1065     if (cycflag) then
1066        if ((cday < 1.) .or. (cday > (daysperyear+1.))) then
1067           write(message,*) 'GETFACTORS:', str, ' bad cday=',cday
1068           call wrf_message( trim(message) )
1069           call wrf_abort  ! ()
1070        end if
1071     else
1072        if (cday < 1.) then
1073           write(message,*) 'GETFACTORS:', str, ' bad cday=',cday
1074           call wrf_message( trim(message) )
1075           call wrf_abort  ! ()
1076        end if
1077     end if
1078     !
1079     ! Determine time interpolation factors.  Account for December-January
1080     ! interpolation if dataset is being cycled yearly.
1081     !
1082     if (cycflag .and. np1 == 1) then                     ! Dec-Jan interpolation
1083        deltat = cdayplus + daysperyear - cdayminus
1084        if (cday > cdayplus) then                         ! We are in December
1085           fact1 = (cdayplus + daysperyear - cday)/deltat
1086           fact2 = (cday - cdayminus)/deltat
1087        else                                              ! We are in January
1088           fact1 = (cdayplus - cday)/deltat
1089           fact2 = (cday + daysperyear - cdayminus)/deltat
1090        end if
1091     else
1092        deltat = cdayplus - cdayminus
1093        fact1 = (cdayplus - cday)/deltat
1094        fact2 = (cday - cdayminus)/deltat
1095     end if
1096
1097     if (.not. valid_timeinterp_factors (fact1, fact2)) then
1098        write(message,*) 'GETFACTORS: ', str, ' bad fact1 and/or fact2=', fact1, fact2
1099        call wrf_message( trim(message) )
1100        call wrf_abort  ! ()
1101     end if
1102
1103     return
1104   end subroutine get_timeinterp_factors
1105
1106   logical function valid_timeinterp_factors (fact1, fact2)
1107     !---------------------------------------------------------------------------
1108     !
1109     ! Purpose: check sanity of time interpolation factors to within 32-bit roundoff
1110     !
1111     !---------------------------------------------------------------------------
1112     implicit none
1113
1114     real(r8), intent(in) :: fact1, fact2           ! time interpolation factors
1115
1116     valid_timeinterp_factors = .true.
1117
1118     ! The fact1 .ne. fact1 and fact2 .ne. fact2 comparisons are to detect NaNs.
1119     if (abs(fact1+fact2-1.) > 1.e-6 .or. &
1120          fact1 > 1.000001 .or. fact1 < -1.e-6 .or. &
1121          fact2 > 1.000001 .or. fact2 < -1.e-6 .or. &
1122          fact1 .ne. fact1 .or. fact2 .ne. fact2) then
1123
1124        valid_timeinterp_factors = .false.
1125     end if
1126
1127     return
1128   end function valid_timeinterp_factors
1129
1130 end module module_interpolate