dp_coupling.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
module dp_coupling
!BOP
!
! !MODULE: dp_coupling --- dynamics-physics coupling module
!
   use precision,     only: r8
   use rgrid,         only: nlon
   use pmgrid,        only: plon, plat, plev, twod_decomp, iam,    &
                            beglat, endlat, beglev, endlev,        &
                            beglatxy, endlatxy, beglonxy, endlonxy
   use ppgrid,        only: pcols, pver
   use phys_grid
   use physics_types, only: physics_state, physics_tend
   use constituents,  only: ppcnst
!
! !PUBLIC MEMBER FUNCTIONS:
      PUBLIC d_p_coupling, p_d_coupling

!
! !DESCRIPTION:
!
!      This module provides 
!
!      \begin{tabular}{|l|l|} \hline \hline
!        d\_p\_coupling    &  dynamics output to physics input \\ \hline
!        p\_d\_coupling    &  physics output to dynamics input \\ \hline 
!                                \hline
!      \end{tabular}
!
! !REVISION HISTORY:
!   00.06.01   Boville    Creation
!   01.10.01   Lin        Various revisions
!   01.03.26   Sawyer     Added ProTeX documentation
!   01.06.27   Mirin      Separate noncoupling coding into new routines
!   01.07.13   Mirin      Some support for multi-2D decompositions
!   02.03.01   Worley     Support for nontrivial physics remapping
!
!EOP
!-----------------------------------------------------------------------

CONTAINS

!-----------------------------------------------------------------------
!BOP
! !IROUTINE: d_p_coupling --- convert dynamics output to physics input
!
! !INTERFACE: 
  subroutine d_p_coupling(ps,   u3s,   v3s,  pt,    coslon,  sinlon,     &
                          t3,   q3, omga, phis,  pe,      peln,    pk,&
                          pkz,  phys_state,  phys_tend, full_phys,       &
                          qtmp, psxy,  u3sxy,v3sxy, ptxy,    t3xy,       &
                          q3xy, omgaxy,phisxy,      pexy,  pelnxy,       &
                          pkxy, pkzxy, qtmpxy,      pe11k,   pe11kln  )

! !USES:
    use physconst, only: zvir
    use dynamics_vars, only: ng_d, ng_s
#if defined (SPMD)
    use spmd_dyn, only : inter_ikj
    use mpishorthand, only : mpicom
    use parutilitiesmodule, only : sumop, parcollective
    use redistributemodule, only : redistributestart, redistributefinish
#endif
!-----------------------------------------------------------------------
    implicit none
!-----------------------------------------------------------------------
! !INPUT PARAMETERS:
!
    real(r8), intent(in) :: ps (plon, beglat:endlat)                      ! surface pressure
    real(r8), intent(inout) :: u3s(plon, beglat-ng_d:endlat+ng_s, beglev:endlev)       ! u-wind on d-grid
    real(r8), intent(in) :: v3s(plon, beglat-ng_s:endlat+ng_d, beglev:endlev)       ! v-wind on d-grid
    real(r8), intent(in) :: pt (plon, beglat-ng_d:endlat+ng_d, beglev:endlev)       ! Virtual pot temp
    real(r8), intent(in) :: t3 (plon, beglev:endlev, beglat:endlat)       ! virtual temperature
    real(r8), intent(in) :: q3 (plon, beglat-ng_d:endlat+ng_d, beglev:endlev, ppcnst) ! constituents
    real(r8), intent(in) :: omga(plon, beglev:endlev, beglat:endlat)      ! vertical velocity
    real(r8), intent(in) :: phis(plon, beglat:endlat)                     ! surface geopotential
    real(r8), intent(in) :: pe  (plon, beglev:endlev+1, beglat:endlat)    ! this fv's pint
    real(r8), intent(in) :: peln(plon, beglev:endlev+1, beglat:endlat)    ! log(pe)
    real(r8), intent(in) :: pk  (plon, beglat:endlat, beglev:endlev+1)    ! pe**cappa
    real(r8), intent(in) :: pkz (plon, beglat:endlat, beglev:endlev)      ! f-v mean of pk
    real(r8), intent(in) :: coslon(plon)                                  ! cosine of longitude
    real(r8), intent(in) :: sinlon(plon)                                  ! sin of longitudes
    logical,  intent(in) :: full_phys

