ecearth3/sources/tm5mp/base/advect_tools.F90

!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5-MP                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:     ADVECT_TOOLS
!
! !DESCRIPTION: 
!\\
!\\
! !INTERFACE: 
!
MODULE ADVECT_TOOLS
  !
  ! !USES:
  !
  USE GO,           ONLY : gol, goPr, goErr   ! General Object messaging
  USE TM5_DISTGRID, ONLY : DGRID, Get_DistGrid

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  PUBLIC :: ADVECT_M       !
  PUBLIC :: M2PHLB1        ! convert air mass to pressure
  !
  ! !REVISION HISTORY: 
  !   15 Nov 2011 - P. Le Sager - adapted to TM5-MP
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5-MP                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    M2PHLB1
  !
  ! !DESCRIPTION: Get pressure at box centers (phlb_dat%data) and at the
  !               surface (sp_dat%data) from air mass (m_dat%data). Only for
  !               region #1!
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE M2PHLB1
    ! 
    ! !USES:
    !
    use binas,       only : grav
    use dims           
    use global_data, only : region_dat
    use meteodata,   only : sp_dat, phlb_dat, m_dat
    use partools,    only : par_reduce, isRoot
    !
    ! !REVISION HISTORY: 
    !   mk december 1998
    !   15 Nov 2011 - P. Le Sager - adapted for TM5-MP
    !   
    ! !REMARKS:
    ! (1) Note: first, we calculate the surface pressure in p(im,jm), then
    !      using the hybrid coordinate system, calculate phlb(im,jm,lmp1)
    ! (2) Different from PressureToMass from meteo.F90, where phlb_dat%data is
    !      filled directly from m_dat%data (and not from at & bt), and where
    !      spm_dat%data (and not sp_dat%data) is filled.
    ! (3) Special routine for region 1. 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    real,dimension(:,:,:), pointer :: m, phlb
    real,dimension(:,:,:), pointer :: p
    real,dimension(:),     pointer :: dxyp
    integer :: i, j, l, i1, i2, j1, j2, halo, lmr, region, st
    real    :: maxdiff, maxval_one, maxdiff_all, maxval_all
    real    :: mnew
    real, allocatable, dimension(:,:) :: pold, phelp

    !--------
    ! start
    !--------
    region = 1
    CALL GET_DISTGRID( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

    halo=2  ! same as p
    allocate(  pold( i1-halo:i2+halo, j1-halo:j2+halo) )
    allocate( phelp( i1-halo:i2+halo, j1-halo:j2+halo) )

    dxyp => region_dat(region)%dxyp
    p    =>   sp_dat(region)%data
    m    =>    m_dat(region)%data
    phlb => phlb_dat(region)%data
    lmr = lm(region)


    pold = p(:,:,1)
    phelp = zero
    p = zero


    do l=  1, lmr
    do j= j1, j2
    do i= i1, i2
             phelp(i,j) = phelp(i,j) + m(i,j,l)*grav/dxyp(j)
    end do
    end do
    end do

    ! surface pressure
    p(:,:,1) = phelp

    
    ! pressure at half level boundaries:
    do  l = 1, lmr+1
       phlb(:,:,l) = at(l) + bt(l)*p(:,:,1)
    end do

    !--------
    ! check
    !--------
    if (okdebug) then

       maxdiff = 0.0

       do l =  1, lmr
       do j = j1, j2
       do i = i1, i2
                mnew=(phlb(i,j,l)-phlb(i,j,l+1))*dxyp(j)/grav
                maxdiff = max(maxdiff,abs(mnew-m(i,j,l))/mnew)
       end do
       end do
       end do

       maxval_one=maxval(abs( p(i1:i2,j1:j2,1)- pold(i1:i2,j1:j2) ) / pold(i1:i2,j1:j2) )

       call Par_Reduce( maxdiff, 'MAX', maxdiff_all, st )
       call Par_Reduce( maxval_one, 'MAX', maxval_all, st )

       if ( isRoot ) then
          write(gol,*) 'm2phlb1: maxdiff (mnew-m)/mnew :',maxdiff_all ; call goPr
          write(gol,*) 'm2phlb1: max pressure difference (%) :',maxval_all ; call goPr
       endif

    end if

    !--------
    ! clean
    !--------
    nullify(p)
    nullify(m)
    nullify(phlb)
    nullify(dxyp)
    deallocate(phelp,pold)

  END SUBROUTINE M2PHLB1
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5-MP                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  ADVECT_M
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine advect_m(region, ntimes)
    !
    ! !USES:
    !
    use binas,       only : grav
    use dims
    use global_data, only : region_dat, wind_dat
    use meteodata,   only : m_dat, sp_dat
    use partools,    only : isRoot
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in) :: region,ntimes
    !
    ! !REVISION HISTORY: 
    !   15 Nov 2011 - P. Le Sager - modified for mpi lat-lon domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(:,:,:), pointer     :: m,am,bm,cm
    real, dimension(:,:,:), pointer     :: p
    real, dimension(:)    , pointer     :: dxyp                                        
    real, dimension(:,:,:), allocatable :: mnew, msave
    real, dimension(:,:)  , allocatable :: pnew
    
    integer :: n,i,j,l, i1, i2, j1, j2, lmr
    real    :: ss,so

    ! start

    if(isRoot) then
       write(gol,*)  'advect_m: called for region ', region ; call goPr
       write(gol,*)  'advect_m: ntimes            ', ntimes ; call goPr
    endif

    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    lmr = lm(region)

    allocate( mnew( i1:i2, j1:j2, lmr))
    allocate(msave( i1:i2, j1:j2, lmr))
    allocate( pnew( i1:i2, j1:j2     ))

    m  => m_dat(region)%data
    am => wind_dat(region)%am
    bm => wind_dat(region)%bm
    cm => wind_dat(region)%cm
    p  => sp_dat(region)%data
    dxyp => region_dat(region)%dxyp

    msave = m(i1:i2, j1:j2, 1:lmr)
    mnew  = m(i1:i2, j1:j2, 1:lmr)
    
    do i=1,ntimes
       mnew = mnew   + revert*am(i1-1:i2-1,   i1:i2  , 1:lmr  ) &   !east   !ADJ re-inserted the revert CMK
                     - revert*am(  i1:i2  ,   i1:i2  , 1:lmr  ) &   !west
                     + revert*bm(  i1:i2  ,   i1:i2  , 1:lmr  ) &   !south
                     - revert*bm(  i1:i2  , i1+1:i2+1, 1:lmr  ) &   !north
                     + revert*cm(  i1:i2  ,   i1:i2  , 0:lmr-1) &   !lower
                     - revert*cm(  i1:i2  ,   i1:i2  , 1:lmr  )     !upper
    end do
    
    pnew = 0.0
    do l= 1, lmr
    do j=j1, j2
    do i=i1, i2
             pnew(i,j) = pnew(i,j) + mnew(i,j,l)*grav/dxyp(j)
    end do
    end do
    end do

    p(i1:i2, j1:j2, 1    ) = pnew   !return new pressure for testing
    m(i1:i2, j1:j2, 1:lmr) = mnew

    nullify(m)
    nullify(am)
    nullify(bm)
    nullify(cm)
    nullify(p)
    nullify(dxyp)

    deallocate(pnew)
    deallocate(mnew)
    deallocate(msave)

  end subroutine advect_m
  !EOC

end module advect_tools