Merge pull request #339 from fvitt/megan_expr

jedwards4b · web-flow · commit 1f06c58ac4ed · 2023-01-20T15:35:17.000-07:00
Updates to the MEGAN specifier string parser
diff --git a/cesm/nuopc_cap_share/shr_expr_parser_mod.F90 b/cesm/nuopc_cap_share/shr_expr_parser_mod.F90
@@ -1,12 +1,12 @@
 !=============================================================================
 ! expression parser utility --
 !   for parsing simple linear mathematical expressions of the form
-!   X = a*Y + b*Z + ...
+!   X = a*R + b*S + c*(X + Y + Z) ...
 !
 !=============================================================================
 module shr_expr_parser_mod
   use shr_kind_mod,only : r8 => shr_kind_r8
-  use shr_kind_mod,only : cx => shr_kind_cx
+  use shr_kind_mod,only : CXX => shr_kind_cxx
 
   implicit none
   private
@@ -35,82 +35,122 @@ function shr_exp_parse( exp_array, nitems ) result(exp_items_list)
     integer, optional, intent(out) :: nitems       ! number of expressions parsed
     type(shr_exp_item_t), pointer  :: exp_items_list ! linked list of items returned
 
-    integer :: i,j, jj, nmax, nterms, n_exp_items
-    character(len=cx) :: tmp_str
+    integer :: i,j, n_exp_items
     type(shr_exp_item_t), pointer :: exp_item, list_item
+    integer :: ndxs(512)
+    integer :: nelem, j1,j2,k
+    character(len=CXX) :: tmp_str, tmp_name
+    character(len=8) :: xchr ! multipler
+    real(r8) :: xdbl
+    real(r8) :: coeff0
+    logical :: more_to_come
+    character(len=CXX), allocatable :: sums_grps(:)
+    character(len=CXX) :: sum_string
+
+    allocate(sums_grps(size(exp_array)))
 
     nullify( exp_items_list )
     nullify( exp_item )
     nullify( list_item )
 
-    n_exp_items = 0
-    nmax = size( exp_array )
+    sums_grps(:) = ' '
 
-    do i = 1,nmax
-       if (len_trim(exp_array(i))>0) then
+    ! combine lines that have a trailing "+" with the next line
+    i=1
+    j=1
+    loop1: do while( len_trim(exp_array(i)) > 0 )
 
-          j = scan( exp_array(i), '=' )
+       k = scan(exp_array(i), '+', back=.true. )
+       more_to_come = k == len_trim(exp_array(i)) ! line ends with "+"
 
-          if ( j>0 ) then
+       if ( more_to_come ) then
+          sums_grps(j) = trim(sums_grps(j)) // trim(adjustl(exp_array(i)))
+       else
+          sums_grps(j) = trim(sums_grps(j)) // trim(adjustl(exp_array(i)))
+          j = j+1
+       endif
+
+       i = i+1
+       if ( i > size(exp_array) ) exit loop1
 
-             n_exp_items = n_exp_items + 1
+    end do loop1
 
-             allocate( exp_item )
-             exp_item%n_terms = 0
-             exp_item%name = trim(adjustl(exp_array(i)(:j-1)))
+    n_exp_items = j-1
 
-             tmp_str = trim(adjustl(exp_array(i)(j+1:)))
+    ! a group is  a summation of terms
 
-             nterms = 1
-             jj = scan( tmp_str, '+' )
-             do while(jj>0)
-                nterms = nterms + 1
-                tmp_str = tmp_str(jj+1:)
-                jj = scan( tmp_str, '+' )
-             enddo
+    ! parse the individual sum strings...  and form the groupings
+    has_grps: if (n_exp_items>0) then
 
-             allocate( exp_item%vars(nterms) )
-             allocate( exp_item%coeffs(nterms) )
+       ! from shr_megan_mod ... should be generalized and shared...
+       grploop: do i = 1,n_exp_items
 
-             tmp_str = trim(adjustl(exp_array(i)(j+1:)))
+          ! parse out the term names
+          ! from first parsing out the terms in the summation equation ("+" separates the terms)
 
+          sum_string = sums_grps(i)
+          j = scan( sum_string, '=' )
+          nelem = 1
+          ndxs(nelem) = j ! ndxs stores the index of each term of the equation
+
+          ! find indices of all the terms in the equation
+          tmp_str = trim( sum_string(j+1:) )
+          j = scan( tmp_str, '+' )
+          do while(j>0)
+             nelem = nelem+1
+             ndxs(nelem) = ndxs(nelem-1) + j
+             tmp_str = tmp_str(j+1:)
              j = scan( tmp_str, '+' )
+          enddo
+          ndxs(nelem+1) = len(sum_string)+1
 
-             if (j>0) then
-                call set_coefvar( tmp_str(:j-1), exp_item )
-                tmp_str = tmp_str(j-1:)
-             else
-                call set_coefvar( tmp_str, exp_item )
-             endif
+          allocate( exp_item )
 
