Changes needed for CAM computing enthalpy -> mediator -> BLOM

cesm/mod_cesm.F90

@@ -33,7 +33,8 @@ module mod_cesm
                              lamult, lasl, ustokes, vstokes, atmco2, atmbrf, &
                              flxdms, flxbrf, &
                              atmn2o, atmnh3, atmnhxdep, atmnoydep, &
-                             use_stream_relaxation
+                             use_stream_relaxation, &
+                             hmat
    use mod_ben02,      only: initai, rdcsic, rdctsf, fnlzai
    use mod_rdcsss,     only: rdcsss
    use mod_idarlx,     only: idarlx
@@ -66,6 +67,7 @@ module mod_cesm
    real(r8), dimension(1 - nbdy:idm + nbdy,1 - nbdy:jdm + nbdy, 2) :: &
       swa_da, &          ! Solar heat flux [W m-2].
       nsf_da, &          ! Non-solar heat flux [W m-2].
+      hmat_da, &         ! surf.mat.enth.flx
       hmlt_da, &         ! Heat flux due to melting [W m-2].
       lip_da, &          ! Liquid water flux [kg m-2 s-1].
       sop_da, &          ! Solid precipitation [kg m-2 s-1].
@@ -106,7 +108,8 @@ module mod_cesm
              slp_da, abswnd_da, ficem_da, lamult_da, lasl_da, flxdms_da, &
              flxbrf_da, ustokes_da, vstokes_da, atmco2_da, atmbrf_da, &
              atmn2o_da, atmnh3_da, atmnhxdep_da,atmnoydep_da, smtfrc, &
-             l1ci, l2ci, inivar_cesm, inicon_cesm, inifrc_cesm, getfrc_cesm
+             l1ci, l2ci, inivar_cesm, inicon_cesm, inifrc_cesm, getfrc_cesm, &
+             hmat_da
 
 contains
 
@@ -229,6 +232,7 @@ subroutine getfrc_cesm
            sfl(i, j)     = w1*sfl_da(i, j, l1ci)     + w2*sfl_da(i, j, l2ci)
            swa(i, j)     = w1*swa_da(i, j, l1ci)     + w2*swa_da(i, j, l2ci)
            nsf(i, j)     = w1*nsf_da(i, j, l1ci)     + w2*nsf_da(i, j, l2ci)
+           hmat(i, j)    = w1*hmat_da(i, j, l1ci)    + w2*hmat_da(i, j, l2ci)
            hmlt(i, j)    = w1*hmlt_da(i, j, l1ci)    + w2*hmlt_da(i, j, l2ci)
            slp(i, j)     = w1*slp_da(i, j, l1ci)     + w2*slp_da(i, j, l2ci)
            abswnd(i, j)  = w1*abswnd_da(i, j, l1ci)  + w2*abswnd_da(i, j, l2ci)
@@ -272,6 +276,7 @@ subroutine getfrc_cesm
         call ncdefvar('sfl_da', 'x y', ndouble, 8)
         call ncdefvar('swa_da', 'x y', ndouble, 8)
         call ncdefvar('nsf_da', 'x y', ndouble, 8)
+        call ncdefvar('hmat_da', 'x y', ndouble, 8)
         call ncdefvar('hmlt_da', 'x y', ndouble, 8)
         call ncdefvar('slp_da', 'x y', ndouble, 8)
         call ncdefvar('abswnd_da', 'x y', ndouble, 8)
@@ -310,6 +315,8 @@ subroutine getfrc_cesm
              ip, 1, 1._r8, 0._r8, 8)
         call ncwrtr('nsf_da', 'x y', nsf_da(1 - nbdy, 1 - nbdy, l2ci), &
              ip, 1, 1._r8, 0._r8, 8)
