LCOV - code coverage report
Current view: | top level - dynamics/se - dp_mapping.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 88 | 342 | 25.7 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 2 | 6 | 33.3 %

              Line data    Source code
           1             : ! Separate dynamics and physics grids
           2             :
           3             : module dp_mapping
           4             :   use dimensions_mod,         only: np, npsq, fv_nphys
           5             :   use shr_kind_mod,           only: r8=>shr_kind_r8, shr_kind_cl
           6             :   use coordinate_systems_mod, only: spherical_polar_t
           7             :   use shr_const_mod,          only: pi => shr_const_pi
           8             :   use fvm_control_volume_mod, only: fvm_struct
           9             :
          10             :   implicit none
          11             :   private
          12             :   save
          13             :
          14             :   public :: dp_init
          15             :   public :: dp_reorder
          16             :   public :: dp_write
          17             :   public :: dp_allocate
          18             :   public :: dp_deallocate
          19             :
          20             :   ! Total number of physics points per spectral element
          21             :   ! no physgrid:   nphys_pts = npsq   (physics on GLL grid)
          22             :   !    physgrid:   nphys_pts = nphys2 (physics on CSLAM grid)
          23             :   ! Value is set when se_fv_nphys namelist variable is read
          24             :   integer,            public :: nphys_pts = npsq
          25             :
          26             :   ! NOTE:  dp_gid() is in space filling curve rank order
          27             :   !        all other global arrays are in block id (global id) order
          28             :   !
          29             :   !        dp_gid() is used to re-order data collected on root via mpi_gatherv
          30             :   !          into block id ordering
          31             :   !
          32             :   ! j=dp_gid(i)   i = element space filling curve rank
          33             :   !               j = element global id = block id = history file ordering
          34             :   !
          35             :   integer,        allocatable,dimension(:)         :: dp_gid  ! NE=240, integer*4 = 1.3MB
          36             :   integer, public,allocatable,dimension(:)         :: dp_owner
          37             :
          38             :   real (r8),public,allocatable,dimension(:,:,:)  :: weights_all_fvm2phys
          39             :   integer  ,public,allocatable,dimension(:,:,:)  :: weights_eul_index_all_fvm2phys,weights_lgr_index_all_fvm2phys
          40             :   real (r8),public,allocatable,dimension(:,:,:)  :: weights_all_phys2fvm
          41             :   integer  ,public,allocatable,dimension(:,:,:)  :: weights_eul_index_all_phys2fvm,weights_lgr_index_all_phys2fvm
          42             :   integer  ,public,allocatable,dimension(:)      :: jall_fvm2phys,jall_phys2fvm
          43             :   integer  ,public                               :: num_weights_fvm2phys,num_weights_phys2fvm
          44             :
          45             :
          46             : contains
          47        1536 :   subroutine dp_init(elem,fvm)
          48             :     use dimensions_mod, only: nelemd, nc, irecons_tracer
          49             :     use element_mod,    only: element_t
          50             :     use spmd_utils,     only: masterproc
          51             :     use cam_logfile,    only: iulog
          52             :
          53             :     type(element_t)  , dimension(nelemd), intent(in) :: elem
          54             :     type (fvm_struct), dimension(nelemd), intent(in) :: fvm
          55             :
          56        1536 :     num_weights_phys2fvm = 0
          57        1536 :     num_weights_fvm2phys = 0
          58        1536 :     if (fv_nphys>0) then
          59        1536 :       num_weights_phys2fvm = (nc+fv_nphys)**2
          60        1536 :       num_weights_fvm2phys = (nc+fv_nphys)**2
          61             :
          62        6144 :        allocate(weights_all_fvm2phys(num_weights_fvm2phys,irecons_tracer,nelemd))
          63        6144 :        allocate(weights_eul_index_all_fvm2phys(num_weights_fvm2phys,2,nelemd))
          64        4608 :        allocate(weights_lgr_index_all_fvm2phys(num_weights_fvm2phys,2,nelemd))
          65             :
          66        4608 :        allocate(weights_all_phys2fvm(num_weights_phys2fvm,irecons_tracer,nelemd))
          67        4608 :        allocate(weights_eul_index_all_phys2fvm(num_weights_phys2fvm,2,nelemd))
          68        4608 :        allocate(weights_lgr_index_all_phys2fvm(num_weights_phys2fvm,2,nelemd))
          69        4608 :        allocate(jall_fvm2phys(nelemd))
          70        3072 :        allocate(jall_phys2fvm(nelemd))
          71             :
          72             :        call fvm2phys_init(elem,fvm,nc,fv_nphys,irecons_tracer,&
          73             :             weights_all_fvm2phys,weights_eul_index_all_fvm2phys,weights_lgr_index_all_fvm2phys,&
          74             :             weights_all_phys2fvm,weights_eul_index_all_phys2fvm,weights_lgr_index_all_phys2fvm,&
          75        1536 :             jall_fvm2phys,jall_phys2fvm)
          76             :
          77        1536 :       if (masterproc) then
          78           2 :         write(iulog, *) 'dp_init: Initialized phys2fvm/fvm2phys mapping vars'
          79             :       end if
          80             :
          81             :     end if
          82        1536 :   end subroutine dp_init
          83             :
          84           0 :   subroutine dp_reorder(before, after)
          85             :     use cam_abortutils, only: endrun
          86             :     use dimensions_mod, only: nelem
          87             :     use cam_logfile,    only: iulog
          88             :     use spmd_utils,     only: masterproc
          89             :     use shr_sys_mod,    only: shr_sys_flush
          90             :
          91             :     real(r8), dimension(fv_nphys*fv_nphys,*), intent(in)  :: before
          92             :     real(r8), dimension(fv_nphys*fv_nphys,*), intent(out) :: after
          93             :     integer :: ie
          94             :
          95             :     ! begin
          96           0 :     do ie = 1,nelem
          97           0 :        if (dp_gid(ie) < 0) then
          98           0 :           if (masterproc) then
          99           0 :              write(iulog,*) 'ie =',ie,', dp_gid(ie) =',dp_gid(ie)
         100           0 :              call shr_sys_flush(iulog)
         101             :           end if
         102           0 :           call endrun('Bad element remap in dp_reorder')
         103             :        end if
         104           0 :        after(:,dp_gid(ie)) = before(:,ie)
         105             :     end do
         106           0 :   end subroutine dp_reorder
         107             :
         108             :   !!!
         109             :
         110           0 :   subroutine dp_allocate(elem)
         111             :     use dimensions_mod, only: nelem, nelemd
         112             :     use element_mod,    only: element_t
         113             :     use spmd_utils,     only: masterproc, masterprocid, npes
         114             :     use spmd_utils,     only: mpicom, mpi_integer
         115             :
         116             :     implicit none
         117             :     type(element_t),dimension(nelemd),intent(in) :: elem
         118             :
         119             :     integer                          :: i,j,ierror
         120           0 :     integer,dimension(nelemd)        :: lgid
         121           0 :     integer,dimension(:),allocatable :: displs,recvcount
         122             :
         123             :     ! begin
         124             :
         125           0 :     allocate(displs(npes))
         126           0 :     allocate(dp_gid(nelem))
         127           0 :     allocate(recvcount(npes))
         128             :     call mpi_gather(nelemd, 1, mpi_integer, recvcount, 1, mpi_integer,        &
         129           0 :          masterprocid, mpicom, ierror)
         130           0 :     lgid(:) = elem(:)%globalid
         131           0 :     if (masterproc) then
         132           0 :       displs(1) = 0
         133           0 :       do i = 2,npes
         134           0 :         displs(i) = displs(i-1)+recvcount(i-1)
         135             :       end do
         136             :     end if
         137             :     call mpi_gatherv(lgid, nelemd, mpi_integer, dp_gid, recvcount, displs,    &
         138           0 :          mpi_integer, masterprocid, mpicom, ierror)
         139           0 :     if (masterproc) then
         140           0 :       allocate(dp_owner(nelem))
         141           0 :       dp_owner(:) = -1
         142           0 :       do i = 1,npes
         143           0 :         do j = displs(i)+1,displs(i)+recvcount(i)
         144           0 :           dp_owner(dp_gid(j)) = i-1
         145             :         end do
         146             :       end do
         147             :     end if
         148           0 :     deallocate(displs)
         149           0 :     deallocate(recvcount)
         150             :     ! minimize global memory use
         151           0 :     call mpi_barrier(mpicom,ierror)
         152           0 :     if (.not.masterproc) then
         153           0 :       allocate(dp_owner(nelem))
         154             :     end if
         155           0 :     call mpi_bcast(dp_gid,nelem,mpi_integer,masterprocid,mpicom,ierror)
         156           0 :     call mpi_bcast(dp_owner,nelem,mpi_integer,masterprocid,mpicom,ierror)
         157           0 :  end subroutine dp_allocate
         158             :
         159             :   !!!
         160             :
         161           0 :   subroutine dp_deallocate()
         162           0 :      deallocate(dp_gid)
         163           0 :      deallocate(dp_owner)
         164           0 :   end subroutine dp_deallocate
         165             :
         166             :   !!!
         167             :
         168           0 :   subroutine dp_write(elem, fvm, grid_format, filename_in)
         169             :     use cam_abortutils,         only: endrun
         170             :     use dimensions_mod,         only: nelem, nelemd
         171             :     use element_mod,            only: element_t
         172             :     use netcdf,                 only: nf90_create, nf90_close, nf90_enddef
         173             :     use netcdf,                 only: nf90_def_dim, nf90_def_var, nf90_put_var
         174             :     use netcdf,                 only: nf90_double, nf90_int, nf90_put_att
         175             :     use netcdf,                 only: nf90_noerr, nf90_strerror, nf90_clobber
         176             :     use spmd_utils,             only: masterproc, masterprocid, mpicom, npes
         177             :     use spmd_utils,             only: mpi_integer, mpi_real8
         178             :     use cam_logfile,            only: iulog
         179             :     use shr_sys_mod,            only: shr_sys_flush
         180             :     use dimensions_mod,         only: ne
         181             :     use coordinate_systems_mod, only: cart2spherical
         182             :
         183             :     ! Inputs
         184             :     type(element_t),   intent(in) :: elem(:)
         185             :     type (fvm_struct), intent(in) :: fvm(:)
         186             :     character(len=*),  intent(in) :: grid_format
         187             :     character(len=*),  intent(in) :: filename_in
         188             :
         189             :     real(r8), parameter :: rad2deg = 180._r8/pi
         190             :
         191             :     ! Local variables
         192             :     integer           :: i, ie, ierror, j, status, ivtx
         193             :     integer           :: grid_corners_id, grid_rank_id, grid_size_id
         194             :     integer           :: ncid
         195             :     integer           :: grid_dims_id, grid_area_id, grid_center_lat_id
         196             :     integer           :: grid_center_lon_id, grid_corner_lat_id
         197             :     integer           :: grid_corner_lon_id, grid_imask_id
         198             :     integer           :: gridsize
         199             :     integer           :: IOrootID
         200             :     logical           :: IOroot
         201             :     character(len=SHR_KIND_CL) :: errormsg
         202             :     character(len=SHR_KIND_CL) :: filename
         203           0 :     integer,allocatable,dimension(:) :: displs,recvcount
         204             :
         205           0 :     real(r8), dimension(fv_nphys, fv_nphys, nelemd, 4, 2)  :: corners
         206           0 :     real(r8), dimension(fv_nphys, fv_nphys, nelemd)        :: lwork
         207           0 :     real(r8), allocatable, dimension(:)                    :: recvbuf
         208           0 :     real(r8), allocatable, dimension(:,:)                  :: gwork
         209             :     real(r8)                                               :: x, y
         210             :     type (spherical_polar_t)                               :: sphere
         211             :
         212             :     ! begin
         213             :
         214             :     !! Check to see if we are doing grid output
         215           0 :     if (trim(grid_format) == "no") then
         216           0 :       if (masterproc) then
         217           0 :         write(iulog, *) 'dp_write: Not writing phys_grid file.'
         218             :       end if
         219           0 :       return
         220           0 :     else if (trim(grid_format) /= 'SCRIP') then
         221           0 :       if (masterproc) then
         222           0 :         write(errormsg, *) 'dp_write: ERROR, bad value for se_write_grid, ',&
         223           0 :              trim(grid_format)
         224           0 :         call endrun(errormsg)
         225             :       end if
         226             :     end if
         227             :
         228             :     ! Create the NetCDF file
         229           0 :     if (len_trim(filename_in) == 0) then
         230           0 :       write(filename, '(3(a,i0),3a)') "ne", ne, "np", np, ".pg", fv_nphys,    &
         231           0 :            "_", trim(grid_format), ".nc"
         232             :     else
         233           0 :       filename = trim(filename_in)
         234             :     end if
         235           0 :     status = nf90_create(trim(filename), nf90_clobber, ncid)
         236           0 :     if (status /= nf90_noerr) then
         237           0 :       call endrun("dp_write: "//trim(nf90_strerror(status)))
         238             :     end if
         239             :
         240             :     ! PIO_put_var puts from its root node, find that (so we do our work there)
         241           0 :     IOrootID = masterprocid
         242           0 :     IOroot = masterproc
         243             :
         244             :     ! Allocate workspace and calculate PE displacement information
         245           0 :     if (IOroot) then
         246           0 :       allocate(displs(npes))
         247           0 :       allocate(recvcount(npes))
         248             :     else
         249           0 :       allocate(displs(0))
         250           0 :       allocate(recvcount(0))
         251             :     end if
         252           0 :     gridsize = nelem * fv_nphys*fv_nphys
         253           0 :     if(masterproc) then
         254           0 :       write(iulog, *) 'Writing physics SCRIP grid file: ', trim(filename)
         255           0 :       write(iulog, '(a,i7,a,i3)') 'nelem = ', nelem, ', fv_nphys = ', fv_nphys
         256           0 :       call shr_sys_flush(iulog)
         257             :     end if
         258             :     call mpi_gather(nelemd*fv_nphys*fv_nphys, 1, mpi_integer, recvcount, 1,           &
         259           0 :          mpi_integer, IOrootID, mpicom, ierror)
         260             :
         261           0 :     if (IOroot) then
         262           0 :       displs(1) = 0
         263           0 :       do i = 2, npes
         264           0 :         displs(i) = displs(i-1)+recvcount(i-1)
         265             :       end do
         266           0 :       allocate(recvbuf(gridsize))
         267             :     else
         268           0 :       allocate(recvbuf(0))
         269             :     end if
         270           0 :     allocate(gwork(4, gridsize))
         271             :
         272           0 :     if (IOroot) then
         273             :       ! Define the horizontal grid dimensions for SCRIP output
         274           0 :       status = nf90_def_dim(ncid, "grid_corners", 4,        grid_corners_id)
         275           0 :       if (status /= nf90_noerr) then
         276           0 :         write(iulog, *) 'dp_write: Error defining dimension, grid_corners'
         277           0 :         call shr_sys_flush(iulog)
         278           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         279             :       end if
         280           0 :       status = nf90_def_dim(ncid, "grid_rank",    1,        grid_rank_id)
         281           0 :       if (status /= nf90_noerr) then
         282           0 :         write(iulog, *) 'dp_write: Error defining dimension, grid_rank'
         283           0 :         call shr_sys_flush(iulog)
         284           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         285             :       end if
         286           0 :       status = nf90_def_dim(ncid, "grid_size",    gridsize, grid_size_id)
         287           0 :       if (status /= nf90_noerr) then
         288           0 :         write(iulog, *) 'dp_write: Error defining dimension, grid_size'
         289           0 :         call shr_sys_flush(iulog)
         290           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         291             :       end if
         292             :
         293             :       ! Define the coordinate variables
         294             :       status = nf90_def_var(ncid, "grid_dims", nf90_int, (/grid_rank_id/),  &
         295           0 :            grid_dims_id)
         296           0 :       if (status /= nf90_noerr) then
         297           0 :         write(iulog, *) 'dp_write: Error defining variable grid_dims'
         298           0 :         call shr_sys_flush(iulog)
         299           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         300             :       end if
         301             :
         302             :       status = nf90_def_var(ncid, "grid_area", nf90_double,                 &
         303           0 :            (/grid_size_id/), grid_area_id)
         304           0 :       if (status /= nf90_noerr) then
         305           0 :         write(iulog, *) 'dp_write: Error defining variable grid_area'
         306           0 :         call shr_sys_flush(iulog)
         307           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         308             :       end if
         309             :
         310           0 :       status = nf90_put_att(ncid, grid_area_id, "units", "radians^2")
         311           0 :       if (status /= nf90_noerr) then
         312           0 :         write(iulog, *) 'dp_write: Error defining attributes for grid_area'
         313           0 :         call shr_sys_flush(iulog)
         314           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         315             :       end if
         316             :
         317           0 :       status = nf90_put_att(ncid, grid_area_id, "long_name", "area weights")
         318           0 :       if (status /= nf90_noerr) then
         319           0 :         write(iulog, *) 'dp_write: Error defining attributes for grid_area'
         320           0 :         call shr_sys_flush(iulog)
         321           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         322             :       end if
         323             :
         324             :       status = nf90_def_var(ncid, "grid_center_lat", nf90_double,           &
         325           0 :            (/grid_size_id/), grid_center_lat_id)
         326           0 :       if (status /= nf90_noerr) then
         327           0 :         write(iulog, *) 'dp_write: Error defining variable grid_center_lat'
         328           0 :         call shr_sys_flush(iulog)
         329           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         330             :       end if
         331             :
         332           0 :       status = nf90_put_att(ncid, grid_center_lat_id, "units", "degrees")
         333           0 :       if (status /= nf90_noerr) then
         334           0 :         write(iulog, *) 'dp_write: Error defining attributes for grid_center_lat'
         335           0 :         call shr_sys_flush(iulog)
         336           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         337             :       end if
         338             :
         339             :       status = nf90_def_var(ncid, "grid_center_lon", nf90_double,           &
         340           0 :            (/grid_size_id/), grid_center_lon_id)
         341           0 :       if (status /= nf90_noerr) then
         342           0 :         write(iulog, *) 'dp_write: Error defining variable grid_center_lon'
         343           0 :         call shr_sys_flush(iulog)
         344           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         345             :       end if
         346             :
         347           0 :       status = nf90_put_att(ncid, grid_center_lon_id, "units", "degrees")
         348           0 :       if (status /= nf90_noerr) then
         349           0 :         write(iulog, *) 'dp_write: Error defining attributes for grid_center_lon'
         350           0 :         call shr_sys_flush(iulog)
         351           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         352             :       end if
         353             :
         354             :       status = nf90_def_var(ncid, "grid_corner_lat", nf90_double,           &
         355           0 :            (/grid_corners_id, grid_size_id/), grid_corner_lat_id)
         356           0 :       if (status /= nf90_noerr) then
         357           0 :         write(iulog, *) 'dp_write: Error defining grid_corner_lat'
         358           0 :         call shr_sys_flush(iulog)
         359           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         360             :       end if
         361             :
         362           0 :       status = nf90_put_att(ncid, grid_corner_lat_id, "units", "degrees")
         363           0 :       if (status /= nf90_noerr) then
         364           0 :         write(iulog, *) 'dp_write: Error defining attributes for grid_corner_lat'
         365           0 :         call shr_sys_flush(iulog)
         366           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         367             :       end if
         368             :
         369             :       status = nf90_def_var(ncid, "grid_corner_lon", nf90_double,           &
         370           0 :            (/grid_corners_id, grid_size_id/), grid_corner_lon_id)
         371           0 :       if (status /= nf90_noerr) then
         372           0 :         write(iulog, *) 'dp_write: Error defining variable grid_corner_lon'
         373           0 :         call shr_sys_flush(iulog)
         374           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         375             :       end if
         376             :
         377           0 :       status = nf90_put_att(ncid, grid_corner_lon_id, "units", "degrees")
         378           0 :       if (status /= nf90_noerr) then
         379           0 :         write(iulog, *) 'dp_write: Error defining attributes for grid_corner_lon'
         380           0 :         call shr_sys_flush(iulog)
         381           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         382             :       end if
         383             :
         384             :       status = nf90_def_var(ncid, "grid_imask", nf90_double,                &
         385           0 :            (/grid_size_id/), grid_imask_id)
         386           0 :       if (status /= nf90_noerr) then
         387           0 :         write(iulog, *) 'dp_write: Error defining variable grid_imask'
         388           0 :         call shr_sys_flush(iulog)
         389           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         390             :       end if
         391             :
         392             :       ! End of NetCDF definitions
         393           0 :       status = nf90_enddef(ncid)
         394           0 :       if (status /= nf90_noerr) then
         395           0 :         write(iulog, *) 'dp_write: Error calling enddef'
         396           0 :         call shr_sys_flush(iulog)
         397           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         398             :       end if
         399             :     end if ! IOroot
         400             :
         401             :     if (IOroot) then
         402           0 :       status = nf90_put_var(ncid, grid_dims_id, (/ gridsize /))
         403           0 :       if (status /= nf90_noerr) then
         404           0 :         write(iulog, *) 'dp_write: Error writing variable grid_dims'
         405           0 :         call shr_sys_flush(iulog)
         406           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         407             :       end if
         408             :     end if
         409             :
         410           0 :     do ie=1,nelemd
         411           0 :       lwork(:,:,ie) = fvm(ie)%area_sphere_physgrid(:,:)
         412             :     end do
         413             :     call mpi_gatherv(lwork, size(lwork), mpi_real8, recvbuf, recvcount, &
         414           0 :          displs, mpi_real8, IOrootID, mpicom, ierror)
         415           0 :     call dp_allocate(elem)
         416           0 :     if (IOroot) then
         417           0 :       call dp_reorder(recvbuf, gwork(1,:))
         418           0 :       status = nf90_put_var(ncid, grid_area_id, gwork(1,:))
         419           0 :       if (status /= nf90_noerr) then
         420           0 :         write(iulog, *) 'dp_write: Error writing variable grid_area'
         421           0 :         call shr_sys_flush(iulog)
         422           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         423             :       end if
         424             :     end if
         425           0 :     do ie=1,nelemd
         426           0 :       lwork(:,:,ie) = rad2deg*fvm(ie)%center_cart_physgrid(:,:)%lat
         427             :     end do
         428             :     call mpi_gatherv(lwork, size(lwork), mpi_real8, recvbuf, recvcount,       &
         429           0 :          displs, mpi_real8, IOrootID, mpicom, ierror)
         430           0 :     if (IOroot) then
         431           0 :       call dp_reorder(recvbuf, gwork(1,:))
         432           0 :       status = nf90_put_var(ncid, grid_center_lat_id, gwork(1,:))
         433           0 :       if (status /= nf90_noerr) then
         434           0 :         write(iulog, *) 'dp_write: Error writing variable grid_center_lat'
         435           0 :         call shr_sys_flush(iulog)
         436           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         437             :       end if
         438             :     end if
         439             :
         440           0 :     do ie=1,nelemd
         441           0 :       lwork(:,:,ie) = rad2deg*fvm(ie)%center_cart_physgrid(:,:)%lon
         442             :     end do
         443             :     call mpi_gatherv(lwork, size(lwork), mpi_real8, recvbuf, recvcount,       &
         444           0 :          displs, mpi_real8, IOrootID, mpicom, ierror)
         445           0 :     if (IOroot) then
         446           0 :       call dp_reorder(recvbuf, gwork(1,:))
         447           0 :       status = nf90_put_var(ncid, grid_center_lon_id, gwork(1,:))
         448           0 :       if (status /= nf90_noerr) then
         449           0 :         write(iulog, *) 'dp_write: Error writing variable grid_center_lon'
         450           0 :         call shr_sys_flush(iulog)
         451           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         452             :       end if
         453             :     end if
         454             :     ! compute physgrid grid corners
         455           0 :     do ie=1,nelemd
         456           0 :       do j=1,fv_nphys
         457           0 :         do i=1,fv_nphys
         458           0 :           do ivtx=1,4
         459           0 :             x      = fvm(ie)%vtx_cart_physgrid(ivtx,1,i,j)
         460           0 :             y      = fvm(ie)%vtx_cart_physgrid(ivtx,2,i,j)
         461           0 :             sphere = cart2spherical(x,y,elem(ie)%FaceNum)
         462           0 :             corners(i,j,ie,ivtx,1) = rad2deg * sphere%lat
         463           0 :             corners(i,j,ie,ivtx,2) = rad2deg * sphere%lon
         464             :           end do
         465             :         end do
         466             :       end do
         467             :     end do
         468             :     ! Collect all information for the grid corner latitude (counter-clockwise)
         469           0 :     do ivtx=1,4
         470             :       call mpi_gatherv(corners(:,:,:,ivtx,1), size(corners(:,:,:,ivtx,1)), mpi_real8, recvbuf, recvcount,     &
         471           0 :            displs, mpi_real8, IOrootID, mpicom, ierror)
         472           0 :       if (IOroot) then
         473           0 :         call dp_reorder(recvbuf, gwork(ivtx,:))
         474             :       end if
         475             :     end do
         476           0 :     if (IOroot) then
         477           0 :       status = nf90_put_var(ncid, grid_corner_lat_id, gwork)
         478           0 :       if (status /= nf90_noerr) then
         479           0 :         write(iulog, *) 'dp_write: Error writing variable grid_corner_lat'
         480           0 :         call shr_sys_flush(iulog)
         481           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         482             :       end if
         483             :     end if
         484             :     ! Collect all information for the grid corner longitudes (counter-clockwise)
         485           0 :     do ivtx=1,4
         486             :       call mpi_gatherv(corners(:,:,:,ivtx,2), size(corners(:,:,:,ivtx,2)), mpi_real8, recvbuf, recvcount,     &
         487           0 :            displs, mpi_real8, IOrootID, mpicom, ierror)
         488           0 :       if (IOroot) then
         489           0 :         call dp_reorder(recvbuf, gwork(ivtx,:))
         490             :       end if
         491             :     end do
         492           0 :     call dp_deallocate()
         493           0 :     if (IOroot) then
         494           0 :       status = nf90_put_var(ncid, grid_corner_lon_id, gwork)
         495           0 :       if (status /= nf90_noerr) then
         496           0 :         write(iulog, *) 'dp_write: Error writing variable grid_corner_lon'
         497           0 :         call shr_sys_flush(iulog)
         498           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         499             :       end if
         500             :     end if
         501             :
         502             :     if (IOroot) then
         503           0 :       gwork(1,:) = 1._r8
         504           0 :       status = nf90_put_var(ncid, grid_imask_id, gwork(1,:))
         505           0 :       if (status /= nf90_noerr) then
         506           0 :         write(iulog, *) 'dp_write: Error writing variable grid_imask'
         507           0 :         call shr_sys_flush(iulog)
         508           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         509             :       end if
         510             :     end if
         511             :
         512             : !    call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
         513             :     ! Close the file
         514           0 :     call mpi_barrier(mpicom, ierror)
         515           0 :     if (IOroot) then
         516           0 :       status = nf90_close(ncid)
         517           0 :       if (status /= nf90_noerr) then
         518           0 :         call endrun("dp_write: "//trim(nf90_strerror(status)))
         519             :       end if
         520             :     end if
         521             :
         522           0 :     call mpi_barrier(mpicom, ierror)
         523           0 :     if(masterproc) then
         524           0 :       write(iulog, *) 'Finished writing physics grid file: ', trim(filename)
         525           0 :       call shr_sys_flush(iulog)
         526             :     end if
         527             :
         528           0 :   end subroutine dp_write
         529             :   !!!
         530             :
         531        1536 :   subroutine fvm2phys_init(elem,fvm,fvm_nc,phys_nc,irecons,&
         532        1536 :             weights_all_fvm2phys,weights_eul_index_all_fvm2phys,weights_lgr_index_all_fvm2phys,&
         533        1536 :             weights_all_phys2fvm,weights_eul_index_all_phys2fvm,weights_lgr_index_all_phys2fvm,&
         534        1536 :             jall_fvm2phys,jall_phys2fvm)
         535             :     use dimensions_mod  , only: ngpc,nelemd
         536             :     use fvm_overlap_mod , only: compute_weights_cell
         537             :     use element_mod     , only: element_t
         538             :     type(element_t)  , dimension(nelemd), intent(in) :: elem
         539             :     type (fvm_struct), dimension(nelemd), intent(in) :: fvm
         540             :     integer        , intent(in) :: fvm_nc, phys_nc, irecons
         541             :     real (kind=r8)       :: dalpha,dbeta
         542        3072 :     real (kind=r8), dimension(0:phys_nc+2):: xgno_phys,ygno_phys
         543        3072 :     real (kind=r8), dimension(0:fvm_nc+2) :: xgno_fvm,ygno_fvm
         544             :
         545             :     real (kind=r8), dimension(ngpc):: gauss_weights, abscissae !dimension(ngauss)
         546             :
         547             :     integer :: i,j,h
         548             :     integer, parameter :: nvertex = 4
         549             :     real (kind=r8), dimension(nvertex)            :: xcell,ycell
         550             :
         551             :     real (kind=r8)   , dimension(num_weights_fvm2phys,irecons,nelemd),intent(out) :: weights_all_fvm2phys
         552             :     integer,  dimension(num_weights_fvm2phys,2,nelemd),intent(out) :: weights_eul_index_all_fvm2phys
         553             :     integer,  dimension(num_weights_fvm2phys,2,nelemd),intent(out) :: weights_lgr_index_all_fvm2phys
         554             :
         555             :     real (kind=r8)   , dimension(num_weights_phys2fvm,irecons,nelemd),intent(out) :: weights_all_phys2fvm
         556             :     integer,  dimension(num_weights_phys2fvm,2,nelemd),intent(out) :: weights_eul_index_all_phys2fvm
         557             :     integer,  dimension(num_weights_phys2fvm,2,nelemd),intent(out) :: weights_lgr_index_all_phys2fvm
         558             :
         559             :     integer                , dimension(nelemd)                       ,intent(out) :: jall_fvm2phys,jall_phys2fvm
         560             :
         561             :     integer, parameter :: jmax_segments_cell = 50
         562        3072 :     real (kind=r8)   , dimension(jmax_segments_cell,irecons)  :: weights_cell
         563             :     integer , dimension(jmax_segments_cell,2)        :: weights_eul_index_cell
         564             :     integer :: jcollect_cell,ie
         565        3072 :     real(kind=r8), dimension(phys_nc,phys_nc) :: phys_area, factor
         566        3072 :     real(kind=r8), dimension(fvm_nc,fvm_nc)   :: fvm_area
         567             :
         568        1536 :     xgno_phys(0) = -1D20; xgno_phys(phys_nc+2) = 1D20
         569        1536 :     xgno_fvm(0)  = -1D20; xgno_fvm(fvm_nc+2)   = 1D20
         570       12336 :     do ie=1,nelemd
         571       10800 :       dalpha         = abs(elem(ie)%corners(1)%x-elem(ie)%corners(2)%x)/phys_nc !in alpha
         572       10800 :       dbeta          = abs(elem(ie)%corners(1)%y-elem(ie)%corners(4)%y)/phys_nc  !in beta
         573       54000 :       do i=1,phys_nc+1
         574       43200 :         xgno_phys(i) = tan(elem(ie)%corners(1)%x+(i-1)*dalpha)
         575       54000 :         ygno_phys(i) = tan(elem(ie)%corners(1)%y+(i-1)*dbeta )
         576             :       end do
         577             :
         578       10800 :       dalpha         = abs(elem(ie)%corners(1)%x-elem(ie)%corners(2)%x)/fvm_nc !in alpha
         579       10800 :       dbeta          = abs(elem(ie)%corners(1)%y-elem(ie)%corners(4)%y)/fvm_nc  !in beta
         580       54000 :       do i=1,fvm_nc+1
         581       43200 :         xgno_fvm(i) = tan(elem(ie)%corners(1)%x+(i-1)*dalpha)
         582       54000 :         ygno_fvm(i) = tan(elem(ie)%corners(1)%y+(i-1)*dbeta )
         583             :       end do
         584             :
         585             :       !
         586             :       ! compute area using line-integrals
         587             :       !
         588             :       !       do j=1,phys_nc
         589             :       !          do i=1,phys_nc
         590             :       !             da_phys(i,j) = (I_00(xgno_phys(i+1),ygno_phys(j+1)) - I_00(xgno_phys(i  ),ygno_phys(j+1)) + &
         591             :       !                  I_00(xgno_phys(i  ),ygno_phys(j  )) - I_00(xgno_phys(i+1),ygno_phys(j  )))
         592             :       !          end do
         593             :       !       end do
         594             :       !
         595             :       !       do j=1,fvm_nc
         596             :       !          do i=1,fvm_nc
         597             :       !             da_fvm(i,j) = (I_00(xgno_fvm(i+1),ygno_fvm(j+1)) - I_00(xgno_fvm(i  ),ygno_fvm(j+1)) + &
         598             :       !                  I_00(xgno_fvm(i  ),ygno_fvm(j  )) - I_00(xgno_fvm(i+1),ygno_fvm(j  )))
         599             :       !          end do
         600             :       !       end do
         601             :
         602       10800 :       gauss_weights = 0.0D0; abscissae=0.0D0!not used since line-segments are parallel to coordinate
         603             :
         604       10800 :       jall_fvm2phys(ie)=1
         605       43200 :       do j=1,phys_nc
         606      140400 :         do i=1,phys_nc
         607       97200 :           xcell(1) = xgno_phys(i)  ; ycell(1) = ygno_phys(j)
         608       97200 :           xcell(2) = xgno_phys(i)  ; ycell(2) = ygno_phys(j+1)
         609       97200 :           xcell(3) = xgno_phys(i+1); ycell(3) = ygno_phys(j+1)
         610       97200 :           xcell(4) = xgno_phys(i+1); ycell(4) = ygno_phys(j)
         611             :
         612             :           call compute_weights_cell(nvertex,.true.,&
         613             :                xcell,ycell,i,j,irecons,xgno_fvm,ygno_fvm,0,fvm_nc+2,&
         614             :                1,fvm_nc+1,1,fvm_nc+1,&
         615             :                ngpc,gauss_weights,abscissae,&
         616       97200 :                weights_cell,weights_eul_index_cell,jcollect_cell,jmax_segments_cell)
         617             :
         618      129600 :           if (jcollect_cell>0) then
         619           0 :             weights_all_fvm2phys(jall_fvm2phys(ie):jall_fvm2phys(ie)+jcollect_cell-1,:,ie) = &
         620     1263600 :                  weights_cell(1:jcollect_cell,:)!/fvm(ie)%area_sphere_physgrid(i,j)!da_phys(i,j)
         621             :
         622           0 :             weights_eul_index_all_fvm2phys(jall_fvm2phys(ie):jall_fvm2phys(ie)+jcollect_cell-1,:,ie) = &
         623      486000 :                  weights_eul_index_cell(1:jcollect_cell,:)
         624      194400 :             weights_lgr_index_all_fvm2phys(jall_fvm2phys(ie):jall_fvm2phys(ie)+jcollect_cell-1,1,ie) = i
         625      194400 :             weights_lgr_index_all_fvm2phys(jall_fvm2phys(ie):jall_fvm2phys(ie)+jcollect_cell-1,2,ie) = j
         626       97200 :             jall_fvm2phys(ie) = jall_fvm2phys(ie)+jcollect_cell
         627             :           endif
         628             :         end do
         629             :       enddo
         630       10800 :       jall_fvm2phys(ie)=jall_fvm2phys(ie)-1
         631             :       !
         632             :       ! make sure sum of area overlap weights exactly match fvm%%area_sphere_physgrid
         633             :       !
         634      140400 :       phys_area = 0.0_r8
         635      108000 :       do h=1,jall_fvm2phys(ie)
         636       97200 :         i  = weights_lgr_index_all_fvm2phys(h,1,ie); j  = weights_lgr_index_all_fvm2phys(h,2,ie)
         637      108000 :         phys_area(i,j) = phys_area(i,j) +weights_all_fvm2phys(h,1,ie)
         638             :       end do
         639      140400 :       factor(:,:) = fvm(ie)%area_sphere_physgrid(:,:)/phys_area(:,:)
         640      108000 :       do h=1,jall_fvm2phys(ie)
         641       97200 :         i  = weights_lgr_index_all_fvm2phys(h,1,ie); j  = weights_lgr_index_all_fvm2phys(h,2,ie)
         642      108000 :         weights_all_fvm2phys(h,1,ie) = weights_all_fvm2phys(h,1,ie)*factor(i,j)
         643             :       end do
         644             :
         645       10800 :       jall_phys2fvm(ie)=1
         646       43200 :       do j=1,fvm_nc
         647      140400 :         do i=1,fvm_nc
         648       97200 :           xcell(1) = xgno_fvm(i)  ; ycell(1) = ygno_fvm(j)
         649       97200 :           xcell(2) = xgno_fvm(i)  ; ycell(2) = ygno_fvm(j+1)
         650       97200 :           xcell(3) = xgno_fvm(i+1); ycell(3) = ygno_fvm(j+1)
         651       97200 :           xcell(4) = xgno_fvm(i+1); ycell(4) = ygno_fvm(j)
         652             :
         653             :           call compute_weights_cell(nvertex,.true.,&
         654             :                xcell,ycell,i,j,irecons,xgno_phys,ygno_phys,0,phys_nc+2,&
         655             :                1,phys_nc+1,1,phys_nc+1,&
         656             :                ngpc,gauss_weights,abscissae,&
         657       97200 :                weights_cell,weights_eul_index_cell,jcollect_cell,jmax_segments_cell)
         658             :
         659      129600 :           if (jcollect_cell>0) then
         660           0 :             weights_all_phys2fvm(jall_phys2fvm(ie):jall_phys2fvm(ie)+jcollect_cell-1,:,ie) &
         661     1263600 :                  = weights_cell(1:jcollect_cell,:)!/fvm(ie)%area_sphere(i,j)!da_fvm(i,j)
         662             :
         663           0 :             weights_eul_index_all_phys2fvm(jall_phys2fvm(ie):jall_phys2fvm(ie)+jcollect_cell-1,:,ie) = &
         664      486000 :                  weights_eul_index_cell(1:jcollect_cell,:)
         665      194400 :             weights_lgr_index_all_phys2fvm(jall_phys2fvm(ie):jall_phys2fvm(ie)+jcollect_cell-1,1,ie) = i
         666      194400 :             weights_lgr_index_all_phys2fvm(jall_phys2fvm(ie):jall_phys2fvm(ie)+jcollect_cell-1,2,ie) = j
         667       97200 :             jall_phys2fvm(ie) = jall_phys2fvm(ie)+jcollect_cell
         668             :           endif
         669             :         end do
         670             :       enddo
         671       10800 :       jall_phys2fvm(ie)=jall_phys2fvm(ie)-1
         672             :       !
         673             :       ! make sure sum of area overlap weights exactly matches fvm%%area_sphere_physgrid
         674             :       !
         675      140400 :       fvm_area = 0.0_r8
         676      108000 :       do h=1,jall_phys2fvm(ie)
         677       97200 :         i  = weights_lgr_index_all_phys2fvm(h,1,ie); j  = weights_lgr_index_all_phys2fvm(h,2,ie)
         678      108000 :         fvm_area(i,j) = fvm_area(i,j) +weights_all_phys2fvm(h,1,ie)
         679             :       end do
         680      140400 :       fvm_area(:,:) = fvm(ie)%area_sphere(:,:)/fvm_area(:,:)
         681      109536 :       do h=1,jall_phys2fvm(ie)
         682       97200 :         i  = weights_lgr_index_all_phys2fvm(h,1,ie); j  = weights_lgr_index_all_phys2fvm(h,2,ie)
         683      108000 :         weights_all_phys2fvm(h,1,ie) = weights_all_phys2fvm(h,1,ie)*fvm_area(i,j)
         684             :       end do
         685             :     end do
         686        1536 :     end subroutine fvm2phys_init
         687             : end module dp_mapping

Generated by: LCOV version 1.14