! xy-decomposed instanciations below:
    real(r8), intent(in) :: psxy (beglonxy:endlonxy, beglatxy:endlatxy)                      ! surface pressure
    real(r8), intent(in) :: u3sxy(beglonxy:endlonxy, beglatxy:endlatxy+1, plev)       ! u-wind on d-grid
    real(r8), intent(in) :: v3sxy(beglonxy:endlonxy, beglatxy:endlatxy, plev)       ! v-wind on d-grid
    real(r8), intent(in) :: ptxy (beglonxy:endlonxy, beglatxy:endlatxy, plev)       ! Virtual pot temp
    real(r8), intent(in) :: t3xy (beglonxy:endlonxy, plev, beglatxy:endlatxy)       ! virtual temperature
    real(r8), intent(in) :: q3xy (beglonxy:endlonxy, beglatxy:endlatxy, plev, ppcnst) ! constituents
    real(r8), intent(in) :: omgaxy(beglonxy:endlonxy, plev, beglatxy:endlatxy)      ! vertical velocity
    real(r8), intent(in) :: phisxy(beglonxy:endlonxy, beglatxy:endlatxy)            ! surface geopotential
    real(r8), intent(in) :: pexy  (beglonxy:endlonxy, plev+1, beglatxy:endlatxy)    ! this fv's pint
    real(r8), intent(in) :: pelnxy(beglonxy:endlonxy, plev+1, beglatxy:endlatxy)    ! log(pe)
    real(r8), intent(in) :: pkxy  (beglonxy:endlonxy, beglatxy:endlatxy, plev+1)    ! pe**cappa
    real(r8), intent(in) :: pkzxy (beglonxy:endlonxy, beglatxy:endlatxy, plev)      ! f-v mean of pk

! !OUTPUT PARAMETERS:

    type(physics_state), intent(out), dimension(begchunk:endchunk) :: phys_state
    type(physics_tend ), intent(out), dimension(begchunk:endchunk) :: phys_tend
    real(r8), intent(out) :: qtmp(plon, beglev:endlev, beglat:endlat)      ! temporary moisture storage
    real(r8), intent(out) :: qtmpxy(beglonxy:endlonxy, plev, beglatxy:endlatxy)      ! temporary moisture storage
    real(r8), intent(out) :: pe11k(plev+1), pe11kln(plev+1)  ! Pres. & log for Rayl. fric

! !DESCRIPTION:
!
!   Coupler for converting dynamics output variables into physics 
!   input variables
!
! !REVISION HISTORY:
!   00.06.01   Boville    Creation
!   01.07.13   AAM        Some support for multi-2D decompositions
!   02.03.01   Worley     Support for nontrivial physics remapping
!
!EOP
!-----------------------------------------------------------------------
!BOC
! !LOCAL VARIABLES:
    integer :: i,ib,j,k,m,lchnk      ! indices
    integer :: ncol                  ! number of columns in current chunk
    integer :: lats(pcols)           ! array of latitude indices
    integer :: lons(pcols)           ! array of longitude indices
    integer :: blksiz                ! number of columns in 2D block
    integer :: tsize                 ! amount of data per grid point passed to physics
    integer, allocatable, dimension(:,:) :: bpter
                                     ! offsets into block buffer for packing data
    integer :: cpter(pcols,0:pver)   ! offsets into chunk buffer for unpacking data

    real(r8) :: pic(pcols)           ! ps**cappa
    real(r8), allocatable :: u3(:, :, :)       ! u-wind on a-grid
    real(r8), allocatable :: v3(:, :, :)       ! v-wind on a-grid
    real(r8), allocatable, dimension(:) :: bbuffer, cbuffer
                                     ! transpose buffers
