Merge pull request #174 from mvertens/feature/ocn_precip_factors

add functionality for ocean precip scaling for CESM
### Description of changes
This change only effects CESM compsets

### Specific notes
This adds the capability currently only activated by CESM C and G compsets (with POP not MOM) to have a precip factor sent by POP to the mediator that is then applied to scale atmosphere rain and snow going back to the ocean and sea-ice (in this case POP and CICE/DICE). In the current testing, the precip factor was only set to 1 by POP - as a result, this PR does not change any baselines.

Contributors other than yourself, if any:

CMEPS Issues Fixed: None

Are changes expected to change answers?
 - [x] bit for bit
 - [ ] different at roundoff level
 - [ ] more substantial 

Any User Interface Changes (namelist or namelist defaults changes)?
 - [ ] Yes
 - [x] No

Testing performed if application target is CESM:(either UFS-S2S or CESM testing is required):
- Verified manually that for SMS_Ld1_Vnuopc.T62_g17.G.cheyenne_intel.pop-default, the scaling factor to the mediator and that was set correctly in med_phases_prep_ocn_mod.F90 and that the results were bfb
- Verified that SMS_Vnuopc_Ld5.T62_t061.GMOM.cheyenne_intel was bfb before and after this change.

cime_config/namelist_definition_drv.xml

@@ -1643,7 +1643,7 @@
   <!-- </entry> -->
 
   <!-- =========================== -->
-  <!-- MED attributes              -->
+  <!-- ALLCOMP attributes          -->
   <!-- =========================== -->
 
   <entry id="ScalarFieldName">
@@ -1666,7 +1666,7 @@
       total number of scalars in the scalar coupling field
     </desc>
     <values>
-      <value>5</value>
+      <value>4</value>
     </values>
   </entry>
 
@@ -1714,8 +1714,8 @@
       index of scalar containing epbal precipitation factor from ocn (only for POP)
     </desc>
     <values>
-      <value flux_epbal='off'>0</value>
-      <value flux_epbal="ocn">5</value>
+      <value COMP_OCN="pop">4</value>
+      <value>0</value>
     </values>
   </entry>
 

mediator/med.F90

@@ -905,7 +905,7 @@ subroutine InitializeIPDv03p1(gcomp, importState, exportState, clock, rc)
     if (isPresent .and. isSet) then
        read(cvalue,*) is_local%wrap%flds_scalar_index_precip_factor
     else
-       is_local%wrap%flds_scalar_index_precip_factor = spval
+       is_local%wrap%flds_scalar_index_precip_factor = 0
     end if
 
     !------------------

mediator/med_internalstate_mod.F90

@@ -62,6 +62,7 @@ module med_internalstate_mod
     integer                :: flds_scalar_index_ny = 0
     integer                :: flds_scalar_index_nextsw_cday = 0
     integer                :: flds_scalar_index_precip_factor = 0
+    real(r8)               :: flds_scalar_precip_factor = 1._r8  ! actual value of precip factor from ocn
 
     ! Import/export States and field bundles (the field bundles have the scalar fields removed)
     type(ESMF_State)       :: NStateImp(ncomps)                  ! Import data from various component, on their grid

mediator/med_io_mod.F90

@@ -1053,7 +1053,7 @@ subroutine med_io_write_FB(filename, iam, FB, whead, wdata, nx, ny, nt, &
        do k = 1,nf
           call FB_getNameN(FB, k, itemc, rc=rc)
           if (chkerr(rc,__LINE__,u_FILE_u)) return
-
+ 
           call FB_getFldPtr(FB, itemc, &
                fldptr1=fldptr1, fldptr2=fldptr2, rank=rank, rc=rc)
           if (chkerr(rc,__LINE__,u_FILE_u)) return

mediator/med_phases_prep_ice_mod.F90

@@ -29,10 +29,12 @@ subroutine med_phases_prep_ice(gcomp, rc)
     use ESMF                  , only : ESMF_FieldBundleGet, ESMF_FieldGet, ESMF_Field
     use ESMF                  , only : ESMF_LOGMSG_ERROR, ESMF_FAILURE
     use ESMF                  , only : ESMF_StateItem_Flag, ESMF_STATEITEM_NOTFOUND
-    use med_utils_mod         , only : chkerr      => med_utils_ChkErr
-    use med_methods_mod       , only : fldchk      => med_methods_FB_FldChk
-    use med_methods_mod       , only : FB_diagnose => med_methods_FB_diagnose
-    use med_constants_mod     , only : dbug_flag   => med_constants_dbug_flag
+    use ESMF                  , only : ESMF_VMBroadCast 
+    use med_utils_mod         , only : chkerr       => med_utils_ChkErr
+    use med_methods_mod       , only : FB_fldchk    => med_methods_FB_FldChk
+    use med_methods_mod       , only : FB_diagnose  => med_methods_FB_diagnose
+    use med_methods_mod       , only : FB_GetFldPtr => med_methods_FB_GetFldPtr
+    use med_constants_mod     , only : dbug_flag    => med_constants_dbug_flag
     use med_merge_mod         , only : med_merge_auto
     use med_internalstate_mod , only : InternalState, logunit, mastertask
     use esmFlds               , only : compatm, compice, compocn, comprof, compglc, ncomps, compname
@@ -45,12 +47,12 @@ subroutine med_phases_prep_ice(gcomp, rc)
     integer, intent(out) :: rc
 
     ! local variables
-    type(ESMF_StateItem_Flag)      :: itemType
     type(InternalState)            :: is_local
     type(ESMF_Field)               :: lfield
     integer                        :: i,n
-    real(R8), pointer              :: dataptr1d(:) => null()
-    real(R8)                       :: precip_fact
+    real(R8), pointer              :: dataptr(:) => null()
+    real(R8), pointer              :: dataptr_scalar_ocn(:,:) => null()
+    real(R8)                       :: precip_fact(1)
     character(len=CS)              :: cvalue
     character(len=64), allocatable :: fldnames(:)
     real(r8)                       :: nextsw_cday
@@ -87,33 +89,50 @@ subroutine med_phases_prep_ice(gcomp, rc)
          fldListTo(compice), rc=rc)
     if (chkerr(rc,__LINE__,u_FILE_u)) return
 
