mksoicol.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

subroutine mksoicol (fsoicol, ndiag, pctgla_o, soil_color_o)

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Make soil color classes for model grid from BATS T42 data
! 
! Method: 
! 
! Author: Gordon Bonan
! 
!-----------------------------------------------------------------------
! $Id: mksoicol.F90,v 1.4.10.2 2002/04/27 15:38:55 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 shr_sys_mod, only : shr_sys_flush 
  implicit none

! ------------------------ arguments ------------------------------
  character(len=*), intent(in) :: fsoicol              !input soicol dataset file name
  integer , intent(in) :: ndiag                        !unit number for diagnostic output
  real(r8), intent(in) :: pctgla_o(lsmlon,lsmlat)      !output grid: percent glacier
  integer , intent(out):: soil_color_o(lsmlon,lsmlat)  !output grid: soil color classes
! -----------------------------------------------------------------

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

  integer :: nlon_i                           !input grid : longitude points (read in)
  integer :: nlat_i                           !input grid : latitude  points (read in)
  integer :: ncid,dimid,varid                 !input netCDF id's
  integer :: ier                              !error status   

  integer :: ii                               !longitude index for BATS grid
  integer :: io                               !longitude index for model grid
  integer :: ji                               !latitude  index for BATS grid
  integer :: jo                               !latitude  index for model grid
  integer :: k                                !temporary BATS or model soil color
  integer :: miss = 99999                     !missing data indicator
  integer :: n                                !loop index
  integer, parameter :: num=2                 !get 1st and 2nd largest areas of overlap
  integer :: wsti(num)                        !index to 1st and 2nd largest values in wst

  integer, parameter :: nsoicol=8             !number of model color classes
  real(r8) :: wst(0:nsoicol)                  !overlap weights, by surface type
  real(r8) :: gast_i(0:nsoicol)               !input grid : global area, by surface type
  real(r8) :: gast_o(0:nsoicol)               !output grid: global area, by surface type
  character(len=35) :: col(0:nsoicol)         !name of each color

  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)
  integer , allocatable :: soil_color_i(:,:)  !input grid: BATS soil color

  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
! -----------------------------------------------------------------

  write (6,*) 'Attempting to make soil color classes .....'
  call shr_sys_flush(6)

! -----------------------------------------------------------------
! Define the model color classes: 0 to nsoicol
! -----------------------------------------------------------------

  col(0) = 'no soil                            '
  col(1) = 'class 1: light                     '
  col(2) = 'class 2:                           '
  col(3) = 'class 3:                           '
  col(4) = 'class 4:                           '
  col(5) = 'class 5:                           '
  col(6) = 'class 6:                           '
  col(7) = 'class 7:                           '
  col(8) = 'class 8: dark                      '
! col(9) = 'class 9: very light North Africa   '

! -----------------------------------------------------------------
! Read input soil colors
! -----------------------------------------------------------------

! Obtain input grid info

  call getfil (fsoicol, 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 (soil_color_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, '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))

! Obtain input data

  call wrap_inq_varid (ncid, 'SOIL_COLOR', varid)
  call wrap_get_var_int (ncid, varid, soil_color_i)

  call wrap_close(ncid)

! -----------------------------------------------------------------
! Map data from input grid to land model grid. Get:
! -----------------------------------------------------------------

! 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,k,mask_o,novr_i2o,iovr_i2o,jovr_i2o,wovr_i2o,fld_i, &
!$OMP & wst, wsti)
  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)                             

! Sum overlap weights by color class

        do k = 0, nsoicol
           wst(k) = 0.
        end do
        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
           k = soil_color_i(ii,ji) !color class (input grid)
           wst(k) = wst(k) + wovr_i2o(n)
        end do

! Rank non-zero weights by color type. wsti(1) is the most extensive
! color type. wsti(2) is the second most extensive color type

        call mkrank (nsoicol, wst, miss, wsti, num)
        soil_color_o(io,jo) = wsti(1)

! If land but no color, set color to 4

        if (landmask(io,jo)==1 .and. soil_color_o(io,jo)==0) soil_color_o(io,jo) = 4

! Set ocean colors to zero

        if (landmask(io,jo) == 0) soil_color_o(io,jo) = 0

! Set color for grid cells that are 100% glacier to zero. Otherwise,
! must have a soil color for the non-glacier portion of grid cell.

        if (abs(pctgla_o(io,jo)-100.)<1.e-06) soil_color_o(io,jo)=0

! Error checks

        if (soil_color_o(io,jo) < 0 .or. soil_color_o(io,jo) > nsoicol) then
           write (6,*) 'MKSOICOL error: land model soil color = ', &
                soil_color_o(io,jo),' is not valid 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,*) 'MKGLACIER 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

! -----------------------------------------------------------------
! Error check2
! Compare global area of each soil color on input and output grids
! -----------------------------------------------------------------

! input grid

  gast_i(:) = 0.
  do ji = 1, nlat_i
     do ii = 1, nlon_i
        k = soil_color_i(ii,ji)
        gast_i(k) = gast_i(k) + area_i(ii,ji)
     end do
  end do

! output grid

  gast_o(:) = 0.
  do jo = 1, lsmlat
     do io = 1, numlon(jo)
        k = soil_color_o(io,jo)
        gast_o(k) = gast_o(k) + area(io,jo)
     end do
  end do

! area comparison

  write (ndiag,*)
  write (ndiag,'(1x,70a1)') ('=',k=1,70)
  write (ndiag,*) 'Soil Color Output'
  write (ndiag,'(1x,70a1)') ('=',k=1,70)

  write (ndiag,*)
  write (ndiag,'(1x,70a1)') ('.',k=1,70)
  write (ndiag,1001)
1001 format (1x,'soil color type',20x,' input grid area output grid area',/ &
             1x,33x,'     10**6 km**2','      10**6 km**2')
  write (ndiag,'(1x,70a1)') ('.',k=1,70)
  write (ndiag,*)

  do k = 0, nsoicol
     write (ndiag,1002) col(k),gast_i(k)*1.e-6,gast_o(k)*1.e-6
1002 format (1x,a35,f16.3,f17.3)
  end do

  if (lsmlat > 1) then
     k = lsmlat/2
     write (ndiag,*)
     write (ndiag,*) 'For reference the area on the output grid of a cell near the equator is: '
     write (ndiag,'(f10.3,a14)')area(1,k)*1.e-06,' x 10**6 km**2'
     write (ndiag,*)
  endif
     
  write (6,*) 'Successfully made soil color classes'
  write (6,*)
  call shr_sys_flush(6)


! Deallocate dynamic memory

  deallocate (latixy_i)
  deallocate (longxy_i)
  deallocate (numlon_i)
  deallocate (lon_i)
  deallocate (lon_i_offset)
  deallocate (lat_i)
  deallocate (area_i)
  deallocate (mask_i)
  deallocate (soil_color_i)

  return
end subroutine mksoicol