mksoitex.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

                           area(io,jo), mask_o      , novr_i2o, iovr_i2o, jovr_i2o, &
                           wovr_i2o)                             

! Process each cell on land grid:
! Find dominant 5 minute x 5 minute IGBP soil mapunit.
! landmask_i=0 means no soil data is available, so assume mapunit=0.
! Map from mapunit values to corresponding %sand and %clay 
! (silt not needed by land model and therefore, not calculated).
! Sum overlap weights by igbp soil mapunit
! landmask_i=0 means no soil data and landmask_i=1 means data exist

        do k = 0, mapunitmax
           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 = mapunit_i(ii,ji) * landmask_i(ii,ji) !mapunit (input grid)
           wst(k) = wst(k) + wovr_i2o(n)
        end do
        
! Rank non-zero weights by soil mapunit. 
! wsti(1) is the most extensive mapunit. 
! wsti(2) is the second most extensive mapunit.

        call mkrank (mapunitmax, wst, miss, wsti, num)

! Set soil texture as follows:
! If land grid cell is ocean or 100% glacier: cell has no soil 
! Otherwise, grid cell needs soil:
!   a. Use dominant igbp soil mapunit based on area of overlap unless
!     'no data' is dominant
!   b. In this case use second most dominant mapunit so long as it has data
!   c. If this has no data or if there isn't a second most dominant
!      mapunit, use loam for soil texture

        if (landmask(io,jo) == 0) then                        !ocean
           mapunit_o(io,jo) = 0.
           do l = 1, nlay
              sand_o(io,jo,l) = 0.
              clay_o(io,jo,l) = 0.
           end do
        else if (abs(pctgla_o(io,jo)-100.) < 1.e-06) then     !glacier
           mapunit_o(io,jo) = 0.
           do l = 1, nlay
              sand_o(io,jo,l) = 0.
              clay_o(io,jo,l) = 0.
           end do
        else                                                  !need soil
           if (wsti(1) /= 0) then                             !not 'no data'
              mapunit_o(io,jo) = wsti(1)
              do l = 1, nlay
                 sand_o(io,jo,l) = sand_i(wsti(1),l)
                 clay_o(io,jo,l) = clay_i(wsti(1),l)
              end do
           else                                               !if (wsti(1) == 0) then
              if (wsti(2) == 0 .or. wsti(2) == miss) then     !no data
                 mapunit_o(io,jo) = wsti(2)
                 do l = 1, nlay
                    sand_o(io,jo,l) = 43.                     !use loam
                    clay_o(io,jo,l) = 18.
                 end do
              else                                            !if (wsti(2) /= 0 and /= miss)
                 mapunit_o(io,jo) = wsti(2)
                 do l = 1, nlay
                    sand_o(io,jo,l) = sand_i(wsti(2),l)
                    clay_o(io,jo,l) = clay_i(wsti(2),l)
                 end do
              end if       !end of wsti(2) if-block
           end if          !end of wsti(1) if-block
        end if             !end of land/ocean if-block

! 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,*) 'MKSOITEX 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 type on input and output grids
! -----------------------------------------------------------------

! input grid: global areas by texture class

  gast_i(:) = 0.
  do l = 1, nlay
     do ji = 1, nlat_i
        do ii = 1, nlon_i
           mapunittemp = nint(mapunit_i(ii,ji))
           if (mapunittemp==0) then
              typ = 'no soil: ocean, glacier, lake, no data'
           else if (clay_i(mapunittemp,l) >= 40.) then
              typ = 'clays'
           else if (sand_i(mapunittemp,l) >= 50.) then
              typ = 'sands'
           else if (clay_i(mapunittemp,l)+sand_i(mapunittemp,l) < 50.) then
              if (landmask_i(ii,ji) /= 0.) then
                 typ = 'silts'
              else            !if (landmask_i(ii,ji) == 0.) then !no data
                 typ = 'no soil: ocean, glacier, lake, no data'
              end if
           else
              typ = 'loams'
           end if
           do m = 0, nlsm
              if (typ == soil(m)) go to 101
           end do
           write (6,*) 'MKSOITEX error: sand = ',sand_i(mapunittemp,l), &
             ' clay = ',clay_i(mapunittemp,l), &
             ' not assigned to soil type for input grid lon,lat,layer = ',ii,ji,l
           call endrun
101        continue
           gast_i(m) = gast_i(m) + area_i(ii,ji)
        end do
     end do
  end do

! output grid: global areas by texture class

  gast_o(:) = 0.
  do l = 1, nlay
     do jo = 1, lsmlat
        do io = 1, numlon(jo)
           if (clay_o(io,jo,l)==0. .and. sand_o(io,jo,l)==0.) then
              typ = 'no soil: ocean, glacier, lake, no data'
           else if (clay_o(io,jo,l) >= 40.) then
              typ = 'clays'
           else if (sand_o(io,jo,l) >= 50.) then
              typ = 'sands'
           else if (clay_o(io,jo,l)+sand_o(io,jo,l) < 50.) then
              typ = 'silts'
           else
              typ = 'loams'
           end if
           do m = 0, nlsm
              if (typ == soil(m)) go to 102
           end do
           write (6,*) 'MKSOITEX error: sand = ',sand_o(io,jo,l), &
             ' clay = ',clay_o(io,jo,l), &
             ' not assigned to soil type for output grid lon,lat,layer = ',io,jo,l
           call endrun
102        continue
           gast_o(m) = gast_o(m) + area(io,jo)
        end do
     end do
  end do

! Diagnostic output

  write (ndiag,*)
  write (ndiag,'(1x,70a1)') ('=',l=1,70)
  write (ndiag,*) 'Soil Texture Output'
  write (ndiag,'(1x,70a1)') ('=',l=1,70)
  write (ndiag,*)

  write (ndiag,*) 'The following table of soil texture classes is for comparison only.'
  write (ndiag,*) 'The actual data is continuous %sand, %silt and %clay not textural classes'
  write (ndiag,*)

  write (ndiag,*)
  write (ndiag,'(1x,70a1)') ('.',l=1,70)
  write (ndiag,1001)
1001 format (1x,'soil texture class',17x,' input grid area output grid area',/ &
             1x,33x,'     10**6 km**2','      10**6 km**2')
  write (ndiag,'(1x,70a1)') ('.',l=1,70)
  write (ndiag,*)

  do l = 0, nlsm
     write (ndiag,1002) soil(l),gast_i(l)*1.e-6,gast_o(l)*1.e-6
1002 format (1x,a38,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 %sand and %clay'
  write (6,*)
  call shr_sys_flush(6)

! Deallocate dynamic memory

  deallocate (mask_i)
  deallocate (latixy_i)
  deallocate (longxy_i)
  deallocate (numlon_i)
  deallocate (lon_i)
  deallocate (lon_i_offset)
  deallocate (lat_i)
  deallocate (area_i)
  deallocate (sand_i)
  deallocate (clay_i)
  deallocate (landmask_i)
  deallocate (mapunit_i)

  return
end subroutine mksoitex