mkgridmod.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <preproc.h>

module mkgridMod

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Routines to create land model grid
! 
! Method: 
!
! Author: Mariana Vertenstein
! 
!-----------------------------------------------------------------------
! $Id: mkgridMod.F90,v 1.9.2.2 2002/04/27 15:38:53 erik Exp $
!-----------------------------------------------------------------------


  use precision
  use clm_varpar          !parameters 
  use clm_varsur          !surface variables 
  use clm_varctl          !run control variables
  use fileutils, only : getfil
  use areaMod             !area averaging routines
  use spmdMod, only: masterproc
  implicit none

!=======================================================================
contains
!=======================================================================

#if (defined OFFLINE) 

  subroutine mkgrid_offline

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Obtain land model grid
! 
! Method: 
! If namelist variable mksrf_offline_fgrid is the empty string, then
! the corresponding file will be used to determine the land model grid
! If namelist variable mksrf_offline is not the empty string, then 
! the land model grid will be generated at run time 
!
! Author: Mariana Vertenstein
! 
!-----------------------------------------------------------------------

    if (mksrf_offline_fgrid /= ' ') then
       call read_grid_offline
       offline_rdgrid = .true.
    else
       call create_grid_offline
       offline_rdgrid = .false.
    endif

  end subroutine mkgrid_offline

!=======================================================================

  subroutine read_grid_offline

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Read land model grid when mode is offline.
! 
! Method: 
! If namelist variable mksrf_offline_fgrid is the empty string, then
! the corresponding file will be used to determine the land model grid
! Assume that the input data file has the land grid in the following 
! form:
!    lon                    => dimension
!    lat                    => dimension
!    lon(lsmlon)            => full grid longitudes
!    nlon(lsmlat)           => reduced grid number of lats per lon
!    rlon(lsmlon,lsmlat)    => reduced grid longitudes
!    lat(lsmlat)            => grid latitudes
!    oro(lsmlon,lsmlat)     => 2d land mask
!    landfrac(lsmlon,lsmlat)=> 2d land fraction
!
! Author: Mariana Vertenstein
! 
!-----------------------------------------------------------------------

    include 'netcdf.inc'

! ------------------------ local variables ------------------------
    integer  :: i,j,k,n                 !indices
    integer  :: nlon_i                  !number of input data longitudes
    integer  :: nlat_i                  !number of input data latitudes
    integer  :: ncid                    !netCDF file id 
    integer  :: dimid                   !netCDF dimension id
    integer  :: varid                   !netCDF variable id
    integer  :: ret                     !netCDF return code
    real(r8) :: lon(lsmlon)             !input longitude array (full grid)
    real(r8) :: lat(lsmlat)             !input latitude array (full grid)
    real(r8) :: oro(lsmlon,lsmlat)      !input oro field 
    character(len=256) :: locfn         !local file name
    logical  :: lndfrac                 !true if landfrac exists in file 
!-----------------------------------------------------------------------

    write (6,*) 'Attempting to read land grid data .....'
    write (6,'(72a1)') ("-",i=1,60)
    
    call getfil (mksrf_offline_fgrid, locfn, 0)
    call wrap_open(locfn, 0, ncid)
    
    call wrap_inq_dimid  (ncid, 'lon', dimid)
    call wrap_inq_dimlen (ncid, dimid, nlon_i)
    if (nlon_i /= lsmlon) then
       write(6,*)'RDGRID_OFFLINE: parameter lsmlon= ',lsmlon, &
            'does not equal input nlon_i= ',nlon_i
       call endrun
    endif
    
    call wrap_inq_dimid  (ncid, 'lat', dimid)
    call wrap_inq_dimlen (ncid, dimid, nlat_i)
    if (nlat_i /= lsmlat) then
       write(6,*)'RDGRID_OFFLINE: parameter lsmlat= ',lsmlat, &
            'does not equal input nlat_i= ',nlat_i
       call endrun
    endif
    
! Determine grid longitudes for either full or reduced grid
! if variable 'rlon' is not on grid file then have full grid
! if variable 'rlon' is on grid file then have reduced grid

    ret = nf_inq_varid (ncid, 'rlon', dimid)
    if (ret == NF_NOERR) then
       fullgrid = .false.
    else
       fullgrid = .true.
    endif

    if (fullgrid) then
       numlon(:) = lsmlon
       call wrap_inq_varid (ncid, 'lon' , varid)
       call wrap_get_var_realx (ncid, varid, lon)
       do j = 1,lsmlat
          do i = 1,lsmlon
             longxy(i,j) = lon(i)
          end do
       end do
    else
       call wrap_inq_varid (ncid, 'nlon' , varid)
       call wrap_get_var_int (ncid, varid, numlon)
       call wrap_inq_varid (ncid, 'rlon' , varid)
       call wrap_get_var_realx (ncid, varid, longxy)
    endif

