ecearth3/sources/tm5mp/base/advect.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"
!
!#################################################################

MODULE ADVECT

  use GO, only : gol, goPr, goErr

  implicit none
  private

  public  ::  Advect_Init, Advect_Done
  public  ::  dynam0

  ! --- const --------------------------------------

  character(len=*), parameter  ::  mname = 'advect'


CONTAINS


  ! ================================================================
  
  
  subroutine Advect_Init( status )
  
    use dims , only : nregions
    use Meteo, only : Set
    use Meteodata, only : sp1_dat, sp2_dat, sp_dat
#ifdef with_prism
    use MeteoData, only : tsp_dat
#endif
    use MeteoData, only : mfu_dat, mfv_dat, mfw_dat
    use MeteoData, only :  pu_dat,  pv_dat,  pw_dat

    ! --- in/out --------------------------------
    
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Advect_Init'
    
    ! --- local --------------------------------
    
    integer          ::  n
    
    ! --- begin --------------------------------
    
    ! loop over all (zoom) regions:
    do n = 1, nregions
    
      ! surface pressures and mass fluxes for advection
      call Set(  sp_dat(n), status, used=.true. )
      call Set( sp1_dat(n), status, used=.true. )
      call Set( sp2_dat(n), status, used=.true. )
#ifdef with_prism
      call Set( tsp_dat(n), status, used=.true. )
