UPP_MATH.f

Go to the documentation of this file.

  module upp_math

  use masks,        only: dx, dy
  use ctlblk_mod,   only: im, jsta_2l, jend_2u, jsta_m, jend_m, spval,&
                              ista_2l, iend_2u, ista_m, iend_m
  use gridspec_mod, only: gridtype
!
  implicit none

  private

  public :: ddvdx, ddudy, uuavg
  public :: dvdxdudy
  public :: h2u, h2v, u2h, v2h

  REAL, ALLOCATABLE :: ddvdx(:,:)
  REAL, ALLOCATABLE :: ddudy(:,:)
  REAL, ALLOCATABLE :: uuavg(:,:)
! 
!-------------------------------------------------------------------------------------
!
  contains
!
!-------------------------------------------------------------------------------------
  subroutine dvdxdudy(uwnd,vwnd)
!
      implicit none

      REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u), intent(in)    :: uwnd, vwnd
!
!-- local variables
!--
      integer i, j
      real r2dx, r2dy
      INTEGER, allocatable ::  ihe(:),ihw(:)
!
      IF(gridtype == 'A')THEN
!$omp parallel do  private(i,j,r2dx,r2dy)
        DO j=jsta_m,jend_m
        DO i=ista_m,iend_m
           IF(vwnd(i+1,j)<spval.AND.vwnd(i-1,j)<spval.AND.              &
     &        uwnd(i,j+1)<spval.AND.uwnd(i,j-1)<spval) THEN
              r2dx   = 1./(2.*dx(i,j))
              r2dy   = 1./(2.*dy(i,j))
              ddvdx(i,j)   = (vwnd(i+1,j)-vwnd(i-1,j))*r2dx
              ddudy(i,j)   = (uwnd(i,j+1)-uwnd(i,j-1))*r2dy
              uuavg(i,j)   = uwnd(i,j)
           END IF
        END DO
        END DO
      ELSE IF (gridtype == 'E')THEN
        allocate(ihw(jsta_2l:jend_2u), ihe(jsta_2l:jend_2u))
!$omp  parallel do private(j)
        DO j=jsta_2l,jend_2u
          ihw(j) = -mod(j,2)
          ihe(j) = ihw(j)+1
        ENDDO
!$omp parallel do  private(i,j,r2dx,r2dy)
        DO j=jsta_m,jend_m
          DO i=ista_m,iend_m
            IF(vwnd(i+ihe(j),j) < spval.AND.vwnd(i+ihw(j),j) < spval .AND.   &
     &         uwnd(i,j+1) < spval     .AND.uwnd(i,j-1) < spval) THEN
               r2dx   = 1./(2.*dx(i,j))
               r2dy   = 1./(2.*dy(i,j))
               ddvdx(i,j)   = (vwnd(i+ihe(j),j)-vwnd(i+ihw(j),j))*r2dx
               ddudy(i,j)   = (uwnd(i,j+1)-uwnd(i,j-1))*r2dy
               uuavg(i,j)   = 0.25*(uwnd(i+ihe(j),j)+uwnd(i+ihw(j),j)               &
     &                      +       uwnd(i,j+1)+uwnd(i,j-1))
            END IF
          END DO
        END DO
        deallocate(ihw, ihe)
      ELSE IF (gridtype == 'B')THEN
!$omp parallel do  private(i,j,r2dx,r2dy)
        DO j=jsta_m,jend_m
          DO i=ista_m,iend_m
            r2dx = 1./dx(i,j)
            r2dy = 1./dy(i,j)
            if(vwnd(i,  j)==spval .or. vwnd(i,  j-1)==spval .or. &
               vwnd(i-1,j)==spval .or. vwnd(i-1,j-1)==spval .or. &
               uwnd(i,  j)==spval .or. uwnd(i-1,j  )==spval .or. &
               uwnd(i,j-1)==spval .or. uwnd(i-1,j-1)==spval) cycle
            ddvdx(i,j) = (0.5*(vwnd(i,j)+vwnd(i,j-1))-0.5*(vwnd(i-1,j) &
     &           +       vwnd(i-1,j-1)))*r2dx
            ddudy(i,j) = (0.5*(uwnd(i,j)+uwnd(i-1,j))-0.5*(uwnd(i,j-1) &
     &           +       uwnd(i-1,j-1)))*r2dy
            uuavg(i,j) = 0.25*(uwnd(i-1,j-1)+uwnd(i-1,j)               &
     &           +       uwnd(i,  j-1)+uwnd(i,  j))
          END DO
        END DO
      END IF

  end subroutine dvdxdudy