-    ! apply precipitation factor from ocean
-    ! TODO (mvertens, 2019-03-18): precip_fact here is not valid if
-    ! the component does not send it - hardwire it to 1 until this is resolved
-    if (trim(coupling_mode) == 'cesm') then
-       precip_fact = 1.0_R8
-       if (precip_fact /= 1.0_R8) then
-          if (first_precip_fact_call .and. mastertask) then
-             write(logunit,'(a)')'(merge_to_ice): Scaling rain, snow, liquid and ice runoff by precip_fact '
-             first_precip_fact_call = .false.
+    ! Apply precipitation factor from ocean (that scales atm rain and snow to ice) if appropriate
+    if (trim(coupling_mode) == 'cesm' .and. is_local%wrap%flds_scalar_index_precip_factor /= 0) then
+
+       ! Note that in med_internal_mod.F90 all is_local%wrap%flds_scalar_index_precip_factor 
+       ! is initialized to 0.
+       ! In addition, in med.F90, if this attribute is not present as a mediator component attribute, 
+       ! it is set to 0. 
+       if (mastertask) then
+          call ESMF_StateGet(is_local%wrap%NstateImp(compocn), & 
+               itemName=trim(is_local%wrap%flds_scalar_name), field=lfield, rc=rc)
+          if (chkerr(rc,__LINE__,u_FILE_u)) return
+          call ESMF_FieldGet(lfield, farrayPtr=dataptr_scalar_ocn, rc=rc)
+          if (chkerr(rc,__LINE__,u_FILE_u)) return
+          scalar_id=is_local%wrap%flds_scalar_index_precip_factor
+          precip_fact(1) = dataptr_scalar_ocn(scalar_id,1)
+          if (first_call) then
+             write(logunit,'(a)')'(merge_to_ice): Scaling rain, snow, liquid and ice runoff by precip_fact from ocn'
+             first_call = .false.
           end if
