LCOV - code coverage report
Current view: | top level - dynamics/se - gravity_waves_sources.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 67 | 74 | 90.5 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 3 | 3 | 100.0 %

              Line data    Source code
           1             : module gravity_waves_sources
           2             :   use derivative_mod, only: derivative_t
           3             :   use dimensions_mod, only: np,nlev
           4             :   use edgetype_mod,   only: EdgeBuffer_t
           5             :   use element_mod,    only: element_t
           6             :   use hybrid_mod,     only: hybrid_t
           7             :   use shr_kind_mod,   only: r8 => shr_kind_r8
           8             :
           9             :   implicit none
          10             :   private
          11             :   save
          12             :
          13             :   !! gravity_waves_sources created by S Santos, 10 Aug 2011
          14             :   !!
          15             :   !! gws_src_fnct starts parallel environment and computes frontogenesis
          16             :   !!   for use by WACCM (via dp_coupling)
          17             :
          18             :   public  :: gws_src_fnct
          19             :   public  :: gws_init
          20             :   private :: compute_frontogenesis
          21             :
          22             :   type (EdgeBuffer_t) :: edge3
          23             :   type (derivative_t)   :: deriv
          24             :   real(r8) :: psurf_ref
          25             :
          26             : !----------------------------------------------------------------------
          27             : CONTAINS
          28             : !----------------------------------------------------------------------
          29             :
          30      372480 :   subroutine gws_init(elem)
          31             :     use parallel_mod, only   : par
          32             :     use edge_mod, only       : initEdgeBuffer
          33             :     use hycoef, only         : hypi
          34             :     use pmgrid, only         : plev
          35             :     use thread_mod, only     : horz_num_threads
          36             :     implicit none
          37             :
          38             :     ! Elem will be needed for future updates to edge code
          39             :     type(element_t), pointer :: elem(:)
          40             :
          41             :     ! Set up variables similar to dyn_comp and prim_driver_mod initializations
          42        1536 :     call initEdgeBuffer(par, edge3, elem, 3*nlev,nthreads=1)
          43             :
          44        1536 :     psurf_ref = hypi(plev+1)
          45             :
          46        1536 :   end subroutine gws_init
          47             :
          48      370944 :   subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga,nphys)
          49        1536 :     use derivative_mod, only  : derivinit
          50             :     use dimensions_mod, only  : npsq, nelemd
          51             :     use dof_mod, only         : UniquePoints
          52             :     use hybrid_mod, only      : config_thread_region, get_loop_ranges
          53             :     use parallel_mod, only    : par
          54             :     use ppgrid, only          : pver
          55             :     use thread_mod, only      : horz_num_threads
          56             :     use dimensions_mod, only  : fv_nphys
          57             :     implicit none
          58             :     type (element_t), intent(inout), dimension(:) :: elem
          59             :     integer, intent(in)          :: tl, nphys, tlq
          60             :     real (kind=r8), intent(out) :: frontgf(nphys*nphys,pver,nelemd)
          61             :     real (kind=r8), intent(out) :: frontga(nphys*nphys,pver,nelemd)
          62             :
          63             :     ! Local variables
          64             :     type (hybrid_t) :: hybrid
          65             :     integer :: nets, nete, ithr, ncols, ie
          66      370944 :     real(kind=r8), allocatable  ::  frontgf_thr(:,:,:,:)
          67      370944 :     real(kind=r8), allocatable  ::  frontga_thr(:,:,:,:)
          68             :
          69             :     ! This does not need to be a thread private data-structure
          70      370944 :     call derivinit(deriv)
          71             :     !!$OMP PARALLEL NUM_THREADS(horz_num_threads),  DEFAULT(SHARED), PRIVATE(nets,nete,hybrid,ie,ncols,frontgf_thr,frontga_thr)
          72             : !    hybrid = config_thread_region(par,'horizontal')
          73      370944 :     hybrid = config_thread_region(par,'serial')
          74      370944 :     call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
          75             :
          76     2225664 :     allocate(frontgf_thr(nphys,nphys,nlev,nets:nete))
          77     1483776 :     allocate(frontga_thr(nphys,nphys,nlev,nets:nete))
          78      370944 :     call compute_frontogenesis(frontgf_thr,frontga_thr,tl,tlq,elem,deriv,hybrid,nets,nete,nphys)
          79      370944 :     if (fv_nphys>0) then
          80     2979144 :       do ie=nets,nete
          81     7824600 :         frontgf(:,:,ie) = RESHAPE(frontgf_thr(:,:,:,ie),(/nphys*nphys,nlev/))
          82     8195544 :         frontga(:,:,ie) = RESHAPE(frontga_thr(:,:,:,ie),(/nphys*nphys,nlev/))
          83             :       end do
          84             :     else
          85           0 :       do ie=nets,nete
          86           0 :         ncols = elem(ie)%idxP%NumUniquePts
          87           0 :         call UniquePoints(elem(ie)%idxP, nlev, frontgf_thr(:,:,:,ie), frontgf(1:ncols,:,ie))
          88           0 :         call UniquePoints(elem(ie)%idxP, nlev, frontga_thr(:,:,:,ie), frontga(1:ncols,:,ie))
          89             :       end do
          90             :     end if
          91      370944 :     deallocate(frontga_thr)
          92      370944 :     deallocate(frontgf_thr)
          93             :     !!$OMP END PARALLEL
          94             :
          95      741888 :   end subroutine gws_src_fnct
          96             :
          97      370944 :   subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets,nete,nphys)
          98             :   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          99             :   ! compute frontogenesis function F
         100             :   !   F =  -gradth dot C
         101             :   ! with:
         102             :   !   theta  = potential temperature
         103             :   !   gradth = grad(theta)
         104             :   !   C = ( gradth dot grad ) U
         105             :   !
         106             :   ! Original by Mark Taylor, July 2011
         107             :   ! Change by Santos, 10 Aug 2011:
         108             :   ! Integrated into gravity_waves_sources module, several arguments made global
         109             :   !  to prevent repeated allocation/initialization
         110             :   !
         111             :   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         112      370944 :     use physconst,      only: cappa
         113             :     use air_composition,only: dry_air_species_num, thermodynamic_active_species_num
         114             :     use air_composition,only: thermodynamic_active_species_idx_dycore
         115             :     use derivative_mod, only: gradient_sphere, ugradv_sphere
         116             :     use edge_mod,       only: edgevpack, edgevunpack
         117             :     use bndry_mod,      only: bndry_exchange
         118             :     use dyn_grid,       only: hvcoord
         119             :     use dimensions_mod, only: fv_nphys,ntrac
         120             :     use fvm_mapping,    only: dyn2phys_vector,dyn2phys
         121             :
         122             :     type(hybrid_t),     intent(in)            :: hybrid
         123             :     type(element_t),    intent(inout), target :: elem(:)
         124             :     type(derivative_t), intent(in)            :: ederiv
         125             :     integer,            intent(in)            :: nets,nete,nphys
         126             :     integer,            intent(in)            :: tl,tlq
         127             :     real(r8),           intent(out)           :: frontgf(nphys,nphys,nlev,nets:nete)
         128             :     real(r8),           intent(out)           :: frontga(nphys,nphys,nlev,nets:nete)
         129             :
         130             :     ! local
         131      741888 :     real(r8) :: area_inv(fv_nphys,fv_nphys), tmp(np,np)
         132      741888 :     real(r8) :: uv_tmp(fv_nphys*fv_nphys,2,nlev)
         133      741888 :     real(r8) :: frontgf_gll(np,np,nlev,nets:nete)
         134             :     real(r8) :: frontga_gll(np,np,nlev,nets:nete)
         135             :     integer  :: k,kptr,i,j,ie,component,h,nq,m_cnst
         136      741888 :     real(r8) :: gradth(np,np,2,nlev,nets:nete) ! grad(theta)
         137             :     real(r8) :: p(np,np)                       ! pressure at mid points
         138             :     real(r8) :: pint(np,np)                    ! pressure at interface points
         139             :     real(r8) :: theta(np,np)                   ! potential temperature at mid points
         140             :     real(r8) :: C(np,np,2), sum_water(np,np)
         141             :
         142     2979144 :     do ie=nets,nete
         143             :       ! pressure at model top
         144    54772200 :       pint(:,:) = hvcoord%hyai(1)*hvcoord%ps0
         145   245170800 :       do k=1,nlev
         146             :         ! moist pressure at mid points
         147  5093814600 :         sum_water(:,:) = 1.0_r8
         148  1697938200 :         do nq=dry_air_species_num+1,thermodynamic_active_species_num
         149  1455375600 :           m_cnst = thermodynamic_active_species_idx_dycore(nq)
         150             :           !
         151             :           ! make sure Q is updated
         152             :           !
         153 30805450200 :           sum_water(:,:) = sum_water(:,:) + elem(ie)%state%Qdp(:,:,k,m_cnst,tlq)/elem(ie)%state%dp3d(:,:,k,tl)
         154             :         end do
         155  5093814600 :         p(:,:) = pint(:,:) + 0.5_r8*sum_water(:,:)*elem(ie)%state%dp3d(:,:,k,tl)
         156             :         ! moist pressure at interface for next iteration
         157  5093814600 :         pint(:,:) = pint(:,:)+elem(ie)%state%dp3d(:,:,k,tl)
         158             :         !
         159  5093814600 :         theta(:,:) = elem(ie)%state%T(:,:,k,tl)*(psurf_ref / p(:,:))**cappa
         160             :         ! gradth(:,:,:,k,ie) = gradient_sphere(theta,ederiv,elem(ie)%Dinv)
         161   242562600 :         call gradient_sphere(theta,ederiv,elem(ie)%Dinv,gradth(:,:,:,k,ie))
         162             :         ! compute C = (grad(theta) dot grad ) u
         163   242562600 :         C(:,:,:) = ugradv_sphere(gradth(:,:,:,k,ie), elem(ie)%state%v(:,:,:,k,tl),ederiv,elem(ie))
         164             :         ! gradth dot C
         165  5093814600 :         frontgf_gll(:,:,k,ie) = -( C(:,:,1)*gradth(:,:,1,k,ie) +  C(:,:,2)*gradth(:,:,2,k,ie)  )
         166             :         ! apply mass matrix
         167  5093814600 :         gradth(:,:,1,k,ie)=gradth(:,:,1,k,ie)*elem(ie)%spheremp(:,:)
         168  5093814600 :         gradth(:,:,2,k,ie)=gradth(:,:,2,k,ie)*elem(ie)%spheremp(:,:)
         169  5096422800 :         frontgf_gll(:,:,k,ie)=frontgf_gll(:,:,k,ie)*elem(ie)%spheremp(:,:)
         170             :       enddo
         171             :       ! pack
         172     2608200 :       call edgeVpack(edge3, frontgf_gll(:,:,:,ie),nlev,0,ie)
         173     2979144 :       call edgeVpack(edge3, gradth(:,:,:,:,ie),2*nlev,nlev,ie)
         174             :     enddo
         175      370944 :     call bndry_exchange(hybrid,edge3,location='compute_frontogenesis')
         176     2979144 :     do ie=nets,nete
         177     2608200 :       call edgeVunpack(edge3, frontgf_gll(:,:,:,ie),nlev,0,ie)
         178     2608200 :       call edgeVunpack(edge3, gradth(:,:,:,:,ie),2*nlev,nlev,ie)
         179             :       ! apply inverse mass matrix,
         180   245170800 :       do k=1,nlev
         181  5093814600 :         gradth(:,:,1,k,ie)=gradth(:,:,1,k,ie)*elem(ie)%rspheremp(:,:)
         182  5093814600 :         gradth(:,:,2,k,ie)=gradth(:,:,2,k,ie)*elem(ie)%rspheremp(:,:)
         183  5096422800 :         frontgf_gll(:,:,k,ie)=frontgf_gll(:,:,k,ie)*elem(ie)%rspheremp(:,:)
         184             :       end do
         185     2979144 :       if (fv_nphys>0) then
         186     2608200 :         uv_tmp(:,:,:) = dyn2phys_vector(gradth(:,:,:,:,ie),elem(ie))
         187   245170800 :         do k=1,nlev
         188   242562600 :           h=0
         189   972858600 :           do j=1,fv_nphys
         190  3153313800 :             do i=1,fv_nphys
         191  2183063400 :               h=h+1
         192  2910751200 :               frontga(i,j,k,ie) = atan2 ( uv_tmp(h,2,k) , uv_tmp(h,1,k) + 1.e-10_r8 )
         193             :             end do
         194             :           end do
         195             :         end do
         196             :         !
         197             :         ! compute inverse physgrid area for mapping of scaler
         198             :         !
         199    54772200 :         tmp = 1.0_r8
         200     2608200 :         area_inv = dyn2phys(tmp,elem(ie)%metdet)
         201    33906600 :         area_inv = 1.0_r8/area_inv
         202   245170800 :         do k=1,nlev
         203   245170800 :           frontgf(:,:,k,ie) = dyn2phys(frontgf_gll(:,:,k,ie),elem(ie)%metdet,area_inv)
         204             :         end do
         205             :       else
         206           0 :         do k=1,nlev
         207           0 :           frontgf(:,:,k,ie)=frontgf_gll(:,:,k,ie)
         208             :           ! Frontogenesis angle
         209           0 :           frontga(:,:,k,ie) = atan2 ( gradth(:,:,2,k,ie) , gradth(:,:,1,k,ie) + 1.e-10_r8 )
         210             :         end do
         211             :       end if
         212             :     enddo
         213      370944 :   end subroutine compute_frontogenesis
         214             :
         215             :
         216             : end module gravity_waves_sources

Generated by: LCOV version 1.14