#endif
      call Set( mfu_dat(n), status, used=.true. )
      call Set( mfv_dat(n), status, used=.true. )
      call Set( mfw_dat(n), status, used=.true. )
      call Set(  pu_dat(n), status, used=.true. )
      call Set(  pv_dat(n), status, used=.true. )
      call Set(  pw_dat(n), status, used=.true. )
    
    end do  ! regions
    
    ! ok
    status = 0
    
  end subroutine Advect_Init
  
  
  ! ***
  
  
  subroutine Advect_Done( status )
  
    ! --- in/out --------------------------------
    
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Advect_Done'

    ! --- begin --------------------------------
    
    ! nothing to be done

    ! ok
    status = 0
    
  end subroutine Advect_Done
  
  
  
  ! ================================================================
  

  !-----------------------------------------------------------------------
  !
  !****  dynam0          - calculates vertical massfluxes        v 9.1
  !
  !      programmed by           mh      mpi HH          1-oct-1991
  !
  !      purpose
  !      -------
  !      This sr is CALLed each time after new massfluxes are read in.
  !      Calculate the new air-masses in each grid box from the
  !      surface pressure.
  !      Calculate the amount of air moved in each advection substep,
  !      also in the vertical direction.
  !
  !      interface
  !      ---------
  !      CALL dynam0
  !
  !      method
  !      ------
  !      integrate air conservation equation in the vertical direction
  !      in order to obtain the vertical massfluxes.
  !
  !      externals
  !      ---------
  !      none
  !
  !      reference
  !      ---------       
  !      see manual
  !
  ! modified by mk (1999) for zoom version.
  !   25 Jan 2012 - P. Le Sager -  adapted for MPI lat-lon domain decomposition
  ! 
  !-----------------------------------------------------------------------

  SUBROUTINE DYNAM0(region, status)

    use dims,         only : im, jm, lm, ndyn, tref, zero, bt, xcyc
    use global_data,  only : wind_dat
    use meteodata  ,  only : pu_dat, pv_dat
    use tm5_distgrid, only : dgrid, Get_DistGrid, update_halo

    ! in/output
    integer,intent(in)  :: region
    integer,intent(out) :: status

    ! local
    real, dimension(:,:,:), pointer     :: pu,pv,am,bm,cm
    real, dimension(:,:),   allocatable :: pit
    real, dimension(:,:,:), allocatable :: sd
    real, dimension(:,:,:), allocatable :: conv_adv

    real    :: dtu, dtv, dtw
    integer :: i, j, l, lmr, i1, i2, j1, j2

    character(len=*), parameter ::  rname = mname//'/dynam0'

    ! start

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2 , J_STRT=j1, J_STOP=j2 )

    allocate(      pit(i1:i2, j1:j2              ) )
    allocate(       sd(i1:i2, j1:j2, lm(region)-1) )
    allocate( conv_adv(i1:i2, j1:j2, lm(region)  ) )

    
    ! leave if fluxes are not in use:
    if ( .not. pu_dat(region)%used ) return
    if ( .not. pv_dat(region)%used ) return
    
    ! set pointers:
    pu => pu_dat(region)%data
    pv => pv_dat(region)%data
    am => wind_dat(region)%am
    bm => wind_dat(region)%bm
    cm => wind_dat(region)%cm

    !---- length of advection substeps (in seconds) in the three directions
    dtu=ndyn/(2.*tref(region))  ! again removed the revert! cmk
    dtv=ndyn/(2.*tref(region))
    dtw=ndyn/(2.*tref(region))

    ! number of levels in this region:
    lmr = lm(region)
    
    ! init to zero:
    am = zero
    bm = zero
    cm = zero

    ! update pu and pv before using them
    CALL UPDATE_HALO( dgrid(region), pu_dat(region)%data, pu_dat(region)%halo, status)
    IF_NOTOK_RETURN(status=1)

    CALL UPDATE_HALO( dgrid(region), pv_dat(region)%data, pv_dat(region)%halo, status)
    IF_NOTOK_RETURN(status=1)

    
    !
    !     compute conv_adv, the horizontal mass convergence
    !      the arrays pu and pv contain the horizontal air mass fluxes crossing 
    !      the grid box boundaries. pu(i,j,l) is the eastward mass flux in kg/sec
    !      at the eastern edge of box i,j,l; pv(i,j,l) is the northward
    !      mass flux at the southern edge of box i,j,l
    !      
    do l=1,lmr
    do j=j1, j2
    do i=i1, i2
             conv_adv(i,j,l)=pu(i-1,j,l)-pu(i,j,l)+pv(i,j,l)-pv(i,j+1,l)
    end do
    end do
    end do
    !
    !     compute pit, the vertically integrated conv_adv
    !
    do j=j1, j2
    do i=i1, i2
          
          pit(i,j)=conv_adv(i,j,1)

          do l=2,lmr
             pit(i,j)=pit(i,j)+conv_adv(i,j,l)
          end do
          !
          !     compute the vertical massflux on the box boundaries (in kg/s)
          !mk  in the hybrid system, the tendency in ps is given by bt. Bt is
          !mk  already normalized. Now distribute pit such that the new mass 
          !mk  distribution after advection exactly matches the 
          !mk  exactly coincides with the distribution that is obtained 
          !mk  from the new surface pressure.
          !
          sd(i,j,lmr-1)=conv_adv(i,j,lmr) - (bt(lmr)-bt(lmr+1))*pit(i,j)
          do  l=lmr-2,1,-1
             sd(i,j,l)=sd(i,j,l+1)+conv_adv(i,j,l+1)-(bt(l+1)-bt(l+2))*pit(i,j)
          end do
          
    end do
    end do
    !
    !     compute amount of air moved each advection substep, am, bm or cm (kg)
    !
    !      am(i,j,l) is the amount of air moved eastward at the eastern
    !      boundary of grid box i,j,l
    !     compute cm
    !      cm(i,j,l) is the amount of air moved in upward direction at the
    !      upper boundary of grid box i,j,l
    !     compute bm
    !      bm(i,j,l) is the amount of air moved northward at the southern
    !      boundary of grid box i,j,l
    !
    do l= 1,    lmr
    do j= j1,   j2
    do i= i1-1, i2
             am(i,j,l)=dtu*pu(i,j,l)
    end do
    end do
    end do

    do l= 1,  lmr
    do j= j1, j2+1
    do i= i1, i2
             bm(i,j,l)=dtv*pv(i,j,l)
    end do
    end do
    end do

    do l= 1,  lmr-1
    do j= j1, j2
    do i= i1, i2
             cm(i,j,l)=-dtw*sd(i,j,l)
    end do
    end do
    end do 

    CALL UPDATE_HALO( dgrid(region), am, wind_dat(region)%halo, status)
    IF_NOTOK_RETURN(status=1)
    
    nullify(pu)
    nullify(pv)
    nullify(am)
    nullify(bm)
    nullify(cm)

    deallocate( pit, sd, conv_adv)
    
  END SUBROUTINE DYNAM0


end module advect