-          else
+          exp_item%n_terms = nelem ! number of terms
 
-             tmp_str = trim(adjustl(exp_array(i))) ! assumed to begin with '+'
+          exp_item%name = trim(adjustl( sum_string(:ndxs(1)-1))) ! thing to the left of the "=" is used as the name of the group
 
-          endif
+          ! now that we have the number of terms in the summation allocate memory for the terms
+          allocate( exp_item%vars(nelem) )
+          allocate( exp_item%coeffs(nelem) )
 
-          ! at this point tmp_str begins with '+'
-          j = scan( tmp_str, '+' )
+          coeff0 = 1._r8 ! default multiplier
 
-          if (j>0) then
+          ! now parse out the multiplier from the terms
+          elmloop: do k = 1,nelem
 
-             ! remove the leading + ...
-             tmp_str = tmp_str(j+1:)
-             j = scan( tmp_str, '+' )
+             exp_item%coeffs(k) = coeff0
 
-             do while(j>0)
+             ! get the term name which follows the '*' operator if the is one
+             tmp_name = adjustl(sum_string(ndxs(k)+1:ndxs(k+1)-1))
 
-                call set_coefvar( tmp_str(:j-1), exp_item )
+             j = scan( tmp_name, '*' )
+             if (j>0) then
 
-                tmp_str = tmp_str(j+1:)
-                j = scan( tmp_str, '+' )
+                xchr = tmp_name(1:j-1) ! get the multipler (left of the '*')
+                read( xchr, * ) xdbl   ! convert the string to a real
+                exp_item%coeffs(k) = xdbl ! store the multiplier
 
-             enddo
+                j1 = scan( tmp_name, '(' )
+                if (j1>0) then
+                   coeff0 = xdbl
+                   tmp_name = trim(adjustl(tmp_name(j1+1:))) ! get the term name (right of the '*')
+                else
+                   coeff0 = 1._r8
+                   tmp_name = trim(adjustl(tmp_name(j+1:))) ! get the term name (right of the '*')
+                endif
 
-             call set_coefvar( tmp_str, exp_item )
+             endif
 
-          endif
+             j2 = scan( tmp_name, ')' )
+             if (j2>0) then
+                coeff0 = 1._r8
+                tmp_name = tmp_name(1:j2-1)
+             endif
 
+             exp_item%vars(k) = trim(tmp_name)
+
+          enddo elmloop
 
           if (associated(exp_item)) then
              if (associated(exp_items_list)) then
@@ -124,13 +164,16 @@ function shr_exp_parse( exp_array, nitems ) result(exp_items_list)
              endif
           endif
 
-       endif
-    enddo
+
+       enddo grploop
+    endif has_grps
 
     if ( present(nitems) ) then
        nitems = n_exp_items
     endif
 
+    deallocate(sums_grps)
+
   end function shr_exp_parse
 
   ! -----------------------------------------------------------------
@@ -157,29 +200,4 @@ subroutine  shr_exp_list_destroy( list )
 
   end subroutine  shr_exp_list_destroy
 
- !==========================
- ! Private Methods
-
-  ! -----------------------------------------------------------------
-  ! -----------------------------------------------------------------
-  subroutine set_coefvar( term, item )
-    character(len=*), intent(in)  :: term
-    type(shr_exp_item_t) , intent(inout) :: item
-
-    integer :: k, n
-
-    item%n_terms = item%n_terms + 1
-    n = item%n_terms
-
-    k = scan( term, '*' )
-    if (k>0) then
-       item%vars(n) = trim(adjustl(term(k+1:)))
-       read( term(:k-1), *) item%coeffs(n)
-    else
-       item%vars(n) = trim(adjustl(term))
-       item%coeffs(n) = 1.0_r8
-    endif
-
-  end subroutine set_coefvar
-
 end module shr_expr_parser_mod
diff --git a/cesm/nuopc_cap_share/shr_megan_mod.F90 b/cesm/nuopc_cap_share/shr_megan_mod.F90
@@ -5,22 +5,22 @@ module shr_megan_mod
   ! MEGAN = Model of Emissions of Gases and Aerosols from Nature
   !
   ! This reads the megan_emis_nl namelist in drv_flds_in and makes the relavent
-  ! information available to CAM, CLM, and driver. 
-  ! - The driver sets up CLM to CAM communication for the  VOC flux fields. 
-  ! - CLM needs to know what specific VOC fluxes need to be passed to the coupler 
+  ! information available to CAM, CLM, and driver.
+  ! - The driver sets up CLM to CAM communication for the  VOC flux fields.
+  ! - CLM needs to know what specific VOC fluxes need to be passed to the coupler
   !   and how to assemble the fluxes.
   ! - CAM needs to know what specific VOC fluxes to expect from CLM.
   !================================================================================
 
   use ESMF                , only : ESMF_VMGetCurrent, ESMF_VM, ESMF_VMGet
