LCOV - code coverage report
Current view: | top level - dynamics/se - dp_coupling.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 253 | 326 | 77.6 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 3 | 3 | 100.0 %

              Line data    Source code
           1             : module dp_coupling
           2             :
           3             : !-------------------------------------------------------------------------------
           4             : ! dynamics - physics coupling module
           5             : !-------------------------------------------------------------------------------
           6             :
           7             : use shr_kind_mod,   only: r8=>shr_kind_r8
           8             : use ppgrid,         only: begchunk, endchunk, pcols, pver, pverp
           9             : use constituents,   only: pcnst, cnst_type
          10             :
          11             : use spmd_dyn,       only: local_dp_map
          12             : use spmd_utils,     only: iam
          13             : use dyn_grid,       only: TimeLevel, edgebuf
          14             : use dyn_comp,       only: dyn_export_t, dyn_import_t
          15             :
          16             : use physics_types,  only: physics_state, physics_tend, physics_cnst_limit
          17             : use phys_grid,      only: get_ncols_p
          18             : use phys_grid,      only: get_dyn_col_p, columns_on_task, get_chunk_info_p, phys_columns_on_task
          19             : use physics_buffer, only: physics_buffer_desc, pbuf_get_chunk, pbuf_get_field
          20             :
          21             : use dp_mapping,     only: nphys_pts
          22             :
          23             : use perf_mod,       only: t_startf, t_stopf
          24             : use cam_abortutils, only: endrun
          25             :
          26             : use parallel_mod,   only: par
          27             : use thread_mod,     only: horz_num_threads, max_num_threads
          28             : use hybrid_mod,     only: config_thread_region, get_loop_ranges, hybrid_t
          29             : use dimensions_mod, only: np, nelemd, nlev, qsize, ntrac, fv_nphys
          30             :
          31             : use dof_mod,        only: UniquePoints, PutUniquePoints
          32             : use element_mod,    only: element_t
          33             :
          34             : implicit none
          35             : private
          36             : save
          37             :
          38             : public :: d_p_coupling, p_d_coupling
          39             :
          40             : real (kind=r8), allocatable :: q_prev(:,:,:,:) ! Previous Q for computing tendencies
          41             :
          42             : !=========================================================================================
          43             : CONTAINS
          44             : !=========================================================================================
          45             :
          46     1112832 : subroutine d_p_coupling(phys_state, phys_tend,  pbuf2d, dyn_out)
          47             :
          48             :    ! Convert the dynamics output state into the physics input state.
          49             :    ! Note that all pressures and tracer mixing ratios coming from the dycore are based on
          50             :    ! dry air mass.
          51             :
          52             :    use gravity_waves_sources,  only: gws_src_fnct
          53             :    use dyn_comp,               only: frontgf_idx, frontga_idx
          54             :    use phys_control,           only: use_gw_front, use_gw_front_igw
          55             :    use hycoef,                 only: hyai, ps0
          56             :    use fvm_mapping,            only: dyn2phys_vector, dyn2phys_all_vars
          57             :    use se_dyn_time_mod,        only: timelevel_qdp
          58             :    use control_mod,            only: qsplit
          59             :    use test_fvm_mapping,       only: test_mapping_overwrite_dyn_state, test_mapping_output_phys_state
          60             :    use prim_advance_mod,       only: tot_energy_dyn
          61             :    ! arguments
          62             :    type(dyn_export_t),  intent(inout)                               :: dyn_out             ! dynamics export
          63             :    type(physics_buffer_desc), pointer                               :: pbuf2d(:,:)
          64             :    type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
          65             :    type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend
          66             :
          67             :
          68             :    ! LOCAL VARIABLES
          69      370944 :    type(element_t), pointer     :: elem(:)             ! pointer to dyn_out element array
          70             :    integer                      :: ie                  ! indices over elements
          71             :    integer                      :: lchnk, icol, ilyr   ! indices over chunks, columns, layers
          72             :
          73      370944 :    real (kind=r8),  allocatable :: ps_tmp(:,:)         ! temp array to hold ps
          74      370944 :    real (kind=r8),  allocatable :: dp3d_tmp(:,:,:)     ! temp array to hold dp3d
          75      370944 :    real (kind=r8),  allocatable :: dp3d_tmp_tmp(:,:)
          76      370944 :    real (kind=r8),  allocatable :: phis_tmp(:,:)       ! temp array to hold phis
          77      370944 :    real (kind=r8),  allocatable :: T_tmp(:,:,:)        ! temp array to hold T
          78      370944 :    real (kind=r8),  allocatable :: uv_tmp(:,:,:,:)     ! temp array to hold u and v
          79      370944 :    real (kind=r8),  allocatable :: q_tmp(:,:,:,:)      ! temp to hold advected constituents
          80      370944 :    real (kind=r8),  allocatable :: omega_tmp(:,:,:)    ! temp array to hold omega
          81             :
          82             :    ! Frontogenesis
          83      370944 :    real (kind=r8),  allocatable :: frontgf(:,:,:)      ! temp arrays to hold frontogenesis
          84      370944 :    real (kind=r8),  allocatable :: frontga(:,:,:)      ! function (frontgf) and angle (frontga)
          85      370944 :    real (kind=r8),  allocatable :: frontgf_phys(:,:,:)
          86      370944 :    real (kind=r8),  allocatable :: frontga_phys(:,:,:)
          87             :                                                         ! Pointers to pbuf
          88      370944 :    real (kind=r8),  pointer     :: pbuf_frontgf(:,:)
          89      370944 :    real (kind=r8),  pointer     :: pbuf_frontga(:,:)
          90             :
          91             :    integer                      :: ncols, ierr
          92             :    integer                      :: col_ind, blk_ind(1), m
          93             :    integer                      :: nphys
          94             :
          95      370944 :    real (kind=r8),  allocatable :: qgll(:,:,:,:)
          96             :    real (kind=r8)               :: inv_dp3d(np,np,nlev)
          97             :    integer                      :: tl_f, tl_qdp_np0, tl_qdp_np1
          98             :
          99      370944 :    type(physics_buffer_desc), pointer :: pbuf_chnk(:)
         100             :    !----------------------------------------------------------------------------
         101             :
         102      370944 :    if (.not. local_dp_map) then
         103           0 :       call endrun('d_p_coupling: Weak scaling does not support load balancing')
         104             :    end if
         105             :
         106      370944 :    elem => dyn_out%elem
         107      370944 :    tl_f = TimeLevel%n0
         108      370944 :    call TimeLevel_Qdp(TimeLevel, qsplit, tl_qdp_np0,tl_qdp_np1)
         109             :
         110      370944 :    nullify(pbuf_chnk)
         111      370944 :    nullify(pbuf_frontgf)
         112      370944 :    nullify(pbuf_frontga)
         113             :
         114      370944 :    if (fv_nphys > 0) then
         115      370944 :       nphys = fv_nphys
         116             :    else
         117           0 :      allocate(qgll(np,np,nlev,pcnst))
         118           0 :      nphys = np
         119             :    end if
         120             :
         121             :    ! Allocate temporary arrays to hold data for physics decomposition
         122     1483776 :    allocate(ps_tmp(nphys_pts,nelemd))
         123     1483776 :    allocate(dp3d_tmp(nphys_pts,pver,nelemd))
         124     1112832 :    allocate(dp3d_tmp_tmp(nphys_pts,pver))
         125     1112832 :    allocate(phis_tmp(nphys_pts,nelemd))
         126     1112832 :    allocate(T_tmp(nphys_pts,pver,nelemd))
         127     1854720 :    allocate(uv_tmp(nphys_pts,2,pver,nelemd))
         128     1854720 :    allocate(q_tmp(nphys_pts,pver,pcnst,nelemd))
         129     1112832 :    allocate(omega_tmp(nphys_pts,pver,nelemd))
         130             :
         131      370944 :    call tot_energy_dyn(elem,dyn_out%fvm, 1, nelemd,tl_f , tl_qdp_np0,'dBF')
         132             :
         133      370944 :    if (use_gw_front .or. use_gw_front_igw) then
         134           0 :       allocate(frontgf(nphys_pts,pver,nelemd), stat=ierr)
         135           0 :       if (ierr /= 0) call endrun("dp_coupling: Allocate of frontgf failed.")
         136           0 :       allocate(frontga(nphys_pts,pver,nelemd), stat=ierr)
         137           0 :       if (ierr /= 0) call endrun("dp_coupling: Allocate of frontga failed.")
         138             :    end if
         139             :
         140      370944 :    if (iam < par%nprocs) then
         141      370944 :       if (use_gw_front .or. use_gw_front_igw) then
         142           0 :          call gws_src_fnct(elem, tl_f, tl_qdp_np0, frontgf, frontga, nphys)
         143             :       end if
         144             :
         145      370944 :       if (fv_nphys > 0) then
         146      370944 :          call test_mapping_overwrite_dyn_state(elem,dyn_out%fvm)
         147             :          !******************************************************************
         148             :          ! physics runs on an FVM grid: map GLL vars to physics grid
         149             :          !******************************************************************
         150      370944 :          call t_startf('dyn2phys')
         151             :          ! note that the fvm halo has been filled in prim_run_subcycle
         152             :          ! if physics grid resolution is not equal to fvm resolution
         153             :          call dyn2phys_all_vars(1,nelemd,elem, dyn_out%fvm,&
         154             :               pcnst,hyai(1)*ps0,tl_f,                      &
         155             :               ! output
         156             :               dp3d_tmp, ps_tmp, q_tmp, T_tmp,              &
         157             :               omega_tmp, phis_tmp                          &
         158      370944 :               )
         159     2979144 :          do ie = 1, nelemd
         160             :             uv_tmp(:,:,:,ie) = &
         161     2979144 :                dyn2phys_vector(elem(ie)%state%v(:,:,:,:,tl_f),elem(ie))
         162             :          end do
         163      370944 :          call t_stopf('dyn2phys')
         164             :       else
         165             :
         166             :          !******************************************************************
         167             :          ! Physics runs on GLL grid: collect unique points before mapping to
         168             :          ! physics decomposition
         169             :          !******************************************************************
         170             :
         171           0 :          if (qsize < pcnst) then
         172           0 :             call endrun('d_p_coupling: Fewer GLL tracers advected than required')
         173             :          end if
         174             :
         175           0 :          call t_startf('UniquePoints')
         176           0 :          do ie = 1, nelemd
         177           0 :            inv_dp3d(:,:,:) = 1.0_r8/elem(ie)%state%dp3d(:,:,:,tl_f)
         178           0 :            do m=1,pcnst
         179           0 :              qgll(:,:,:,m) = elem(ie)%state%Qdp(:,:,:,m,tl_qdp_np0)*inv_dp3d(:,:,:)
         180             :            end do
         181           0 :             ncols = elem(ie)%idxP%NumUniquePts
         182           0 :             call UniquePoints(elem(ie)%idxP, elem(ie)%state%psdry(:,:), ps_tmp(1:ncols,ie))
         183           0 :             call UniquePoints(elem(ie)%idxP, nlev, elem(ie)%state%dp3d(:,:,:,tl_f), dp3d_tmp(1:ncols,:,ie))
         184           0 :             call UniquePoints(elem(ie)%idxP, nlev, elem(ie)%state%T(:,:,:,tl_f), T_tmp(1:ncols,:,ie))
         185           0 :             call UniquePoints(elem(ie)%idxV, 2, nlev, elem(ie)%state%V(:,:,:,:,tl_f), uv_tmp(1:ncols,:,:,ie))
         186           0 :             call UniquePoints(elem(ie)%idxV, nlev, elem(ie)%derived%omega, omega_tmp(1:ncols,:,ie))
         187             :
         188           0 :             call UniquePoints(elem(ie)%idxP, elem(ie)%state%phis, phis_tmp(1:ncols,ie))
         189           0 :             call UniquePoints(elem(ie)%idxP, nlev, pcnst, qgll,Q_tmp(1:ncols,:,:,ie))
         190             :          end do
         191           0 :          call t_stopf('UniquePoints')
         192             :
         193             :       end if ! if fv_nphys>0
         194             :
         195             :    else
         196             :
         197           0 :       ps_tmp(:,:)      = 0._r8
         198           0 :       T_tmp(:,:,:)     = 0._r8
         199           0 :       uv_tmp(:,:,:,:)  = 0._r8
         200           0 :       omega_tmp(:,:,:) = 0._r8
         201           0 :       phis_tmp(:,:)    = 0._r8
         202           0 :       Q_tmp(:,:,:,:)   = 0._r8
         203             :
         204           0 :       if (use_gw_front .or. use_gw_front_igw) then
         205           0 :          frontgf(:,:,:) = 0._r8
         206           0 :          frontga(:,:,:) = 0._r8
         207             :       end if
         208             :
         209             :    endif ! iam < par%nprocs
         210             :
         211      370944 :    if (fv_nphys < 1) then
         212           0 :       deallocate(qgll)
         213             :    end if
         214             :
         215             :    ! q_prev is for saving the tracer fields for calculating tendencies
         216      370944 :    if (.not. allocated(q_prev)) then
         217        4608 :       allocate(q_prev(pcols,pver,pcnst,begchunk:endchunk))
         218             :    end if
         219  1989210384 :    q_prev = 0.0_R8
         220             :
         221      370944 :    call t_startf('dpcopy')
         222      370944 :    if (use_gw_front .or. use_gw_front_igw) then
         223           0 :       allocate(frontgf_phys(pcols, pver, begchunk:endchunk))
         224           0 :       allocate(frontga_phys(pcols, pver, begchunk:endchunk))
         225             :    end if
         226             :    !$omp parallel do num_threads(max_num_threads) private (col_ind, lchnk, icol, ie, blk_ind, ilyr, m)
         227    23844744 :    do col_ind = 1, phys_columns_on_task
         228    23473800 :       call get_dyn_col_p(col_ind, ie, blk_ind)
         229    23473800 :       call get_chunk_info_p(col_ind, lchnk, icol)
         230    23473800 :       phys_state(lchnk)%ps(icol)   = ps_tmp(blk_ind(1), ie)
         231    23473800 :       phys_state(lchnk)%phis(icol) = phis_tmp(blk_ind(1), ie)
         232   633792600 :       do ilyr = 1, pver
         233   610318800 :          phys_state(lchnk)%pdel(icol, ilyr)  = dp3d_tmp(blk_ind(1), ilyr, ie)
         234   610318800 :          phys_state(lchnk)%t(icol, ilyr)     = T_tmp(blk_ind(1), ilyr, ie)
         235   610318800 :          phys_state(lchnk)%u(icol, ilyr)     = uv_tmp(blk_ind(1), 1, ilyr, ie)
         236   610318800 :          phys_state(lchnk)%v(icol, ilyr)     = uv_tmp(blk_ind(1), 2, ilyr, ie)
         237   610318800 :          phys_state(lchnk)%omega(icol, ilyr) = omega_tmp(blk_ind(1), ilyr, ie)
         238             :
         239   633792600 :          if (use_gw_front .or. use_gw_front_igw) then
         240           0 :             frontgf_phys(icol, ilyr, lchnk) = frontgf(blk_ind(1), ilyr, ie)
         241           0 :             frontga_phys(icol, ilyr, lchnk) = frontga(blk_ind(1), ilyr, ie)
         242             :          end if
         243             :       end do
         244             :
         245   117739944 :       do m = 1, pcnst
         246  1924851600 :          do ilyr = 1, pver
         247  1901377800 :             phys_state(lchnk)%q(icol, ilyr,m) = Q_tmp(blk_ind(1), ilyr,m, ie)
         248             :          end do
         249             :       end do
         250             :    end do
         251      370944 :    if (use_gw_front .or. use_gw_front_igw) then
         252             :       !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, icol, ilyr, pbuf_chnk, pbuf_frontgf, pbuf_frontga)
         253           0 :       do lchnk = begchunk, endchunk
         254           0 :          ncols = get_ncols_p(lchnk)
         255           0 :          pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
         256           0 :          call pbuf_get_field(pbuf_chnk, frontgf_idx, pbuf_frontgf)
         257           0 :          call pbuf_get_field(pbuf_chnk, frontga_idx, pbuf_frontga)
         258           0 :          do icol = 1, ncols
         259           0 :             do ilyr = 1, pver
         260           0 :                pbuf_frontgf(icol, ilyr) = frontgf_phys(icol, ilyr, lchnk)
         261           0 :                pbuf_frontga(icol, ilyr) = frontga_phys(icol, ilyr, lchnk)
         262             :             end do
         263             :          end do
         264             :       end do
         265           0 :       deallocate(frontgf_phys)
         266           0 :       deallocate(frontga_phys)
         267             :    end if
         268             :
         269      370944 :    call t_stopf('dpcopy')
         270             :
         271             :    ! Save the tracer fields input to physics package for calculating tendencies
         272             :    ! The mixing ratios are all dry at this point.
         273     1866312 :    do lchnk = begchunk, endchunk
         274     1495368 :       ncols = phys_state(lchnk)%ncol
         275  1953947520 :       q_prev(1:ncols,1:pver,1:pcnst,lchnk) = phys_state(lchnk)%q(1:ncols,1:pver,1:pcnst)
         276             :    end do
         277      370944 :    call test_mapping_output_phys_state(phys_state,dyn_out%fvm)
         278             :
         279             :    ! Deallocate the temporary arrays
         280      370944 :    deallocate(ps_tmp)
         281      370944 :    deallocate(dp3d_tmp)
         282      370944 :    deallocate(phis_tmp)
         283      370944 :    deallocate(T_tmp)
         284      370944 :    deallocate(uv_tmp)
         285      370944 :    deallocate(q_tmp)
         286      370944 :    deallocate(omega_tmp)
         287             :
         288             :    ! ps, pdel, and q in phys_state are all dry at this point.  After return from derived_phys_dry
         289             :    ! ps and pdel include water vapor only, and the 'wet' constituents have been converted to wet mmr.
         290      370944 :    call t_startf('derived_phys')
         291      370944 :    call derived_phys_dry(phys_state, phys_tend, pbuf2d)
         292      370944 :    call t_stopf('derived_phys')
         293             :
         294             :   !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, ilyr, icol)
         295     1866312 :    do lchnk = begchunk, endchunk
         296     1495368 :       ncols=get_ncols_p(lchnk)
         297     1866312 :       if (pcols > ncols) then
         298     1193976 :          phys_state(lchnk)%phis(ncols+1:) = 0.0_r8
         299             :       end if
         300             :    end do
         301      741888 : end subroutine d_p_coupling
         302             :
         303             : !=========================================================================================
         304             :
         305      738816 : subroutine p_d_coupling(phys_state, phys_tend, dyn_in, tl_f, tl_qdp)
         306             :
         307             :    ! Convert the physics output state into the dynamics input state.
         308             :
         309      370944 :    use phys_grid,        only: get_dyn_col_p, columns_on_task, get_chunk_info_p, phys_columns_on_task
         310             :    use bndry_mod,        only: bndry_exchange
         311             :    use edge_mod,         only: edgeVpack, edgeVunpack
         312             :    use fvm_mapping,      only: phys2dyn_forcings_fvm
         313             :    use test_fvm_mapping, only: test_mapping_overwrite_tendencies
         314             :    use test_fvm_mapping, only: test_mapping_output_mapped_tendencies
         315             :    use dimensions_mod,   only: use_cslam
         316             :    ! arguments
         317             :    type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
         318             :    type(physics_tend),  intent(inout), dimension(begchunk:endchunk) :: phys_tend
         319             :    integer,             intent(in)                                  :: tl_qdp, tl_f
         320             :    type(dyn_import_t),  intent(inout)                               :: dyn_in
         321             :    type(hybrid_t)                                                   :: hybrid
         322             :
         323             :    ! LOCAL VARIABLES
         324             :    integer                      :: ncols             ! index
         325      369408 :    type(element_t), pointer     :: elem(:)           ! pointer to dyn_in element array
         326             :    integer                      :: ie                ! index for elements
         327             :    integer                      :: col_ind           ! index over columns
         328             :    integer                      :: blk_ind(1)        ! element offset
         329             :    integer                      :: lchnk, icol, ilyr ! indices for chunk, column, layer
         330             :
         331      369408 :    real (kind=r8),  allocatable :: dp_phys(:,:,:)    ! temp array to hold dp on physics grid
         332      369408 :    real (kind=r8),  allocatable :: T_tmp(:,:,:)      ! temp array to hold T
         333      369408 :    real (kind=r8),  allocatable :: dq_tmp(:,:,:,:)   ! temp array to hold q
         334      369408 :    real (kind=r8),  allocatable :: uv_tmp(:,:,:,:)   ! temp array to hold uv
         335             :    integer                      :: m, i, j, k
         336             :
         337             :    real (kind=r8)               :: factor
         338             :    integer                      :: num_trac
         339             :    integer                      :: nets, nete
         340             :    integer                      :: kptr, ii
         341             :    !----------------------------------------------------------------------------
         342             :
         343      369408 :    if (.not. local_dp_map) then
         344           0 :       call endrun('p_d_coupling: Weak scaling does not support load balancing')
         345             :    end if
         346             :
         347      369408 :    if (iam < par%nprocs) then
         348      369408 :       elem => dyn_in%elem
         349             :    else
         350           0 :       nullify(elem)
         351             :    end if
         352             :
         353     1477632 :    allocate(T_tmp(nphys_pts,pver,nelemd))
         354     1847040 :    allocate(uv_tmp(nphys_pts,2,pver,nelemd))
         355     1847040 :    allocate(dq_tmp(nphys_pts,pver,pcnst,nelemd))
         356     1108224 :    allocate(dp_phys(nphys_pts,pver,nelemd))
         357             :
         358   678290808 :    T_tmp  = 0.0_r8
         359  1421147208 :    uv_tmp = 0.0_r8
         360  2036731008 :    dq_tmp = 0.0_r8
         361             :
         362      369408 :    if (.not. allocated(q_prev)) then
         363           0 :       call endrun('p_d_coupling: q_prev not allocated')
         364             :    end if
         365             :
         366             :    ! Convert wet to dry mixing ratios and modify the physics temperature
         367             :    ! tendency to be thermodynamically consistent with the dycore.
         368             :    !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, icol, ilyr, m, factor)
         369     1858584 :    do lchnk = begchunk, endchunk
         370     1489176 :       ncols = get_ncols_p(lchnk)
         371    25235184 :       do icol = 1, ncols
         372   632657376 :          do ilyr = 1, pver
         373             :             ! convert wet mixing ratios to dry
         374   607791600 :             factor = phys_state(lchnk)%pdel(icol,ilyr)/phys_state(lchnk)%pdeldry(icol,ilyr)
         375  2454543000 :             do m = 1, pcnst
         376  2431166400 :                if (cnst_type(m) == 'wet') then
         377  1823374800 :                   phys_state(lchnk)%q(icol,ilyr,m) = factor*phys_state(lchnk)%q(icol,ilyr,m)
         378             :                end if
         379             :             end do
         380             :          end do
         381             :       end do
         382             :     end do
         383             :
         384      369408 :    call t_startf('pd_copy')
         385             :    !$omp parallel do num_threads(max_num_threads) private (col_ind, lchnk, icol, ie, blk_ind, ilyr, m)
         386    23746008 :    do col_ind = 1, phys_columns_on_task
         387    23376600 :       call get_dyn_col_p(col_ind, ie, blk_ind)
         388    23376600 :       call get_chunk_info_p(col_ind, lchnk, icol)
         389             :
         390             :       ! test code -- does nothing unless cpp macro debug_coupling is defined.
         391    23376600 :       call test_mapping_overwrite_tendencies(phys_state(lchnk),            &
         392    23376600 :            phys_tend(lchnk), ncols, lchnk, q_prev(1:ncols,:,:,lchnk),      &
         393    70129800 :            dyn_in%fvm)
         394             :
         395   654914208 :       do ilyr = 1, pver
         396   607791600 :          dp_phys(blk_ind(1),ilyr,ie)  = phys_state(lchnk)%pdeldry(icol,ilyr)
         397   607791600 :          T_tmp(blk_ind(1),ilyr,ie)    = phys_tend(lchnk)%dtdt(icol,ilyr)
         398   607791600 :          uv_tmp(blk_ind(1),1,ilyr,ie) = phys_tend(lchnk)%dudt(icol,ilyr)
         399   607791600 :          uv_tmp(blk_ind(1),2,ilyr,ie) = phys_tend(lchnk)%dvdt(icol,ilyr)
         400  2454543000 :          do m = 1, pcnst
         401  1823374800 :             dq_tmp(blk_ind(1),ilyr,m,ie) =                                    &
         402  4254541200 :                  (phys_state(lchnk)%q(icol,ilyr,m) - q_prev(icol,ilyr,m,lchnk))
         403             :          end do
         404             :       end do
         405             :    end do
         406      369408 :    call t_stopf('pd_copy')
         407             :
         408      369408 :    if (iam < par%nprocs) then
         409             :
         410      369408 :       if (fv_nphys > 0) then
         411             :
         412             :          ! put forcings into fvm structure
         413      369408 :          num_trac = max(qsize,ntrac)
         414     2966808 :          do ie = 1, nelemd
         415    10759008 :             do j = 1, fv_nphys
         416    33766200 :                do i = 1, fv_nphys
         417    23376600 :                   ii = i + (j-1)*fv_nphys
         418   631168200 :                   dyn_in%fvm(ie)%ft(i,j,1:pver)                 = T_tmp(ii,1:pver,ie)
         419  1846751400 :                   dyn_in%fvm(ie)%fm(i,j,1:2,1:pver)             = uv_tmp(ii,1:2,1:pver,ie)
         420  1916881200 :                   dyn_in%fvm(ie)%fc_phys(i,j,1:pver,1:num_trac) = dq_tmp(ii,1:pver,1:num_trac,ie)
         421   638960400 :                   dyn_in%fvm(ie)%dp_phys(i,j,1:pver)            = dp_phys(ii,1:pver,ie)
         422             :                end do
         423             :             end do
         424             :          end do
         425             :
         426             :          !JMD $OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,n)
         427             :          !JMD        hybrid = config_thread_region(par,'horizontal')
         428      369408 :          hybrid = config_thread_region(par,'serial')
         429      369408 :          call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
         430             :          !
         431             :          ! high-order mapping of ft and fm using fvm technology
         432             :          !
         433      369408 :          call t_startf('phys2dyn')
         434      369408 :          call phys2dyn_forcings_fvm(elem, dyn_in%fvm, hybrid,nets,nete,ntrac==0, tl_f, tl_qdp)
         435      369408 :          call t_stopf('phys2dyn')
         436             :       else
         437             :
         438           0 :          call t_startf('putUniquePoints')
         439             :
         440             :          !$omp parallel do num_threads(max_num_threads) private(ie,ncols)
         441           0 :          do ie = 1, nelemd
         442           0 :             ncols = elem(ie)%idxP%NumUniquePts
         443           0 :             call putUniquePoints(elem(ie)%idxP, nlev, T_tmp(1:ncols,:,ie),       &
         444           0 :                elem(ie)%derived%fT(:,:,:))
         445           0 :             call putUniquePoints(elem(ie)%idxV, 2, nlev, uv_tmp(1:ncols,:,:,ie), &
         446           0 :                elem(ie)%derived%fM(:,:,:,:))
         447           0 :             call putUniquePoints(elem(ie)%idxV, nlev,pcnst, dq_tmp(1:ncols,:,:,ie), &
         448           0 :                elem(ie)%derived%fQ(:,:,:,:))
         449             :          end do
         450           0 :          call t_stopf('putUniquePoints')
         451             :       end if
         452             :    end if
         453             :
         454      369408 :    deallocate(T_tmp)
         455      369408 :    deallocate(uv_tmp)
         456      369408 :    deallocate(dq_tmp)
         457             :
         458             :    ! Boundary exchange for physics forcing terms.
         459             :    ! For physics on GLL grid, for points with duplicate degrees of freedom,
         460             :    ! putuniquepoints() set one of the element values and set the others to zero,
         461             :    ! so do a simple sum (boundary exchange with no weights).
         462             :    ! For physics grid, we interpolated into all points, so do weighted average.
         463             :
         464      369408 :    call t_startf('p_d_coupling:bndry_exchange')
         465             :
         466     2966808 :    do ie = 1, nelemd
         467     2597400 :       if (fv_nphys > 0) then
         468    70129800 :          do k = 1, nlev
         469             :             dyn_in%elem(ie)%derived%FM(:,:,1,k) =                          &
         470    67532400 :                  dyn_in%elem(ie)%derived%FM(:,:,1,k) *                     &
         471  1485712800 :                  dyn_in%elem(ie)%spheremp(:,:)
         472             :             dyn_in%elem(ie)%derived%FM(:,:,2,k) =                          &
         473           0 :                  dyn_in%elem(ie)%derived%FM(:,:,2,k) *                     &
         474  1418180400 :                  dyn_in%elem(ie)%spheremp(:,:)
         475             :             dyn_in%elem(ie)%derived%FT(:,:,k) =                            &
         476           0 :                  dyn_in%elem(ie)%derived%FT(:,:,k) *                       &
         477  1420777800 :                  dyn_in%elem(ie)%spheremp(:,:)
         478             :          end do
         479             :       end if
         480     2597400 :       kptr = 0
         481     2597400 :       call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FM(:,:,:,:), 2*nlev, kptr, ie)
         482     2597400 :       kptr = kptr + 2*nlev
         483     2597400 :       call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FT(:,:,:), nlev, kptr, ie)
         484     2966808 :       if (.not. use_cslam) then
         485             :          !
         486             :          ! if using CSLAM qdp is being overwritten with CSLAM values in the dynamics
         487             :          ! so no need to do boundary exchange of tracer tendency on GLL grid here
         488             :          !
         489           0 :          kptr = kptr + nlev
         490           0 :          call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FQ(:,:,:,:), nlev*qsize, kptr, ie)
         491             :       end if
         492             :    end do
         493             :
         494      369408 :    if (iam < par%nprocs) then
         495      369408 :      call bndry_exchange(par, edgebuf, location='p_d_coupling')
         496             :    end if
         497             :
         498     2966808 :    do ie = 1, nelemd
         499     2597400 :       kptr = 0
         500     2597400 :       call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FM(:,:,:,:), 2*nlev, kptr, ie)
         501     2597400 :       kptr = kptr + 2*nlev
         502     2597400 :       call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FT(:,:,:), nlev, kptr, ie)
         503     2597400 :       kptr = kptr + nlev
         504     2597400 :       if (.not. use_cslam) then
         505           0 :          call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FQ(:,:,:,:), nlev*qsize, kptr, ie)
         506             :       end if
         507     2966808 :       if (fv_nphys > 0) then
         508    70129800 :          do k = 1, nlev
         509             :             dyn_in%elem(ie)%derived%FM(:,:,1,k) =                             &
         510    67532400 :                  dyn_in%elem(ie)%derived%FM(:,:,1,k) *                        &
         511  1485712800 :                  dyn_in%elem(ie)%rspheremp(:,:)
         512             :             dyn_in%elem(ie)%derived%FM(:,:,2,k) =                             &
         513           0 :                  dyn_in%elem(ie)%derived%FM(:,:,2,k) *                        &
         514  1418180400 :                  dyn_in%elem(ie)%rspheremp(:,:)
         515             :             dyn_in%elem(ie)%derived%FT(:,:,k) =                               &
         516           0 :                  dyn_in%elem(ie)%derived%FT(:,:,k) *                          &
         517  1420777800 :                  dyn_in%elem(ie)%rspheremp(:,:)
         518             :          end do
         519             :       end if
         520             :    end do
         521      369408 :    call t_stopf('p_d_coupling:bndry_exchange')
         522             :
         523      369408 :    if (iam < par%nprocs .and. fv_nphys > 0) then
         524      738816 :       call test_mapping_output_mapped_tendencies(dyn_in%fvm(1:nelemd), elem(1:nelemd), &
         525     1108224 :                                                  1, nelemd, tl_f, tl_qdp)
         526             :    end if
         527      738816 : end subroutine p_d_coupling
         528             :
         529             : !=========================================================================================
         530             :
         531      370944 : subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d)
         532             :
         533             :    ! The ps, pdel, and q components of phys_state are all dry on input.
         534             :    ! On output the psdry and pdeldry components are initialized; ps and pdel are
         535             :    ! updated to contain contribution from water vapor only; the 'wet' constituent
         536             :    ! mixing ratios are converted to a wet basis.  Initialize geopotential heights.
         537             :    ! Finally compute energy and water column integrals of the physics input state.
         538             :
         539      369408 :    use constituents,    only: qmin
         540             :    use physconst,       only: gravit, zvir
         541             :    use cam_thermo,      only: cam_thermo_dry_air_update, cam_thermo_water_update
         542             :    use air_composition, only: thermodynamic_active_species_num
         543             :    use air_composition, only: thermodynamic_active_species_idx
         544             :    use air_composition, only: cpairv, rairv, cappav, dry_air_species_num
         545             :    use shr_const_mod,   only: shr_const_rwv
         546             :    use phys_control,    only: waccmx_is
         547             :    use geopotential,    only: geopotential_t
         548             :    use static_energy,   only: update_dry_static_energy_run
         549             :    use check_energy,    only: check_energy_timestep_init
         550             :    use hycoef,          only: hyai, ps0
         551             :    use shr_vmath_mod,   only: shr_vmath_log
         552             :    use qneg_module,     only: qneg3
         553             :    use dyn_tests_utils, only: vc_dry_pressure
         554             :    use shr_kind_mod,    only: shr_kind_cx
         555             :    ! arguments
         556             :    type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
         557             :    type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend
         558             :    type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         559             :
         560             :    ! local variables
         561             :    integer :: lchnk
         562             :
         563             :    real(r8) :: zvirv(pcols,pver)    ! Local zvir array pointer
         564             :    real(r8) :: factor_array(pcols,nlev)
         565             :
         566             :    integer :: m, i, k, ncol, m_cnst
         567      370944 :    type(physics_buffer_desc), pointer :: pbuf_chnk(:)
         568             :
         569             :    !Needed for "update_dry_static_energy" CCPP scheme
         570             :    integer :: errflg
         571             :    character(len=shr_kind_cx) :: errmsg
         572             :    !----------------------------------------------------------------------------
         573             :
         574             :    ! Evaluate derived quantities
         575             :
         576             :    !!$omp parallel do num_threads(horz_num_threads) private (lchnk, ncol, k, i, m , zvirv, pbuf_chnk, factor_array)
         577     1866312 :    do lchnk = begchunk,endchunk
         578             :
         579     1495368 :       ncol = get_ncols_p(lchnk)
         580             :
         581             :       ! dry pressure variables
         582             :
         583    24969168 :       do i = 1, ncol
         584   635287968 :          phys_state(lchnk)%psdry(i) = hyai(1)*ps0 + sum(phys_state(lchnk)%pdel(i,:))
         585             :       end do
         586    24969168 :       do i = 1, ncol
         587    24969168 :          phys_state(lchnk)%pintdry(i,1) = hyai(1)*ps0
         588             :       end do
         589     2990736 :       call shr_vmath_log(phys_state(lchnk)%pintdry(1:ncol,1), &
         590     4486104 :                          phys_state(lchnk)%lnpintdry(1:ncol,1),ncol)
         591    40374936 :       do k = 1, nlev
         592   649198368 :          do i = 1, ncol
         593  1220637600 :             phys_state(lchnk)%pintdry(i,k+1) = phys_state(lchnk)%pintdry(i,k) + &
         594  1869835968 :                                                phys_state(lchnk)%pdel(i,k)
         595             :          end do
         596    77759136 :          call shr_vmath_log(phys_state(lchnk)%pintdry(1:ncol,k+1),&
         597   118134072 :                             phys_state(lchnk)%lnpintdry(1:ncol,k+1),ncol)
         598             :       end do
         599             :
         600    40374936 :       do k=1,nlev
         601   649198368 :          do i=1,ncol
         602   610318800 :             phys_state(lchnk)%pdeldry (i,k) = phys_state(lchnk)%pdel(i,k)
         603   610318800 :             phys_state(lchnk)%rpdeldry(i,k) = 1._r8/phys_state(lchnk)%pdeldry(i,k)
         604           0 :             phys_state(lchnk)%pmiddry (i,k) = 0.5D0*(phys_state(lchnk)%pintdry(i,k+1) + &
         605   649198368 :                                                      phys_state(lchnk)%pintdry(i,k))
         606             :          end do
         607    77759136 :          call shr_vmath_log(phys_state(lchnk)%pmiddry(1:ncol,k), &
         608   118134072 :                             phys_state(lchnk)%lnpmiddry(1:ncol,k),ncol)
         609             :       end do
         610             :
         611             :       ! wet pressure variables (should be removed from physics!)
         612   662448024 :       factor_array(:,:) = 1.0_r8
         613     5981472 :       do m_cnst=dry_air_species_num + 1,thermodynamic_active_species_num
         614     4486104 :         m = thermodynamic_active_species_idx(m_cnst)
         615   122620176 :         do k=1,nlev
         616  1952081208 :           do i=1,ncol
         617             :             ! at this point all q's are dry
         618  1947595104 :             factor_array(i,k) = factor_array(i,k)+phys_state(lchnk)%q(i,k,m)
         619             :           end do
         620             :         end do
         621             :       end do
         622             :
         623    40374936 :       do k=1,nlev
         624   650693736 :          do i=1,ncol
         625   649198368 :             phys_state(lchnk)%pdel (i,k) = phys_state(lchnk)%pdeldry(i,k)*factor_array(i,k)
         626             :          end do
         627             :       end do
         628             :
         629             :       ! initialize vertical loop - model top pressure
         630             :
         631    24969168 :       do i=1,ncol
         632    23473800 :          phys_state(lchnk)%ps(i)     = phys_state(lchnk)%pintdry(i,1)
         633    24969168 :          phys_state(lchnk)%pint(i,1) = phys_state(lchnk)%pintdry(i,1)
         634             :       end do
         635    40374936 :       do k = 1, nlev
         636   649198368 :          do i=1,ncol
         637   610318800 :             phys_state(lchnk)%pint(i,k+1) =  phys_state(lchnk)%pint(i,k)+phys_state(lchnk)%pdel(i,k)
         638   610318800 :             phys_state(lchnk)%pmid(i,k)   = (phys_state(lchnk)%pint(i,k+1)+phys_state(lchnk)%pint(i,k))/2._r8
         639   649198368 :             phys_state(lchnk)%ps  (i)     =  phys_state(lchnk)%ps(i) + phys_state(lchnk)%pdel(i,k)
         640             :          end do
         641    38879568 :          call shr_vmath_log(phys_state(lchnk)%pint(1:ncol,k),phys_state(lchnk)%lnpint(1:ncol,k),ncol)
         642    40374936 :          call shr_vmath_log(phys_state(lchnk)%pmid(1:ncol,k),phys_state(lchnk)%lnpmid(1:ncol,k),ncol)
         643             :       end do
         644     1495368 :       call shr_vmath_log(phys_state(lchnk)%pint(1:ncol,pverp),phys_state(lchnk)%lnpint(1:ncol,pverp),ncol)
         645             :
         646    40374936 :       do k = 1, nlev
         647   650693736 :          do i = 1, ncol
         648   649198368 :             phys_state(lchnk)%rpdel(i,k)  = 1._r8/phys_state(lchnk)%pdel(i,k)
         649             :          end do
         650             :       end do
         651             :
         652             :       !------------------------------------------------------------
         653             :       ! Apply limiters to mixing ratios of major species (waccmx)
         654             :       !------------------------------------------------------------
         655     1495368 :       if (dry_air_species_num>0) then
         656           0 :         call physics_cnst_limit( phys_state(lchnk) )
         657             :         !-----------------------------------------------------------------------------
         658             :         ! Call cam_thermo_dry_air_update to compute cpairv, rairv, mbarv, and cappav as
         659             :         ! constituent dependent variables.
         660             :         ! Compute molecular viscosity(kmvis) and conductivity(kmcnd).
         661             :         ! Fill local zvirv variable; calculated for WACCM-X.
         662             :         !-----------------------------------------------------------------------------
         663           0 :         call cam_thermo_dry_air_update(phys_state(lchnk)%q, phys_state(lchnk)%t, lchnk, ncol)
         664           0 :         zvirv(:,:) = shr_const_rwv / rairv(:,:,lchnk) -1._r8
         665             :       else
         666   662448024 :         zvirv(:,:) = zvir
         667             :       end if
         668             :       !
         669             :       ! update cp_dycore in module air_composition.
         670             :       ! (note: at this point q is dry)
         671             :       !
         672     1495368 :       call cam_thermo_water_update(phys_state(lchnk)%q(1:ncol,:,:), lchnk, ncol, vc_dry_pressure)
         673    40374936 :       do k = 1, nlev
         674   650693736 :          do i = 1, ncol
         675   610318800 :             phys_state(lchnk)%exner(i,k) = (phys_state(lchnk)%pint(i,pver+1) &
         676  1259517168 :                                             / phys_state(lchnk)%pmid(i,k))**cappav(i,k,lchnk)
         677             :          end do
         678             :       end do
         679             :       !
         680             :       ! CAM physics: water tracers are moist; the rest dry
         681             :       !
         682   650693736 :       factor_array(1:ncol,1:nlev) = 1._r8/factor_array(1:ncol,1:nlev)
         683     5981472 :       do m = 1,pcnst
         684     5981472 :          if (cnst_type(m) == 'wet') then
         685   121124808 :             do k = 1, nlev
         686  1952081208 :                do i = 1, ncol
         687  1947595104 :                   phys_state(lchnk)%q(i,k,m) = factor_array(i,k)*phys_state(lchnk)%q(i,k,m)
         688             :                end do
         689             :             end do
         690             :          end if
         691             :       end do
         692             :
         693             :       ! Ensure tracers are all positive
         694             :       call qneg3('D_P_COUPLING',lchnk  ,ncol    ,pcols   ,pver    , &
         695     1495368 :            1, pcnst, qmin  ,phys_state(lchnk)%q)
         696             :
         697             :       ! Compute initial geopotential heights - based on full pressure
         698     1495368 :       call geopotential_t(phys_state(lchnk)%lnpint, phys_state(lchnk)%lnpmid  , phys_state(lchnk)%pint, &
         699             :          phys_state(lchnk)%pmid  , phys_state(lchnk)%pdel    , phys_state(lchnk)%rpdel                , &
         700             :          phys_state(lchnk)%t     , phys_state(lchnk)%q(:,:,:), rairv(:,:,lchnk), gravit, zvirv        , &
         701     2990736 :          phys_state(lchnk)%zi    , phys_state(lchnk)%zm      , ncol)
         702             :       ! Compute initial dry static energy, include surface geopotential
         703     2990736 :       call update_dry_static_energy_run(pver, gravit, phys_state(lchnk)%t(1:ncol,:),  &
         704     1495368 :                                         phys_state(lchnk)%zm(1:ncol,:),               &
         705     1495368 :                                         phys_state(lchnk)%phis(1:ncol),               &
         706     1495368 :                                         phys_state(lchnk)%s(1:ncol,:),                &
         707     8972208 :                                         cpairv(1:ncol,:,lchnk), errflg, errmsg)
         708             :       ! Compute energy and water integrals of input state
         709     1495368 :       pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
         710     1866312 :       call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk)
         711             :
         712             :
         713             :    end do  ! lchnk
         714             :
         715      741888 : end subroutine derived_phys_dry
         716             : end module dp_coupling

Generated by: LCOV version 1.14