mkgridmod.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

    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

#endif

end module mkgridMod