dy = (mksrf_offline_edgen - mksrf_offline_edges) / lsmlat    do j = 1, lsmlat       do i = 1, lsmlon          longxy(i,j) = mksrf_offline_edgew + (2*i-1)*dx / 2.          latixy(i,j) = mksrf_offline_edges + (2*j-1)*dy / 2       end do    end do! Define edges and area of land model grid cells    call celledge (lsmlat     , lsmlon     , numlon     , longxy     ,&                   latixy     , lsmedge(1) , lsmedge(2) , lsmedge(3) ,&                   lsmedge(4) , lats       , lonw       )    call cellarea (lsmlat     , lsmlon     , numlon     , lats       ,&                   lonw       , lsmedge(1) , lsmedge(2) , lsmedge(3) ,&                   lsmedge(4) , area       )! read navy orography data and obtain fractional land    call getfil (mksrf_offline_fnavyoro, locfn, 0)    call wrap_open(locfn, 0, ncid)    call wrap_inq_dimid  (ncid, 'lon', dimid)    call wrap_inq_dimlen (ncid, dimid, nlon_i)    call wrap_inq_dimid  (ncid, 'lat', dimid)    call wrap_inq_dimlen (ncid, dimid, nlat_i)    allocate (latixy_i(nlon_i,nlat_i), stat=ier)       if (ier/=0) call endrun    allocate (longxy_i(nlon_i,nlat_i), stat=ier)       if (ier/=0) call endrun    allocate (numlon_i(nlat_i), stat=ier)    if (ier/=0) call endrun    allocate (lon_i(nlon_i+1,nlat_i), stat=ier)    if (ier/=0) call endrun    allocate (lon_i_offset(nlon_i+1,nlat_i), stat=ier)    if (ier/=0) call endrun    allocate (lat_i(nlat_i+1), stat=ier)            if (ier/=0) call endrun    allocate (area_i(nlon_i,nlat_i), stat=ier)      if (ier/=0) call endrun    allocate (mask_i(nlon_i,nlat_i), stat=ier)         if (ier/=0) call endrun    allocate (fland_i(nlon_i,nlat_i), stat=ier)         if (ier/=0) call endrun    call wrap_inq_varid (ncid, 'LATIXY', varid)    call wrap_get_var_realx (ncid, varid, latixy_i)    call wrap_inq_varid (ncid, 'LONGXY', varid)    call wrap_get_var_realx (ncid, varid, longxy_i)    call wrap_inq_varid (ncid, 'NUMLON', varid)    call wrap_get_var_int (ncid, varid, numlon_i)    call wrap_inq_varid (ncid, 'EDGEN', varid)    call wrap_get_var_realx (ncid, varid, edge_i(1))    call wrap_inq_varid (ncid, 'EDGEE', varid)    call wrap_get_var_realx (ncid, varid, edge_i(2))    call wrap_inq_varid (ncid, 'EDGES', varid)    call wrap_get_var_realx (ncid, varid, edge_i(3))    call wrap_inq_varid (ncid, 'EDGEW', varid)    call wrap_get_var_realx (ncid, varid, edge_i(4))    call wrap_inq_varid (ncid, 'LANDFRAC', varid)    call wrap_get_var_realx (ncid, varid, fland_i)    call wrap_close(ncid)! determine maximim overlap for height resolution and model grids and ! and allocate dynamic memory for overlap arrays! first determine input grid cell and cell areas    numlon_i(:) = nlon_i    call celledge (nlat_i    , nlon_i    , numlon_i  , longxy_i  , &                   latixy_i  , edge_i(1) , edge_i(2) , edge_i(3) , &                   edge_i(4) , lat_i     , lon_i     )    call cellarea (nlat_i    , nlon_i    , numlon_i  , lat_i     , &                   lon_i     , edge_i(1) , edge_i(2) , edge_i(3) , &                   edge_i(4) , area_i    )    do ji = 1, nlat_i       do ii = 1, numlon_i(ji)          mask_i(ii,ji) = 1.       end do    end do! Shift x-grid to locate periodic grid intersections. This! assumes that all lon_i(1,j) have the same value for all! latitudes j and that the same holds for lon_o(1,j)  if (lon_i(1,1) < lonw(1,1)) then     offset = 360.0  else     offset = -360.0  end if    do ji = 1, nlat_i     do ii = 1, numlon_i(ji) + 1        lon_i_offset(ii,ji) = lon_i(ii,ji) + offset     end do  end do! Process each cell on land model grid! novr_i2o - number of input grid cells that overlap each land grid cell! iovr_i2o - longitude index of overlapping input grid cell! jovr_i2o - latitude  index of overlapping input grid cell! wovr_i2o - fraction of land grid cell overlapped by input grid cell!$OMP PARALLEL DO PRIVATE (io,jo,ii,ji,n,mask_o,novr_i2o,iovr_i2o,jovr_i2o,wovr_i2o,fld_i)  do jo = 1, lsmlat     do io = 1, numlon(jo)! Determine areas of overlap and indices        mask_o = 1.        call areaini_point (io        , jo          , nlon_i  , nlat_i  , numlon_i , &                           lon_i      , lon_i_offset, lat_i   , area_i  , mask_i   , &                           lsmlon     , lsmlat      , numlon  , lonw    , lats     , &                           area(io,jo), mask_o      , novr_i2o, iovr_i2o, jovr_i2o , &                           wovr_i2o)                             ! Make area average        landfrac(io,jo) = 0.        do n = 1, novr_i2o   !overlap cell index           ii = iovr_i2o(n)  !lon index (input grid) of overlap cell           ji = jovr_i2o(n)  !lat index (input grid) of overlap cell           landfrac(io,jo) = landfrac(io,jo) + fland_i(ii,ji) * wovr_i2o(n)        end do        if (landfrac(io,jo) > 100.000001) then           write (6,*) 'MKGRID error: fland = ',landfrac(io,jo), &                ' is greater than 100 for lon,lat = ',io,jo           call endrun        end if! Global sum of output field -- must multiply by fraction of! output grid that is land as determined by input grid        fld_o(io,jo) = 0.        do n = 1, novr_i2o           ii = iovr_i2o(n)           ji = jovr_i2o(n)           fld_i = ((ji-1)*nlon_i + ii)            fld_o(io,jo) = fld_o(io,jo) + wovr_i2o(n) * fld_i         end do     end do  !end of output longitude loop  end do     !end of output latitude  loop!$OMP END PARALLEL DO! -----------------------------------------------------------------! Error check1! Compare global sum fld_o to global sum fld_i. ! -----------------------------------------------------------------! This check is true only if both grids span the same domain. ! To obtain global sum of input field must multiply by ! fraction of input grid that is land as determined by input grid  sum_fldo = 0.  do jo = 1,lsmlat     do io = 1,numlon(jo)        sum_fldo = sum_fldo + area(io,jo) * fld_o(io,jo)      end do  end do  sum_fldi = 0.  do ji = 1, nlat_i           do ii = 1, numlon_i(ji)        fld_i = ((ji-1)*nlon_i + ii)         sum_fldi = sum_fldi + area_i(ii,ji) * fld_i     end do  end do  if ( abs(mksrf_offline_edgen - mksrf_offline_edges) == 180. .and. &       abs(mksrf_offline_edgee - mksrf_offline_edgew) == 360. ) then     if ( abs(sum_fldo/sum_fldi-1.) > relerr ) then        write (6,*) 'MKGRID error: input field not conserved'        write (6,'(a30,e20.10)') 'global sum output field = ',sum_fldo        write (6,'(a30,e20.10)') 'global sum input  field = ',sum_fldi        call endrun     end if  end if! Determine land mask     where (landfrac(:,:) < flandmin)       landmask(:,:) = 0     !ocean    elsewhere       landmask(:,:) = 1     !land    endwhere! deallocate dynamic memory    deallocate (numlon_i)    deallocate (latixy_i)    deallocate (longxy_i)    deallocate (lon_i)    deallocate (lon_i_offset)    deallocate (lat_i)    deallocate (area_i)    deallocate (mask_i)    deallocate (fland_i)    write (6,*) 'Successfully makde land grid data'    write (6,*)  end subroutine create_grid_offline#endif!=======================================================================#if (defined COUP_CAM || defined COUP_CSM)  subroutine mkgrid_cam(cam_longxy, cam_latixy, cam_numlon, cam_landfrac, cam_landmask)!----------------------------------------------------------------------- ! ! Purpose: ! Generate land model grid when mode is CAM or CSM! For CAM mode get grid AND fractional land and land mask from cam model! For CSM mode get grid AND fractional land and land mask from coupler! ! Method: ! ! Author: Mariana Vertenstein! !-----------------------------------------------------------------------! ------------------------ arguments -----------------------------------    real(r8), intent(in) :: cam_longxy(:,:)      !cam longitudes     integer , intent(in) :: cam_numlon(:)        !number of cam longitudes per latitude    real(r8), intent(in) :: cam_latixy(:,:)      !cam latitudes      real(r8), intent(in) :: cam_landfrac(:,:)    !cam land fraction    integer , intent(in) :: cam_landmask(:,:)    !cam land mask! -----------------------------------------------------------------! ------------------------ local variables ------------------------    integer i,j     ! loop indices! -----------------------------------------------------------------! Determine if grid has pole points - if so, make sure that north pole! is non-land and south pole is land    if (abs((cam_latixy(1,lsmlat) - 90.)) < 1.e-6) then       pole_points = .true.       if ( masterproc ) write(6,*)'MKGRIDMOD: model has pole_points'        do i = 1,cam_numlon(1)          if (cam_landmask(i,1) /= 1) then             write(6,*)'cam grid with pole points has non-land at south pole'             write(6,*)'longitude index= ',i             call endrun          endif       end do       do i = 1,cam_numlon(lsmlat)          if (cam_landmask(i,lsmlat) == 1) then             write(6,*)'cam grid with pole points has land at north pole'             write(6,*)'longitude index= ',i             call endrun          endif       end do    else       pole_points = .false.       if ( masterproc ) write(6,*)'MKGRIDMOD: model does not have pole_points'     endif          ! Determine land grid, land mask and land fraction    numlon(:) = cam_numlon(:)      fullgrid = .true.    do j = 1,lsmlat       if (cam_numlon(j) < lsmlon) fullgrid = .false.    end do    do j = 1,lsmlat       do i = 1,numlon(j)          longxy(i,j)   = cam_longxy(i,j)           latixy(i,j)   = cam_latixy(i,j)          landmask(i,j) = cam_landmask(i,j)          landfrac(i,j) = cam_landfrac(i,j)       end do    end do! Define land grid edges and grid cell areas    call celledge (lsmlat, lsmlon, numlon, longxy, latixy, &                   lats  , lonw  )    call cellarea (lsmlat, lsmlon, numlon, lats, lonw, &                   area   )! Determine land fraction and land mask    return  end subroutine mkgrid_cam#endifend module mkgridMod