LCOV - code coverage report
Current view: | top level - dynamics/se - native_mapping.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 19 | 254 | 7.5 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 2 | 50.0 %

              Line data    Source code
           1             : module native_mapping
           2             : !
           3             : !  Create mapping files using the SE basis functions.  This module looks for the namelist 'native_mapping' in
           4             : !  file NLFileName (usually atm_in) and reads from it a list of up to maxoutgrids grid description files
           5             : !  It then creates a grid mapping file from the currently defined SE grid to the grid described in each file
           6             : !  using the SE basis functions.   The output mapping file name is generated based on the SE model resolution
           7             : !  and the input grid file name and ends in '_date_native.nc'
           8             : !
           9             :   use cam_logfile,       only : iulog
          10             :   use shr_kind_mod,      only : r8 => shr_kind_r8, shr_kind_cl
          11             :   use shr_const_mod,     only : pi=>shr_const_pi
          12             :   use cam_abortutils,    only : endrun
          13             :   use spmd_utils,        only : iam, masterproc, mpi_character, mpi_logical, mpi_integer, mpi_max, &
          14             :                                 mpicom, mstrid=>masterprocid
          15             :
          16             :   implicit none
          17             :   private
          18             :   public :: native_mapping_readnl, create_native_mapping_files, do_native_mapping
          19             :
          20             :   integer, parameter :: maxoutgrids=5
          21             :   character(len=shr_kind_cl) :: native_mapping_outgrids(maxoutgrids)
          22             :   logical, protected :: do_native_mapping
          23             :
          24             : !=============================================================================================
          25             : contains
          26             : !=============================================================================================
          27             :
          28        1536 : subroutine native_mapping_readnl(NLFileName)
          29             :
          30             :    use units,          only : getunit, freeunit
          31             :    use namelist_utils, only : find_group_name
          32             :
          33             :    character(len=*), intent(in) :: NLFileName
          34             :
          35             :    character(len=shr_kind_cl) :: mappingfile, fname
          36             :
          37             :    namelist /native_mapping_nl/  native_mapping_outgrids
          38             :    integer :: nf, unitn, ierr
          39             :    logical :: exist
          40             :    character(len=*), parameter ::  sub="native_mapping_readnl"
          41             :    !-----------------------------------------------------------------------------
          42             :
          43        1536 :    do_native_mapping=.false.
          44             :
          45        9216 :    do nf=1,maxoutgrids
          46        9216 :       native_mapping_outgrids(nf)=''
          47             :    enddo
          48             :
          49        1536 :    if(masterproc) then
          50           2 :       exist=.true.
          51           2 :       write(iulog,*) sub//': Check for native_mapping_nl namelist in ',trim(nlfilename)
          52           2 :       unitn = getunit()
          53           2 :       open( unitn, file=trim(nlfilename), status='old' )
          54             :
          55           2 :       call find_group_name(unitn, 'native_mapping_nl', status=ierr)
          56           2 :       if(ierr/=0) then
          57           2 :          write(iulog,*) sub//': No native_mapping_nl namelist found'
          58           2 :          exist=.false.
          59             :       end if
          60           2 :       if(exist) then
          61           0 :          read(unitn, native_mapping_nl, iostat=ierr)
          62           0 :          if(ierr/=0) then
          63           0 :             call endrun(sub//': namelist read returns an error condition for native_mapping_nl')
          64             :          end if
          65           0 :          if(len_trim(native_mapping_outgrids(1))==0) exist=.false.
          66             :       end if
          67           2 :       close(unitn)
          68           2 :       call freeunit(unitn)
          69             :    end if
          70             :
          71        1536 :    call mpi_bcast(exist, 1, mpi_logical, mstrid, mpicom, ierr)
          72        1536 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: exist")
          73             :
          74        1536 :    if(.not. exist) return
          75             :
          76           0 :    call mpi_bcast(native_mapping_outgrids, maxoutgrids*shr_kind_cl, mpi_character, mstrid, mpicom, ierr)
          77           0 :    if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: native_mapping_outgrids")
          78             :
          79           0 :    do_native_mapping=.true.
          80             :
          81             : end subroutine native_mapping_readnl
          82             :
          83             : !=============================================================================================
          84             :
          85           0 : subroutine create_native_mapping_files(par, elem, maptype, ncol, clat, clon, areaa)
          86             :
          87             :     use parallel_mod, only : parallel_t, global_shared_buf, global_shared_sum
          88             :     use global_norms_mod, only: wrap_repro_sum
          89             :     use cam_pio_utils, only : cam_pio_openfile, cam_pio_createfile
          90             :     use element_mod, only : element_t
          91             :     use hybrid_mod, only : hybrid_t, config_thread_region
          92             :     use pio, only : pio_noerr, pio_openfile, pio_createfile, pio_closefile, &
          93             :          pio_get_var, pio_put_var, pio_write_darray,pio_int, pio_double, &
          94             :          pio_def_var, pio_put_att, pio_global, file_desc_t, var_desc_t, &
          95             :          io_desc_t, pio_internal_error,pio_inq_dimlen, pio_inq_varid, &
          96             :          pio_get_att, pio_enddef, pio_bcast_error,pio_internal_error, &
          97             :          pio_def_dim, pio_inq_dimid, pio_seterrorhandling, pio_initdecomp
          98             :     use quadrature_mod, only : quadrature_t, gauss, gausslobatto
          99             :     use interpolate_mod, only : interpdata_t, cube_facepoint_ne, interpolate_scalar, set_interp_parameter, interp_init, &
         100             :          get_interp_parameter
         101             :     use coordinate_systems_mod, only : spherical_polar_t, cartesian2d_t
         102             :     use dimensions_mod, only : nelemd, ne, np, npsq, nelem
         103             :     use reduction_mod, only : ParallelMin,ParallelMax
         104             :     use cube_mod, only : convert_gbl_index
         105             :     use infnan, only : isnan
         106             :     use dof_mod, only : CreateMetaData
         107             :     use thread_mod,     only: omp_get_thread_num
         108             :     use datetime_mod, only: datetime
         109             :
         110             :
         111             :     use cam_history_support, only : fillvalue
         112             :
         113             :
         114             :     type(parallel_t), intent(in) :: par
         115             :     type(element_t),  intent(in) :: elem(:)
         116             :     character(len=*), intent(in) :: maptype
         117             :     integer,          intent(in) :: ncol
         118             :     real(r8),         intent(in) :: clat(ncol)
         119             :     real(r8),         intent(in) :: clon(ncol)
         120             :     real(r8),         intent(in) :: areaa(ncol)
         121             :
         122             :     character(len=shr_kind_cl) :: mappingfile, fname
         123             :
         124             :
         125             :     type(hybrid_t) :: hybrid
         126             :     logical :: exist
         127             :
         128             :     type (spherical_polar_t) :: sphere
         129             :     type(file_desc_t) :: ogfile, agfile
         130           0 :     type (interpdata_t)  :: interpdata(nelemd)
         131             :     integer :: ierr, dimid, npts, vid
         132           0 :     real(r8), allocatable :: lat(:), lon(:)
         133             :     integer :: i, ii, ie2, je2, ie, je, face_no, face_no2, k, j, n, ngrid, tpts, nf, number
         134             :     real(r8) :: countx, count_max, count_total
         135           0 :     integer :: fdofp(np,np,nelemd)
         136             :     type (cartesian2D_t) :: cart
         137             :     real(r8) :: f(np,np)
         138           0 :     real(r8), allocatable :: h(:), h1d(:)
         139           0 :     integer, allocatable :: grid_imask(:), row(:), col(:), ldof(:), dg_dims(:)
         140             :     integer :: ns_dim, cnt, na_dim, nb_dim, sg_dim, dg_dim
         141             :     type(var_desc_t) :: rowid, colid, sid, xca_id, yca_id, xcb_id, ycb_id, maskb_id, maska_id
         142             :     type(var_desc_t) :: areaA_id, areaB_id, dg_id, sg_id
         143             :     type(io_desc_t) :: iodesci, iodescd
         144             :     character(len=12) :: unit_str
         145           0 :     real(r8), allocatable :: areaB(:)
         146           0 :     integer :: cntperelem_in(nelem),  cntperelem_out(nelem)
         147             :     integer :: ithr, dg_rank, substr1, substr2
         148             :
         149             :     type(interpdata_t), pointer :: mapping_interpolate(:)
         150             :     character(len=8) :: cdate, ctime
         151             :     integer :: olditype, oldnlat, oldnlon, itype
         152             :
         153             :
         154             :
         155           0 :     if(.not. do_native_mapping) return
         156             :
         157           0 :     if (maptype=='native') then
         158           0 :        itype=0
         159           0 :     else if (maptype=='bilin') then
         160           0 :        itype=1
         161             :     else
         162           0 :        call endrun('bad interp_type')
         163             :     endif
         164             :
         165             :
         166             :
         167             :
         168           0 :     if(iam > par%nprocs) then
         169             :        ! The special case of npes_se < npes_cam is not worth dealing with here
         170           0 :        call endrun('Native mapping code requires npes_se==npes_cam')
         171             :     end if
         172             :
         173             :
         174           0 :     call interp_init()
         175             :
         176             :
         177           0 :     oldnlon = get_interp_parameter('nlon')
         178           0 :     oldnlat = get_interp_parameter('nlat')
         179           0 :     olditype = get_interp_parameter('itype')
         180             :
         181           0 :     call datetime(cdate, ctime)
         182             :
         183           0 :     do nf=1,maxoutgrids
         184           0 :        fname = native_mapping_outgrids(nf)
         185           0 :        if(masterproc) then
         186           0 :           write(iulog,*) 'looking for target grid = ',trim(fname)
         187             :        endif
         188           0 :        if(len_trim(fname)==0) cycle
         189           0 :        inquire(file=fname,exist=exist)
         190           0 :        if(.not. exist) then
         191           0 :           write(iulog,*) 'WARNING: Could not find or open grid file ',fname
         192           0 :           cycle
         193             :        end if
         194           0 :        if(masterproc) then
         195           0 :           write(iulog,*) 'Creating ',trim(maptype),' mapping to grid ',fname
         196             :        endif
         197           0 :        call cam_pio_openfile( ogfile, fname, 0)
         198             :
         199           0 :        ierr = pio_inq_dimid( ogfile, 'grid_size', dimid)
         200           0 :        ierr = pio_inq_dimlen( ogfile, dimid, npts)
         201           0 :        allocate(lat(npts), lon(npts), grid_imask(npts), areab(npts))
         202             :
         203           0 :        ierr = pio_inq_dimid( ogfile, 'grid_rank', dimid)
         204           0 :        ierr = pio_inq_dimlen(ogfile, dimid, dg_rank)
         205           0 :        allocate(dg_dims(dg_rank))
         206           0 :        ierr = pio_inq_varid( ogfile, 'grid_dims', vid)
         207           0 :        ierr = pio_get_var( ogfile, vid, dg_dims)
         208             :
         209             :
         210           0 :        ierr = pio_inq_varid( ogfile, 'grid_center_lat', vid)
         211           0 :        ierr = pio_get_var(ogfile, vid, lat)
         212           0 :        ierr = pio_get_att(ogfile, vid, 'units', unit_str)
         213             :
         214           0 :        ierr = pio_inq_varid( ogfile, 'grid_center_lon', vid)
         215           0 :        ierr = pio_get_var(ogfile, vid, lon)
         216             :
         217           0 :        call pio_seterrorhandling(ogfile, PIO_BCAST_ERROR)
         218           0 :        ierr = pio_inq_varid( ogfile, 'grid_area', vid)
         219           0 :        call pio_seterrorhandling(ogfile, PIO_INTERNAL_ERROR)
         220           0 :        if(ierr == PIO_NOERR) then
         221           0 :           ierr = pio_get_var(ogfile, vid, areaB)
         222             :        else
         223           0 :           areaB=fillvalue
         224             :        end if
         225             :
         226           0 :        if(unit_str .eq. 'degrees') then
         227           0 :           lat = lat * pi/180_r8
         228           0 :           lon = lon * pi/180_r8
         229             :        end if
         230             :
         231           0 :        ierr = pio_inq_varid( ogfile, 'grid_imask', vid)
         232           0 :        ierr = pio_get_var(ogfile, vid, grid_imask)
         233           0 :        call pio_closefile(ogfile)
         234             :
         235           0 :        do ie=1,nelemd
         236           0 :           interpdata(ie)%n_interp=0
         237             :        end do
         238             :
         239           0 :        call set_interp_parameter('itype',itype)   ! itype=0 native, 1 for bilinear
         240           0 :        if(lon(1)==lon(2)) then
         241           0 :           call set_interp_parameter('nlon',dg_dims(1))
         242           0 :           call set_interp_parameter('nlat',dg_dims(2))
         243             :        else
         244           0 :           call set_interp_parameter('nlon',dg_dims(2))
         245           0 :           call set_interp_parameter('nlat',dg_dims(1))
         246             :        end if
         247             :
         248             :
         249             :
         250             :
         251             :
         252             : !       call setup_latlon_interp(elem, cam_interpolate, hybrid, 1, nelemd)
         253             :        ! go through once, counting the number of points on each element
         254             :
         255           0 :        sphere%r=1
         256           0 :        do i=1,npts
         257           0 :           if(grid_imask(i)==1) then
         258           0 :              sphere%lat=lat(i)
         259           0 :              sphere%lon=lon(i)
         260           0 :              call cube_facepoint_ne(sphere, ne, cart, number)  ! new interface
         261           0 :              if (number /= -1) then
         262           0 :                 do ii=1,nelemd
         263           0 :                    if (number == elem(ii)%vertex%number) then
         264           0 :                       interpdata(ii)%n_interp = interpdata(ii)%n_interp + 1
         265           0 :                       exit
         266             :                    endif
         267             :                 enddo
         268             :              endif
         269             :
         270             :
         271           0 :              if(masterproc) then
         272           0 :                 if(mod(i,npts/10).eq.1) then
         273           0 :                    print *,'finished point ',i,' of ',npts
         274             :                 endif
         275             :              end if
         276             :           end if
         277             :        enddo
         278             :
         279           0 :        hybrid = config_thread_region(par,'serial')
         280             : !       ithr=omp_get_thread_num()
         281             : !       hybrid = hybrid_create(par,ithr,1)
         282             :
         283             :
         284             :
         285             :        ! check if every point in interpolation grid was claimed by an element:
         286           0 :        countx=sum(interpdata(1:nelemd)%n_interp)
         287           0 :        global_shared_buf(1,1) = countx
         288           0 :        call wrap_repro_sum(nvars=1, comm=hybrid%par%comm, nsize=1)
         289           0 :        count_total = global_shared_sum(1)
         290           0 :        tpts = sum(grid_imask)
         291           0 :        if (count_total /= tpts ) then
         292           0 :           write(iulog,*)__FILE__,__LINE__,iam, count_total, tpts, npts
         293           0 :           call endrun('Error setting up interpolation grid count_total<>npts')
         294             :        endif
         295             :
         296           0 :        countx=maxval(interpdata(1:nelemd)%n_interp)
         297           0 :        count_max = ParallelMax(countx,hybrid)
         298             :
         299           0 :        if (masterproc) then
         300           0 :           write(iulog,'(a,f8.1)') 'Average number of interpolation points per element: ',count_total/real(6*ne*ne)
         301           0 :           write(iulog,'(a,f8.0)') 'Maximum number of interpolation points on any element: ',count_max
         302             :        endif
         303             :
         304             :
         305             :        ! allocate storage
         306           0 :        do ii=1,nelemd
         307           0 :           ngrid = interpdata(ii)%n_interp
         308           0 :           allocate(interpdata(ii)%interp_xy( ngrid ) )
         309           0 :           allocate(interpdata(ii)%ilat( ngrid ) )
         310           0 :           allocate(interpdata(ii)%ilon( ngrid ) )
         311           0 :           interpdata(ii)%n_interp=0  ! reset counter
         312             :        enddo
         313             :
         314             :        ! now go through the list again, adding the coordinates
         315             :        ! if this turns out to be slow, then it can be done in the loop above
         316             :        ! but we have to allocate and possibly resize the interp_xy() array.
         317           0 :        do i=1,npts
         318           0 :           if(grid_imask(i)==1) then
         319           0 :              sphere%lat=lat(i)
         320           0 :              sphere%lon=lon(i)
         321           0 :              call cube_facepoint_ne(sphere, ne, cart, number)  ! new interface
         322           0 :              if (number /= -1) then
         323           0 :                 do ii=1,nelemd
         324           0 :                    if (number == elem(ii)%vertex%number) then
         325           0 :                       ngrid = interpdata(ii)%n_interp + 1
         326           0 :                       interpdata(ii)%n_interp = ngrid
         327           0 :                       interpdata(ii)%interp_xy( ngrid ) = cart
         328           0 :                       interpdata(ii)%ilon( ngrid ) = i
         329           0 :                       interpdata(ii)%ilat( ngrid ) = i
         330             :                    endif
         331             :                 enddo
         332             :              endif
         333             :           end if
         334             :        end do
         335             :
         336             :
         337           0 :        allocate(h(int(countx)))
         338           0 :        allocate(h1d(int(countx)*npsq*nelemd))
         339           0 :        allocate(row(int(countx)*npsq*nelemd))
         340           0 :        allocate(col(int(countx)*npsq*nelemd))
         341             :
         342           0 :        row = 0
         343           0 :        col = 0
         344             :
         345           0 :        ngrid=0
         346           0 :        cntperelem_in=0
         347           0 :        call CreateMetaData(hybrid%par, elem, fdofp=fdofp)
         348             :
         349           0 :        do ie=1,nelemd
         350             :           ii=0
         351           0 :           do j=1,np
         352           0 :              do i=1,np
         353           0 :                 ii=ii+1
         354           0 :                 f = 0.0_R8
         355           0 :                 f(i,j) = 1.0_R8
         356           0 :                 h = 0
         357           0 :                 call interpolate_scalar(interpdata(ie), f, np, 0, h(:))
         358             :
         359           0 :                 do n=1,interpdata(ie)%n_interp
         360           0 :                    if(any(isnan(h ))) then
         361             :
         362           0 :                       call endrun('nan generated')
         363             :                    end if
         364           0 :                    if(h(n)/=0) then
         365           0 :                       ngrid=ngrid+1
         366           0 :                       h1d(ngrid) = h(n)
         367           0 :                       row(ngrid) = interpdata(ie)%ilon(n)
         368           0 :                       col(ngrid) =  fdofp(i,j,ie)
         369           0 :                       cntperelem_in(elem(ie)%Globalid)=cntperelem_in(elem(ie)%Globalid)+1
         370             :                    end if
         371             :                 enddo
         372             :
         373             :              enddo
         374             :           end do
         375             :        end do
         376             :
         377           0 :        countx=ngrid
         378           0 :        global_shared_buf(1,1) = countx
         379           0 :        call wrap_repro_sum(nvars=1, comm=hybrid%par%comm, nsize=1)
         380           0 :        count_total = global_shared_sum(1)
         381             :
         382             :
         383           0 :        call mpi_allreduce(cntperelem_in, cntperelem_out, nelem, MPI_INTEGER, MPI_MAX, par%comm, ierr)
         384             :
         385             :
         386           0 :        allocate(ldof(ngrid))
         387           0 :        ldof = 0
         388           0 :        ii=1
         389           0 :        do ie=1,nelemd
         390           0 :           if(elem(ie)%GlobalID==1) then
         391             :              cnt = 0
         392             :           else
         393           0 :              cnt = sum(cntperelem_out(1:elem(ie)%globalid-1))
         394             :           endif
         395           0 :           do i=1,cntperelem_out(elem(ie)%globalid)
         396           0 :              ldof(ii) = cnt+i
         397           0 :              ii=ii+1
         398             :           end do
         399             :        end do
         400             :
         401           0 :        deallocate(h)
         402             :
         403           0 :        ngrid = int(count_total)
         404             :
         405           0 :        substr1 = index(fname,'/',BACK=.true.)
         406           0 :        substr2 = index(fname,'.nc',BACK=.true.)
         407             :
         408           0 :        if(ne<100) then
         409           0 :           write(mappingfile,113) ne,np,fname(substr1+1:substr2-1),trim(maptype),cdate(7:8),cdate(1:2),cdate(4:5)
         410           0 :        else if(ne<1000) then
         411           0 :           write(mappingfile,114) ne,np,fname(substr1+1:substr2-1),trim(maptype),cdate(7:8),cdate(1:2),cdate(4:5)
         412             :        else
         413           0 :           write(mappingfile,115) ne,np,fname(substr1+1:substr2-1),trim(maptype),cdate(7:8),cdate(1:2),cdate(4:5)
         414             :        end if
         415             :
         416             : 113    format('map_ne',i2.2,'np',i1,'_to_',a,'_',a,'_',3a2,'.nc')
         417             : 114    format('map_ne',i3.3,'np',i1,'_to_',a,'_',a,'_',3a2,'.nc')
         418             : 115    format('map_ne',i4.4,'np',i1,'_to_',a,'_',a,'_',3a2,'.nc')
         419             :
         420           0 :        call cam_pio_createfile( ogfile,mappingfile , 0)
         421             :
         422           0 :        ierr = pio_def_dim( ogfile, 'n_a', ncol, na_dim)
         423           0 :        ierr = pio_def_dim( ogfile, 'n_b', npts, nb_dim)
         424           0 :        ierr = pio_def_dim( ogfile, 'n_s', ngrid, ns_dim)
         425             :
         426           0 :        ierr = pio_def_dim( ogfile, 'src_grid_rank', 1, sg_dim)
         427           0 :        ierr = pio_def_var( ogfile, 'src_grid_dims',pio_int, (/sg_dim/),sg_id)
         428             :
         429           0 :        ierr = pio_def_dim( ogfile, 'dst_grid_rank',dg_rank, dg_dim)
         430           0 :        ierr = pio_def_var( ogfile, 'dst_grid_dims',pio_int, (/dg_dim/),dg_id)
         431             :
         432             :
         433             :
         434             :
         435             :
         436           0 :        ierr = pio_def_var( ogfile, 'col', pio_int, (/ns_dim/), colid)
         437           0 :        ierr = pio_def_var( ogfile, 'row', pio_int, (/ns_dim/), rowid)
         438           0 :        ierr = pio_def_var( ogfile, 'S', pio_double, (/ns_dim/), sid)
         439             :
         440           0 :        ierr = pio_def_var( ogfile, 'xc_a', pio_double, (/na_dim/), xca_id)
         441           0 :        ierr = pio_def_var( ogfile, 'yc_a', pio_double, (/na_dim/), yca_id)
         442             :
         443           0 :        ierr = pio_def_var( ogfile, 'xc_b', pio_double, (/nb_dim/), xcb_id)
         444           0 :        ierr = pio_def_var( ogfile, 'yc_b', pio_double, (/nb_dim/), ycb_id)
         445             :
         446           0 :        ierr = pio_def_var( ogfile, 'area_a', pio_double, (/na_dim/), areaA_id)
         447           0 :        ierr = pio_def_var( ogfile, 'area_b', pio_double, (/nb_dim/), areaB_id)
         448           0 :        ierr = pio_put_att( ogfile, areaB_id, '_FillValue',fillvalue)
         449             :
         450           0 :        ierr = pio_def_var( ogfile, 'mask_a', pio_int, (/na_dim/), maska_id)
         451           0 :        ierr = pio_def_var( ogfile, 'mask_b', pio_int, (/nb_dim/), maskb_id)
         452             :
         453             :
         454             :
         455           0 :        ierr = pio_put_att( ogfile, xca_id, 'units','radians')
         456           0 :        ierr = pio_put_att( ogfile, yca_id, 'units','radians')
         457           0 :        ierr = pio_put_att( ogfile, xcb_id, 'units','radians')
         458           0 :        ierr = pio_put_att( ogfile, ycb_id, 'units','radians')
         459             :
         460           0 :        ierr = pio_put_att( ogfile, PIO_GLOBAL, 'title', 'SE NATIVE Regridding Weights')
         461           0 :        ierr = pio_put_att( ogfile, PIO_GLOBAL, 'normalization', 'none')
         462           0 :        if (itype==0 ) then
         463           0 :           ierr = pio_put_att( ogfile, PIO_GLOBAL, 'map_method', 'Spectral-Element remapping')
         464           0 :        else if (itype==1) then
         465           0 :           ierr = pio_put_att( ogfile, PIO_GLOBAL, 'map_method', 'Bilinear remapping')
         466             :        endif
         467           0 :        ierr = pio_put_att( ogfile, PIO_GLOBAL, 'conventions', 'NCAR-CSM')
         468             :
         469           0 :        ierr = pio_put_att( ogfile, PIO_GLOBAL, 'grid_file_out', fname  )
         470           0 :        ierr = pio_put_att( ogfile, PIO_GLOBAL, 'grid_file_atm', 'none - model generated')
         471             :
         472             :
         473           0 :        ierr = pio_enddef ( ogfile )
         474             :
         475           0 :        ierr = pio_put_var(ogfile, sg_id, ncol)
         476           0 :        ierr = pio_put_var(ogfile, dg_id, dg_dims(1:dg_rank))
         477             :
         478             :
         479           0 :        call pio_initdecomp( ogfile%iosystem, pio_int, (/ngrid/), ldof, iodesci)
         480           0 :        call pio_initdecomp( ogfile%iosystem, pio_double, (/ngrid/), ldof, iodescd)
         481             :
         482           0 :        call pio_write_darray(ogfile, colid, iodesci, col, ierr)
         483           0 :        call pio_write_darray(ogfile, rowid, iodesci, row, ierr)
         484           0 :        call pio_write_darray(ogfile, sid, iodescd, h1d, ierr)
         485             :
         486             :
         487           0 :        ierr = pio_put_var(ogfile, xcb_id, lon)
         488           0 :        ierr = pio_put_var(ogfile, ycb_id, lat)
         489             :
         490           0 :        ierr = pio_put_var(ogfile, xca_id, clon)
         491           0 :        ierr = pio_put_var(ogfile, yca_id, clat)
         492             :
         493           0 :        ierr = pio_put_var(ogfile, maskb_id, grid_imask)
         494           0 :        deallocate(grid_imask)
         495             :
         496           0 :        ierr = pio_put_var(ogfile, areaA_id, areaA)
         497           0 :        ierr = pio_put_var(ogfile, areaB_id, areaB)
         498           0 :        deallocate(areaB)
         499             :
         500           0 :        allocate(grid_imask(ncol))
         501           0 :        grid_imask=1
         502             :
         503           0 :        ierr = pio_put_var(ogfile, maska_id, grid_imask)
         504             :
         505           0 :        call pio_closefile(ogfile)
         506             :
         507           0 :        deallocate(grid_imask, lat,lon, h1d, col, row, dg_dims, ldof)
         508           0 :        do ii=1,nelemd
         509           0 :           if(associated(interpdata(ii)%interp_xy))then
         510           0 :              deallocate(interpdata(ii)%interp_xy)
         511             :           endif
         512           0 :           if(associated(interpdata(ii)%ilat))then
         513           0 :              deallocate(interpdata(ii)%ilat)
         514             :           endif
         515           0 :           if (associated(interpdata(ii)%ilon))then
         516           0 :              deallocate(interpdata(ii)%ilon)
         517             :           endif
         518             :        end do
         519             :
         520             :
         521             :     end do
         522             :
         523           0 :     call set_interp_parameter('itype',olditype)
         524           0 :     call set_interp_parameter('nlon',oldnlon)
         525           0 :     call set_interp_parameter('nlat',oldnlat)
         526             :
         527             :
         528           0 :   end subroutine create_native_mapping_files
         529             :
         530             :
         531             :
         532             :
         533             :
         534             : end module native_mapping

Generated by: LCOV version 1.14