#include <misc.h>#include <preproc.h>module surfFileMod!=======================================================================contains!=======================================================================  subroutine surfrd(veg, wt,  &       cam_longxy, cam_latixy, cam_numlon, cam_landfrac, cam_landmask)!----------------------------------------------------------------------- ! ! Purpose: ! Read surface data file and make subgrid patches! ! Method: ! The model's surface dataset recognizes 5 basic land cover types within ! a grid cell: lake, wetland, urban, glacier, and vegetated. The vegetated ! portion of the grid cell is comprised of up to [maxpatch_pft] PFTs ! (patches). These subgrid patches are read in explicitly from the surface! dataset at the model resolution for each grid cell. This is in contrast ! to LSMv1, where the PFTs were built implicitly from biome types.!! The monthly surface dataset contains the following (on the ! [lsmlon] x [lsmlat] grid). !    o real edges of grid!    o real latitude  of grid cell (degrees)!    o real longitude of grid cell (degrees)!    o integer number of longitudes per latitude!    o integer surface type: 0 = ocean or 1 = land!    o integer soil color (1 to 9) for use with soil albedos!    o real soil texture, %sand, for thermal and hydraulic properties!    o real soil texture, %clay, for thermal and hydraulic properties!    o real % of cell covered by lake    for use as subgrid patch!    o real % of cell covered by wetland for use as subgrid patch!    o real % of cell that is urban      for use as subgrid patch!    o real % of cell that is glacier    for use as subgrid patch!    o integer PFTs!    o real % abundance PFTs (as a percent of vegetated area)!! OFFLINE MODE ONLY:! Surface grid edges -- Grids do not have to be global. ! If the model grid is read in from the surface dataset, the grid is ! assumed to be global although it does not have to be regular (it! can be a gaussian grid, for example). If the model grid is generated ! by clm2, the grid does not have to be global but it but must be ! regular and the user must define the north, east, south, and west edges! as noted below.!!    o lsmedge(1) = northern edge of grid (degrees): >  -90 and <= 90!    o lsmedge(2) = eastern edge of grid (degrees) : see following notes!    o lsmedge(3) = southern edge of grid (degrees): >= -90 and <  90!    o lsmedge(4) = western edge of grid (degrees) : see following notes!!    For partial grids, the northern and southern edges are any latitude!    between 90 (North Pole) and -90 (South Pole). Western and eastern!    edges are any longitude between -180 and 180, with longitudes!    west of Greenwich negative. That is, western edge >= -180 and < 180;!    eastern edge > western edge and <= 180.!!    For global grids, the northern and southern edges are 90 (North Pole)!    and -90 (South Pole). The western and eastern edges depend on!    whether the grid starts at Dateline or Greenwich. Regardless,!    these edges must span 360 degrees. !! Author: Gordon Bonan! !-----------------------------------------------------------------------! $Id: surfFileMod.F90,v 1.13.2.4.6.1 2002/05/13 19:25:08 erik Exp $ !-----------------------------------------------------------------------    use precision    use clm_varpar                      !parameters    use clm_varctl                      !control variables     use clm_varsur                      !surface data     use pft_varcon                      !vegetation type (PFT)     use fileutils, only : getfil    use spmdMod                         !spmd routines and variables           use areaMod                         !area avg routines    implicit none    include 'netcdf.inc'! ------------------------ arguments------------------------------------    integer , intent(out) :: veg(lsmlon,lsmlat,maxpatch) !PFT     real(r8), intent(out) :: wt(lsmlon,lsmlat,maxpatch)  !subgrid weights    real(r8), optional, intent(in) :: cam_longxy(:,:)    !cam lon values    real(r8), optional, intent(in) :: cam_latixy(:,:)    !cam lat values     integer , optional, intent(in) :: cam_numlon(:)      !cam number of longitudes     real(r8), optional, intent(in) :: cam_landfrac(:,:)  !cam fractional land    integer , optional, intent(in) :: cam_landmask(:,:)  !cam land mask! ----------------------------------------------------------------------! ------------------------ local variables -----------------------------    character(len=256) :: locfn                    !local file name    integer  :: i,j,k,m                            !indices    integer  :: ier                                !error status     integer  :: pft(lsmlon,lsmlat,maxpatch_pft)    !PFT     integer  :: ncid,dimid,varid                   !input netCDF id's    integer  :: beg4d(4),len4d(4)                  !netCDF variable edges    integer  :: nlon_i                             !number of input data longitudes    integer  :: nlat_i                             !number of input data latitudes    integer  :: nlev_i                             !number of input data levels    integer  :: maxpatch_pft_i                     !number of input data pft types    real(r8) :: pctpft(lsmlon,lsmlat,maxpatch_pft) !percent of vegetated area for PFTs    real(r8) :: sumscl                             !temporory scalar sum     real(r8) :: sumvec(lsmlon,lsmlat)              !temporary vector sum ! ----------------------------------------------------------------------    if (masterproc) then! Initialize surface data to unused value       landmask(:,:) = -999       landfrac(:,:) = -999.       latixy(:,:)   = -999.       longxy(:,:)   = -999.       soic2d(:,:)   = -999       sand3d(:,:,:) = -999.       clay3d(:,:,:) = -999.       pctlak(:,:)   = -999.       pctwet(:,:)   = -999.       pcturb(:,:)   = -999.       pctgla(:,:)   = -999.! Obtain netcdf file and read surface data       write (6,*) 'Attempting to read surface boundary data .....'       call getfil (fsurdat, locfn, 0)       call wrap_open(locfn, 0, ncid)       call wrap_inq_dimid  (ncid, 'lsmlon', dimid)       call wrap_inq_dimlen (ncid, dimid, nlon_i)       if (nlon_i /= lsmlon) then          write(6,*)'SURFRD: parameter lsmlon= ',lsmlon, &               'does not equal input nlon_i= ',nlon_i          call endrun       endif       call wrap_inq_dimid  (ncid, 'lsmlat', dimid)       call wrap_inq_dimlen (ncid, dimid, nlat_i)       if (nlat_i /= lsmlat) then          write(6,*)'SURFRD: parameter lsmlat= ',lsmlat, &               'does not equal input nlat_i= ',nlat_i          call endrun       endif       call wrap_inq_dimid  (ncid, 'nlevsoi', dimid)       call wrap_inq_dimlen (ncid, dimid, nlev_i)       if (nlev_i /= nlevsoi) then          write(6,*)'SURFRD: parameter nlevsoi= ',nlevsoi, &               'does not equal input nlev_i= ',nlev_i          call endrun       endif       call wrap_inq_dimid  (ncid, 'lsmpft', dimid)       call wrap_inq_dimlen (ncid, dimid, maxpatch_pft_i)       if (maxpatch_pft_i /= maxpatch_pft) then          write(6,*)'SURFRD: parameter maxpatch_pft',maxpatch_pft, &               'does not equal input maxpatch_pft_i= ',maxpatch_pft_i          call endrun       endif#if (defined OFFLINE)       ier = nf_inq_varid (ncid, 'EDGEN', dimid)       if (ier == NF_NOERR) then            offline_rdgrid = .false. ! variable is on history tape - regional grid       else          offline_rdgrid = .true.  ! variable is not on history tape       endif       if (.not. offline_rdgrid) then          call wrap_inq_varid (ncid, 'EDGEN', varid)          call wrap_get_var_realx (ncid, varid, lsmedge(1))                 call wrap_inq_varid (ncid, 'EDGEE', varid)          call wrap_get_var_realx (ncid, varid, lsmedge(2))                    call wrap_inq_varid (ncid, 'EDGES', varid)          call wrap_get_var_realx (ncid, varid, lsmedge(3))                    call wrap_inq_varid (ncid, 'EDGEW', varid)          call wrap_get_var_realx (ncid, varid, lsmedge(4))       endif#endif       call wrap_inq_varid (ncid, 'NUMLON', varid)       call wrap_get_var_int   (ncid, varid  , numlon)       fullgrid = .true.       do j = 1,lsmlat          if (numlon(j) < lsmlon) fullgrid = .false.       end do       if (fullgrid) then          call wrap_inq_varid (ncid, 'LONGXY' , varid)       else          call wrap_inq_varid (ncid, 'RLONGXY', varid)       endif       call wrap_get_var_realx (ncid, varid, longxy)       call wrap_inq_varid (ncid, 'LATIXY', varid)       call wrap_get_var_realx (ncid, varid, latixy)       call wrap_inq_varid (ncid, 'LANDMASK', varid)        call wrap_get_var_int (ncid, varid, landmask)       call wrap_inq_varid (ncid, 'LANDFRAC', varid)        call wrap_get_var_realx (ncid, varid, landfrac)       call wrap_inq_varid (ncid, 'SOIL_COLOR', varid)       call wrap_get_var_int   (ncid, varid, soic2d)       call wrap_inq_varid (ncid, 'PCT_SAND', varid)       call wrap_get_var_realx (ncid, varid, sand3d)       call wrap_inq_varid (ncid, 'PCT_CLAY', varid)       call wrap_get_var_realx (ncid, varid, clay3d)        call wrap_inq_varid (ncid, 'PCT_WETLAND', varid)       call wrap_get_var_realx (ncid, varid, pctwet)       call wrap_inq_varid (ncid, 'PCT_LAKE', varid)       call wrap_get_var_realx (ncid, varid, pctlak)       call wrap_inq_varid (ncid, 'PCT_GLACIER', varid)       call wrap_get_var_realx (ncid, varid, pctgla)       call wrap_inq_varid (ncid, 'PCT_URBAN', varid)       call wrap_get_var_realx (ncid, varid, pcturb)       call wrap_inq_varid (ncid, 'PFT', varid)       call wrap_get_var_int (ncid, varid, pft)       call wrap_inq_varid (ncid, 'PCT_PFT', varid)       call wrap_get_var_realx (ncid, varid, pctpft)       call wrap_close(ncid)! Define edges and area of grid cells#if (defined OFFLINE)        if (offline_rdgrid) then          call celledge (lsmlat, lsmlon, numlon, longxy, latixy, &                         lats  , lonw  )          call cellarea (lsmlat, lsmlon, numlon, lats, lonw, &                         area   )       else	          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       )       endif#else       ! CLM2 is run as part of CAM model or through the CCSM flux coupler      ! The input grid, land mask and fractional land must match the corresponding      ! CAM values. Consequently, the following check must be done on input       ! latitudes and longitudes since these are computed each time by the CAM        ! model and might be different to roundoff depending on the platform       do j = 1,lsmlat          do i = 1,numlon(j)             if ( abs(cam_latixy(i,j)-latixy(i,j)) > 1.e-6 ) then                write(6,*)'CAM latitude ',cam_latixy(i,j),' and clm2 input latitude ', &                     latixy(i,j),' has difference too large at i,j= ',i,j                call endrun             else                latixy(i,j) = cam_latixy(i,j)             endif             if ( abs(cam_longxy(i,j)-longxy(i,j)) > 1.e-6 ) then                write(6,*)'CAM longitude ',cam_longxy(i,j),' and clm2 input longitude ', &                     longxy(i,j),' has difference too large at i,j= ',i,j