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           0 :   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           0 :     call initEdgeBuffer(par, edge3, elem, 3*nlev,nthreads=1)
      43             : 
      44           0 :     psurf_ref = hypi(plev+1)
      45             : 
      46           0 :   end subroutine gws_init
      47             : 
      48           0 :   subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga,nphys)
      49           0 :     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           0 :     real(kind=r8), allocatable  ::  frontgf_thr(:,:,:,:)
      67           0 :     real(kind=r8), allocatable  ::  frontga_thr(:,:,:,:)
      68             : 
      69             :     ! This does not need to be a thread private data-structure
      70           0 :     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           0 :     hybrid = config_thread_region(par,'serial')
      74           0 :     call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
      75             : 
      76           0 :     allocate(frontgf_thr(nphys,nphys,nlev,nets:nete))
      77           0 :     allocate(frontga_thr(nphys,nphys,nlev,nets:nete))
      78           0 :     call compute_frontogenesis(frontgf_thr,frontga_thr,tl,tlq,elem,deriv,hybrid,nets,nete,nphys)
      79           0 :     if (fv_nphys>0) then
      80           0 :       do ie=nets,nete
      81           0 :         frontgf(:,:,ie) = RESHAPE(frontgf_thr(:,:,:,ie),(/nphys*nphys,nlev/))
      82           0 :         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           0 :     deallocate(frontga_thr)
      92           0 :     deallocate(frontgf_thr)
      93             :     !!$OMP END PARALLEL
      94             : 
      95           0 :   end subroutine gws_src_fnct
      96             : 
      97           0 :   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           0 :     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           0 :     real(r8) :: area_inv(fv_nphys,fv_nphys), tmp(np,np)
     132           0 :     real(r8) :: uv_tmp(fv_nphys*fv_nphys,2,nlev)
     133           0 :     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           0 :     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           0 :     do ie=nets,nete
     143             :       ! pressure at model top
     144           0 :       pint(:,:) = hvcoord%hyai(1)*hvcoord%ps0
     145           0 :       do k=1,nlev
     146             :         ! moist pressure at mid points
     147           0 :         sum_water(:,:) = 1.0_r8
     148           0 :         do nq=dry_air_species_num+1,thermodynamic_active_species_num
     149           0 :           m_cnst = thermodynamic_active_species_idx_dycore(nq)
     150             :           !
     151             :           ! make sure Q is updated
     152             :           !
     153           0 :           sum_water(:,:) = sum_water(:,:) + elem(ie)%state%Qdp(:,:,k,m_cnst,tlq)/elem(ie)%state%dp3d(:,:,k,tl)
     154             :         end do
     155           0 :         p(:,:) = pint(:,:) + 0.5_r8*sum_water(:,:)*elem(ie)%state%dp3d(:,:,k,tl)
     156             :         ! moist pressure at interface for next iteration
     157           0 :         pint(:,:) = pint(:,:)+elem(ie)%state%dp3d(:,:,k,tl)
     158             :         !
     159           0 :         theta(:,:) = elem(ie)%state%T(:,:,k,tl)*(psurf_ref / p(:,:))**cappa
     160             :         ! gradth(:,:,:,k,ie) = gradient_sphere(theta,ederiv,elem(ie)%Dinv)
     161           0 :         call gradient_sphere(theta,ederiv,elem(ie)%Dinv,gradth(:,:,:,k,ie))
     162             :         ! compute C = (grad(theta) dot grad ) u
     163           0 :         C(:,:,:) = ugradv_sphere(gradth(:,:,:,k,ie), elem(ie)%state%v(:,:,:,k,tl),ederiv,elem(ie))
     164             :         ! gradth dot C
     165           0 :         frontgf_gll(:,:,k,ie) = -( C(:,:,1)*gradth(:,:,1,k,ie) +  C(:,:,2)*gradth(:,:,2,k,ie)  )
     166             :         ! apply mass matrix
     167           0 :         gradth(:,:,1,k,ie)=gradth(:,:,1,k,ie)*elem(ie)%spheremp(:,:)
     168           0 :         gradth(:,:,2,k,ie)=gradth(:,:,2,k,ie)*elem(ie)%spheremp(:,:)
     169           0 :         frontgf_gll(:,:,k,ie)=frontgf_gll(:,:,k,ie)*elem(ie)%spheremp(:,:)
     170             :       enddo
     171             :       ! pack
     172           0 :       call edgeVpack(edge3, frontgf_gll(:,:,:,ie),nlev,0,ie)
     173           0 :       call edgeVpack(edge3, gradth(:,:,:,:,ie),2*nlev,nlev,ie)
     174             :     enddo
     175           0 :     call bndry_exchange(hybrid,edge3,location='compute_frontogenesis')
     176           0 :     do ie=nets,nete
     177           0 :       call edgeVunpack(edge3, frontgf_gll(:,:,:,ie),nlev,0,ie)
     178           0 :       call edgeVunpack(edge3, gradth(:,:,:,:,ie),2*nlev,nlev,ie)
     179             :       ! apply inverse mass matrix,
     180           0 :       do k=1,nlev
     181           0 :         gradth(:,:,1,k,ie)=gradth(:,:,1,k,ie)*elem(ie)%rspheremp(:,:)
     182           0 :         gradth(:,:,2,k,ie)=gradth(:,:,2,k,ie)*elem(ie)%rspheremp(:,:)
     183           0 :         frontgf_gll(:,:,k,ie)=frontgf_gll(:,:,k,ie)*elem(ie)%rspheremp(:,:)
     184             :       end do
     185           0 :       if (fv_nphys>0) then
     186           0 :         uv_tmp(:,:,:) = dyn2phys_vector(gradth(:,:,:,:,ie),elem(ie))
     187           0 :         do k=1,nlev
     188           0 :           h=0
     189           0 :           do j=1,fv_nphys
     190           0 :             do i=1,fv_nphys
     191           0 :               h=h+1
     192           0 :               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           0 :         tmp = 1.0_r8
     200           0 :         area_inv = dyn2phys(tmp,elem(ie)%metdet)
     201           0 :         area_inv = 1.0_r8/area_inv
     202           0 :         do k=1,nlev
     203           0 :           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           0 :   end subroutine compute_frontogenesis
     214             : 
     215             : 
     216             : end module gravity_waves_sources