mkpft.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

           pctlnd_o(io,jo) = pctlnd_o(io,jo) + 100.*landmask_i(ii,ji)*wovr_i2o(n)
        end do
        mask_o = pctlnd_o(io,jo)/100.

        call areaini_point (io        , jo          , nlon_i  , nlat_i  , numlon_i  , &
                           lon_i      , lon_i_offset, lat_i   , area_i  , landmask_i, &
                           lsmlon     , lsmlat      , numlon  , lonw    , lats      , &
                           area(io,jo), mask_o      , novr_i2o, iovr_i2o, jovr_i2o  , &
                           wovr_i2o   ) 


! Area-average percent cover on input grid [pctpft_i] to output grid [pctpft_o] 
! and correct [pctpft_o] according to land landmask
! Note that percent cover is in terms of total grid area. 

        wst_sum = 0.
        do m = 0, numpft
           pctpft_o(io,jo,m) = 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
              pctpft_o(io,jo,m) = pctpft_o(io,jo,m) + pctpft_i(ii,ji,m)*wovr_i2o(n)
           end do
           if (pctlnd_o(io,jo) > 0.) then
              if (landmask(io,jo) == 0) pctpft_o(io,jo,m) = 0.
           else
              pctpft_o(io,jo,m) = 0.
              if (landmask(io,jo) == 1) pctpft_o(io,jo,0) = 100.
           end if
           wst_sum = wst_sum + pctpft_o(io,jo,m) 
        end do

! Error check: percents should sum to 100 for land grid cells

        if (landmask(io,jo) == 1) then
           if (abs(wst_sum-100.) > 0.000001) then
              write (6,*) 'MKPFT error: pft = ', &
                   (pctpft_o(io,jo,m), m = 0, numpft), &
                   ' do not sum to 100. at column, row = ',io,jo, &
                   ' but to ', wst_sum
              call endrun
           end if
        end if

! Find the output pft and pct arrays
! Save percent cover by PFT [wst] and total percent cover [wst_sum]

        wst_sum = 0. 
        sumpct = 0
        do m = 0, numpft
           wst(m) = pctpft_o(io,jo,m)       
           wst_sum = wst_sum + pctpft_o(io,jo,m) 
        end do

! Rank [wst] in ascending order to obtain the top [maxpatch_pft] PFTs 

        if (landmask(io,jo) == 1) then
           call mkrank (numpft, wst, miss, wsti, maxpatch_pft)
        end if

! Fill in [pft] and [pctpft] with data for top [maxpatch_pft] PFTs. 
! If land model grid cell is ocean, set to no PFTs. 
! If land model grid cell is land, there are three possibilities: 
!  1. If [pctlnd_o] = 0, there is no PFT data from the input grid. 
!     Since need land data, use bare ground.
!  2. If [pctlnd_o] > 0, there is some PFT data from the input grid but:
!     a. use the chosen PFT so long as it is not a missing value
!     b. missing value means no more PFTs with cover > 0

        if (landmask(io,jo) == 1) then       !model grid wants land
           do m = 1, maxpatch_pft
              if (wsti(m) .ne. miss) then
                 pft(io,jo,m) = wsti(m)          
                 pctpft(io,jo,m) = wst(wsti(m))
              else
                 pft(io,jo,m) = noveg
                 pctpft(io,jo,m) = 0.
              end if
              sumpct = sumpct + pctpft(io,jo,m)
           end do
        else                                     !model grid wants ocean
           do m = 1, maxpatch_pft
              pft(io,jo,m) = 0
              pctpft(io,jo,m) = 0.
           end do
        end if

! Correct for the case of more than [maxpatch_pft] PFTs present 

        if (sumpct < wst_sum) then
           diff  = wst_sum - sumpct 
           sumpct = 0.
           do m = 1, maxpatch_pft
              pctpft(io,jo,m) = pctpft(io,jo,m) + diff/maxpatch_pft
              sumpct = sumpct + pctpft(io,jo,m)
           end do
        end if

! Error check: make sure have a valid PFT

        do m = 1, maxpatch_pft
           if (pft(io,jo,m)<0 .or. pft(io,jo,m)>numpft) then
              write (6,*) 'MKPFT error: invalid PFT for i,j=',io,jo,pft(io,jo,m)
              call endrun
           end if
        end do

! Error check: make sure PFTs sum to 100% cover

        if (landmask(io,jo) == 1) then    
           if (abs(sumpct - 100.) > 0.000001) then
              write(6,*) 'MKPFT error: sum(pct) over maxpatch_pft'
              write(6,*) '             is not = 100.'
              write(6,*) sumpct, io,jo
              call endrun
           end if
           if (sumpct < -0.000001) then
              write(6,*) 'MKPFT error: sum(pct) over maxpatch_pft'
              write(6,*) '             is < 0.'
              write(6,*) sumpct, io,jo
              call endrun
           end if
        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) * landmask_i(ii,ji)
           fld_o(io,jo) = fld_o(io,jo) + wovr_i2o(n) * fld_i * mask_o
        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) * landmask_i(ii,ji)
        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 areas on input and output grids
! -----------------------------------------------------------------

! input grid

  gpft_i(:) = 0.
  garea_i = 0.
  do ji = 1, nlat_i
     do ii = 1, nlon_i
        garea_i = garea_i + area_i(ii,ji)
        do m = 0, numpft
           gpft_i(m) = gpft_i(m) + pctpft_i(ii,ji,m)*area_i(ii,ji)
        end do
     end do
  end do

! output grid

  gpft_o(:) = 0.
  garea_o = 0.
  do jo = 1, lsmlat
     do io = 1, numlon(jo)
        garea_o = garea_o + area(io,jo)
        do m = 0, numpft
           gpft_o(m) = gpft_o(m) + pctpft_o(io,jo,m)*area(io,jo)
        end do
     end do
  end do

! comparison

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

  write (ndiag,*)
  write (ndiag,'(1x,70a1)') ('.',k=1,70)
  write (ndiag,1001)
1001 format (1x,'plant 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 m = 0, numpft
     write (ndiag,1002) veg(m), gpft_i(m)*1.e-06/100.,gpft_o(m)*1.e-06/100.
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 PFTs'
  write (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 (landmask_i)
  deallocate (pctpft_i)

  return
end subroutine mkpft