LCOV - code coverage report
Current view: | top level - dynamics/se - test_fvm_mapping.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 12 | 12 | 100.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 5 | 5 | 100.0 %

              Line data    Source code
           1             : module test_fvm_mapping
           2             :   use shr_kind_mod,           only: r8=>shr_kind_r8
           3             :   use fvm_control_volume_mod, only: fvm_struct
           4             :   use cam_history,            only: outfld
           5             :   use physconst,              only: pi
           6             :   use dimensions_mod,         only: np, nelemd, nlev, npsq, ntrac, use_cslam
           7             :   use element_mod,            only: element_t
           8             :   implicit none
           9             :   private
          10             :
          11             :   real(r8), parameter, private :: deg2rad = pi/180.0_r8
          12             :   real(r8), parameter, private :: psurf_moist = 100000.0_r8 !moist surface pressure
          13             :   integer,  parameter, private :: cl_idx  = 12
          14             :   integer,  parameter, private :: cl2_idx = 13
          15             :   real(r8), parameter, private :: cly_constant = 4.e-6_r8
          16             :   integer,  parameter, private :: num_fnc = 26
          17             :   integer,  parameter, private :: offset = 15
          18             :
          19             :   public :: test_mapping_overwrite_tendencies, test_mapping_addfld
          20             :   public :: test_mapping_output_mapped_tendencies, test_mapping_overwrite_dyn_state
          21             :   public :: test_mapping_output_phys_state
          22             : contains
          23             :
          24      372480 :   subroutine test_mapping_addfld
          25             : #ifdef debug_coupling
          26             :     use cam_history,        only: addfld, add_default, horiz_only, register_vector_field
          27             :     use constituents,       only: cnst_get_ind,cnst_name
          28             :     character(LEN=128) :: name
          29             :     integer :: nq,m_cnst
          30             :
          31             :     name = 'd2p_u_gll'
          32             :     call addfld(trim(name),   (/ 'lev' /),  'I','m/2','Exact zonal wind on GLL grid',gridname='GLL')
          33             :     !call add_default (trim(name), 1, ' ')
          34             :
          35             :     name = 'd2p_v_gll'
          36             :     call addfld(trim(name),   (/ 'lev' /),  'I','m/2','Exact meridional wind on GLL grid',gridname='GLL')
          37             :     !call add_default (trim(name), 1, ' ')
          38             :
          39             :     name = 'd2p_scalar_gll'
          40             :     call addfld(trim(name),   (/ 'lev' /),  'I','','Exact scalar on GLL grid',gridname='GLL')
          41             :     !call add_default (trim(name), 1, ' ')
          42             :
          43             :     name = 'd2p_u'
          44             :     call addfld(trim(name),   (/ 'lev' /),  'I','m/2','Zonal wind mapped to physics grid')
          45             :     !call add_default (trim(name), 1, ' ')
          46             :
          47             :     name = 'd2p_u_err'
          48             :     call addfld(trim(name),   (/ 'lev' /),  'I','m/2','Error in zonal wind mapped to physics grid')
          49             :     !call add_default (trim(name), 1, ' ')
          50             :
          51             :     name = 'd2p_v_err'
          52             :     call addfld(trim(name),   (/ 'lev' /),  'I','m/2','Error in meridional wind mapped to physics grid')
          53             :     !call add_default (trim(name), 1, ' ')
          54             :
          55             :     name = 'd2p_v'
          56             :     call addfld(trim(name),   (/ 'lev' /),  'I','m/s','Meridional wind mapped to physics grid')
          57             :     !call add_default (trim(name), 1, ' ')
          58             :
          59             :     name = 'd2p_scalar'
          60             :     call addfld(trim(name),   (/ 'lev' /),  'I','','Scalar mapped to physics grid')
          61             :     call add_default (trim(name), 1, ' ')
          62             :
          63             :     name = 'd2p_scalar_err'
          64             :     call addfld(trim(name),   (/ 'lev' /),  'I','','Error in scalar mapped to physics grid')
          65             :     call add_default (trim(name), 1, ' ')
          66             :
          67             :     do nq=ntrac,ntrac
          68             :       m_cnst = nq
          69             :       name = 'f2p_'//trim(cnst_name(m_cnst))//'_fvm'
          70             :       call addfld(trim(name),   (/ 'lev' /),  'I','','Exact water tracer on fvm grid',gridname='FVM')
          71             :       call add_default (trim(name), 1, ' ')
          72             :       name = 'f2p_'//trim(cnst_name(m_cnst))//'_err'
          73             :       call addfld(trim(name),   (/ 'lev' /),  'I','','Error in water tracer on physics grid (mapped from fvm grid)')
          74             :       call add_default (trim(name), 1, ' ')
          75             :       name = 'f2p_'//trim(cnst_name(m_cnst))//''
          76             :       call addfld(trim(name),   (/ 'lev' /),  'I','','Water tracer on physics grid (mapped from fvm grid')
          77             :       call add_default (trim(name), 1, ' ')
          78             :       !
          79             :       ! physgrid to gll (condensate loading tracers)
          80             :       !
          81             :       name = 'p2d_'//trim(cnst_name(m_cnst))//''
          82             :       call addfld(trim(name),   (/ 'lev' /),  'I','','Water tracer on physics grid')
          83             :       !call add_default (trim(name), 1, ' ')
          84             :       name = 'p2d_'//trim(cnst_name(m_cnst))//'_gll'
          85             :       call addfld(trim(name),   (/ 'lev' /),  'I','','Water tracer on GLL grid',gridname='GLL')
          86             :       !call add_default (trim(name), 1, ' ')
          87             :       name = 'p2d_'//trim(cnst_name(m_cnst))//'_err_gll'
          88             :       call addfld(trim(name),   (/ 'lev' /),  'I','','Error in water tracer mapped to GLL grid',gridname='GLL')
          89             :       !call add_default (trim(name), 1, ' ')
          90             :       !
          91             :       ! physgrid to fvm (condensate loading tracers)
          92             :       !
          93             :       name = 'p2f_'//trim(cnst_name(m_cnst))//''
          94             :       call addfld(trim(name),   (/ 'lev' /),  'I','','Water tracer on physics grid')
          95             :       call add_default (trim(name), 1, ' ')
          96             :       name = 'p2f_'//trim(cnst_name(m_cnst))//'_fvm'
          97             :       call addfld(trim(name),   (/ 'lev' /),  'I','','Water tracer on FVM grid',gridname='FVM')
          98             :       call add_default (trim(name), 1, ' ')
          99             :       name = 'p2f_'//trim(cnst_name(m_cnst))//'_err_fvm'
         100             :       call addfld(trim(name),   (/ 'lev' /),  'I','','Error in water tracer mapped to FVM grid',gridname='FVM')
         101             :       call add_default (trim(name), 1, ' ')
         102             :     end do
         103             :     !
         104             :     ! temperature tendency
         105             :     !
         106             :     name = 'p2d_ptend'
         107             :     call addfld(trim(name),   (/ 'lev' /),  'I','','T tendency on physics grid')
         108             :     !call add_default (trim(name), 1, ' ')
         109             :     name = 'p2d_ptend_gll'
         110             :     call addfld(trim(name),   (/ 'lev' /),  'I','','T tendency on GLL grid',gridname='GLL')
         111             :     !call add_default (trim(name), 1, ' ')
         112             :     name = 'p2d_ptend_err_gll'
         113             :     call addfld(trim(name),   (/ 'lev' /),  'I','','Error in T tendency mapped to GLL grid',gridname='GLL')
         114             :     !call add_default (trim(name), 1, ' ')
         115             :
         116             :     call addfld('p2d_u',   (/ 'lev' /),  'I','m/2','Zonal wind on physics grid')
         117             :     !call add_default ('p2d_u', 1, ' ')
         118             :     call addfld('p2d_v',   (/ 'lev' /),  'I','m/2','Meridional wind on physics grid')
         119             :     !call add_default ('p2d_v', 1, ' ')
         120             :     call addfld('p2d_u_gll',   (/ 'lev' /),  'I','m/2','Zonal wind on physics grid',gridname='GLL')
         121             :     !call add_default ('p2d_u_gll', 1, ' ')
         122             :     call addfld('p2d_v_gll',   (/ 'lev' /),  'I','m/2','Meridional wind on physics grid',gridname='GLL')
         123             :     !call add_default ('p2d_v_gll', 1, ' ')
         124             :     call addfld('p2d_u_gll_err',   (/ 'lev' /),  'I','m/2','Error in zonal wind interpolation to GLL grid',gridname='GLL')
         125             :     !call add_default ('p2d_u_gll_err', 1, ' ')
         126             :     call addfld('p2d_v_gll_err',   (/ 'lev' /),  'I','m/2','Error in meridional wind interpolation to GLL grid',&
         127             :          gridname='GLL')
         128             :     !call add_default ('p2d_v_gll_err', 1, ' ')
         129             :
         130             : !      name = 'phys2dyn_'//trim(cnst_name(m_cnst))//'_physgrid'
         131             : !      call outfld(trim(name),phys_state%q(:ncols,:,m_cnst),ncols,lchnk)
         132             : #endif
         133        1536 :   end subroutine test_mapping_addfld
         134             :
         135    23376600 :   subroutine test_mapping_overwrite_tendencies(phys_state,phys_tend,ncols,lchnk,q_prev,fvm)
         136             :     use dimensions_mod,         only: fv_nphys
         137             :     use constituents,           only: cnst_get_ind,pcnst,cnst_name
         138             :     use physics_types,  only: physics_state, physics_tend
         139             :     type(physics_state), intent(inout) :: phys_state
         140             :     type(physics_tend),  intent(inout) :: phys_tend
         141             :     real(r8), dimension(:,:,:), intent(inout) :: q_prev
         142             :     type(fvm_struct), intent(inout):: fvm(:)
         143             :     integer,          intent(in)   :: ncols,lchnk
         144             : #ifdef debug_coupling
         145             :     integer :: icol,k
         146             :     character(LEN=128) :: name
         147             :     integer :: m_cnst, nq, ie
         148             :
         149             :     q_prev(:,:,ntrac) = 0.0_r8
         150             :     phys_state%pdel(1:ncols,:) = phys_state%pdeldry(1:ncols,:) !make sure there is no conversion from wet to dry
         151             :     do nq=ntrac,ntrac
         152             :       m_cnst = nq
         153             :       do icol=1,ncols
         154             :         do k=1,num_fnc
         155             :           phys_state%q(icol,k,m_cnst)   = test_func(phys_state%lat(icol), phys_state%lon(icol), k, k)
         156             :         end do
         157             :       enddo
         158             :       name = 'p2f_'//trim(cnst_name(m_cnst))//''
         159             :       call outfld(trim(name),phys_state%q(:ncols,:,m_cnst),ncols,lchnk)
         160             :       name = 'p2d_'//trim(cnst_name(m_cnst))//''
         161             :       call outfld(trim(name),phys_state%q(:ncols,:,m_cnst),ncols,lchnk)
         162             :     end do
         163             :
         164             :     do icol=1,ncols
         165             :       do k=ntrac,ntrac
         166             :         phys_tend%dudt(icol,k) = test_func(phys_state%lat(icol), phys_state%lon(icol), k, k)
         167             :         phys_tend%dvdt(icol,k) = test_func(phys_state%lat(icol), phys_state%lon(icol), k, k)
         168             :         phys_tend%dtdt(icol,k) = test_func(phys_state%lat(icol), phys_state%lon(icol), k, k)
         169             :       end do
         170             :     enddo
         171             :     name = 'p2d_u'
         172             :     call outfld(trim(name),phys_tend%dudt(:ncols,:),ncols,lchnk)
         173             :     name = 'p2d_v'
         174             :     call outfld(trim(name),phys_tend%dvdt(:ncols,:),ncols,lchnk)
         175             :     name = 'p2d_ptend'
         176             :     call outfld(trim(name),phys_tend%dtdt(:ncols,:),ncols,lchnk)
         177             :
         178             :
         179             :     do icol=1,ncols
         180             :       do k=1,nlev
         181             : !        phys_tend%dudt(icol,k) = 0.0_r8
         182             : !        phys_tend%dvdt(icol,k) = 0.0_r8
         183             : !        phys_tend%dtdt(icol,k) = 0.0_r8
         184             :       end do
         185             :     enddo
         186             : #endif
         187    23376600 :   end subroutine test_mapping_overwrite_tendencies
         188             :
         189      369408 :   subroutine test_mapping_output_mapped_tendencies(fvm,elem,nets,nete,tl_f,tl_qdp)
         190    23376600 :     use dimensions_mod,         only: fv_nphys,nlev,nc
         191             :     use constituents,           only: cnst_get_ind,cnst_name
         192             :     integer,          intent(in)   :: nets,nete,tl_f,tl_qdp
         193             :     type(fvm_struct), intent(inout):: fvm(nets:nete)
         194             :     type(element_t),  intent(inout):: elem(nets:nete)             ! pointer to dyn_out element array
         195             : #ifdef debug_coupling
         196             :     integer :: ie,i,j,k
         197             :     character(LEN=128) :: name
         198             :     integer :: nq,m_cnst
         199             :     real(r8) :: diff(nc,nc,nlev,ntrac)
         200             :     diff = 0.0_r8
         201             :     do ie = nets,nete
         202             :       call outfld('p2d_u_gll', RESHAPE(elem(ie)%derived%fm(:,:,1,:),(/npsq,nlev/)), npsq, ie)
         203             :       call outfld('p2d_v_gll', RESHAPE(elem(ie)%derived%fm(:,:,2,:),(/npsq,nlev/)), npsq, ie)
         204             :       call outfld('p2d_ptend_gll', RESHAPE(elem(ie)%derived%ft(:,:,:),(/npsq,nlev/)), npsq, ie)
         205             :       do k=ntrac,ntrac
         206             :         do j=1,np
         207             :           do i=1,np
         208             :             elem(ie)%derived%fm(i,j,1,k)   = elem(ie)%derived%fm(i,j,1,k) -&
         209             :                  test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k,k)
         210             :             elem(ie)%derived%fm(i,j,2,k)   = elem(ie)%derived%fm(i,j,2,k) - &
         211             :                  test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k,k)
         212             :             elem(ie)%derived%ft(i,j,k)   = elem(ie)%derived%ft(i,j,k) - &
         213             :                  test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k,k)
         214             :           end do
         215             :         end do
         216             :       end do
         217             :       call outfld('p2d_u_gll_err'    , RESHAPE(elem(ie)%derived%fm(:,:,1,:),(/npsq,nlev/)), npsq, ie)
         218             :       call outfld('p2d_v_gll_err'    , RESHAPE(elem(ie)%derived%fm(:,:,2,:),(/npsq,nlev/)), npsq, ie)
         219             :       call outfld('p2d_ptend_err_gll', RESHAPE(elem(ie)%derived%ft(:,:,:),(/npsq,nlev/)), npsq, ie)
         220             :       elem(ie)%derived%ft(:,:,:)   = 0.0_r8
         221             :     end do
         222             :
         223             :     do ie = nets,nete
         224             :       do nq=ntrac,ntrac
         225             :         m_cnst = nq
         226             :         name = 'p2d_'//trim(cnst_name(m_cnst))//'_gll'
         227             :         call outfld(TRIM(name), RESHAPE(elem(ie)%derived%fq(:,:,:,nq),(/npsq,nlev/)), npsq, ie)
         228             :         !        call outfld(trim(name),&
         229             :         !             RESHAPE(fvm(ie)%fc(1:nc,1:nc,:,m_cnst),&
         230             :         !             (/nc*nc,nlev/)),nc*nc,ie)
         231             :         do k=1,num_fnc
         232             :           do j=1,np
         233             :             do i=1,np
         234             :               elem(ie)%derived%fq(i,j,k,nq) = elem(ie)%derived%fq(i,j,k,nq)-&
         235             :                    test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k, k)
         236             :             end do
         237             :           end do
         238             :         end do
         239             :         name = 'p2d_'//trim(cnst_name(m_cnst))//'_err_gll'
         240             :         call outfld(TRIM(name), RESHAPE(elem(ie)%derived%fq(:,:,:,nq),(/npsq,nlev/)), npsq, ie)
         241             :       end do
         242             :       if (use_cslam) then
         243             :         do nq=ntrac,ntrac
         244             :           m_cnst = nq
         245             :           name = 'p2f_'//trim(cnst_name(m_cnst))//'_fvm'
         246             :           !
         247             :           ! cly
         248             :           !
         249             : !          k=num_tracer+1
         250             : !          fvm(ie)%fc(1:nc,1:nc,k,:) = fvm(ie)%fc(1:nc,1:nc,cl_idx,:)+&
         251             : !                                    2.0_r8*fvm(ie)%fc(1:nc,1:nc,cl2_idx,:)
         252             :           call outfld(trim(name),&
         253             :                RESHAPE(fvm(ie)%fc(1:nc,1:nc,:,m_cnst)/fvm(ie)%dp_fvm(1:nc,1:nc,:),&
         254             :                (/nc*nc,nlev/)),nc*nc,ie)
         255             :           do k=1,num_fnc
         256             :             do j=1,nc
         257             :               do i=1,nc
         258             :                 diff(i,j,k,m_cnst) = fvm(ie)%fc(i,j,k,m_cnst)/fvm(ie)%dp_fvm(i,j,k)-&
         259             :                      test_func(fvm(ie)%center_cart(i,j)%lat,fvm(ie)%center_cart(i,j)%lon, k, k)
         260             :               end do
         261             :             end do
         262             :           end do
         263             :           name = 'p2f_'//trim(cnst_name(m_cnst))//'_err_fvm'
         264             :           call outfld(TRIM(name), RESHAPE(diff(:,:,:,m_cnst),(/nc*nc,nlev/)), nc*nc, ie)
         265             :
         266             :         end do
         267             :       endif
         268             :     end do
         269             : #endif
         270      369408 :   end subroutine test_mapping_output_mapped_tendencies
         271             :
         272      370944 :   subroutine test_mapping_overwrite_dyn_state(elem,fvm)
         273             :     use fvm_control_volume_mod, only: fvm_struct
         274             :     use constituents,           only: cnst_name
         275             :     use dimensions_mod,         only: nc,nhc
         276             :     use hybrid_mod,             only: get_loop_ranges, hybrid_t,config_thread_region
         277             :     use control_mod,            only: north, south, east, west, neast, nwest, seast, swest
         278             :     !    use fvm_mod,                only: fill_halo_fvm_noprealloc
         279             :     use fvm_mod,                only: fill_halo_fvm,ghostBufQnhc_h
         280             :     use parallel_mod,           only: par
         281             :     type (fvm_struct), intent(inout)    :: fvm(:)
         282             :     type(element_t), intent(inout) :: elem(:)             ! pointer to dyn_out element array
         283             : #ifdef debug_coupling
         284             :     integer            :: i,j,k,ie,nq,m_cnst
         285             :     character(LEN=128) :: name
         286             :     integer            :: nets,nete
         287             :     type(hybrid_t)                                                   :: hybrid
         288             :
         289             :     hybrid = config_thread_region(par,'serial')
         290             :     call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
         291             :     do ie=nets,nete
         292             :       do nq=ntrac,ntrac
         293             :         m_cnst = nq
         294             :         name = 'f2p_'//trim(cnst_name(m_cnst))//'_fvm'
         295             :         do k=1,num_fnc
         296             :           do j=1,nc
         297             :             do i=1,nc
         298             :               fvm(ie)%c(i,j,k,m_cnst) = test_func(fvm(ie)%center_cart(i,j)%lat,fvm(ie)%center_cart(i,j)%lon, k, k)
         299             :             end do
         300             :           end do
         301             :         end do
         302             :         !
         303             :         ! cly
         304             :         !
         305             : !        k=num_tracer+1
         306             : !        do j=1,nc
         307             : !          do i=1,nc
         308             : !            fvm(ie)%c(i,j,k,m_cnst) = test_func(fvm(ie)%center_cart(i,j)%lat,fvm(ie)%center_cart(i,j)%lon, k,cl_idx)+&
         309             : !                                 2.0_r8*test_func(fvm(ie)%center_cart(i,j)%lat,fvm(ie)%center_cart(i,j)%lon, k,cl2_idx)
         310             : !          end do
         311             : !        end do
         312             :         call outfld(TRIM(name), RESHAPE(fvm(ie)%c(1:nc,1:nc,:,m_cnst),(/nc*nc,nlev/)), nc*nc, ie)
         313             :       end do
         314             :
         315             :       elem(ie)%state%Qdp(:,:,:,:,:)   = 0.0_r8 !for testing the p2d map
         316             :       do k=1,num_fnc
         317             :         do j=1,np
         318             :           do i=1,np
         319             :             elem(ie)%state%v(i,j,1,k,:)   = test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k, k )
         320             :             elem(ie)%state%v(i,j,2,k,:)   = test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k, k)
         321             :           end do
         322             :         end do
         323             :       end do
         324             :       do k=1,num_fnc
         325             :         do j=1,np
         326             :           do i=1,np
         327             :             elem(ie)%derived%omega(i,j,k) = test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k, k)
         328             :           end do
         329             :         end do
         330             :       end do
         331             :       call outfld('d2p_scalar_gll', RESHAPE(elem(ie)%derived%omega(:,:,:)  ,(/npsq,nlev/)), npsq, ie)
         332             :       call outfld('d2p_u_gll', RESHAPE(elem(ie)%state%v(:,:,1,:,1),(/npsq,nlev/)), npsq, ie)
         333             :       call outfld('d2p_v_gll', RESHAPE(elem(ie)%state%v(:,:,2,:,1),(/npsq,nlev/)), npsq, ie)
         334             :     end do
         335             :     !
         336             :     ! do boundary exchange (this call should be indentical to call in prim_driver)
         337             :     !
         338             :     call fill_halo_fvm(ghostBufQnhc_h,elem,fvm,hybrid,nets,nete,nhc,1,nlev,nlev)
         339             :     do ie=nets,nete
         340             :       if (fvm(ie)%cubeboundary>4) then
         341             :         do k=ntrac,ntrac
         342             :           select case(fvm(ie)%cubeboundary)
         343             :           case (nwest)
         344             :             fvm(ie)%c(0,nc+1,:,k) = fvm(ie)%c(1,nc+1,:,k)
         345             :           case (swest)
         346             :             fvm(ie)%c(0,0,:,k) = fvm(ie)%c(0,1,:,k)
         347             :           case (seast)
         348             :             fvm(ie)%c(nc+1,0,:,k) = fvm(ie)%c(0,nc,:,k)
         349             :           case (neast)
         350             :             fvm(ie)%c(nc+1,nc+1,:,k) = fvm(ie)%c(nc,nc+1,:,k)
         351             :           end select
         352             :         end do
         353             :       end if
         354             :     end do
         355             : #endif
         356      370944 :   end subroutine test_mapping_overwrite_dyn_state
         357             :
         358      370944 :   subroutine test_mapping_output_phys_state(phys_state,fvm)
         359      370944 :     use physics_types, only: physics_state
         360             :     use ppgrid,        only: begchunk, endchunk, pver, pcols
         361             :     use constituents,  only: cnst_get_ind,cnst_name
         362             :     type(physics_state), intent(inout) :: phys_state(begchunk:endchunk)
         363             :     type(fvm_struct), pointer:: fvm(:)
         364             : #ifdef debug_coupling
         365             :     integer            :: lchnk, ncol,k,icol,m_cnst,nq,ie
         366             :     character(LEN=128) :: name
         367             :
         368             :     do lchnk = begchunk, endchunk
         369             :       call outfld('d2p_scalar', phys_state(lchnk)%omega(1:pcols,1:pver), pcols, lchnk)
         370             :       call outfld('d2p_u', phys_state(lchnk)%U(1:pcols,1:pver), pcols, lchnk)
         371             :       call outfld('d2p_v', phys_state(lchnk)%V(1:pcols,1:pver), pcols, lchnk)
         372             :       if (use_cslam) then
         373             :         do nq=ntrac,ntrac
         374             :           m_cnst = nq
         375             :           name = 'f2p_'//trim(cnst_name(m_cnst))
         376             :           !
         377             :           ! cly
         378             :           !
         379             :           !phys_state(lchnk)%q(1:pcols,num_tracer+1,m_cnst)=phys_state(lchnk)%q(1:pcols,cl_idx,m_cnst)+&
         380             :           !      2.0_r8*phys_state(lchnk)%q(1:pcols,12,m_cnst)
         381             :           call outfld(TRIM(name), phys_state(lchnk)%q(1:pcols,1:pver,m_cnst), pcols, lchnk)
         382             : !          k=num_tracer+1
         383             : !          do icol=1,phys_state(lchnk)%ncol
         384             : !            phys_state(lchnk)%q(icol,k,m_cnst) = phys_state(lchnk)%q(icol,cl_idx,m_cnst)+&
         385             : !                                          2.0_r8*phys_state(lchnk)%q(icol,cl2_idx,m_cnst)-&
         386             : !                                                 cly_constant
         387             : !          end do
         388             :           do k=1,num_fnc
         389             :             do icol=1,phys_state(lchnk)%ncol
         390             :               phys_state(lchnk)%q(icol,k,m_cnst) = phys_state(lchnk)%q(icol,k,m_cnst)&
         391             :                    -test_func(phys_state(lchnk)%lat(icol), phys_state(lchnk)%lon(icol), k,k)
         392             :             end do
         393             :           enddo
         394             :           name = 'f2p_'//trim(cnst_name(m_cnst))//'_err'
         395             :           call outfld(TRIM(name), phys_state(lchnk)%q(1:pcols,1:pver,m_cnst), pcols, lchnk)
         396             :           phys_state(lchnk)%q(1:pcols,1:pver,m_cnst) = 0.0_r8
         397             :         end do
         398             :       end if
         399             :     end do
         400             :
         401             :
         402             :     do lchnk = begchunk, endchunk
         403             :       do k=1,nlev
         404             :         do icol=1,phys_state(lchnk)%ncol
         405             :           phys_state(lchnk)%U(icol,k) = phys_state(lchnk)%U(icol,k)&
         406             :                -test_func(phys_state(lchnk)%lat(icol), phys_state(lchnk)%lon(icol), k, 9)
         407             :           phys_state(lchnk)%V(icol,k) = phys_state(lchnk)%V(icol,k)&
         408             :                -test_func(phys_state(lchnk)%lat(icol), phys_state(lchnk)%lon(icol), k,10)
         409             :         end do
         410             :       enddo
         411             :       name = 'd2p_u_err'
         412             :       call outfld(trim(name),phys_state(lchnk)%U(:pcols,:),pcols,lchnk)
         413             :       name = 'd2p_v_err'
         414             :       call outfld(trim(name),phys_state(lchnk)%V(:pcols,:),pcols,lchnk)
         415             :       do k=1,num_fnc
         416             :         do icol=1,phys_state(lchnk)%ncol
         417             :           phys_state(lchnk)%omega(icol,k) = phys_state(lchnk)%omega(icol,k)&
         418             :                -test_func(phys_state(lchnk)%lat(icol), phys_state(lchnk)%lon(icol), k,k)
         419             :         end do
         420             :       end do
         421             :       name = 'd2p_scalar_err'
         422             :       call outfld(trim(name),phys_state(lchnk)%omega(:pcols,:),pcols,lchnk)
         423             :     end do
         424             : #endif
         425      370944 :   end subroutine test_mapping_output_phys_state
         426             :
         427             : !  subroutine test_mapping_overwrite_state(phys_tend,nets,nete)
         428             : !#ifdef debug_coupling
         429             : !    phys_tend(lchnk)%dtdt(icol,ilyr)   = 0!test_func(phys_state(lchnk)%lat(icol),phys_state(lchnk)%lon(icol),ilyr,9)
         430             : !    phys_tend(lchnk)%dudt(icol,ilyr)   = 0!test_func(phys_state(lchnk)%lat(icol),phys_state(lchnk)%lon(icol),ilyr,12)
         431             : !    phys_tend(lchnk)%dvdt(icol,ilyr)   = 0!test_func(phys_state(lchnk)%lat(icol),phys_state(lchnk)%lon(icol),ilyr,13)
         432             : !    q_prev(icol, ilyr, 2:pcnst, lchnk) = 0.0D0
         433             : !    do m=2,pcnst
         434             : !      phys_state(lchnk)%Q(icol,ilyr,m)=0!test_func(phys_state(lchnk)%lat(icol),phys_state(lchnk)%lon(icol),ilyr,m)
         435             : !    end do
         436             : !#endif
         437             : !  end subroutine test_mapping_overwrite_state
         438             :
         439             : #ifdef debug_coupling
         440             :   function test_func(lat_in, lon_in, k, funcnum) result(fout)
         441             :     use hycoef,        only: hyai, hybi, hyam, hybm, ps0
         442             :     use shr_sys_mod,   only: shr_sys_flush
         443             :     use cam_abortutils, only: endrun
         444             :     real(r8), intent(in) :: lon_in
         445             :     real(r8), intent(in) :: lat_in
         446             :     integer,  intent(in) :: k
         447             :     integer,  intent(in) :: funcnum
         448             :     real(r8)             :: fout
         449             :     real(r8)             :: lon1,lat1,R0,Rg1,Rg2,lon2,lat2,cl,cl2
         450             :     real(r8)             :: eta_c
         451             :
         452             :     real(r8)             :: radius      = 10.0_r8 ! radius of the perturbation
         453             :     real(r8)             :: perturb_lon = 20.0_r8 ! longitudinal position, 20E
         454             :     real(r8)             :: perturb_lat = 40.0_r8 ! latitudinal position, 40N
         455             :     real(r8)             :: cos_tmp, sin_tmp, eta
         456             :     real(r8)             :: u_wind, v_wind, lat, lon, u_tmp, v_tmp
         457             :     real(r8)             :: rotation_angle
         458             :     real(r8)             :: det,r,k1,k2
         459             :     real(r8), parameter :: pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164_r8
         460             :     real(r8), parameter :: half_pi = pi*0.5_r8
         461             :     real(r8), parameter :: degrees_to_radians = pi/180.0_r8
         462             :     real(r8), parameter :: k1_lat_center =   20.d0*degrees_to_radians
         463             :     real(r8), parameter :: k1_lon_center =  300.d0*degrees_to_radians
         464             :
         465             :     lon = lon_in
         466             :     lat = lat_in
         467             :
         468             :
         469             :     select case(MOD(funcnum,8)+1)
         470             :     case(1)
         471             :       !
         472             :       !   Non-smooth scalar field (slotted cylinder)
         473             :       !
         474             :       R0 = 0.5_r8
         475             :       lon1 = 5.0_r8 * PI / 6.0_r8
         476             :       lat1 = 0.0_r8
         477             :       Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
         478             :       lon2 = 7.0_r8 * PI / 6.0_r8
         479             :       lat2 = 0.0_r8
         480             :       Rg2 = acos(sin(lat2)*sin(lat)+cos(lat2)*cos(lat)*cos(lon-lon2))
         481             :
         482             :       if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
         483             :         fout = 2.0_r8
         484             :       elseif ((Rg2 <= R0) .AND. (abs(lon-lon2) >= R0/6)) then
         485             :         fout = 2.0_r8
         486             :       elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
         487             :            .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
         488             :         fout = 2.0_r8
         489             :       elseif ((Rg2 <= R0) .AND. (abs(lon-lon2) < R0/6) &
         490             :            .AND. (lat-lat2 > 5.0_r8*R0/12.0_r8)) then
         491             :         fout = 2.0_r8
         492             :       else
         493             :         fout = 1.0_r8
         494             :       endif
         495             :     case(2)
         496             :       !
         497             :       ! Smooth Gaussian "ball"
         498             :       !
         499             :       R0    = 10.0_r8           ! radius of the perturbation
         500             :       lon1  = 20.0_r8*deg2rad   ! longitudinal position, 20E
         501             :       lat1  = 40.0_r8 *deg2rad  ! latitudinal position, 40N
         502             :       eta_c = 0.6_r8
         503             :       sin_tmp = SIN(lat1)*SIN(lat)
         504             :       cos_tmp = COS(lat1)*COS(lat)
         505             :       Rg1 = ACOS( sin_tmp + cos_tmp*COS(lon-lon1) )    ! great circle distance
         506             :       eta =  (hyam(k)*ps0 + hybm(k)*psurf_moist)/psurf_moist
         507             :       fout = EXP(- ((Rg1*R0)**2 + ((eta-eta_c)/0.1_r8)**2))
         508             :       if (ABS(fout) < 1.0E-8_r8) then
         509             :         fout = 0.0_r8
         510             :       end IF
         511             :     case(3)
         512             :       !
         513             :       !
         514             :       !
         515             :       fout = 0.5_r8 * ( tanh( 3.0_r8*abs(lat)-pi ) + 1.0_r8)
         516             :     case(4)
         517             :       fout = 2.0_r8+cos(5.0_r8+40*lon)!1.0e-8_r8
         518             :       fout = -0.5_r8-0.5_r8*(cos(16*lon)*(sin(2_r8*lat)**16))
         519             :     case(5)
         520             :       !
         521             :       ! approximately Y^2_2 spherical harmonic
         522             :       !
         523             :       fout = sin(lon)*cos(40*lat)!1.0e-8_r8
         524             :       fout = 0.5_r8*(cos(16*lon)*(sin(2_r8*lat)**16))
         525             :     case(6)
         526             :       !
         527             :       ! approximately Y32_16 spherical harmonic
         528             :       !
         529             :       fout = 0.5_r8 + 0.5_r8*(cos(16*lon)*(sin(2_r8*lat)**16))
         530             :     case(7)
         531             :       fout = 2.0_r8 + lat
         532             :     case(8)
         533             :       fout = 2.0_r8 + cos(lon)
         534             :     case(9)
         535             :       rotation_angle = 45.0_r8*pi/180.0_r8
         536             :       CALL regrot(lon_in,lat_in,lon,lat,0.0_r8,-0.5_r8*pi+rotation_angle,1)
         537             :       call Rossby_Haurwitz (lon, lat,u_wind, v_wind)
         538             :       CALL turnwi(u_wind,v_wind,u_tmp,v_tmp,lon_in,lat_in,lon,lat,0.0_r8,-0.5_r8*pi+rotation_angle,-1)
         539             :       fout = u_tmp
         540             :     case(10)
         541             :       rotation_angle = 45.0_r8*pi/180.0_r8
         542             :       CALL regrot(lon_in,lat_in,lon,lat,0.0_r8,-0.5_r8*pi+rotation_angle,1)
         543             :       call Rossby_Haurwitz (lon, lat,u_wind, v_wind)
         544             :       CALL turnwi(u_wind,v_wind,u_tmp,v_tmp,lon_in,lat_in,lon,lat,0.0_r8,-0.5_r8*pi+rotation_angle,-1)
         545             :       fout = v_tmp
         546             :     case(11)
         547             :       fout = 1.0E-8_r8
         548             :     case(12)
         549             :       !
         550             :       ! Terminator chemistry initial condition
         551             :       !
         552             :       k1 = 1.0_r8*max(0.d0,sin(lat)*sin(k1_lat_center) + cos(lat)*cos(k1_lat_center)*cos(lon-k1_lon_center))
         553             :       k2 = 1._r8
         554             :
         555             :       r = k1 / (4._r8*k2)
         556             :       det = sqrt(r*r + 2._r8*cly_constant*r)
         557             :
         558             :       fout  = (det-r)
         559             : !      fout = cly_constant/2._r8 - (det-r)/2._r8
         560             :     case(13)
         561             :       !
         562             :       ! Terminator chemistry initial condition
         563             :       !
         564             :       k1 = 1.0_r8*max(0.d0,sin(lat)*sin(k1_lat_center) + cos(lat)*cos(k1_lat_center)*cos(lon-k1_lon_center))
         565             :       k2 = 1._r8
         566             :
         567             :       r = k1 / (4._r8*k2)
         568             :       det = sqrt(r*r + 2._r8*cly_constant*r)
         569             :
         570             : !      fout  = (det-r)
         571             :       fout = cly_constant/2._r8 - (det-r)/2._r8
         572             :     case default
         573             : !      call endrun("Illegal funcnum_arg in test_func")
         574             :       fout = 1.0_r8
         575             :     end select
         576             :   end function test_func
         577             :
         578             :   function test_wind(lat, lon, iwind) result(fout)
         579             :    use cam_abortutils, only: endrun
         580             :     real(r8), intent(in) :: lon
         581             :     real(r8), intent(in) :: lat
         582             :     integer,  intent(in) :: iwind
         583             :
         584             :     real(r8)             :: fout
         585             :
         586             :
         587             :     fout = 0
         588             :   end function test_wind
         589             :
         590             :
         591             :   SUBROUTINE regrot(pxreg,pyreg,pxrot,pyrot,pxcen,pycen,kcall)
         592             :     use physconst, only: pi
         593             : !
         594             : !----------------------------------------------------------------------
         595             : !
         596             : !*    conversion between regular and rotated spherical coordinates.
         597             : !*
         598             : !*    pxreg     longitudes of the regular coordinates
         599             : !*    pyreg     latitudes of the regular coordinates
         600             : !*    pxrot     longitudes of the rotated coordinates
         601             : !*    pyrot     latitudes of the rotated coordinates
         602             : !*              all coordinates given in degrees n (negative for s)
         603             : !*              and degrees e (negative values for w)
         604             : !*    pxcen     regular longitude of the south pole of the rotated grid
         605             : !*    pycen     regular latitude of the south pole of the rotated grid
         606             : !*
         607             : !*    kcall=-1: find regular as functions of rotated coordinates.
         608             : !*    kcall= 1: find rotated as functions of regular coordinates.
         609             : !
         610             : !-----------------------------------------------------------------------
         611             : !
         612             :       integer kxdim,kydim,kx,ky,kcall
         613             :       real(r8) :: pxreg,pyreg,&
         614             :                   pxrot,pyrot,&
         615             :                   pxcen,pycen
         616             : !
         617             : !-----------------------------------------------------------------------
         618             : !
         619             :       real(r8) zsycen,zcycen,zxmxc,zsxmxc,zcxmxc,zsyreg,zcyreg, &
         620             :                zsyrot,zcyrot,zcxrot,zsxrot,zpi,zpih
         621             :       integer jy,jx
         622             :
         623             :       zpih = pi*0.5_r8
         624             : !
         625             :       !----------------------------------------------------------------------
         626             : !
         627             :       zsycen = SIN((pycen+zpih))
         628             :       zcycen = COS((pycen+zpih))
         629             : !
         630             :       IF (kcall.eq.1) then
         631             : !
         632             :          zxmxc  = pxreg - pxcen
         633             :          zsxmxc = SIN(zxmxc)
         634             :          zcxmxc = COS(zxmxc)
         635             :          zsyreg = SIN(pyreg)
         636             :          zcyreg = COS(pyreg)
         637             :          zsyrot = zcycen*zsyreg - zsycen*zcyreg*zcxmxc
         638             :          zsyrot = max(zsyrot,-1.0_r8)
         639             :          zsyrot = min(zsyrot,+1.0_r8)
         640             :          !
         641             :          pyrot = ASIN(zsyrot)
         642             :          !
         643             :          zcyrot = COS(pyrot)
         644             :          zcxrot = (zcycen*zcyreg*zcxmxc +zsycen*zsyreg)/zcyrot
         645             :          zcxrot = max(zcxrot,-1.0_r8)
         646             :          zcxrot = min(zcxrot,+1.0_r8)
         647             :          zsxrot = zcyreg*zsxmxc/zcyrot
         648             :          !
         649             :          pxrot = ACOS(zcxrot)
         650             :          !
         651             :          IF (zsxrot < 0.0_r8) then
         652             :            pxrot = -pxrot
         653             :          end IF
         654             :                !
         655             :       ELSEIF (kcall.eq.-1) then
         656             :          !
         657             :          zsxrot = SIN(pxrot)
         658             :          zcxrot = COS(pxrot)
         659             :          zsyrot = SIN(pyrot)
         660             :          zcyrot = COS(pyrot)
         661             :          zsyreg = zcycen*zsyrot + zsycen*zcyrot*zcxrot
         662             :          zsyreg = max(zsyreg,-1.0_r8)
         663             :          zsyreg = min(zsyreg,+1.0_r8)
         664             :          !
         665             :          pyreg = ASIN(zsyreg)
         666             :          !
         667             :          zcyreg = COS(pyreg)
         668             :          zcxmxc = (zcycen*zcyrot*zcxrot -&
         669             :               zsycen*zsyrot)/zcyreg
         670             :          zcxmxc = max(zcxmxc,-1.0_r8)
         671             :          zcxmxc = min(zcxmxc,+1.0_r8)
         672             :          zsxmxc = zcyrot*zsxrot/zcyreg
         673             :          zxmxc  = ACOS(zcxmxc)
         674             :          IF (zsxmxc < 0.0_r8) then
         675             :            zxmxc = -zxmxc
         676             :          end IF
         677             :          !
         678             :          pxreg = zxmxc + pxcen
         679             :          !
         680             :       ELSE
         681             :          WRITE(6,'(1x,''invalid kcall in regrot'')')
         682             :          STOP
         683             :       ENDIF
         684             :     END SUBROUTINE regrot
         685             :
         686             :   SUBROUTINE turnwi(puarg,pvarg,pures,pvres,pxreg,pyreg,pxrot,pyrot,pxcen,pycen,kcall)
         687             :     use physconst, only: pi
         688             :     !
         689             :     !-----------------------------------------------------------------------
         690             :     !
         691             :     !*    turn horizontal velocity components between regular and
         692             :     !*    rotated spherical coordinates.
         693             :     !
         694             :     !*    puarg : input u components
         695             :     !*    pvarg : input v components
         696             :     !*    pures : output u components
         697             :     !*    pvres : output v components
         698             :     !*    pa    : transformation coefficients
         699             :     !*    pb    :    -"-
         700             :     !*    pc    :    -"-
         701             :     !*    pd    :    -"-
         702             :     !*    pxreg : regular longitudes
         703             :     !*    pyreg : regular latitudes
         704             :     !*    pxrot : rotated longitudes
         705             :     !*    pyrot : rotated latitudes
         706             :     !*    kxdim              : dimension in the x (longitude) direction
         707             :     !*    kydim              : dimension in the y (latitude) direction
         708             :     !*    kx                 : number of gridpoints in the x direction
         709             :     !*    ky                 : number of gridpoints in the y direction
         710             :     !*    pxcen              : regular longitude of the south pole of the
         711             :     !*                         transformed grid
         712             :     !*    pycen              : regular latitude of the south pole of the
         713             :     !*                         transformed grid
         714             :     !*
         715             :     !*    kcall < 0          : find wind components in regular coordinates
         716             :     !*                         from wind components in rotated coordinates
         717             :     !*    kcall > 0          : find wind components in rotated coordinates
         718             :     !*                         from wind components in regular coordinates
         719             :     !*    note that all coordinates are given in degrees n and degrees e.
         720             :     !*       (negative values for s and w)
         721             :     !
         722             :     !-----------------------------------------------------------------------
         723             :
         724             :     integer kxdim,kydim,kx,ky,kcall
         725             :     real(r8) puarg,pvarg,    &
         726             :          pures,pvres,    &
         727             :          pa,   pb,       &
         728             :          pc,   pd,       &
         729             :          pxreg,pyreg,    &
         730             :          pxrot,pyrot
         731             :       real(r8) pxcen,pycen
         732             :       !
         733             :       !-----------------------------------------------------------------------
         734             :       !
         735             :       integer jy,jx
         736             :       real(r8) zpih,zsyc,zcyc,zsxreg,zcxreg,zsyreg,zcyreg,zxmxc,&
         737             :            zsxmxc,zcxmxc,zsxrot,zcxrot,zsyrot,zcyrot
         738             :       !
         739             :       !-----------------------------------------------------------------------
         740             :       !
         741             :       IF (kcall.eq.1) then
         742             :         zpih = pi*0.5_r8
         743             :         zsyc = SIN(pycen+zpih)
         744             :         zcyc = COS(pycen+zpih)
         745             :         !
         746             :         zsxreg = SIN(pxreg)
         747             :         zcxreg = COS(pxreg)
         748             :         zsyreg = SIN(pyreg)
         749             :         zcyreg = COS(pyreg)
         750             :         !
         751             :         zxmxc  = pxreg - pxcen
         752             :         zsxmxc = SIN(zxmxc)
         753             :         zcxmxc = COS(zxmxc)
         754             :         !
         755             :         zsxrot = SIN(pxrot)
         756             :         zcxrot = COS(pxrot)
         757             :         zsyrot = SIN(pyrot)
         758             :         zcyrot = COS(pyrot)
         759             :         !
         760             :         pa = zcyc*zsxmxc*zsxrot + zcxmxc*zcxrot
         761             :         pb = zcyc*zcxmxc*zsyreg*zsxrot - zsyc*zcyreg*zsxrot - &
         762             :              zsxmxc*zsyreg*zcxrot
         763             :         pc = zsyc*zsxmxc/zcyrot
         764             :         pd = (zsyc*zcxmxc*zsyreg + zcyc*zcyreg)/zcyrot
         765             :         !
         766             :         pures = pa*puarg + pb*pvarg
         767             :         pvres = pc*puarg + pd*pvarg
         768             :       ELSEIF (kcall.eq.-1) then
         769             :         zpih = pi*0.5_r8
         770             :         zsyc = SIN(pycen+zpih)
         771             :         zcyc = COS(pycen+zpih)
         772             :         !
         773             :         zsxreg = SIN(pxreg)
         774             :         zcxreg = COS(pxreg)
         775             :         zsyreg = SIN(pyreg)
         776             :         zcyreg = COS(pyreg)
         777             :         !
         778             :         zxmxc  = pxreg - pxcen
         779             :         zsxmxc = SIN(zxmxc)
         780             :         zcxmxc = COS(zxmxc)
         781             :         !
         782             :         zsxrot = SIN(pxrot)
         783             :         zcxrot = COS(pxrot)
         784             :         zsyrot = SIN(pyrot)
         785             :         zcyrot = COS(pyrot)
         786             :         !
         787             :         pa = zcxmxc*zcxrot + zcyc*zsxmxc*zsxrot
         788             :         pb = zcyc*zsxmxc*zcxrot*zsyrot + zsyc*zsxmxc*zcyrot -&
         789             :              zcxmxc*zsxrot*zsyrot
         790             :         pc =-zsyc*zsxrot/zcyreg
         791             :         pd = (zcyc*zcyrot - zsyc*zcxrot*zsyrot)/zcyreg
         792             :         !
         793             :         pures = pa*puarg + pb*pvarg
         794             :         pvres = pc*puarg + pd*pvarg
         795             :       ELSE
         796             :         write(6,'(1x,''invalid kcall in turnwi'')')
         797             :         STOP
         798             :       ENDIF
         799             :     END SUBROUTINE turnwi
         800             :
         801             :   SUBROUTINE Rossby_Haurwitz (lon, lat,u_wind, v_wind)
         802             :     use physconst, only: rearth
         803             : !-----------------------------------------------------------------------
         804             : !     input parameters
         805             : !-----------------------------------------------------------------------
         806             :       real(r8), intent(in)  :: lon,                                     & ! longitude in radians
         807             :                                lat                                        ! latitude in radians
         808             :                                                                           ! both coefficients 'a' and 'b' are needed at the full model level
         809             : !-----------------------------------------------------------------------
         810             : !     input parameters
         811             : !-----------------------------------------------------------------------
         812             :       real(r8), intent(out) :: u_wind,                                  & ! zonal wind in m/s
         813             :                                v_wind                                     ! meridional wind in m/s
         814             :
         815             : !-----------------------------------------------------------------------
         816             : !     test case parameters
         817             : !-----------------------------------------------------------------------
         818             :       real(r8),parameter :: u0      = 50._r8,                          &   ! reference wind
         819             :                             n       = 4._r8                                ! wavenumber
         820             :
         821             : !-----------------------------------------------------------------------
         822             : !     local
         823             : !-----------------------------------------------------------------------
         824             :       real(r8) :: tmp1, tmp2, tmp3, KK, MM
         825             :       real(r8) :: sin_lat, cos_lat, sin_slat, cos_slat
         826             :
         827             : !-----------------------------------------------------------------------
         828             : !     initialize the wind components
         829             : !-----------------------------------------------------------------------
         830             :       MM      = u0/(n*rearth)   ! parameter M
         831             :       KK      = u0/(n*rearth)   ! parameter K
         832             :
         833             :
         834             :       cos_lat = cos(lat)
         835             :       sin_lat = sin(lat)
         836             :       tmp1 = rearth * MM * cos_lat
         837             :       tmp2 = rearth * KK * cos_lat**(n-1._r8)*(n*sin_lat**2 - cos_lat**2)
         838             :       tmp3 = -rearth * KK * n * cos_lat**(n-1._r8) * sin_lat
         839             :       u_wind = tmp1 + tmp2 * cos(n*lon)
         840             :       v_wind = tmp3 * sin(n*lon)
         841             :   end subroutine Rossby_Haurwitz
         842             :
         843             : #endif
         844             : end module test_fvm_mapping

Generated by: LCOV version 1.14