LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - prim_init.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 168 | 195 | 86.2 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 1 | 1 | 100.0 %

              Line data    Source code
           1             : module prim_init
           2             :
           3             :   use shr_kind_mod,   only: r8=>shr_kind_r8
           4             :   use dimensions_mod, only: nc, use_cslam
           5             :   use reduction_mod,  only: reductionbuffer_ordered_1d_t
           6             :   use quadrature_mod, only: quadrature_t, gausslobatto
           7             :
           8             :   implicit none
           9             :   private
          10             :   save
          11             :
          12             :   public :: prim_init1
          13             :
          14             :   real(r8), public :: fvm_corners(nc+1) ! fvm cell corners on reference element
          15             :   real(r8), public :: fvm_points(nc)    ! fvm cell centers on reference element
          16             :
          17             :   type (quadrature_t), public               :: gp    ! element GLL points
          18             :   type (ReductionBuffer_ordered_1d_t)       :: red   ! reduction buffer (shared)
          19             :
          20             : contains
          21        1536 :   subroutine prim_init1(elem, fvm, par, Tl)
          22             :     use cam_logfile,            only: iulog
          23             :     use shr_sys_mod,            only: shr_sys_flush
          24             :     use thread_mod,             only: max_num_threads
          25             :     use dimensions_mod,         only: np, nlev, nelem, nelemd, nelemdmax, qsize_d
          26             :     use dimensions_mod,         only: GlobalUniqueCols, fv_nphys,irecons_tracer
          27             :     use control_mod,            only: topology, partmethod
          28             :     use element_mod,            only: element_t, allocate_element_desc
          29             :     use fvm_mod,                only: fvm_init1
          30             :     use mesh_mod,               only: MeshUseMeshFile
          31             :     use se_dyn_time_mod,        only: timelevel_init, timelevel_t
          32             :     use mass_matrix_mod,        only: mass_matrix
          33             :     use derivative_mod,         only: allocate_subcell_integration_matrix_cslam
          34             :     use derivative_mod,         only: allocate_subcell_integration_matrix_physgrid
          35             :     use cube_mod,               only: cubeedgecount , cubeelemcount, cubetopology
          36             :     use cube_mod,               only: cube_init_atomic, rotation_init_atomic, set_corner_coordinates
          37             :     use cube_mod,               only: assign_node_numbers_to_elem
          38             :     use mesh_mod,               only: MeshSetCoordinates, MeshUseMeshFile, MeshCubeTopology
          39             :     use mesh_mod,               only: MeshCubeElemCount, MeshCubeEdgeCount
          40             :     use metagraph_mod,          only: metavertex_t, localelemcount, initmetagraph, printmetavertex
          41             :     use gridgraph_mod,          only: gridvertex_t, gridedge_t
          42             :     use gridgraph_mod,          only: allocate_gridvertex_nbrs, deallocate_gridvertex_nbrs
          43             :     use schedtype_mod,          only: schedule
          44             :     use schedule_mod,           only: genEdgeSched
          45             :     use prim_advection_mod,     only: prim_advec_init1
          46             :     use cam_abortutils,         only: endrun
          47             :     use spmd_utils,             only: mpi_integer, mpi_max
          48             :     use parallel_mod,           only: parallel_t, syncmp, global_shared_buf, nrepro_vars
          49             :     use spacecurve_mod,         only: genspacepart
          50             :     use dof_mod,                only: global_dof, CreateUniqueIndex, SetElemOffset
          51             :     use params_mod,             only: SFCURVE
          52             :     use physconst,              only: pi
          53             :     use reduction_mod,          only: red_min, red_max, red_max_int, red_flops
          54             :     use reduction_mod,          only: red_sum, red_sum_int, initreductionbuffer
          55             :     use infnan,                 only: nan, assignment(=)
          56             :     use shr_reprosum_mod,       only: repro_sum => shr_reprosum_calc
          57             :     use fvm_analytic_mod,       only: compute_basic_coordinate_vars
          58             :     use fvm_control_volume_mod, only: fvm_struct, allocate_physgrid_vars
          59             :     use air_composition,        only: thermodynamic_active_species_num
          60             :
          61             :     type(element_t),  pointer        :: elem(:)
          62             :     type(fvm_struct), pointer        :: fvm(:)
          63             :     type(parallel_t),  intent(inout) :: par
          64             :     type(timelevel_t), intent(out)   :: Tl
          65             :
          66             :     ! Local Variables
          67        1536 :     type (GridVertex_t), target,allocatable :: GridVertex(:)
          68        1536 :     type (GridEdge_t),   target,allocatable :: Gridedge(:)
          69        1536 :     type (MetaVertex_t), target,allocatable :: MetaVertex(:)
          70             :
          71             :     integer                        :: ie
          72             :     integer                        :: nets, nete
          73             :     integer                        :: nelem_edge
          74             :     integer                        :: ierr=0, j
          75             :     logical,           parameter   :: Debug = .FALSE.
          76             :
          77        1536 :     real(r8),          allocatable :: aratio(:,:)
          78             :     real(r8)                       :: area(1), xtmp
          79             :     character(len=80)              :: rot_type ! cube edge rotation type
          80             :
          81             :     integer                        :: i
          82             :
          83             :     character(len=128)             :: errmsg
          84             :     character(len=*),  parameter   :: subname = 'PRIM_INIT1: '
          85             :
          86             :     ! ====================================
          87             :     ! Set cube edge rotation type for model
          88             :     ! unnecessary complication here: all should
          89             :     ! be on the same footing. RDL
          90             :     ! =====================================
          91        1536 :     rot_type = "contravariant"
          92             :
          93             :     ! ===============================================================
          94             :     ! Allocate and initialize the graph (array of GridVertex_t types)
          95             :     ! ===============================================================
          96             :
          97        1536 :     if (topology=="cube") then
          98             :
          99        1536 :       if (par%masterproc) then
         100           2 :         write(iulog,*) subname, "creating cube topology..."
         101           2 :         call shr_sys_flush(iulog)
         102             :       end if
         103             :
         104        1536 :       if (MeshUseMeshFile) then
         105           0 :         nelem = MeshCubeElemCount()
         106           0 :         nelem_edge = MeshCubeEdgeCount()
         107             :       else
         108        1536 :         nelem      = CubeElemCount()
         109        1536 :         nelem_edge = CubeEdgeCount()
         110             :       end if
         111             :
         112     8299008 :       allocate(GridVertex(nelem))
         113    66322944 :       allocate(GridEdge(nelem_edge))
         114             :
         115     8295936 :       do j = 1, nelem
         116     8295936 :         call allocate_gridvertex_nbrs(GridVertex(j))
         117             :       end do
         118             :
         119        1536 :       if (MeshUseMeshFile) then
         120           0 :         if (par%masterproc) then
         121           0 :           write(iulog,*) subname, "Set up grid vertex from mesh..."
         122             :         end if
         123           0 :         call MeshCubeTopology(GridEdge, GridVertex)
         124             :       else
         125        1536 :         call CubeTopology(GridEdge,GridVertex)
         126             :       end if
         127             :
         128        1536 :       if (par%masterproc) then
         129           2 :         write(iulog,*)"...done."
         130             :       end if
         131             :     end if
         132        1536 :     if(par%masterproc) then
         133           2 :       write(iulog,*) subname, "total number of elements nelem = ",nelem
         134             :     end if
         135             :
         136        1536 :     if(partmethod == SFCURVE) then
         137        1536 :       if(par%masterproc) then
         138           2 :         write(iulog,*) subname, "partitioning graph using SF Curve..."
         139             :       end if
         140        1536 :       call genspacepart(GridVertex)
         141             :     else
         142           0 :       write(errmsg, *) 'Unsupported partition method, ',partmethod
         143           0 :       call endrun(subname//trim(errmsg))
         144             :     end if
         145             :
         146             :     ! ===========================================================
         147             :     ! given partition, count number of local element descriptors
         148             :     ! ===========================================================
         149        1536 :     allocate(MetaVertex(1))
         150        1536 :     allocate(Schedule(1))
         151             :
         152        1536 :     nelem_edge = SIZE(GridEdge)
         153             :
         154             :     ! ====================================================
         155             :     !  Generate the communication graph
         156             :     ! ====================================================
         157        1536 :     call initMetaGraph(par%rank+1,MetaVertex(1),GridVertex,GridEdge)
         158             :
         159        1536 :     nelemd = LocalElemCount(MetaVertex(1))
         160             :     if (par%masterproc .and. Debug) then
         161             :       call PrintMetaVertex(MetaVertex(1))
         162             :     endif
         163             :
         164        1536 :     if(nelemd <= 0) then
         165           0 :       call endrun(subname//'Not yet ready to handle nelemd = 0 yet' )
         166             :     end if
         167        1536 :     call mpi_allreduce(nelemd, nelemdmax, 1, MPI_INTEGER, MPI_MAX, par%comm, ierr)
         168             :
         169             :     !Allocate elements:
         170        1536 :     if (nelemd > 0) then
         171       15408 :        allocate(elem(nelemd))
         172        1536 :        call allocate_element_desc(elem)
         173             :        !Allocate Qdp and derived FQ arrays:
         174        1536 :        if(fv_nphys > 0) then !SE-CSLAM
         175       12336 :           do ie=1,nelemd
         176       32400 :              allocate(elem(ie)%state%Qdp(np,np,nlev,thermodynamic_active_species_num,1), stat=ierr)
         177       10800 :              if( ierr /= 0 ) then
         178           0 :                 call endrun('prim_init1: failed to allocate Qdp array')
         179             :              end if
         180       32400 :              allocate(elem(ie)%derived%FQ(np,np,nlev,thermodynamic_active_species_num), stat=ierr)
         181       12336 :              if( ierr /= 0 ) then
         182           0 :                 call endrun('prim_init1: failed to allocate fq array')
         183             :              end if
         184             :           end do
         185             :        else !Regular SE
         186           0 :           do ie=1,nelemd
         187           0 :              allocate(elem(ie)%state%Qdp(np,np,nlev,qsize_d,2), stat=ierr)
         188           0 :              if( ierr /= 0 ) then
         189           0 :                 call endrun('prim_init1: failed to allocate Qdp array')
         190             :              end if
         191           0 :              allocate(elem(ie)%derived%FQ(np,np,nlev,qsize_d), stat=ierr)
         192           0 :              if( ierr /= 0 ) then
         193           0 :                 call endrun('prim_init1: failed to allocate fq array')
         194             :              end if
         195             :           end do
         196             :        end if
         197             :        !Allocate remaining derived quantity arrays:
         198       12336 :        do ie=1,nelemd
         199       10800 :           allocate(elem(ie)%derived%FDP(np,np,nlev), stat=ierr)
         200       10800 :           if( ierr /= 0 ) then
         201           0 :              call endrun('prim_init1: failed to allocate fdp array')
         202             :           end if
         203       10800 :           allocate(elem(ie)%derived%divdp(np,np,nlev), stat=ierr)
         204       10800 :           if( ierr /= 0 ) then
         205           0 :              call endrun('prim_init1: failed to allocate divdp array')
         206             :           end if
         207       10800 :           allocate(elem(ie)%derived%divdp_proj(np,np,nlev), stat=ierr)
         208       12336 :           if( ierr /= 0 ) then
         209           0 :              call endrun('prim_init1: failed to allocate divdp_proj array')
         210             :           end if
         211             :        end do
         212             :     end if
         213             :
         214        1536 :     if (fv_nphys > 0) then
         215       15408 :       allocate(fvm(nelemd))
         216        1536 :       call allocate_physgrid_vars(fvm,par)
         217             :     else
         218             :       ! Even if fvm not needed, still desirable to allocate it as empty
         219             :       ! so it can be passed as a (size zero) array rather than pointer.
         220           0 :       allocate(fvm(0))
         221             :     end if
         222             :
         223             :     ! ====================================================
         224             :     !  Generate the communication schedule
         225             :     ! ====================================================
         226             :
         227        1536 :     call genEdgeSched(par, elem, par%rank+1, Schedule(1), MetaVertex(1))
         228             :
         229        4608 :     allocate(global_shared_buf(nelemd, nrepro_vars))
         230    12621264 :     global_shared_buf = 0.0_r8
         231             :
         232        1536 :     call syncmp(par)
         233             :
         234             :     ! =================================================================
         235             :     ! Set number of domains (for 'decompose') equal to number of threads
         236             :     !  for OpenMP across elements, equal to 1 for OpenMP within element
         237             :     ! =================================================================
         238             :
         239             :     ! =================================================================
         240             :     ! Initialize shared boundary_exchange and reduction buffers
         241             :     ! =================================================================
         242        1536 :     if(par%masterproc) then
         243           2 :       write(iulog,*) subname, 'init shared boundary_exchange buffers'
         244           2 :       call shr_sys_flush(iulog)
         245             :     end if
         246        1536 :     call InitReductionBuffer(red,3*nlev,max_num_threads)
         247        1536 :     call InitReductionBuffer(red_sum,5)
         248        1536 :     call InitReductionBuffer(red_sum_int,1)
         249        1536 :     call InitReductionBuffer(red_max,1)
         250        1536 :     call InitReductionBuffer(red_max_int,1)
         251        1536 :     call InitReductionBuffer(red_min,1)
         252        1536 :     call initReductionBuffer(red_flops,1)
         253             :
         254        1536 :     gp = gausslobatto(np)  ! GLL points
         255             :
         256             :     ! fvm nodes are equally spaced in alpha/beta
         257             :     ! HOMME with equ-angular gnomonic projection maps alpha/beta space
         258             :     ! to the reference element via simple scale + translation
         259             :     ! thus, fvm nodes in reference element [-1,1] are a tensor product of
         260             :     ! array 'fvm_corners(:)' computed below:
         261        1536 :     xtmp = nc
         262        7680 :     do i = 1, nc+1
         263        7680 :       fvm_corners(i)= 2*(i-1)/xtmp - 1  ! [-1,1] including end points
         264             :     end do
         265        6144 :     do i = 1, nc
         266        6144 :       fvm_points(i)= ( fvm_corners(i)+fvm_corners(i+1) ) /2
         267             :     end do
         268             :
         269        1536 :     if (topology == "cube") then
         270        1536 :       if(par%masterproc) then
         271           2 :         write(iulog,*) subname, "initializing cube elements..."
         272           2 :         call shr_sys_flush(iulog)
         273             :       end if
         274        1536 :       if (MeshUseMeshFile) then
         275           0 :         call MeshSetCoordinates(elem)
         276             :       else
         277       12336 :         do ie = 1, nelemd
         278       12336 :           call set_corner_coordinates(elem(ie))
         279             :         end do
         280        1536 :         call assign_node_numbers_to_elem(elem, GridVertex)
         281             :       end if
         282       12336 :       do ie = 1, nelemd
         283       12336 :         call cube_init_atomic(elem(ie),gp%points)
         284             :       end do
         285             :     end if
         286             :
         287             :     ! =================================================================
         288             :     ! Initialize mass_matrix
         289             :     ! =================================================================
         290        1536 :     if(par%masterproc) then
         291           2 :       write(iulog,*) subname, 'running mass_matrix'
         292           2 :       call shr_sys_flush(iulog)
         293             :     end if
         294        1536 :     call mass_matrix(par, elem)
         295        4608 :     allocate(aratio(nelemd,1))
         296             :
         297        1536 :     if (topology == "cube") then
         298        1536 :       area = 0
         299       12336 :       do ie = 1, nelemd
         300      228336 :         aratio(ie,1) = sum(elem(ie)%mp(:,:)*elem(ie)%metdet(:,:))
         301             :       end do
         302        1536 :       call repro_sum(aratio, area, nelemd, nelemd, 1, commid=par%comm)
         303        1536 :       area(1) = 4.0_r8*pi/area(1)  ! ratio correction
         304        1536 :       deallocate(aratio)
         305        1536 :       if (par%masterproc) then
         306           2 :         write(iulog,'(2a,f20.17)') subname, "re-initializing cube elements: area correction=", area(1)
         307           2 :         call shr_sys_flush(iulog)
         308             :       end if
         309             :
         310       12336 :       do ie = 1, nelemd
         311       10800 :         call cube_init_atomic(elem(ie),gp%points,area(1))
         312       12336 :         call rotation_init_atomic(elem(ie),rot_type)
         313             :       end do
         314             :     end if
         315             :
         316        1536 :     if(par%masterproc) then
         317           2 :       write(iulog,*) subname, 're-running mass_matrix'
         318           2 :       call shr_sys_flush(iulog)
         319             :     end if
         320        1536 :     call mass_matrix(par, elem)
         321             :
         322             :     ! =================================================================
         323             :     ! Determine the global degree of freedome for each gridpoint
         324             :     ! =================================================================
         325        1536 :     if(par%masterproc) then
         326           2 :       write(iulog,*) subname, 'running global_dof'
         327           2 :       call shr_sys_flush(iulog)
         328             :     end if
         329        1536 :     call global_dof(par, elem)
         330             :
         331             :     ! =================================================================
         332             :     ! Create Unique Indices
         333             :     ! =================================================================
         334             :
         335       12336 :     do ie = 1, nelemd
         336       12336 :       call CreateUniqueIndex(elem(ie)%GlobalId,elem(ie)%gdofP,elem(ie)%idxP)
         337             :     end do
         338             :
         339        1536 :     call SetElemOffset(par,elem, GlobalUniqueCols)
         340             :
         341       12336 :     do ie = 1, nelemd
         342       12336 :       elem(ie)%idxV=>elem(ie)%idxP
         343             :     end do
         344             :
         345             :     ! initialize flux terms to 0
         346       12336 :     do ie = 1, nelemd
         347    43200000 :       elem(ie)%derived%FM=0.0_r8
         348   126630000 :       elem(ie)%derived%FQ=0.0_r8
         349    21103200 :       elem(ie)%derived%FT=0.0_r8
         350    21103200 :       elem(ie)%derived%FDP=0.0_r8
         351    21103200 :       elem(ie)%derived%pecnd=0.0_r8
         352             :
         353    21103200 :       elem(ie)%derived%Omega=0
         354    63320400 :       elem(ie)%state%dp3d=0
         355             :
         356       10800 :       elem(ie)%derived%etadot_prescribed = nan
         357       10800 :       elem(ie)%derived%u_met = nan
         358       10800 :       elem(ie)%derived%v_met = nan
         359       10800 :       elem(ie)%derived%dudt_met = nan
         360       10800 :       elem(ie)%derived%dvdt_met = nan
         361       10800 :       elem(ie)%derived%T_met = nan
         362       10800 :       elem(ie)%derived%dTdt_met = nan
         363       10800 :       elem(ie)%derived%ps_met = nan
         364       10800 :       elem(ie)%derived%dpsdt_met = nan
         365       10800 :       elem(ie)%derived%nudge_factor = nan
         366             :
         367    17085600 :       elem(ie)%derived%Utnd=0._r8
         368    17085600 :       elem(ie)%derived%Vtnd=0._r8
         369    17087136 :       elem(ie)%derived%Ttnd=0._r8
         370             :     end do
         371             :
         372             :     ! ==========================================================
         373             :     !  This routines initalizes a Restart file.  This involves:
         374             :     !      I)  Setting up the MPI datastructures
         375             :     ! ==========================================================
         376        1536 :     deallocate(GridEdge)
         377     8295936 :     do j = 1, nelem
         378     8295936 :       call deallocate_gridvertex_nbrs(GridVertex(j))
         379             :     end do
         380        1536 :     deallocate(GridVertex)
         381             :
         382       12336 :     do j = 1, MetaVertex(1)%nmembers
         383       12336 :       call deallocate_gridvertex_nbrs(MetaVertex(1)%members(j))
         384             :     end do
         385        1536 :     deallocate(MetaVertex)
         386             :
         387             :     ! =====================================
         388             :     ! Set number of threads...
         389             :     ! =====================================
         390        1536 :     if(par%masterproc) then
         391           2 :       write(iulog,*) subname, "max_num_threads=",max_num_threads
         392           2 :       call shr_sys_flush(iulog)
         393             :     end if
         394             :
         395        1536 :     nets = 1
         396        1536 :     nete = nelemd
         397        1536 :     call Prim_Advec_Init1(par, elem)
         398        1536 :     if (fv_nphys > 0) then
         399        1536 :       call fvm_init1(par,elem)
         400             :     end if
         401             :
         402             :     ! =======================================================
         403             :     ! Allocate memory for subcell flux calculations.
         404             :     ! =======================================================
         405        1536 :     call allocate_subcell_integration_matrix_cslam(np, nc)
         406        1536 :     if (fv_nphys > 0) then
         407        1536 :       call allocate_subcell_integration_matrix_physgrid(np, fv_nphys)
         408             :     end if
         409             :
         410        1536 :     call TimeLevel_init(tl)
         411             :
         412        1536 :     if (fv_nphys > 0) then
         413        1536 :       if(par%masterproc) then
         414           2 :         write(iulog,*) subname, 'initialize basic fvm coordinate variables'
         415           2 :         call shr_sys_flush(iulog)
         416             :       end if
         417       12336 :       do ie = 1, nelemd
         418           0 :         call compute_basic_coordinate_vars(elem(ie), nc, irecons_tracer,      &
         419           0 :              fvm(ie)%dalpha, fvm(ie)%dbeta, fvm(ie)%vtx_cart(:,:,1:nc,1:nc),  &
         420           0 :              fvm(ie)%center_cart(1:nc,1:nc), fvm(ie)%area_sphere(1:nc,1:nc),  &
         421       10800 :              fvm(ie)%spherecentroid(:,1:nc,1:nc))
         422           0 :         call compute_basic_coordinate_vars(elem(ie), fv_nphys, irecons_tracer,&
         423           0 :              fvm(ie)%dalpha_physgrid, fvm(ie)%dbeta_physgrid,                 &
         424       21600 :              fvm(ie)%vtx_cart_physgrid   (:,:,1:fv_nphys,1:fv_nphys),         &
         425       21600 :              fvm(ie)%center_cart_physgrid(1:fv_nphys,1:fv_nphys),             &
         426       21600 :              fvm(ie)%area_sphere_physgrid(1:fv_nphys,1:fv_nphys),             &
         427       77136 :              fvm(ie)%spherecentroid_physgrid(:,1:fv_nphys,1:fv_nphys))
         428             :       end do
         429             :     end if
         430             :
         431        1536 :     if(par%masterproc) then
         432           2 :       write(iulog,*) subname, 'end of prim_init'
         433           2 :       call shr_sys_flush(iulog)
         434             :     end if
         435        1536 :   end subroutine prim_init1
         436             : end module prim_init

Generated by: LCOV version 1.14