oas_defineMod.F90 Source File


Source Code

module oas_defineMod
  use mod_oasis
  implicit none
  save
  private

  public  :: oas_definitions_init
  integer :: ierror
contains

  subroutine oas_definitions_init(bounds)
    use spmdMod      , only : masterproc
    use domainMod    , only : ldomain
    use clm_varpar   , only : nlevsoi, nlevgrnd
    use decompMod    , only : ldecomp, bounds_type
    use oas_vardefMod

    type(bounds_type) , intent(in)  :: bounds ! start and end gridcell indices for this MPI task

    ! oasis_def_partition
    integer, allocatable :: partition(:)      ! partition descriptor; input to oasis_def_partition
    integer              :: gcell_start       ! starting gridcell index
    integer              :: gcell_previous    ! gridcell index from previous loop iteration
    integer              :: num_local_gcells  ! total number of gridcells for this process
    integer              :: k, g              ! array/loop indices
    integer              :: grid_id           ! id returned after call to oasis_def_partition
#ifdef COUP_OAS_ICON
    integer              :: jg                ! loop counter
#endif

    ! oasis_def_var
    integer              :: var_nodims(2)     ! var dimension parameters


    ! TODO: Fix grids.nc and masks.nc generation for 1D grid representation
    !if (masterproc) then
    ! call define_grid()
    !end if

    ! -----------------------------------------------------------------
    ! ... Define partition
    ! -----------------------------------------------------------------
    if (ldomain%nj == 1) then
      ! POINTS partitioning scheme
      ! partition length = 2 + total number of gridpoints allocated to this MPI process
      num_local_gcells = bounds%endg - bounds%begg + 1
      allocate(partition(num_local_gcells + 2))
      partition(:) = 0; k = 0

      ! Use POINTS partitioning scheme. This partition is a list of global indices associated with each process.
      partition(1) = 4
      partition(2) = num_local_gcells

      do g = bounds%begg, bounds%endg
        partition(3+k) = ldecomp%gdc2glo(g)
        k = k + 1
      enddo
    else
      ! ORANGE partitioning scheme
      ! partition length = (# elements for partition info) + (max segments ORANGE partition) x (# elements per segment info)
      !                  = 2 + 200*2 = 402
      allocate(partition(402))
      partition(:) = 0; k = 0

      ! Use ORANGE partitioning scheme. This scheme defines an ensemble
      ! of gridcell segments. See OASIS3-MCT User's guide for more info.
      partition(1) = 3

      ! Mark 1st segment
      gcell_start = ldecomp%gdc2glo(bounds%begg)
      partition(2) = 1
      gcell_previous = gcell_start

      ! Capture segments by detecting segment boundaries. A boundary is 
      ! detected when the current and previous gridcells are not consecutive.
      do g = bounds%begg+1, bounds%endg
        if (ldecomp%gdc2glo(g) - gcell_previous /= 1) then
          ! Previous segment complete; its partition params could now be defined
          partition(3+k) = gcell_start - 1                  ! segment global offset (0-based)
          partition(4+k) = gcell_previous - gcell_start + 1 ! segment length
          k = k + 2
          gcell_start  = ldecomp%gdc2glo(g) ! current gridcell marks the start of a new segment 
          partition(2) = partition(2) + 1   ! increment number of segments (limited to 200 based from OASIS3-MCT User's guide)
        end if
        gcell_previous = ldecomp%gdc2glo(g)
      enddo

      ! Define partition params for last segment
      partition(3+k) = gcell_start - 1
      partition(4+k) = gcell_previous - gcell_start + 1
      partition(2) = partition(2) + 1
    end if

    call oasis_def_partition(grid_id, partition, ierror, ldomain%ns)
    deallocate(partition)

    ! -----------------------------------------------------------------
    ! ... Define coupling fields
    ! -----------------------------------------------------------------

#ifdef COUP_OAS_PFL
    var_nodims(1) = 1         ! unused
    var_nodims(2) = nlevsoi   ! number of fields in a bundle

    call oasis_def_var(oas_et_loss_id, "ECLM_ET", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) 
    
    var_nodims(2) = nlevgrnd         ! number of fields in a bundle
    call oasis_def_var(oas_sat_id, "ECLM_SOILLIQ", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror)
    call oasis_def_var(oas_psi_id, "ECLM_PSI", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror)
#endif 

#ifdef COUP_OAS_ICON

    var_nodims(1) = 1         ! unused
    var_nodims(2) = 1         ! number of fields in a bundle

    ! send to ICON
    CALL oasis_def_var(oas_id_t,  "CLMTEMPE", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) ! 1 
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMTEMPE.')
    CALL oasis_def_var(oas_id_u,  "CLMUWIND", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) ! 2
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMUWIND.')
    CALL oasis_def_var(oas_id_v,  "CLMVWIND", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) ! 3
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMVWIND.')
    CALL oasis_def_var(oas_id_qv, "CLMSPWAT", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) ! 4
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMSPWAT.')
    CALL oasis_def_var(oas_id_ht, "CLMTHICK", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) ! 5
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMTHICK.')
    CALL oasis_def_var(oas_id_pr, "CLMPRESS", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) ! 6
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMPRESS.')
    CALL oasis_def_var(oas_id_rs, "CLMDIRSW", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) ! 7
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMDIRSW.')
    CALL oasis_def_var(oas_id_fs, "CLMDIFSW", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) ! 8
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMDIFSW.')
    CALL oasis_def_var(oas_id_lw, "CLMLONGW", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) ! 9
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMLONGW.')
    CALL oasis_def_var(oas_id_cr, "CLMCVPRE", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) !10
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMCVPRE.')
    CALL oasis_def_var(oas_id_gr, "CLMGSPRE", grid_id, var_nodims, OASIS_In, OASIS_Real, ierror) !11
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMGSPRE.')

    ! receive from ICON
    CALL oasis_def_var(oas_id_it, "CLMINFRA", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) !12
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMINFRA.')
    CALL oasis_def_var(oas_id_ad, "CLMALBED", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) !13
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMALBED.')
    CALL oasis_def_var(oas_id_ai, "CLMALBEI", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) !14
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMALBEI.')
    CALL oasis_def_var(oas_id_tx, "CLMTAUX" , grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) !15
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMTAUX.')
    CALL oasis_def_var(oas_id_ty, "CLMTAUY" , grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) !16
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMTAUY.')
    CALL oasis_def_var(oas_id_sh, "CLMSHFLX", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) !17
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMSHFLX.')
    CALL oasis_def_var(oas_id_lh, "CLMLHFLX", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) !18
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMLHFLX.')
    CALL oasis_def_var(oas_id_st, "CLMTGRND", grid_id, var_nodims, OASIS_Out, OASIS_Real, ierror) !19
    IF (ierror /= 0) CALL oasis_abort(oas_comp_id, oas_comp_name, 'Failure in oasis_def_var for CLMTGRND.')

#endif

    ! End definition phase
    call oasis_enddef(ierror)
  end subroutine oas_definitions_init

  subroutine define_grid()
    use shr_kind_mod , only : r8 => shr_kind_r8
    use domainMod    , only : ldomain

    character(len=4), parameter   :: grid_name='gclm'
    integer,          parameter   :: SOUTH = 1
    integer,          parameter   :: NORTH = 2
    integer,          parameter   :: WEST  = 1
    integer,          parameter   :: EAST  = 2   
    real(kind=r8),    allocatable :: corner_lon(:,:), corner_lat(:,:)
    real(kind=r8),    allocatable :: oas_lon(:,:), oas_lat(:,:)
    real(kind=r8),    allocatable :: oas_corner_lon(:,:,:), oas_corner_lat(:,:,:)
    integer,          allocatable :: oas_mask(:,:)
    real(kind=r8)                 :: center, offset, neighbor
    integer                       :: write_grid_files
    integer                       :: i, j
    integer                       :: n_lons, n_lats, n_gridcells

    call oasis_start_grids_writing(write_grid_files)

    if (write_grid_files == 1) then
      n_lons = ldomain%ni
      n_lats = ldomain%nj
      n_gridcells = ldomain%ns

      ! -----------------------------------------------------------------
      ! ... Define centers
      ! -----------------------------------------------------------------
      allocate(oas_lon(n_gridcells, 1))
      allocate(oas_lat(n_gridcells, 1))
      do j = 1, n_lats
        i = (j-1)*n_lons + 1
        oas_lon(i:n_lons, 1) = ldomain%lonc(:)
        oas_lat(i:n_lons, 1) = ldomain%latc(j)
      enddo
      call oasis_write_grid(grid_name, n_gridcells, 1, oas_lon, oas_lat)

      ! -----------------------------------------------------------------
      ! ... Define corners
      ! -----------------------------------------------------------------
      allocate(corner_lon(n_lons, WEST:EAST))
      allocate(corner_lat(n_lats, SOUTH:NORTH))
      allocate(oas_corner_lon(n_gridcells, 1, 4))
      allocate(oas_corner_lat(n_gridcells, 1, 4))

      do j = 1, n_lats
        if (j == 1) then
          neighbor = ldomain%latc(j+1)
        else
          neighbor = ldomain%latc(j-1)
        end if
        center = ldomain%latc(j)
        offset = abs(center - neighbor) / 2.0_r8
        corner_lat(j, SOUTH:NORTH) = [center-offset, center+offset]
      enddo
    
      do i = 1, n_lons
        if (i == 1) then
          neighbor = ldomain%lonc(i+1)
        else
          neighbor = ldomain%lonc(i-1)
        end if
        center = ldomain%lonc(i)
        offset = abs(center - neighbor) / 2.0_r8
        corner_lon(i, WEST:EAST) = [center-offset, center+offset]
      enddo

      !  oas_corner indexing scheme                       
      !       4 +-----+ 3  NORTH    
      !         |     |             
      !         |     |             
      !       1 +-----+ 2  SOUTH   
      !       WEST   EAST                                    
      do j = 1, n_lats
        i = (j-1)*n_lons + 1
        oas_corner_lon(i:n_lons,1,1:2) = corner_lon(:, WEST:EAST) ! bottom side
        oas_corner_lat(i:n_lons,1,2)   = corner_lat(j, SOUTH)     ! right side
        oas_corner_lat(i:n_lons,1,3)   = corner_lat(j, NORTH)     ! right side
      enddo
    
      ! Fill missing longitudes and latitudes
      oas_corner_lon(:,1,[4,3]) = oas_corner_lon(:,1,[1,2]) ! West-East for top side
      oas_corner_lat(:,1,[1,4]) = oas_corner_lat(:,1,[2,3]) ! South-North for left side
      call oasis_write_corner(grid_name, n_gridcells, 1, 4, oas_corner_lon, oas_corner_lat)

      ! -----------------------------------------------------------------
      ! ... Define mask
      ! -----------------------------------------------------------------
      ! CLM5 landmask convention: 0 = ocean, 1 = land
      allocate(oas_mask(n_gridcells, 1))
      do j = 1, n_lats
        i = (j-1)*n_lons + 1
        oas_mask(i:n_lons,1) = ldomain%mask(:)
      enddo

      ! Invert mask to conform to OASIS convention: 0 = not masked, 1 = masked
      where (oas_mask(:,1) == 0)
        oas_mask(:,1) = 1
      else where
        oas_mask(:,1) = 0
      end where
      call oasis_write_mask(grid_name, n_gridcells, 1, oas_mask)
      call oasis_terminate_grids_writing()
      deallocate(oas_lon, oas_lat, oas_corner_lon, oas_corner_lat, oas_mask, corner_lon, corner_lat)
    end if

  end subroutine define_grid
end module oas_defineMod