!---------------------------End Local workspace-------------------------

    if (twod_decomp .eq. 1) then

!-----------------------------------------------------------------------
! Store moisture in temporary array, to be used after physics update
!-----------------------------------------------------------------------
       if (full_phys) then
!$omp parallel do private(i, j, k)
          do j=beglatxy,endlatxy
             do k=1,plev
                do i=beglonxy,endlonxy
                   qtmpxy(i,k,j) = q3xy(i,j,k,1)
                enddo
             enddo
          enddo
!
! Transpose temporary moisture array back to yz decomposition
!
#if defined (SPMD)
          call redistributestart (inter_ikj, .false., qtmpxy)
          call redistributefinish(inter_ikj, .false., qtmp)
#endif

       endif

!-----------------------------------------------------------------------
! Transform dynamics staggered winds to physics grid (D=>A)
!-----------------------------------------------------------------------

       allocate (u3(beglonxy:endlonxy, plev, beglatxy:endlatxy))
       allocate (v3(beglonxy:endlonxy, plev, beglatxy:endlatxy))

       call d2a3dikj(u3sxy,  v3sxy,   u3,    v3,   plon,  plat, plev,    &
                     beglatxy, endlatxy,  0, 0, 1, 0, 0,                 &
                     beglonxy, endlonxy,    coslon,      sinlon)

!-----------------------------------------------------------------------
! Copy data from dynamics data structure to physics data structure
!-----------------------------------------------------------------------

       if (local_dp_map) then

!$OMP PARALLEL DO PRIVATE (LCHNK, NCOL, I, K, M, LONS, LATS, PIC)

          do lchnk = begchunk,endchunk
             ncol = get_ncols_p(lchnk)
             call get_lon_all_p(lchnk, ncol, lons)
             call get_lat_all_p(lchnk, ncol, lats)
             phys_state(lchnk)%ncol  = ncol
             phys_state(lchnk)%lchnk = lchnk

             do i=1,ncol
                phys_state(lchnk)%ps(i)   = psxy(lons(i),lats(i))
                phys_state(lchnk)%phis(i) = phisxy(lons(i),lats(i))
                pic(i) = pkxy(lons(i),lats(i),pver+1)
             enddo

             do k=1,plev
                do i=1,ncol
                   phys_state(lchnk)%u    (i,k) = u3(lons(i),k,lats(i))
                   phys_state(lchnk)%v    (i,k) = v3(lons(i),k,lats(i))
                   phys_state(lchnk)%omega(i,k) = omgaxy(lons(i),k,lats(i))
                   if (full_phys) then
                      phys_state(lchnk)%t    (i,k) = t3xy(lons(i),k,lats(i)) / (1. + zvir*q3xy(lons(i),lats(i),k,1))
                      phys_state(lchnk)%exner(i,k) = pic(i) / pkzxy(lons(i),lats(i),k) 
                   else
                      phys_state(lchnk)%t    (i,k) = ptxy(lons(i),lats(i),k) * pkzxy(lons(i),lats(i),k)
                   end if
                end do
             end do

             do k=1,plev+1
                do i=1,ncol
!
! edge-level pressure arrays: copy from the arrays computed by dynpkg
!
                   phys_state(lchnk)%pint  (i,k) = pexy  (lons(i),k,lats(i))
                   phys_state(lchnk)%lnpint(i,k) = pelnxy(lons(i),k,lats(i))
                end do
             end do

