!### 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