!
!-------------------------------------------------------------------------------------
!
      subroutine h2u(ingrid,outgrid)
! This subroutine interpolates from H points onto U points
! Author: CHUANG, EMC, Dec. 2010

      use ctlblk_mod, only: spval, jsta, jend, jsta_m, jend_m, me, num_procs, jm,&
              im, jsta_2l, jend_2u, ista, iend, ista_m, iend_m, ista_2l, iend_2u 
      use gridspec_mod, only: gridtype

      implicit none

      include "mpif.h"
      integer:: i,j,ie,iw
      real,dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),intent(in)::ingrid
      real,dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),intent(out)::outgrid
      outgrid=spval
      if(gridtype == 'A')THEN
       do j=jsta,jend
        do i=ista,iend
     outgrid(i,j)=ingrid(i,j)
    end do
       end do
      else IF(gridtype == 'E')THEN
       call exch(ingrid(ista_2l,jsta_2l))
       DO j=jsta_m,jend_m
        DO i=ista_m,iend_m
     ie=i+mod(j,2)
         iw=ie-1
         outgrid(i,j)=(ingrid(iw,j)+ingrid(ie,j)+ingrid(i,j+1)+ingrid(i,j-1))/4.0
    end do
       end do
      ELSE IF(gridtype == 'B')THEN
       call exch(ingrid(ista_2l,jsta_2l))
       DO j=jsta,jend_m
        DO i=ista,iend_m
     outgrid(i,j)=(ingrid(i,j)+ingrid(i,j+1)+ingrid(i+1,j)+ingrid(i+1,j+1))/4.0
    end do
       end do
! Fill in boundary points because hysplit fails when 10 m wind has bitmaps
       do j=jsta,jend_m
        outgrid(iend,j)=outgrid(iend-1,j)
       end do
       IF(me == (num_procs-1) .and. jend_2u >= jm) then
        DO i=ista,iend
         outgrid(i,jend) = outgrid(i,jend-1)
        END DO
       END IF      
      ELSE IF(gridtype == 'C')THEN
       DO j=jsta,jend
        DO i=ista,iend_m
          outgrid(i,j)=(ingrid(i,j)+ingrid(i+1,j))/2.0
        end do
       end do
      end if 
         
      end subroutine h2u
!
!-------------------------------------------------------------------------------------
!   
      subroutine h2v(ingrid,outgrid)
! This subroutine interpolates from H points onto V points
! Author: CHUANG, EMC, Dec. 2010
      use ctlblk_mod, only: spval, jsta, jend, jsta_m, jend_m, im, jsta_2l, jend_2u,&
                                   ista, iend, ista_m, iend_m,     ista_2l, iend_2u
      use gridspec_mod, only: gridtype
      implicit none
      include "mpif.h"
      integer:: i,j,ie,iw
      real,dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),intent(in)::ingrid
      real,dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),intent(out)::outgrid
      outgrid=spval
      if(gridtype == 'A')THEN
       do j=jsta,jend
        do i=ista,iend
          outgrid(i,j)=ingrid(i,j)
        end do
       end do
      else IF(gridtype == 'E')THEN
       call exch(ingrid(ista_2l,jsta_2l))
       DO j=jsta_m,jend_m
        DO i=ista_m,iend_m
     ie=i+mod(j,2)
         iw=ie-1
         outgrid(i,j)=(ingrid(iw,j)+ingrid(ie,j)+ingrid(i,j+1)+ingrid(i,j-1))/4.0
    end do
       end do
      ELSE IF(gridtype == 'B')THEN
       call exch(ingrid(ista_2l,jsta_2l))
       DO j=jsta,jend_m
        DO i=ista,iend_m
     outgrid(i,j)=(ingrid(i,j)+ingrid(i,j+1)+ingrid(i+1,j)+ingrid(i+1,j+1))/4.0
    end do
       end do
      ELSE IF(gridtype == 'C')THEN
       call exch(ingrid(ista_2l,jsta_2l))
       DO j=jsta,jend_m
        DO i=ista,iend
     outgrid(i,j)=(ingrid(i,j)+ingrid(i,j+1))/2.0
    end do
       end do 
      end if 
      
      end subroutine h2v