!
! Copy constituents
!
             do m=1,ppcnst
                do k=1,plev
                   do i=1,ncol
                      phys_state(lchnk)%q(i,k,m) = q3xy(lons(i),lats(i),k,m)
                   end do
                end do
             end do
 
          end do   ! begchunk:endchunk loop

       else

          tsize = 7 + ppcnst

          blksiz = (endlatxy-beglatxy+1)*(endlonxy-beglonxy+1)
          allocate(bpter(blksiz,0:plev))
          allocate(bbuffer(tsize*block_buf_nrecs))
          allocate(cbuffer(tsize*chunk_buf_nrecs))

          call block_to_chunk_send_pters(iam+1,blksiz,plev+1,tsize,bpter)

          ib = 0
          do j=beglatxy,endlatxy
             do i=beglonxy,endlonxy
                ib = ib + 1

                bbuffer(bpter(ib,0))   = pexy(i,plev+1,j)
                bbuffer(bpter(ib,0)+1) = pelnxy(i,plev+1,j)
                bbuffer(bpter(ib,0)+2) = psxy(i,j)
                bbuffer(bpter(ib,0)+3) = phisxy(i,j)

                do k=1,plev

                   bbuffer(bpter(ib,k))   = pexy(i,k,j)
                   bbuffer(bpter(ib,k)+1) = pelnxy(i,k,j)
                   bbuffer(bpter(ib,k)+2) = u3    (i,k,j)
                   bbuffer(bpter(ib,k)+3) = v3    (i,k,j)
                   bbuffer(bpter(ib,k)+4) = omgaxy(i,k,j)
                   if (full_phys) then
                      bbuffer(bpter(ib,k)+5) = t3xy(i,k,j) / (1. + zvir*q3xy(i,j,k,1))
                      bbuffer(bpter(ib,k)+6) = pkxy(i,j,pver+1) / pkzxy(i,j,k) 
                   else
                      bbuffer(bpter(ib,k)+6) = ptxy(i,j,k) * pkzxy(i,j,k)
                   end if

                   do m=1,ppcnst
                      bbuffer(bpter(ib,k)+6+m) = q3xy(i,j,k,m)
                   end do

                end do
             end do
          end do

          call transpose_block_to_chunk(tsize, bbuffer, cbuffer)

          do lchnk = begchunk,endchunk
             ncol = get_ncols_p(lchnk)
             phys_state(lchnk)%ncol  = ncol
             phys_state(lchnk)%lchnk = lchnk

             call block_to_chunk_recv_pters(lchnk,pcols,pver+1,tsize,cpter)

             do i=1,ncol

                phys_state(lchnk)%pint  (i,pver+1) = cbuffer(cpter(i,0))
                phys_state(lchnk)%lnpint(i,pver+1) = cbuffer(cpter(i,0)+1)
                phys_state(lchnk)%ps(i)            = cbuffer(cpter(i,0)+2)
                phys_state(lchnk)%phis(i)          = cbuffer(cpter(i,0)+3)

                do k=1,plev

                   phys_state(lchnk)%pint  (i,k) = cbuffer(cpter(i,k))
                   phys_state(lchnk)%lnpint(i,k) = cbuffer(cpter(i,k)+1)
                   phys_state(lchnk)%u     (i,k) = cbuffer(cpter(i,k)+2)
                   phys_state(lchnk)%v     (i,k) = cbuffer(cpter(i,k)+3)
                   phys_state(lchnk)%omega (i,k) = cbuffer(cpter(i,k)+4)
                   if (full_phys) then
                      phys_state(lchnk)%t    (i,k) = cbuffer(cpter(i,k)+5)
                      phys_state(lchnk)%exner(i,k) = cbuffer(cpter(i,k)+6)
                   else
                      phys_state(lchnk)%t    (i,k) = cbuffer(cpter(i,k)+6)
                   end if

                   do m=1,ppcnst
                      phys_state(lchnk)%q(i,k,m) = cbuffer(cpter(i,k)+6+m)
                   end do

                end do
             end do

          end do   ! begchunk:endchunk loop

          deallocate(bpter)
          deallocate(bbuffer)
          deallocate(cbuffer)

       endif

    else

!-----------------------------------------------------------------------
! Store moisture in temporary array, to be used after physics update
!-----------------------------------------------------------------------