Merge branch 'iulian787/data_land_moab' (PR #6567)

Add data land capability for moab driver;
case tested: --res r05_r05 --compset RMOSGPCC
standard data moab model, mainly adapt from data ocean implementation

major change: Initialize all factors set by MOAB area correction init routine to 1.0. It affects all models,
but it seems that some values were not initialized for data land case.
On chrysalis, with intel compiler, it actually complained about non initialized values when applying correction factors.

The compressing file errors seen during development were caused by non-ascii characters in the log files.
Those were caused by an uninitialized mapper, for which we were still creating a MOAB map. It was never used,
it is fixed properly now, the moab map is not computed anymore.

[BFB]

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,15 @@ subroutine dlnd_comp_init(Eclock, x2l, l2x, &
        scmMode, scmlat, scmlon)
 
     ! !DESCRIPTION: initialize dlnd model
+#ifdef HAVE_MOAB
+       use iMOAB, only: iMOAB_DefineTagStorage, &
+                        iMOAB_SetIntTagStorage, iMOAB_SetDoubleTagStorage, &
+                        iMOAB_ResolveSharedEntities, iMOAB_CreateVertices, &
+                        iMOAB_UpdateMeshInfo
+#ifdef MOABDEBUG
+       use iMOAB, only: iMOAB_WriteMesh
+#endif
+#endif
     implicit none
 
     ! !INPUT/OUTPUT PARAMETERS:
@@ -135,6 +148,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 
+    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 +281,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 +477,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 +511,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 +620,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/prep_lnd_mod.F90

@@ -1,6 +1,6 @@
 module prep_lnd_mod
 
-  use shr_kind_mod    , only: r8 => SHR_KIND_R8
+  use shr_kind_mod    , only: R8 => SHR_KIND_R8
   use shr_kind_mod    , only: cs => SHR_KIND_CS
   use shr_kind_mod    , only: cl => SHR_KIND_CL
   use shr_kind_mod    , only: cxx => SHR_KIND_CXX
@@ -162,7 +162,7 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln
     integer  nrflds  ! number of rof fields projected on land
     integer arrsize  ! for setting the r2x fields on land to 0
     integer ent_type ! for setting tags
-    real (kind=r8) , allocatable :: tmparray (:) ! used to set the r2x fields to 0
+    real (kind=R8) , allocatable :: tmparray (:) ! used to set the r2x fields to 0
 
 #endif
     character(*), parameter  :: subname = '(prep_lnd_init)'
@@ -222,7 +222,6 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln
           call seq_map_init_rcfile(mapper_Fr2l, rof(1), lnd(1), &
                'seq_maps.rc','rof2lnd_fmapname:','rof2lnd_fmaptype:',samegrid_lr, &
                string='mapper_Fr2l initialization',esmf_map=esmf_map_flag)
-       end if
 ! symmetric of l2r, from prep_rof
 #ifdef HAVE_MOAB
           ! Call moab intx only if land and river are init in moab
@@ -370,7 +369,7 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln
             arrsize = nrflds*mlsize
             allocate (tmparray(arrsize)) ! mlsize is the size of local land
             ! do we need to zero out others or just river ?
-            tmparray = 0._r8
+            tmparray = 0._R8
             ierr = iMOAB_SetDoubleTagStorage(mblxid, tagname, arrsize , ent_type, tmparray)
             if (ierr .ne. 0) then
                write(logunit,*) subname,' cant zero out r2x tags on land'
@@ -381,6 +380,7 @@ subroutine prep_lnd_init(infodata, atm_c2_lnd, rof_c2_lnd, glc_c2_lnd, iac_c2_ln
          end if ! if ((mbrxid .ge. 0) .and.  (mblxid .ge. 0))
 ! endif HAVE_MOAB
 #endif
+       end if
        call shr_sys_flush(logunit)
 
        if (atm_c2_lnd) then
@@ -693,7 +693,7 @@ subroutine prep_lnd_mrg_moab (infodata)
 #endif
 #ifdef MOABCOMP
     character(CXX)           :: tagname, mct_field
-    real(r8)                 :: difference
+    real(R8)                 :: difference
     type(mct_list) :: temp_list
     integer :: size_list, index_list, ent_type
     type(mct_string)    :: mctOStr  !