Fix lcompute_nhx_surface_emis

Only call marbl_comp_nhx_surface_emis if lcompute_nhx_surface_emis ==
.true.. Also, set all components of marbl_surface_forcing_internal_type
to 0 at construction (if lcompute_nhx_surface_emis == .false. then
nhx_surface_emis is never set; the same may be true for other components
as well).

All POP tests run with lcompute_nhx_surface_emis = .true., so this
should not change answers in our test suite... but I think we should add
a lcompute_nhx_surface_emis = .false. test to POP.

src/marbl_interface_private_types.F90

@@ -466,27 +466,47 @@ end subroutine marbl_PAR_destructor
 
   subroutine marbl_surface_forcing_internal_constructor(this, num_elements)
 
+    use marbl_constants_mod, only : c0
+
     class(marbl_surface_forcing_internal_type), intent(out) :: this
     integer (int_kind),                         intent(in)  :: num_elements
 
     allocate(this%piston_velocity (num_elements))
+    this%piston_velocity  = c0
     allocate(this%flux_co2        (num_elements))
+    this%flux_co2         = c0
     allocate(this%flux_alt_co2    (num_elements))
+    this%flux_alt_co2     = c0
     allocate(this%co2star         (num_elements))
+    this%co2star          = c0
     allocate(this%dco2star        (num_elements))
+    this%dco2star         = c0
     allocate(this%pco2surf        (num_elements))
+    this%pco2surf         = c0
     allocate(this%dpco2           (num_elements))
+    this%dpco2            = c0
     allocate(this%co3             (num_elements))
+    this%co3              = c0
     allocate(this%co2star_alt     (num_elements))
+    this%co2star_alt      = c0
     allocate(this%dco2star_alt    (num_elements))
+    this%dco2star_alt     = c0
     allocate(this%pco2surf_alt    (num_elements))
+    this%pco2surf_alt     = c0
     allocate(this%dpco2_alt       (num_elements))
+    this%dpco2_alt        = c0
     allocate(this%schmidt_co2     (num_elements))
+    this%schmidt_co2      = c0
     allocate(this%schmidt_o2      (num_elements))
+    this%schmidt_o2       = c0
     allocate(this%pv_o2           (num_elements))
+    this%pv_o2            = c0
     allocate(this%pv_co2          (num_elements))
+    this%pv_co2           = c0
     allocate(this%o2sat           (num_elements))
+    this%o2sat            = c0
     allocate(this%nhx_surface_emis(num_elements))
+    this%nhx_surface_emis = c0
 
   end subroutine marbl_surface_forcing_internal_constructor
 

src/marbl_mod.F90

@@ -2089,22 +2089,22 @@ subroutine marbl_set_surface_forcing( &
     !  FIXME: ensure that ph is computed, even if lflux_gas_co2=.false.
     !-----------------------------------------------------------------------
 
-    call marbl_comp_nhx_surface_emis(                &
-         num_elements     = num_elements,            &
-         nh4              = surface_vals(:,nh4_ind), &
-         ph               = ph_prev_surf,            &
-         sst              = sst,                     &
-         sss              = sss,                     &
-         u10_sqr          = u10_sqr,                 &
-         atmpres          = ap_used,                 &
-         ifrac            = ifrac,                   &
-         nhx_surface_emis = nhx_surface_emis)
-
-    if (sfo_ind%flux_nhx_id.ne.0) then
-       surface_forcing_output%sfo(sfo_ind%flux_nhx_id)%forcing_field = nhx_surface_emis
-    end if
-
     if (lcompute_nhx_surface_emis) then
+      call marbl_comp_nhx_surface_emis(                &
+           num_elements     = num_elements,            &
+           nh4              = surface_vals(:,nh4_ind), &
+           ph               = ph_prev_surf,            &
+           sst              = sst,                     &
+           sss              = sss,                     &
+           u10_sqr          = u10_sqr,                 &
+           atmpres          = ap_used,                 &
+           ifrac            = ifrac,                   &
+           nhx_surface_emis = nhx_surface_emis)
+
+      if (sfo_ind%flux_nhx_id.ne.0) then
+         surface_forcing_output%sfo(sfo_ind%flux_nhx_id)%forcing_field = nhx_surface_emis
+      end if
+
       stf(:, nh4_ind) = stf(:, nh4_ind) - nhx_surface_emis(:)
     endif