-  use ESMF                , only : ESMF_LogFoundError, ESMF_LOGERR_PASSTHRU, ESMF_SUCCESS
-  use shr_kind_mod        , only : r8 => shr_kind_r8, cl=>shr_kind_cl, cx=>shr_kind_cx, cs=>shr_kind_cs
+  use ESMF                , only : ESMF_LogFoundError, ESMF_LOGERR_PASSTHRU
+  use shr_kind_mod        , only : r8 => shr_kind_r8, cl=>shr_kind_cl, cx=>shr_kind_cx
   use shr_sys_mod         , only : shr_sys_abort
   use shr_log_mod         , only : shr_log_getLogUnit
   use shr_mpi_mod         , only : shr_mpi_bcast
   use shr_nl_mod          , only : shr_nl_find_group_name
   use shr_expr_parser_mod , only : shr_exp_parse, shr_exp_item_t, shr_exp_list_destroy
-  
+
   implicit none
   private
 
@@ -68,6 +68,9 @@ module shr_megan_mod
   ! switch to use mapped emission factors
   logical :: shr_megan_mapped_emisfctrs = .false.
 
+  integer :: localPet = -huge(1)
+  integer :: logunit = -huge(1)
+
 !--------------------------------------------------------
 contains
 !--------------------------------------------------------
@@ -100,7 +103,8 @@ subroutine shr_megan_readnl( NLFileName, megan_nflds)
     ! Example:
     ! &megan_emis_nl
     !  megan_specifier = 'ISOP = isoprene',
-    !     'C10H16 = myrcene + sabinene + limonene + carene_3 + ocimene_t_b + pinene_b + ...',
+    !     'C10H16 = myrcene + sabinene + limonene + carene_3 + ocimene_t_b + pinene_b + ',
+    !     ' thujene_a + bornene + 0.5*(terpineol_4 + terpineol_a + terpinyl_ACT_a + myrtenal) + ...',
     !     'CH3OH = methanol',
     !     'C2H5OH = ethanol',
     !     'CH2O = formaldehyde',
@@ -109,25 +113,22 @@ subroutine shr_megan_readnl( NLFileName, megan_nflds)
     !  megan_factors_file = '$datapath/megan_emis_factors.nc'
     ! /
     !-------------------------------------------------------------------------
-    
+
     ! input/output variables
     character(len=*), intent(in)  :: NLFileName
     integer,          intent(out) :: megan_nflds
 
     ! local variables
     type(ESMF_VM)       :: vm
-    integer             :: localPet
     integer             :: mpicom
     integer             :: unitn            ! namelist unit number
     integer             :: ierr             ! error code
     logical             :: exists           ! if file exists or not
-    integer, parameter  :: maxspc = 100
-    character(len=2*CX) :: megan_specifier(maxspc) = ' '
+    integer, parameter  :: maxspc = 200
+    character(len=CX)   :: megan_specifier(maxspc) = ' '
     logical             :: megan_mapped_emisfctrs = .false.
     character(len=CL)   :: megan_factors_file = ' '
     integer             :: rc
-    integer             :: logunit
-    integer             :: i, tmp(1)
     character(*), parameter :: F00   = "('(shr_megan_readnl) ',2a)"
     character(len=*), parameter :: subname='(shr_megan_readnl)'
     !--------------------------------------------------------------
@@ -140,12 +141,12 @@ subroutine shr_megan_readnl( NLFileName, megan_nflds)
     end if
 
     call ESMF_VMGetCurrent(vm, rc=rc)
-    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return 
+    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
 
     call ESMF_VMGet(vm, localPet=localPet, mpiCommunicator=mpicom, rc=rc)
-    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return 
+    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
     call shr_log_getLogUnit(logunit)
-    ! Note the following still needs to be called on all processors since the mpi_bcast is a collective 
+    ! Note the following still needs to be called on all processors since the mpi_bcast is a collective
     ! call on all the pes of mpicom
     if (localPet==0) then
        inquire( file=trim(NLFileName), exist=exists)
@@ -204,6 +205,8 @@ subroutine shr_megan_init( specifier)
     allocate(shr_megan_mechcomps(n_entries))
     shr_megan_mechcomps(:)%n_megan_comps = 0
 
+    if (localPet==0) write(logunit,*) 'MEGAN entries:'
+
     item => items_list
     i = 1
     do while(associated(item))
@@ -221,7 +224,9 @@ subroutine shr_megan_init( specifier)
        shr_megan_mechcomps(i)%n_megan_comps = item%n_terms
        allocate(shr_megan_mechcomps(i)%megan_comps(item%n_terms))
 
+       if (localPet==0) write(logunit,*) ' species : ', item%name
        do j = 1,item%n_terms
+          if (localPet==0) write(logunit,'(f12.4,a,a)')  item%coeffs(j),' * ', item%vars(j)
           shr_megan_mechcomps(i)%megan_comps(j)%ptr => add_megan_comp( item%vars(j), item%coeffs(j) )
        enddo
        shr_megan_mechcomps_n = shr_megan_mechcomps_n+1