-          write(cvalue,*) precip_fact
+          if (precip_fact(1) /= 1._r8) then
+             write(logunit,'(a,f21.13)')&
+                  '(merge_to_ice): Scaling rain, snow, liquid and ice runoff by non-unity precip_fact ',&
+                  precip_fact(1)
+          end if
+       end if
+       call ESMF_VMBroadCast(is_local%wrap%vm, precip_fact, 1, 0, rc=rc)
+       if (chkerr(rc,__LINE__,u_FILE_u)) return
+       is_local%wrap%flds_scalar_precip_factor = precip_fact(1)      
+       if (dbug_flag > 5) then
+          write(cvalue,*) precip_fact(1)
           call ESMF_LogWrite(trim(subname)//" precip_fact is "//trim(cvalue), ESMF_LOGMSG_INFO)
-
-          allocate(fldnames(3))
-          fldnames = (/'Faxa_rain', 'Faxa_snow', 'Fixx_rofi'/)
-          do n = 1,size(fldnames)
-             if (fldchk(is_local%wrap%FBExp(compice), trim(fldnames(n)), rc=rc)) then
-                call ESMF_FieldBundleGet(is_local%wrap%FBExp(compice), fieldname=trim(fldnames(n)), &
-                     field=lfield, rc=rc)
-                if (chkerr(rc,__LINE__,u_FILE_u)) return
-                call ESMF_FieldGet(lfield, farrayptr=dataptr1d, rc=rc)
-                if (chkerr(rc,__LINE__,u_FILE_u)) return
-                dataptr1d(:) = dataptr1d(:) * precip_fact
-             end if
-          end do
-          deallocate(fldnames)
        end if
+
+       ! Scale rain and snow to ice from atm by the precipitation factor received from the ocean
+       allocate(fldnames(3))
+       fldnames = (/'Faxa_rain', 'Faxa_snow', 'Fixx_rofi'/)
+       do n = 1,size(fldnames)
+          if (FB_fldchk(is_local%wrap%FBExp(compice), trim(fldnames(n)), rc=rc)) then
+             call FB_GetFldPtr(is_local%wrap%FBExp(compice), trim(fldnames(n)), dataptr, rc=rc)
+             if (ChkErr(rc,__LINE__,u_FILE_u)) return
+             dataptr(:) = dataptr(:) * is_local%wrap%flds_scalar_precip_factor
+          end if
+       end do
+       deallocate(fldnames)
     end if
 
     ! obtain nextsw_cday from atm if it is in the import state and send it to ice

mediator/med_phases_prep_ocn_mod.F90

@@ -69,7 +69,7 @@ subroutine med_phases_prep_ocn_accum(gcomp, rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
     ! auto merges to ocn
-    if (trim(coupling_mode) == 'cesm' .or. &
+    if ( trim(coupling_mode) == 'cesm' .or. &
          trim(coupling_mode) == 'nems_orig_data' .or. &
          trim(coupling_mode) == 'hafs') then
        call med_merge_auto(compocn, &
@@ -193,7 +193,8 @@ subroutine med_phases_prep_ocn_custom_cesm(gcomp, rc)
     ! custom calculations for cesm
     !---------------------------------------
 
-    use ESMF , only : ESMF_GridComp
+    use ESMF , only : ESMF_GridComp, ESMF_StateGet, ESMF_Field, ESMF_FieldGet
+    use ESMF , only : ESMF_VMBroadCast
     use ESMF , only : ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_SUCCESS
     use ESMF , only : ESMF_FAILURE,  ESMF_LOGMSG_ERROR
 
@@ -203,6 +204,7 @@ subroutine med_phases_prep_ocn_custom_cesm(gcomp, rc)
 
     ! local variables
     type(InternalState) :: is_local
+    type(ESMF_Field)    :: lfield
     real(R8), pointer   :: ifrac(:) => null()
     real(R8), pointer   :: ofrac(:) => null()
     real(R8), pointer   :: ifracr(:) => null()
@@ -227,20 +229,21 @@ subroutine med_phases_prep_ocn_custom_cesm(gcomp, rc)
     real(R8), pointer   :: Fioi_swpen_idf(:) => null()
     real(R8), pointer   :: Fioi_swpen(:) => null()
     real(R8), pointer   :: dataptr(:) => null()
-    real(R8), pointer   :: dataptr_o(:) => null()
+    real(R8), pointer   :: dataptr_scalar_ocn(:,:) => null()
     real(R8)            :: frac_sum
     real(R8)            :: ifrac_scaled, ofrac_scaled
     real(R8)            :: ifracr_scaled, ofracr_scaled
     logical             :: export_swnet_by_bands
     logical             :: import_swpen_by_bands
     logical             :: export_swnet_afracr
-    logical             :: first_precip_fact_call = .true.
-    real(R8)            :: precip_fact
+    real(R8)            :: precip_fact(1)
     character(CS)       :: cvalue
     real(R8)            :: fswabsv, fswabsi
+    integer             :: scalar_id
     integer             :: n
     integer             :: lsize
     real(R8)            :: c1,c2,c3,c4
+    logical             :: first_call = .true.
     character(len=64), allocatable :: fldnames(:)
     character(len=*), parameter    :: subname='(med_phases_prep_ocn_custom_cesm)'
     !---------------------------------------
@@ -359,8 +362,8 @@ subroutine med_phases_prep_ocn_custom_cesm(gcomp, rc)
           import_swpen_by_bands = .false.
        end if
 
-       ! Swnet without swpen from sea-ice
        if ( FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_swnet_afracr',rc=rc)) then
+          ! Swnet without swpen from sea-ice
           call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_swnet_afracr', Foxx_swnet_afracr, rc=rc)
           if (ChkErr(rc,__LINE__,u_FILE_u)) return
           export_swnet_afracr = .true.
@@ -416,14 +419,14 @@ subroutine med_phases_prep_ocn_custom_cesm(gcomp, rc)
 
        ! Output to ocean per ice thickness fraction and sw penetrating into ocean
        if ( FB_fldchk(is_local%wrap%FBExp(compocn), 'Sf_afrac', rc=rc)) then
-          call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Sf_afrac', fldptr1=dataptr_o, rc=rc)
+          call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Sf_afrac', fldptr1=dataptr, rc=rc)
           if (ChkErr(rc,__LINE__,u_FILE_u)) return
-          dataptr_o(:) = ofrac(:)
+          dataptr(:) = ofrac(:)
        end if
        if ( FB_fldchk(is_local%wrap%FBExp(compocn), 'Sf_afracr', rc=rc)) then
-          call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Sf_afracr', fldptr1=dataptr_o, rc=rc)
+          call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Sf_afracr', fldptr1=dataptr, rc=rc)
           if (ChkErr(rc,__LINE__,u_FILE_u)) return
-          dataptr_o(:) = ofracr(:)
+          dataptr(:) = ofracr(:)
        end if
 
     end if  ! if sea-ice is present
@@ -433,25 +436,47 @@ subroutine med_phases_prep_ocn_custom_cesm(gcomp, rc)
        deallocate(Foxx_swnet)
     end if
 
-    !---------------------------------------
-    ! application of precipitation factor from ocean
-    !---------------------------------------
-    precip_fact = 1.0_R8
-    if (precip_fact /= 1.0_R8) then
-       if (first_precip_fact_call .and. mastertask) then
-          write(logunit,'(a)')'(merge_to_ocn): Scaling rain, snow, liquid and ice runoff by precip_fact '
-          first_precip_fact_call = .false.
+    ! Apply precipitation factor from ocean (that scales atm rain and snow back to ocn ) if appropriate
+    if (trim(coupling_mode) == 'cesm' .and. is_local%wrap%flds_scalar_index_precip_factor /= 0) then
+
+       ! Note that in med_internal_mod.F90 all is_local%wrap%flds_scalar_index_precip_factor 
+       ! is initialized to 0.
+       ! In addition, in med.F90, if this attribute is not present as a mediator component attribute, 
+       ! it is set to 0. 
+       if (mastertask) then
+          call ESMF_StateGet(is_local%wrap%NstateImp(compocn), & 
+               itemName=trim(is_local%wrap%flds_scalar_name), field=lfield, rc=rc)
+          if (chkerr(rc,__LINE__,u_FILE_u)) return
+          call ESMF_FieldGet(lfield, farrayPtr=dataptr_scalar_ocn, rc=rc)
+          if (chkerr(rc,__LINE__,u_FILE_u)) return
+          scalar_id=is_local%wrap%flds_scalar_index_precip_factor
+          precip_fact(1) = dataptr_scalar_ocn(scalar_id,1)
+          if (first_call) then
+             write(logunit,'(a)')'(merge_to_ocn): Scaling rain, snow, liquid and ice runoff by precip_fact from ocn'
+             first_call = .false.
+          end if
+          if (precip_fact(1) /= 1._r8) then
+             write(logunit,'(a,f21.13)')&
+                  '(merge_to_ocn): Scaling rain, snow, liquid and ice runoff by non-unity precip_fact ',&
+                  precip_fact(1)
+          end if
+       end if
+       call ESMF_VMBroadCast(is_local%wrap%vm, precip_fact, 1, 0, rc=rc)
+       if (chkerr(rc,__LINE__,u_FILE_u)) return
+       is_local%wrap%flds_scalar_precip_factor = precip_fact(1)      
+       if (dbug_flag > 5) then
+          write(cvalue,*) precip_fact(1)
+          call ESMF_LogWrite(trim(subname)//" precip_fact is "//trim(cvalue), ESMF_LOGMSG_INFO)
        end if
-       write(cvalue,*) precip_fact
-       call ESMF_LogWrite(trim(subname)//" precip_fact is "//trim(cvalue), ESMF_LOGMSG_INFO)
 
+       ! Scale rain and snow to ocn from atm by the precipitation factor received from the ocean
        allocate(fldnames(4))
-       fldnames = (/'Faxa_rain','Faxa_snow', 'Foxx_rofl', 'Foxx_rofi'/)
+       fldnames = (/'Faxa_rain', 'Faxa_snow', 'Foxx_rofl', 'Foxx_rofi'/)
        do n = 1,size(fldnames)
           if (FB_fldchk(is_local%wrap%FBExp(compocn), trim(fldnames(n)), rc=rc)) then
              call FB_GetFldPtr(is_local%wrap%FBExp(compocn), trim(fldnames(n)) , dataptr, rc=rc)
              if (ChkErr(rc,__LINE__,u_FILE_u)) return
-             dataptr(:) = dataptr(:) * precip_fact
+             dataptr(:) = dataptr(:) * is_local%wrap%flds_scalar_precip_factor
           end if
        end do
        deallocate(fldnames)