Processing MARBL Output

Surface Flux Output

After calling surface_flux_compute(), the GCM needs to manage all the MARBL output. The block of code below shows how POP accesses surface_flux_saved_state, surface_flux_output, surface_fluxes, surface_flux_diags, and glo_avg_fields_surface_flux. Note that POP aborts if any values in surface_fluxes are NaN.

!-----------------------------------------------------------------------
! Copy data from marbl output column to pop slab data structure
!-----------------------------------------------------------------------

do index_marbl = 1, marbl_col_cnt(iblock)
   i = marbl_col_to_pop_i(index_marbl,iblock)
   j = marbl_col_to_pop_j(index_marbl,iblock)

   do n=1,size(surface_flux_saved_state)
     surface_flux_saved_state(n)%field_2d(i,j,iblock) = &
       marbl_instances(iblock)%surface_flux_saved_state%state(n)%field_2d(index_marbl)
   end do

   do n=1,sfo_cnt
     surface_flux_outputs(i,j,iblock,n) = &
        marbl_instances(iblock)%surface_flux_output%sfo(n)%forcing_field(index_marbl)
   end do

   !-----------------------------------------------------------
   ! before copying surface fluxes, check to see if any are NaNs
   !-----------------------------------------------------------

   if (any(shr_infnan_isnan(marbl_instances(iblock)%surface_fluxes(index_marbl,:)))) then
      write(stdout, *) subname, ': NaN in stf_module, (i,j)=(', &
         this_block%i_glob(i), ',', this_block%j_glob(j), ')'
      write(stdout, *) '(lon,lat)=(', TLOND(i,j,iblock), ',', TLATD(i,j,iblock), ')'
      do n = 1, ecosys_tracer_cnt
         write(stdout, *) trim(marbl_instances(1)%tracer_metadata(n)%short_name), ' ', &
            marbl_instances(iblock)%tracers_at_surface(index_marbl,n), ' ', &
            marbl_instances(iblock)%surface_fluxes(index_marbl,n)
      end do
      do n = 1, size(surface_flux_forcings)
         associate (forcing_field => surface_flux_forcings(n))
            write(stdout, *) trim(forcing_field%metadata%marbl_varname)
            if (forcing_field%rank == 2) then
               write(stdout, *) forcing_field%field_0d(i,j,iblock)
            else
               write(stdout, *) forcing_field%field_1d(i,j,:,iblock)
            end if
         end associate
      end do
      call exit_POP(sigAbort, 'Stopping in ' // subname)
   end if

   do n = 1,ecosys_tracer_cnt
      stf_module(i,j,n) = &
           marbl_instances(iblock)%surface_fluxes(index_marbl,n)
   end do

   do n=1,size(marbl_instances(1)%surface_flux_diags%diags)
      surface_flux_diags(i,j,n,iblock) = &
           marbl_instances(iblock)%surface_flux_diags%diags(n)%field_2d(index_marbl)
   end do

   ! copy values to be used in computing requested global averages
   ! arrays have zero extent if none are requested
   glo_avg_fields_surface(i,j,iblock,:) = marbl_instances(iblock)%glo_avg_fields_surface_flux(index_marbl,:)
end do

Interior Tendency Output

After calling interior_tendency_compute(), the GCM needs to manage all the MARBL output. The block of code below shows how POP accesses interior_tendency_saved_state, interior_tendencies and glo_avg_fields_interior_tendency. Diagnostics are accumulated in ecosys_tavg_accumulate_interior(), which is a wrapper to handle column-wise accumulation (POP typically accumulates level by level). Note that POP aborts if any values in interior_tendencies are NaN.

do c = this_block%jb,this_block%je
   do i = this_block%ib,this_block%ie

      if (land_mask(i,c,bid)) then

         ! [snipped prestage calls]

         !-----------------------------------------------------------
         ! copy marbl column data back to slab
         !-----------------------------------------------------------

         call timer_start(ecosys_interior_marbl_to_pop, block_id=bid)

         do n=1,size(interior_tendency_saved_state)
           interior_tendency_saved_state(n)%field_3d(:,i,c,bid) =               &
             marbl_instances(bid)%interior_tendency_saved_state%state(n)%field_3d(:,1)
         end do

         !-----------------------------------------------------------
         ! before copying tendencies, check to see if any are NaNs
         !-----------------------------------------------------------

         do k = 1, KMT(i, c, bid)
            if (any(shr_infnan_isnan(marbl_instances(bid)%interior_tendencies(:, k)))) then
               write(stdout, *) subname, ': NaN in dtracer_module, (i,j,k)=(', &
                  this_block%i_glob(i), ',', this_block%j_glob(c), ',', k, ')'
               write(stdout, *) '(lon,lat)=(', TLOND(i,c,bid), ',', TLATD(i,c,bid), ')'
               do n = 1, ecosys_tracer_cnt
                  write(stdout, *) trim(marbl_instances(1)%tracer_metadata(n)%short_name), ' ', &
                     marbl_instances(bid)%tracers(n, k), ' ', &
                     marbl_instances(bid)%interior_tendencies(n, k)
               end do
               do n = 1, size(interior_tendency_forcings)
                  associate (forcing_field => interior_tendency_forcings(n))
                     write(stdout, *) trim(forcing_field%metadata%marbl_varname)
                     if (forcing_field%rank == 2) then
                        write(stdout, *) forcing_field%field_0d(i,c,bid)
                     else
                        if (forcing_field%ldim3_is_depth) then
                           write(stdout, *) forcing_field%field_1d(i,c,k,bid)
                        else
                           write(stdout, *) forcing_field%field_1d(i,c,:,bid)
                        end if
                     end if
                  end associate
               end do
               call exit_POP(sigAbort, 'Stopping in ' // subname)
            end if
         end do

         do n = 1, ecosys_tracer_cnt
            dtracer_module(i, c, 1:KMT(i, c, bid), n) = marbl_instances(bid)%interior_tendencies(n, 1:KMT(i, c, bid))
         end do

         ! copy values to be used in computing requested global averages
         ! arrays have zero extent if none are requested
         glo_avg_fields_interior(i, c, bid, :) = marbl_instances(bid)%glo_avg_fields_interior_tendency(:)
         call timer_stop(ecosys_interior_marbl_to_pop, block_id=bid)

         !-----------------------------------------------------------
         ! Update pop tavg diags
         !-----------------------------------------------------------

         call timer_start(ecosys_interior_marbl_tavg, block_id=bid)

         call ecosys_tavg_accumulate_interior(i, c, marbl_instances(bid), bid)

         call timer_stop(ecosys_interior_marbl_tavg, block_id=bid)

      end if ! end if land_mask > 0

   end do ! do i
end do ! do c