Add DLND to moab coupler

components/data_comps/dlnd/src/dlnd_comp_mod.F90

@@ -29,6 +29,10 @@ module dlnd_comp_mod
   use dlnd_shr_mod   , only: domain_fracname ! namelist input
   use dlnd_shr_mod   , only: nullstr
 
+#ifdef HAVE_MOAB
+  use seq_comm_mct,     only : mlnid ! id of moab lnd app
+  use iso_c_binding
+#endif
   ! !PUBLIC TYPES:
   implicit none
   private ! except
@@ -100,6 +104,12 @@ subroutine dlnd_comp_init(Eclock, x2l, l2x, &
        scmMode, scmlat, scmlon)
 
     ! !DESCRIPTION: initialize dlnd model
+#ifdef HAVE_MOAB
+       use iMOAB, only: iMOAB_DefineTagStorage, iMOAB_GetDoubleTagStorage, &
+                        iMOAB_SetIntTagStorage, iMOAB_SetDoubleTagStorage, &
+                        iMOAB_ResolveSharedEntities, iMOAB_CreateVertices, &
+                        iMOAB_GetMeshInfo, iMOAB_UpdateMeshInfo, iMOAB_WriteMesh
+#endif
     implicit none
 
     ! !INPUT/OUTPUT PARAMETERS:
@@ -135,6 +145,18 @@ subroutine dlnd_comp_init(Eclock, x2l, l2x, &
     character(nec_len) :: nec_str   ! elevation class, as character string
     character(*), parameter :: domain_fracname_unset = 'null'
 
+#ifdef HAVE_MOAB
+    character*400  tagname
+    real(R8) latv, lonv
+    integer iv, tagindex, ilat, ilon  !, arrsize, nfields ! ierr is already defined as integer above
+    real(R8), allocatable, target :: data(:)
+    integer(IN), pointer :: idata(:)   ! temporary
+    real(r8), dimension(:), allocatable :: moab_vert_coords  ! temporary
+#ifdef MOABDEBUG
+    character*100 outfile, wopts
+#endif
+#endif
+
     !--- formats ---
     character(*), parameter :: F00   = "('(dlnd_comp_init) ',8a)"
     character(*), parameter :: F0L   = "('(dlnd_comp_init) ',a, l2)"
@@ -256,6 +278,119 @@ subroutine dlnd_comp_init(Eclock, x2l, l2x, &
 
     call t_stopf('dlnd_initmctdom')
 
+#ifdef HAVE_MOAB
+    ilat = mct_aVect_indexRA(ggrid%data,'lat')
+    ilon = mct_aVect_indexRA(ggrid%data,'lon')
+    allocate(moab_vert_coords(lsize*3))
+    do iv = 1, lsize
+       lonv = ggrid%data%rAttr(ilon, iv) * SHR_CONST_PI/180.
+       latv = ggrid%data%rAttr(ilat, iv) * SHR_CONST_PI/180.
+       moab_vert_coords(3*iv-2)=COS(latv)*COS(lonv)
+       moab_vert_coords(3*iv-1)=COS(latv)*SIN(lonv)
+       moab_vert_coords(3*iv  )=SIN(latv)
+    enddo
+
+    ! create the vertices with coordinates from MCT domain
+    ierr = iMOAB_CreateVertices(mlnid, lsize*3, 3, moab_vert_coords)
+    if (ierr .ne. 0)  &
+       call shr_sys_abort('Error: fail to create MOAB vertices in data lnd model')
+
+    tagname='GLOBAL_ID'//C_NULL_CHAR
+    ierr = iMOAB_DefineTagStorage(mlnid, tagname, &
+                                  0, & ! dense, integer
+                                  1, & ! number of components
+                                  tagindex )
+    if (ierr .ne. 0)  &
+       call shr_sys_abort('Error: fail to retrieve GLOBAL_ID tag ')
+
+    ! get list of global IDs for Dofs
+    call mct_gsMap_orderedPoints(gsMap, my_task, idata)
+
+    ierr = iMOAB_SetIntTagStorage ( mlnid, tagname, lsize, &
+                                     0, & ! vertex type
+                                     idata)
+    if (ierr .ne. 0)  &
+       call shr_sys_abort('Error: fail to set GLOBAL_ID tag ')
+
+    ierr = iMOAB_ResolveSharedEntities( mlnid, lsize, idata );
+    if (ierr .ne. 0)  &
+       call shr_sys_abort('Error: fail to resolve shared entities')
+
+    deallocate(moab_vert_coords)
+    deallocate(idata)
+
+    ierr = iMOAB_UpdateMeshInfo( mlnid )
+    if (ierr .ne. 0)  &
+       call shr_sys_abort('Error: fail to update mesh info ')
+
+    allocate(data(lsize))
+    ierr = iMOAB_DefineTagStorage( mlnid, "area:aream:frac:mask"//C_NULL_CHAR, &
+                                     1, & ! dense, double
+                                     1, & ! number of components
+                                     tagindex )
+    if (ierr > 0 )  &
+       call shr_sys_abort('Error: fail to create tag: area:aream:frac:mask' )
+
+    data(:) = ggrid%data%rAttr(mct_aVect_indexRA(ggrid%data,'area'),:)
+    tagname='area'//C_NULL_CHAR
+    ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsize, &
+                                     0, & ! set data on vertices
+                                     data)
+    if (ierr > 0 )  &
+       call shr_sys_abort('Error: fail to get area tag ')
+
+    ! set the same data for aream (model area) as area
+    ! data(:) = ggrid%data%rAttr(mct_aVect_indexRA(ggrid%data,'aream'),:)
+    tagname='aream'//C_NULL_CHAR
+    ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsize, &
+                                     0, & ! set data on vertices
+                                     data)
+    if (ierr > 0 )  &
+       call shr_sys_abort('Error: fail to set aream tag ')
+
+    data(:) = ggrid%data%rAttr(mct_aVect_indexRA(ggrid%data,'mask'),:)
+    tagname='mask'//C_NULL_CHAR
+    ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsize, &
+                                     0, & ! set data on vertices
+                                     data)
+    if (ierr > 0 )  &
+       call shr_sys_abort('Error: fail to set mask tag ')
+
+    data(:) = ggrid%data%rAttr(mct_aVect_indexRA(ggrid%data,'frac'),:)
+    tagname='frac'//C_NULL_CHAR
+    ierr = iMOAB_SetDoubleTagStorage ( mlnid, tagname, lsize, &
+                                     0, & ! set data on vertices
+                                     data)
+    if (ierr > 0 )  &
+       call shr_sys_abort('Error: fail to set frac tag ')
+
+    deallocate(data)
+
+    ! define tags
+    ierr = iMOAB_DefineTagStorage( mlnid, trim(seq_flds_x2l_fields)//C_NULL_CHAR, &
+                                     1, & ! dense, double
+                                     1, & ! number of components
+                                     tagindex )
+    if (ierr > 0 )  &
+       call shr_sys_abort('Error: fail to create seq_flds_x2l_fields tags ')
+
+    ierr = iMOAB_DefineTagStorage( mlnid, trim(seq_flds_l2x_fields)//C_NULL_CHAR, &
+                                     1, & ! dense, double
+                                     1, & ! number of components
+                                     tagindex )
+    if (ierr > 0 )  &
+       call shr_sys_abort('Error: fail to create seq_flds_l2x_fields tags ')
+#ifdef MOABDEBUG
+       !      debug test
+    outfile = 'LndDataMesh.h5m'//C_NULL_CHAR
+    wopts   = ';PARALLEL=WRITE_PART'//C_NULL_CHAR !
+       !      write out the mesh file to disk
+    ierr = iMOAB_WriteMesh(mlnid, trim(outfile), trim(wopts))
+    if (ierr .ne. 0) then
+       call shr_sys_abort(subname//' ERROR in writing data mesh lnd ')
+    endif
+#endif
+#endif
     !----------------------------------------------------------------------------
     ! Initialize MCT attribute vectors
     !----------------------------------------------------------------------------
@@ -339,8 +474,15 @@ subroutine dlnd_comp_run(EClock, x2l, l2x, &
        inst_suffix, logunit, case_name)
 
     ! !DESCRIPTION:  run method for dlnd model
-    implicit none
+#ifdef HAVE_MOAB
+#ifdef MOABDEBUG
+    use iMOAB, only: iMOAB_WriteMesh
+#endif
+    use seq_flds_mod    , only: seq_flds_l2x_fields 
+    use seq_flds_mod    , only: moab_set_tag_from_av
+#endif
 
+    implicit none
     ! !INPUT/OUTPUT PARAMETERS:
     type(ESMF_Clock)       , intent(in)    :: EClock
     type(mct_aVect)        , intent(inout) :: x2l
@@ -366,6 +508,17 @@ subroutine dlnd_comp_run(EClock, x2l, l2x, &
     integer(IN)   :: nu                    ! unit number
     logical       :: write_restart         ! restart now
     character(len=18) :: date_str
+#ifdef HAVE_MOAB
+    real(R8), allocatable, target :: datam(:)
+    type(mct_list) :: temp_list
+    integer :: size_list, index_list, lsize
+    type(mct_string)    :: mctOStr  !
+    character*400  tagname, mct_field
+#ifdef MOABDEBUG
+    integer  :: cur_dlnd_stepno, ierr
+    character*100 outfile, wopts, lnum
+#endif
+#endif
 
     character(*), parameter :: F00   = "('(dlnd_comp_run) ',8a)"
     character(*), parameter :: F04   = "('(dlnd_comp_run) ',2a,2i8,'s')"
@@ -464,6 +617,32 @@ subroutine dlnd_comp_run(EClock, x2l, l2x, &
 
     call t_stopf('DLND_RUN')
 
+#ifdef HAVE_MOAB
+    lsize = mct_avect_lsize(l2x) ! is it the same as mct_avect_lsize(avstrm) ?
+    allocate(datam(lsize)) ! 
+    call mct_list_init(temp_list ,seq_flds_l2x_fields)
+    size_list=mct_list_nitem (temp_list)
+    do index_list = 1, size_list
+      call mct_list_get(mctOStr,index_list,temp_list)
+      mct_field = mct_string_toChar(mctOStr)
+      tagname= trim(mct_field)//C_NULL_CHAR
+      call moab_set_tag_from_av(tagname, l2x, index_list, mlnid, datam, lsize) ! loop over all a2x fields, not just a few
+    enddo
+    call mct_list_clean(temp_list)
+    deallocate(datam) ! maybe we should keep it around, deallocate at the final only?
+
+#ifdef MOABDEBUG
+    call seq_timemgr_EClockGetData( EClock, stepno=cur_dlnd_stepno )
+    write(lnum,"(I0.2)")cur_dlnd_stepno
+    outfile = 'dlnd_comp_run_'//trim(lnum)//'.h5m'//C_NULL_CHAR
+    wopts   = 'PARALLEL=WRITE_PART'//C_NULL_CHAR
+    ierr = iMOAB_WriteMesh(mlnid, outfile, wopts)
+    if (ierr > 0 )  then
+       write(logunit,*) 'Failed to write data lnd component state '
+    endif
+#endif
+#endif
+
   end subroutine dlnd_comp_run
 
   !===============================================================================

components/data_comps/dlnd/src/lnd_comp_mct.F90

@@ -16,7 +16,11 @@ module lnd_comp_mct
   use dlnd_comp_mod   , only: dlnd_comp_init, dlnd_comp_run, dlnd_comp_final
   use dlnd_shr_mod    , only: dlnd_shr_read_namelists
   use seq_flds_mod    , only: seq_flds_x2l_fields, seq_flds_l2x_fields
-
+#ifdef HAVE_MOAB
+  use seq_comm_mct, only : mlnid !            iMOAB app id for lnd
+  use iso_c_binding
+  use iMOAB           , only: iMOAB_RegisterApplication
+#endif
   ! !PUBLIC TYPES:
   implicit none
   private ! except
@@ -52,7 +56,9 @@ module lnd_comp_mct
 
   !===============================================================================
   subroutine lnd_init_mct( EClock, cdata, x2l, l2x, NLFilename )
-
+#ifdef HAVE_MOAB
+    use shr_stream_mod, only: shr_stream_getDomainInfo, shr_stream_getFile
+#endif
     ! !DESCRIPTION:  initialize dlnd model
     implicit none
 
@@ -146,13 +152,25 @@ subroutine lnd_init_mct( EClock, cdata, x2l, l2x, NLFilename )
     !----------------------------------------------------------------------------
     ! Initialize dlnd
     !----------------------------------------------------------------------------
-
+#ifdef HAVE_MOAB
+    ierr = iMOAB_RegisterApplication(trim("DLND")//C_NULL_CHAR, mpicom, compid, mlnid)
+    if (ierr .ne. 0) then
+      write(logunit,*) subname,' error in registering data lnd comp'
+      call shr_sys_abort(subname//' ERROR in registering data lnd comp')
+    endif
+#endif
     call dlnd_comp_init(Eclock, x2l, l2x, &
          seq_flds_x2l_fields, seq_flds_l2x_fields, &
          SDLND, gsmap, ggrid, mpicom, compid, my_task, master_task, &
          inst_suffix, inst_name, logunit, read_restart, &
          scmMode, scmlat, scmlon)
-
+#ifdef HAVE_MOAB
+     if (my_task == master_task) then
+          call seq_infodata_PutData( infodata, lnd_domain=SDLND%domainFile) ! we use the same one for regular case
+          ! in regular case, it is copied from fatmlndfrc ; so we don't know if it is data land or not
+          write(logunit,*), ' use this land domain file: ', SDLND%domainFile
+     endif
+#endif
     !----------------------------------------------------------------------------
     ! Fill infodata that needs to be returned from dlnd
     !----------------------------------------------------------------------------

driver-moab/main/component_mod.F90

@@ -748,6 +748,7 @@ subroutine component_init_areacor_moab (comp, mbccid, mbcxid, seq_flds_c2x_fluxe
          lsize = comp(1)%mblsize
          allocate(areas (lsize, 3)) ! lsize is along grid; read mask too
          allocate(factors (lsize, 2))
+         factors = 1.0 ! initialize with 1.0 all factors; then maybe correct them
          ! get areas
          tagname='area:aream:mask'//C_NULL_CHAR
          arrsize = 3 * lsize

driver-moab/main/seq_frac_mct.F90

@@ -180,7 +180,7 @@ module seq_frac_mct
 
   use iMOAB, only : iMOAB_DefineTagStorage, iMOAB_SetDoubleTagStorage, &
         iMOAB_GetMeshInfo, iMOAB_SetDoubleTagStorageWithGid, iMOAB_WriteMesh, &
-        iMOAB_ApplyScalarProjectionWeights, iMOAB_SendElementTag, iMOAB_ReceiveElementTag, &
+        iMOAB_SendElementTag, iMOAB_ReceiveElementTag, &
         iMOAB_FreeSenderBuffers, iMOAB_GetDoubleTagStorage
 #ifdef MOABDEBUG
    use seq_comm_mct,     only : num_moab_exports