! Determine grid latitudes

    call wrap_inq_varid (ncid, 'lat' , varid)
    call wrap_get_var_realx (ncid, varid, lat)
    do j = 1,lsmlat
       do i =1,lsmlon
          latixy(i,j) = lat(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 mask and land fraction

    call wrap_inq_varid (ncid, 'ORO' , varid)
    call wrap_get_var_realx (ncid, varid, oro)

! Get land fraction if it exists, otherwise set according to land mask

    ret = nf_inq_varid (ncid, 'LANDFRAC', dimid)
    if (ret == NF_NOERR) then
       lndfrac = .true.
    else
       lndfrac = .false.
    endif

    if (lndfrac) then
      call wrap_inq_varid (ncid, 'LANDFRAC' , varid)
      call wrap_get_var_realx (ncid, varid, landfrac)
      do j = 1,lsmlat
         do i = 1,numlon(j)
            if (nint(oro(i,j)) == 1) then
               landmask(i,j) = 1
            else
               landmask(i,j) = 0
            endif
         end do
      end do
    else
      do j = 1,lsmlat
         do i = 1,numlon(j)
            if (nint(oro(i,j)) == 1) then
               landmask(i,j) = 1
               landfrac(i,j) = 1.0
            else
               landmask(i,j) = 0
               landfrac(i,j) = 0.
            endif
         end do
      end do
    endif

    write (6,'(72a1)') ("-",i=1,60)
    write (6,*) 'Successfully read land grid data'
    write (6,*)

  end subroutine read_grid_offline

!=======================================================================

  subroutine create_grid_offline

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Generate land model grid when mode is offline.
! 
! Method: 
! Surface grid edges -- Grids do not have to be global. To allow this, grids 
! must define the north, east, south, and west edges:
!
! If namelist variable mksrf_offline is not the empty string, then the land model
! grid will be generated at run time using the settings of the
! namelist variables
!    o mksrf_offline_fnavyoro : 20 min navy orography dataset
!    o mksrf_offline_edgen (edge(1)) : northern edge of grid (degrees): >  -90 and <= 90
!    o mksrf_offline_edgee (edge(2)) : eastern edge of grid (degrees) : see following notes
!    o mksrf_offline_edges (edge(3)) : southern edge of grid (degrees): >= -90 and <  90
!    o mksrf_offline_edgew (edge(4)) : western edge of grid (degrees) : see following notes
!
! For partial grids, 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, northern and southern edges are 90 (North Pole)
! and -90 (South Pole). The western edge of the longitude grid starts 
! at the dateline if the grid is generated (the longitudes for each grid 
! cell correspond with the edges (i.e., range from -180 to 180)). 
! 
! Author: Mariana Vertenstein
! 
!-----------------------------------------------------------------------

! ------------------------ local variables ------------------------
    character(len=256) :: locfn                !local file name

    integer  :: i,j,k,n                        !indices
    integer  :: ii,ji,io,jo                    !indices
    integer  :: ncid                           !netCDF file id 
    integer  :: dimid                          !netCDF dimension id
    integer  :: varid                          !netCDF variable id
    integer  :: ier                            !error status   

    integer  :: nlon_i                         !input number of longitudes
    integer  :: nlat_i                         !input number of latitudes
    real(r8) :: dx                             !land model cell width
    real(r8) :: dy                             !land model cell length
    real(r8) :: edge_i(4)                      !input grid: N,E,S,W edges (degrees)
    real(r8), allocatable :: latixy_i(:,:)     !input grid: latitude (degrees)
    real(r8), allocatable :: longxy_i(:,:)     !input grid: longitude (degrees)
    integer , allocatable :: numlon_i(:)       !input grid: number longitude points by lat
    real(r8), allocatable :: lon_i(:,:)        !input grid: longitude, west edge (degrees)
    real(r8), allocatable :: lon_i_offset(:,:) !input grid: longitude, west edge (degrees)
    real(r8), allocatable :: lat_i(:)          !input grid: latitude, south edge (degrees)
    real(r8), allocatable :: area_i(:,:)       !input grid: cell area
    real(r8), allocatable :: mask_i(:,:)       !input grid: mask (0, 1)
    real(r8), allocatable :: fland_i(:,:)      !input grid: fractional land

    real(r8) :: mask_o                         !output grid: mask (0, 1)
    integer  :: novr_i2o                       !number of overlapping input cells
    integer  :: iovr_i2o(maxovr)               !lon index of overlap input cell
    integer  :: jovr_i2o(maxovr)               !lat index of overlap input cell
    real(r8) :: wovr_i2o(maxovr)               !weight    of overlap input cell
    real(r8) :: offset                         !used to shift x-grid 360 degrees
    
    real(r8) :: fld_o(lsmlon,lsmlat)           !output grid: dummy field 
    real(r8) :: fld_i                          !input grid: dummy field 
    real(r8) :: sum_fldo                       !global sum of dummy output field
    real(r8) :: sum_fldi                       !global sum of dummy input field
    real(r8) :: relerr = 0.00001               !max error: sum overlap weights ne 1
    
    real(r8) :: flandmin = 0.50                !minimum land fraction for grid cell to be called land  
! -----------------------------------------------------------------

! Set numlon to uniform grid (offline ASSUMES that never have a reduced grid)

    numlon(:) = lsmlon

! Determine model grid edges

    lsmedge(1) = mksrf_offline_edgen
    lsmedge(2) = mksrf_offline_edgee
    lsmedge(3) = mksrf_offline_edges
    lsmedge(4) = mksrf_offline_edgew

! Determine grid longitudes and latitudes in increments of dx and dy
! Global latitude grid goes from south pole to north pole
! Global longitude grid starts at Dateline with western edge on Dateline

    dx = (mksrf_offline_edgee - mksrf_offline_edgew) / lsmlon