mklai.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

     call wrap_get_vara_realx (ncid, varid , beg4d, len4d, mlai_i )

     call wrap_inq_varid (ncid, 'MONTHLY_SAI', varid)
     call wrap_get_vara_realx (ncid, varid , beg4d, len4d, msai_i )

     call wrap_inq_varid (ncid, 'MONTHLY_HEIGHT_TOP', varid)
     call wrap_get_vara_realx (ncid, varid, beg4d, len4d, mhgtt_i)

     call wrap_inq_varid (ncid, 'MONTHLY_HEIGHT_BOT', varid)
     call wrap_get_vara_realx (ncid, varid, beg4d, len4d, mhgtb_i)

! 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

     mlai_o(:,:,:)  = 0.
     msai_o(:,:,:)  = 0.
     mhgtt_o(:,:,:) = 0.
     mhgtb_o(:,:,:) = 0.

     mlai(:,:,:,m)  = 0.
     msai(:,:,:,m)  = 0.
     mhgtt(:,:,:,m) = 0.
     mhgtb(:,:,:,m) = 0.

!$OMP PARALLEL DO PRIVATE (io,jo,ii,ji,l,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)                             

           mask_o = 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
              mask_o = mask_o + landmask_i(ii,ji) * wovr_i2o(n)
           end do
              
           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   ) 

! Make area average and set oceans to zero

           if (landmask(io,jo) == 0) then
              do l = 0, numpft
                 mlai_o(io,jo,l)  = 0.
                 msai_o(io,jo,l)  = 0.
                 mhgtt_o(io,jo,l) = 0.
                 mhgtb_o(io,jo,l) = 0.
              end do
           else
              do l = 0, numpft
                 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
                    mlai_o(io,jo,l)  = mlai_o(io,jo,l)  + mlai_i(ii,ji,l)  * wovr_i2o(n)
                    msai_o(io,jo,l)  = msai_o(io,jo,l)  + msai_i(ii,ji,l)  * wovr_i2o(n)
                    mhgtt_o(io,jo,l) = mhgtt_o(io,jo,l) + mhgtt_i(ii,ji,l) * wovr_i2o(n)
                    mhgtb_o(io,jo,l) = mhgtb_o(io,jo,l) + mhgtb_i(ii,ji,l) * wovr_i2o(n)
                 end do
              end do
           endif  

! Assign lai/sai/hgtt/hgtb to the top [maxpatch_pft] PFTS
! as determined by mkpft.F

           do l = 0, numpft
              do n = 1, maxpatch_pft
                 if (l == pft(io,jo,n)) then
                    mlai(io,jo,n,m) = mlai_o(io,jo,l)
                    msai(io,jo,n,m) = msai_o(io,jo,l)
                    mhgtt(io,jo,n,m) = mhgtt_o(io,jo,l)
                    mhgtb(io,jo,n,m) = mhgtb_o(io,jo,l)
                 end if
              end do
           end do

! 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 for month m
! -----------------------------------------------------------------

! 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,*) 'MKLAI 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 global area

     glai_i(:)  = 0.
     gsai_i(:)  = 0.
     ghgtt_i(:) = 0.
     ghgtb_i(:) = 0.
     garea_i    = 0.

     do ji = 1, nlat_i
        do ii = 1, nlon_i
           garea_i = garea_i + area_i(ii,ji)
           do l = 0, numpft
              glai_i(l) = glai_i(l)+mlai_i(ii,ji,l)*area_i(ii,ji)
              gsai_i(l) = gsai_i(l)+msai_i(ii,ji,l)*area_i(ii,ji)
              ghgtt_i(l) = ghgtt_i(l)+mhgtt_i(ii,ji,l)*area_i(ii,ji)
              ghgtb_i(l) = ghgtb_i(l)+mhgtb_i(ii,ji,l)*area_i(ii,ji)
           end do
        end do
     end do

! Output grid global area

     glai_o(:)  = 0.
     gsai_o(:)  = 0.
     ghgtt_o(:) = 0.
     ghgtb_o(:) = 0.
     garea_o    = 0.

     do jo = 1, lsmlat
        do io = 1, numlon(jo)
           garea_o = garea_o + area(io,jo)
           do l = 0, numpft
              glai_o(l) = glai_o(l)+mlai_o(io,jo,l)*area(io,jo)
              gsai_o(l) = gsai_o(l)+msai_o(io,jo,l)*area(io,jo)
              ghgtt_o(l) = ghgtt_o(l)+mhgtt_o(io,jo,l)*area(io,jo)
              ghgtb_o(l) = ghgtb_o(l)+mhgtb_o(io,jo,l)*area(io,jo)
           end do
        end do
     end do

! Comparison

     write (ndiag,*)
     write (ndiag,'(1x,70a1)') ('=',k=1,70)
     write (ndiag,*) 'LAI Output for month ',m
     write (ndiag,'(1x,70a1)') ('=',k=1,70)

     write (ndiag,*)
     write (ndiag,'(1x,70a1)') ('.',k=1,70)
     write (ndiag,1001)
1001 format (1x,'PFT input grid area output grid area',/ &
             1x,3x,'     10**6 km**2','      10**6 km**2')
     write (ndiag,'(1x,70a1)') ('.',k=1,70)
     write (ndiag,*)

     do l = 0, numpft
        write (ndiag,1002) l, glai_i(l)*1.e-06*1.e-02,glai_o(l)*1.e-06*1.e-02
1002    format (1x,i3,f16.3,f17.3)
     end do

     write (6,*) 'Successfully made LAIs/SAIs/heights for month ', m
     write (6,*)
     call shr_sys_flush(6)

  end do  ! end loop over time

! Close input file

  call wrap_close(ncid)

! 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 (mlai_i)
  deallocate (msai_i)
  deallocate (mhgtt_i)
  deallocate (mhgtb_i)

  return
end subroutine mklai