!
!-------------------------------------------------------------------------------------
!
      subroutine u2h(ingrid,outgrid)
! This subroutine interpolates from U points onto H points
! Author: CHUANG, EMC, Dec. 2010
      use ctlblk_mod, only: spval, jsta, jend, jsta_m, jend_m, im, jsta_2l, jend_2u,&
                                   ista, iend, ista_m, iend_m,     ista_2l, iend_2u
      use gridspec_mod, only: gridtype
      implicit none
      include "mpif.h"
      integer:: i,j,ie,iw
      real,dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),intent(in)::ingrid
      real,dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),intent(out)::outgrid
      outgrid=spval
      if(gridtype == 'A')THEN
       do j=jsta,jend
        do i=ista,iend
     outgrid(i,j)=ingrid(i,j)
    end do
       end do
      else IF(gridtype == 'E')THEN
       call exch(ingrid(ista_2l,jsta_2l))
       DO j=jsta_m,jend_m
        DO i=ista_m,iend_m
     ie=i+mod(j+1,2)
         iw=ie-1
         outgrid(i,j)=(ingrid(iw,j)+ingrid(ie,j)+ingrid(i,j+1)+ingrid(i,j-1))/4.0
    end do
       end do
      ELSE IF(gridtype == 'B')THEN
       call exch(ingrid(ista_2l,jsta_2l))
       DO j=jsta_m,jend_m
        DO i=ista_m,iend_m
     outgrid(i,j)=(ingrid(i-1,j-1)+ingrid(i,j-1)+ingrid(i-1,j)+ingrid(i,j))/4.0
    end do
       end do
      ELSE IF(gridtype == 'C')THEN
       DO j=jsta,jend
        DO i=ista_m,iend
     outgrid(i,j)=(ingrid(i-1,j)+ingrid(i,j))/2.0
    end do
       end do       
      end if 
         
      end subroutine u2h       
!
!-------------------------------------------------------------------------------------
!        
      subroutine v2h(ingrid,outgrid)
! This subroutine interpolates from V points onto H points
! Author: CHUANG, EMC, Dec. 2010
      use ctlblk_mod, only: spval, jsta, jend, jsta_m, jend_m, im, jsta_2l, jend_2u,&
                                   ista, iend, ista_m, iend_m,     ista_2l, iend_2u
      use gridspec_mod, only: gridtype
      implicit none
      include "mpif.h"
      integer:: i,j,ie,iw
      real,dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),intent(in)::ingrid
      real,dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),intent(out)::outgrid
      outgrid=spval
      if(gridtype == 'A')THEN
       do j=jsta,jend
        do i=ista,iend
     outgrid(i,j)=ingrid(i,j)
    end do
       end do
      else IF(gridtype == 'E')THEN
       call exch(ingrid(ista_2l,jsta_2l))
       DO j=jsta_m,jend_m
        DO i=ista_m,iend_m
     ie=i+mod(j,2)
         iw=ie-1
         outgrid(i,j)=(ingrid(iw,j)+ingrid(ie,j)+ingrid(i,j+1)+ingrid(i,j-1))/4.0
    end do
       end do
      ELSE IF(gridtype == 'B')THEN
       call exch(ingrid(ista_2l,jsta_2l))
       DO j=jsta_m,jend_m
        DO i=ista_m,iend_m
     outgrid(i,j)=(ingrid(i-1,j-1)+ingrid(i,j-1)+ingrid(i-1,j)+ingrid(i,j))/4.0
    end do
       end do
      ELSE IF(gridtype == 'C')THEN
       call exch(ingrid(ista_2l,jsta_2l))
       DO j=jsta_m,jend
        DO i=ista,iend
     outgrid(i,j)=(ingrid(i,j-1)+ingrid(i,j))/2.0
    end do
       end do 
      end if 
         
      end subroutine v2h                 
!
!-------------------------------------------------------------------------------------
!        
  end module upp_math