+        call ncwrtr('hmat_da', 'x y', hmat_da(1 - nbdy, 1 - nbdy, l2ci), &
+             ip, 1, 1._r8, 0._r8, 8)
         call ncwrtr('hmlt_da', 'x y', hmlt_da(1 - nbdy, 1 - nbdy, l2ci), &
              ip, 1, 1._r8, 0._r8, 8)
         call ncwrtr('slp_da', 'x y', slp_da(1 - nbdy, 1 - nbdy, l2ci), &
@@ -362,6 +369,7 @@ subroutine getfrc_cesm
          call chksum(swa, 1, halo_ps, 'swa')
          call chksum(nsf, 1, halo_ps, 'nsf')
          call chksum(hmlt, 1, halo_ps, 'hmlt')
+         call chksum(hmat, 1, halo_ps, 'hmat')
          call chksum(slp, 1, halo_ps, 'slp')
          call chksum(abswnd, 1, halo_ps, 'abswnd')
          call chksum(ficem, 1, halo_ps, 'ficem')

cime_config/testdefs/ExpectedTestFails.xml

@@ -44,13 +44,6 @@
     </phase>
   </test>
 
-  <test name="ERS_Ld3.TL319_tn05.NOIIAJRA.betzy_intel">
-    <phase name="COMPARE_base_rest">
-      <status>FAIL</status>
-      <issue>BLOM#664</issue>
-    </phase>
-  </test>
-
   <test name="SMS_D_Ld1.T62_tn14.NOIIAOC20TR.betzy_intel">
     <phase name="RUN">
       <status>FAIL</status>

drivers/nuopc/ocn_comp_nuopc.F90

@@ -390,6 +390,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
       character(len=cllen) :: msg, cvalue
       logical :: isPresent, isSet
       logical :: ocn2glc_coupling
+      logical :: atm_computes_enthalpy_flux
       logical :: flds_co2a, flds_co2c, flds_dms, flds_brf
       logical :: hamocc_defined
 #ifndef HAMOCC
@@ -538,6 +539,13 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
       read(cvalue,*) flds_co2c
       call blom_logwrite(subname//': flds_co2c = '//trim(cvalue))
 
+      ! Determine if atm computes enthalpy fluxes
+      call NUOPC_CompAttributeGet(gcomp, name="atm_computes_enthalpy_flux", value=cvalue, rc=rc)
+      if (ChkErr(rc, __LINE__, u_FILE_u)) return
+      read(cvalue,*) atm_computes_enthalpy_flux
+      write(msg,'(a,l1)') subname//': atm_computes_enthalpy_flux is ', atm_computes_enthalpy_flux
+      call blom_logwrite(msg)
+
       ! Determine if ocn is sending temperature and salinity data to glc
       ! If data is sent to glc will need to determine number of ocean
       ! levels and ocean level indices
@@ -594,12 +602,14 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
       write(msg,'(a,l1)') subname//': export brf ', flds_brf
       call blom_logwrite(msg)
 
+
+
       ! ------------------------------------------------------------------------
       ! Advertise import fields.
       ! ------------------------------------------------------------------------
 
       call blom_advertise_imports(flds_scalar_name, fldsToOcn_num, fldsToOcn, &
-           flds_co2a, flds_co2c)
+           flds_co2a, flds_co2c, atm_computes_enthalpy_flux)
 
       do n = 1,fldsToOcn_num
          call NUOPC_Advertise(importState, standardName=fldsToOcn(n)%stdname, &

drivers/nuopc/ocn_import_export.F90

@@ -32,7 +32,7 @@ module ocn_import_export
    use mod_constants,  only: rearth, onem
    use mod_time,       only: nstep, baclin, delt1, dlt
    use mod_xc
-   use mod_grid,       only: scuy, scvx, scp2, scuxi, scvyi, plon, plat, cosang, sinang
+   use mod_grid,       only: scuy, scvx, scp2, scuxi, scvyi, plon, plat, cosang, sinang, area
    use mod_state,      only: u, v, dp, temp, saln, pbu, pbv, ubflxs, vbflxs, sealv
    use mod_forcing,    only: wavsrc_opt, wavsrc_extern, sprfac, prfac, &
                              flxco2, flxdms, flxbrf, flxn2o, flxnh3
@@ -44,8 +44,8 @@ module ocn_import_export
                              ustarw_da, slp_da, abswnd_da, ficem_da, lamult_da, &
                              lasl_da, ustokes_da, vstokes_da, atmco2_da, &
                              atmnhxdep_da, atmnoydep_da, &
-                             l1ci, l2ci
-   use mod_utility,    only: util1, util2
+                             l1ci, l2ci, hmat_da
+   use mod_utility,    only: util1, util2, util3, util4
    use mod_checksum,   only: csdiag, chksum
 #ifdef HAMOCC
    use mo_control_bgc, only: use_BROMO, ocn_co2_type
@@ -143,6 +143,9 @@ module ocn_import_export
         index_Faxa_snow   = -1, &
         index_Faxa_rain   = -1, &
         index_Faxa_ndep   = -1, &
+        index_Faxa_hmat   = -1, &
+        index_Faxa_hlat   = -1, &
+        index_Faxa_hmoa   = -1, &
         index_Sa_pslv     = -1, &
         index_Sa_co2diag  = -1, &
         index_Sa_co2prog  = -1, &
@@ -216,7 +219,7 @@ end subroutine fldlist_add
    ! ---------------------------------------------------------------------------
 
    subroutine blom_advertise_imports(flds_scalar_name, fldsToOcn_num, fldsToOcn, &
-        flds_co2a, flds_co2c)
+        flds_co2a, flds_co2c, atm_computes_enthalpy)
 
      ! -------------------------------------------------------------------
      ! Determine fldsToOcn for import fields
@@ -227,6 +230,7 @@ subroutine blom_advertise_imports(flds_scalar_name, fldsToOcn_num, fldsToOcn, &
      type(fldlist_type) , intent(inout) :: fldsToOcn(:)
      logical            , intent(in)    :: flds_co2a
      logical            , intent(in)    :: flds_co2c
+     logical            , intent(in)    :: atm_computes_enthalpy
 
      integer :: index_scalar
 
@@ -272,6 +276,13 @@ subroutine blom_advertise_imports(flds_scalar_name, fldsToOcn_num, fldsToOcn, &
      call fldlist_add(fldsToOcn_num, fldsToOcn, 'Faxa_rain' , index_Faxa_rain)
      call fldlist_add(fldsToOcn_num, fldsToOcn, 'Faxa_ndep' , index_Faxa_ndep, &
           ungridded_lbound=1, ungridded_ubound=2)
+     if (atm_computes_enthalpy) then
+        call fldlist_add(fldsToOcn_num, fldsToOcn, 'Faxa_hmat'    , index_Faxa_hmat)
+        call fldlist_add(fldsToOcn_num, fldsToOcn, 'Faxa_hmat_oa' , index_Faxa_hmoa)
+        ! Note the following was added to avoid a mapping in the mediator of
+        ! Faxa_hlat from the atm to the ocn grid - it is not used in BLOM at the moment
+        call fldlist_add(fldsToOcn_num, fldsToOcn, 'Faxa_hlat'    , index_Faxa_hlat)
+     end if
      if (flds_co2a .or. flds_co2c) then
         call fldlist_add(fldsToOcn_num, fldsToOcn, 'Sa_co2diag' ,index_Sa_co2diag)
         call fldlist_add(fldsToOcn_num, fldsToOcn, 'Sa_co2prog', index_Sa_co2prog)
@@ -725,13 +736,14 @@ subroutine blom_importflds(fldlist_num, fldlist)
       real(r8), parameter :: &
          mval = - 1.e12_r8, &
          fval = - 1.e13_r8
+
       logical :: first_call = .true.
+      integer :: hmat_method = 2
 
       ! Local variables.
-      real(r8) :: afac, utmp, vtmp
-      integer :: n, i, j, l
-      integer :: index_co2
-      real(r8):: rofi_heat_flx, snow_heat_flx
+      real(r8) :: afac, utmp, vtmp, rofi_heat_flx, snow_heat_flx, &
+                  hmat_oa_asum, oocn_asum, hmat_asum, hmat_oa_avg, hmat_avg
+      integer :: n, i, j, l, index_co2
 
       ! Update time level indices.
       if (l1ci == 1 .and. l2ci == 1) then
@@ -808,7 +820,6 @@ subroutine blom_importflds(fldlist_num, fldlist)
                sfl_da(i,j,l2ci) = mval
                swa_da(i,j,l2ci) = mval
                nsf_da(i,j,l2ci) = mval
-               hmlt_da(i,j,l2ci) = mval
                slp_da(i,j,l2ci) = mval
                abswnd_da(i,j,l2ci) = mval
                ficem_da(i,j,l2ci) = mval
@@ -824,7 +835,6 @@ subroutine blom_importflds(fldlist_num, fldlist)
                sfl_da(i,j,l2ci) = 0._r8
                swa_da(i,j,l2ci) = 0._r8
                nsf_da(i,j,l2ci) = 0._r8
-               hmlt_da(i,j,l2ci) = 0._r8
                slp_da(i,j,l2ci) = fval
                abswnd_da(i,j,l2ci) = fval
                ficem_da(i,j,l2ci) = fval
@@ -900,6 +910,112 @@ subroutine blom_importflds(fldlist_num, fldlist)
       enddo
       !$omp end parallel do
 
+      if (index_Faxa_hmat > 0 .and. index_Faxa_hmoa > 0) then
+         !$omp parallel do private(i, n, afac)
+         do j = 1, jjcpl
+            do i = 1, ii
+               if (ip(i,j) == 0) then
+                  hmat_da(i,j,l1ci)= mval
+                  hmlt_da(i,j,l2ci) = mval
+               elseif (cplmsk(i,j) == 0) then
+                  hmat_da(i,j,l1ci) = 0._r8
+                  hmlt_da(i,j,l2ci) = 0._r8
+               else
+                  n = (j - 1)*ii + i
+                  afac = med2mod_areacor(n)
+                  ! Heat flux components due to material material enthalpy flux of
+                  ! water exchange with atmosphere [W m-2]. Here, index_Faxa_hmat
+                  ! is related to enthalpy flux of evaporation and index_Faxa_hmoa
+                  ! related to the ocean average of all other enthalpy flux
+                  ! components.
+                  util1(i,j) = fldlist(index_Faxa_hmat)%dataptr(n)*afac
+                  util2(i,j) = fldlist(index_Faxa_hmoa)%dataptr(n)*afac
+               end if
+            end do
+         end do
+         select case (hmat_method)
+         case (1)
+            ! Apply enthalpy flux components directly.
+            do j = 1, jjcpl
+               do l = 1, isp(j)
+                  do i = max(1, ifp(j,l)), min(ii, ilp(j,l))
+                     hmat_da(i,j,l2ci) = util1(i,j) + util2(i,j)
+                     nsf_da(i,j,l2ci) = nsf_da(i,j,l2ci) + hmat_da(i,j,l2ci)
+                  enddo
+               enddo
+            enddo
+         case (2)
+            ! Redistribute 'hmat_oa' to open ocean area.
+            do j = 1, jjcpl
+               do l = 1, isp(j)
+                  do i = max(1, ifp(j,l)), min(ii, ilp(j,l))
+                     util3(i,j) = util2(i,j)*scp2(i,j)
+                     util4(i,j) = (1._r8 - ficem_da(i,j,l2ci))*scp2(i,j)
+                  enddo
+               enddo
+            enddo
+            call xcsum(hmat_oa_asum, util3, ips)
+            call xcsum(oocn_asum   , util4, ips)
+            hmat_oa_avg = hmat_oa_asum/oocn_asum
+            do j = 1, jjcpl
+               do l = 1, isp(j)
+                  do i = max(1, ifp(j,l)), min(ii, ilp(j,l))
+                     hmat_da(i,j,l2ci) = &
+                          util1(i,j) + hmat_oa_avg*(1._r8 - ficem_da(i,j,l2ci))
+                     nsf_da(i,j,l2ci) = nsf_da(i,j,l2ci) + hmat_da(i,j,l2ci)
+                  enddo
+               enddo
+            enddo
+         case (3)
+            ! Apply global average enthalpy flux.
+            do j = 1, jjcpl
+               do l = 1, isp(j)
+                  do i = max(1, ifp(j,l)), min(ii, ilp(j,l))
+                     util3(i,j) = (util1(i,j) + util2(i,j))*scp2(i,j)
+                  enddo
+               enddo
+            enddo
+            call xcsum(hmat_asum, util3, ips)
+            hmat_avg = hmat_asum/area
+            do j = 1, jjcpl
+               do l = 1, isp(j)
+                  do i = max(1, ifp(j,l)), min(ii, ilp(j,l))
+                     hmat_da(i,j,l2ci) = hmat_avg
+                     nsf_da(i,j,l2ci) = nsf_da(i,j,l2ci) + hmat_da(i,j,l2ci)
+                  enddo
+               enddo
+            enddo
+         case (4)
+            ! Apply global average enthalpy flux over open ocean.
+            do j = 1, jjcpl
+               do l = 1, isp(j)
+                  do i = max(1, ifp(j,l)), min(ii, ilp(j,l))
+                     util3(i,j) = (util1(i,j) + util2(i,j))*scp2(i,j)
+                     util4(i,j) = (1._r8 - ficem_da(i,j,l2ci))*scp2(i,j)
+                  enddo
+               enddo
+            enddo
+            call xcsum(hmat_asum, util3, ips)
+            call xcsum(oocn_asum, util4, ips)
+            hmat_avg = hmat_asum/oocn_asum
+            do j = 1, jjcpl
+               do l = 1, isp(j)
+                  do i = max(1, ifp(j,l)), min(ii, ilp(j,l))
+                     hmat_da(i,j,l2ci) = hmat_avg*(1._r8 - ficem_da(i,j,l2ci))
+                     nsf_da(i,j,l2ci) = nsf_da(i,j,l2ci) + hmat_da(i,j,l2ci)
+                  enddo
+               enddo
+            enddo
+         case default
+            write(lp,*) subname//': BLOM ERROR: Unsupported hmat_method'
+            call xcstop(subname)
+            stop(subname)
+         end select
+      else
+         hmat_da(:,:,:) = mval
+         hmlt_da(:,:,:) = mval
+      end if
+
       if (nreg == 2) then
          call xctilr(lip_da(1-nbdy,1-nbdy,l2ci), 1,1, 0,0, halo_ps)
          call xctilr(sop_da(1-nbdy,1-nbdy,l2ci), 1,1, 0,0, halo_ps)
@@ -911,6 +1027,9 @@ subroutine blom_importflds(fldlist_num, fldlist)
          call xctilr(swa_da(1-nbdy,1-nbdy,l2ci), 1,1, 0,0, halo_ps)
          call xctilr(nsf_da(1-nbdy,1-nbdy,l2ci), 1,1, 0,0, halo_ps)
          call xctilr(hmlt_da(1-nbdy,1-nbdy,l2ci), 1,1, 0,0, halo_ps)
+         if (index_Faxa_hmat > 0) then
+            call xctilr(hmat_da(1-nbdy,1-nbdy,l2ci), 1,1, 0,0, halo_ps)
+         end if
          call xctilr(atmnhxdep_da(1-nbdy,1-nbdy,l2ci), 1,1, 0,0, halo_ps)
          call xctilr(atmnoydep_da(1-nbdy,1-nbdy,l2ci), 1,1, 0,0, halo_ps)
       endif
@@ -1069,6 +1188,9 @@ subroutine blom_importflds(fldlist_num, fldlist)
          call chksum(atmco2_da   (1-nbdy,1-nbdy,l2ci), 1, halo_ps, 'atmco2'   )
          call chksum(atmnhxdep_da(1-nbdy,1-nbdy,l2ci), 1, halo_ps, 'atmnhxdep')
          call chksum(atmnoydep_da(1-nbdy,1-nbdy,l2ci), 1, halo_ps, 'atmnoydep')
+         if (index_Faxa_hmat > 0) then
+            call chksummsk(hmat_da(1-nbdy,1-nbdy,l2ci),1,halo_ps,'hmat')
+         end if
       endif
 
       if (first_call) then