macrozooplankton mortality modified by light and sea ice

[kristenkrumhardt]

We would like to some zooplankton mortality rates to be modified by surface shortwave radiation and sea ice fraction. This would be helpful for simulating the sea ice dependence and diapause during winter of antarctic krill (macrozooplankton).

[mnlevy1981]

@kristenkrumhardt and I talked about this a bit earlier today and the plan we came up with is

1. Add mort_multiplier_opt to zooplankton_settings_type with valid values of 'constant_one' and 'sw_and_ice_dependent'; there will also be an iopt to avoid string comparisons in logic, but I will continue to use string comparisons in this comment for clarity. The idea is that we will replace z_mort_0 and z_mort2_0 in the zoo_loss computation with z_mort_multiplier * z_mort_0 and z_mort2_multiplier * z_mort2_0, respectively; the two _multiplier variables will be 1 if zooplankton_settings(zoo_ind)%mort_multiplier_opt == 'constant_one' and some function of SW and ice frac if zooplankton_settings(zoo_ind)%mort_multiplier_opt == 'sw_and_ice_dependent'.
2. If zooplankton_settings(zoo_ind)%mort_multiplier_opt == 'sw_and_ice_dependent' for any zoo_ind, then we'll need to add ice fraction in interior_tendency_forcings (I believe we finalize all settings values before determining what forcings are necessary so this shouldn't be too difficult)
3. Currently we are computing zoo_loss in compute_Zprime(), though we could address #335 and create compute_zoo_loss(); regardless, the subroutine responsible for computing zoo_loss will need interior_tendency_forcings and PAR%col_frac as arguments
4. Once the subroutine computing zoo_loss can access shortwave and ice fraction, it will should be easy to compute the _multiplier terms as described in (1) above

That's just a rough roadmap we came up with from walking through the existing code, I'm certainly open to suggestions. The variable names / values are not set in stone, but that's what I plan to use in the first pass. And it's a minor point, but I'm on the fence about whether or not we should roll in the #335 fix now that there will be more logic behind computing zoo_loss.

[mnlevy1981]

Slight change of plans from what's outlined above -- rather than make local variables for the _multipliers, I just make local copies of z_mort_0 and z_mort2_0 that get defined properly based on zooplankton_settings(zoo_ind)%mort_coeff_opt (which can be always_the_same or sw_and_ice_dependent).

The code snippet below comes from my enhancement/new_zoo_mortality_func branch, which currently keeps z_mort_0 and z_mort2_0 unchanged even with the sw_and_ice_dependent option but at least has placeholders to show where changes need to be made.

        select case (zooplankton_settings(zoo_ind)%mort_coeff_iopt)
          case (mort_coeff_iopt_sw_and_ice)
            associate(&
              surf_sw => interior_tendency_forcings(interior_tendency_forcing_ind%surf_shortwave_id)%field_1d(1,:), &
              ice_frac => interior_tendency_forcings(interior_tendency_forcing_ind%ifrac_id)%field_0d(:) &
              )
              ! TODO: Determine relationship between surf_sw, ice_frac, and mort rates
              ! if ((sum(surf_sw).le.1._r8) .and. (ice_frac(1).gt.0.9_r8)) then
              !   z_mort_0 = 0.85_r8*zooplankton_settings(zoo_ind)%z_mort_0
              !   z_mort2_0 = 0.85_r8*zooplankton_settings(zoo_ind)%z_mort2_0
              ! else
              !   z_mort_0 = zooplankton_settings(zoo_ind)%z_mort_0
              !   z_mort2_0 = zooplankton_settings(zoo_ind)%z_mort2_0
              ! end if
              z_mort_0 = zooplankton_settings(zoo_ind)%z_mort_0
              z_mort2_0 = zooplankton_settings(zoo_ind)%z_mort2_0
            end associate
          case DEFAULT
            z_mort_0 = zooplankton_settings(zoo_ind)%z_mort_0
            z_mort2_0 = zooplankton_settings(zoo_ind)%z_mort2_0
        end select

        zoo_loss(zoo_ind,:) = (z_mort2_0 * Zprime(zoo_ind,:)**zoo_mort2_exp &
                               + z_mort_0 * Zprime(zoo_ind,:)) * Tfunc_zoo(zoo_ind,:)

@kristenkrumhardt after Cheyenne comes back up, let's talk about about how to set up a sandbox for you to to make changes in. I think adding ice fraction to the list of interior tendency forcings will be handled fine by POP since it already needs to know what to do with that forcing for the surface fluxes, so hopefully you just need to uncomment and modify the if block (and delete the two lines below it) and then you can open a PR. If POP does need to be updated, I can tackle that pretty quickly.

[matt-long]

Just a word of caution here: I think it makes sense to explore the effect of this parameterization before totally dialing in the technical details of the implementation. In that sense, I would proceed with augmenting the interface to pass the requisite fields (i.e., PAR, IFRAC), then punt on everything else until we think we know what actually works.

[mnlevy1981]

@matt-long we are passing interior_tendency_forcings(:) so PAR is already available, we just need to update the associate block to add easy access to it.

Another thing to add to the associate statement is the PAR column fractions -- in my commented out example if statement, sum(surf_sw) is definitely NOT what we want to base our conditional on. If shortwave is a proxy for time-of-year, then surf_sw(1) [surf_sw is an array with dimension of num_PAR_subcols; I believe subcolumn 1 is "open ocean"] is correct otherwise we want sum(surf_sw * col_frac) to get the total incoming shortwave... I'll update the comments on my branch and then @kristenkrumhardt can explore and we can clean everything up at the end.

[matt-long]

I think PAR/SW may not be the best thing to use: what we really want is day length or zenith angle.

[kristenkrumhardt]

Agreed. So we actually need 1) day of year and 2) latitude (which we have already read into MARBL from POP for the great calcite belt/alk project) . After doing more reading on the subject, studies that report on krill diapause during winter say that it's triggered by day length (not a mean daily light level).

[mnlevy1981]

Just to make sure I'm understanding this right - are the two options

(1) Adding zenith angle as a forcing field (if this option is turned on), or (2) Add latitude to marbl_domain_type and also pass day of year to MARBL

If so, I think (1) would be easier -- for (2) I could see having "day of year" as an optional argument to marbl_instances%interior_tendency_compute(), but if we would just be calculating something similar to zenith angle based on day of year & latitude anyway, I'd just as soon let the forcing type manage it.

I think this forcing update would be tricky to figure out in the stand-alone driver (I might need to recreate the initial condition / forcing netCDF file to include zenith angle for each column or have the driver itself compute zenith angle assuming we are running on Jan 1), but it should be straightforward to have POP do.

[kristenkrumhardt]

I think option 1 would be completely fine. I just to figure out at what zenith angle to trigger the diapause.. I'll work on getting some studies together that have latitude and day length when krill diapause is triggered, then we can calculate it...

[matt-long]

Can we assume that zenith angle is readily available to pass to MARBL in other models? I think if it's straightforward in POP, we should start there.

[mnlevy1981]

POP has a compute_cosz() function inforcing_coupled.F90, though I think we might want to call shr_orb_decl() and shr_orb_cosz() directly from the MARBL driver? I could be way offbase, but my thought is "we're using cosz as a stand-in for "number of hours per day with sunlight", so we probably want to use the value at local noon rather than compute it in a manner that is dependent on longitude?

So there are details to work out, but it is definitely available in POP already.

[matt-long]

What about MOM?

We might want to compute day length in the driver. cosz computed at the timestep level is not useful; rather we want "hours of day when cosz > 0" or similar. Not sure what the best way to proceed here is.

[kristenkrumhardt]

This linear relationship between day length and krill metabolic activity (using O2 consumption as an indicator of metabolic activity) is the clearest relationship that I've found in the literature so far (from a review article by B. Meyer (2012) "The overwintering of Antarctic krill, Euphausia superba, from an ecophysiological perspective", in Polar Biology).

[mnlevy1981]

I haven't looked for similar functions in MOM yet, but once we settle on a formula I'll either search through the code or ask @ashao . From Wikipedia, the most trustworthy of sources, it looks like the length of the day can be computed using one of the share-code routines I mentioned above:

delta = shr_orb_decl(calday, orb_eccen, orb_mvelpp, orb_lambm0,orb_obliqr, delta, eccf)
hours_of_sunlight = acos(max(min(-tan(lat)*tan(delta), 1), -1))*24/pi

Obviously I'll play with that a little bit to make sure I haven't messed up the sign or the min / max, but I could see adding a compute_hours_of_sunlight() function in forcing_coupled.F90 (where we already have a call to shr_orb_decl so presumably it would be easy to access the necessary arguments) that would be used to populate an array (SUN_HOURS?) of dimension (nx_block, ny_block, max_blocks_clinic) in the ecosystem driver based on TLAT.

[matt-long]

Here's gist with a purported day-length computation: https://gist.github.com/anttilipp/ed3ab35258c7636d87de6499475301ce

@kristenkrumhardt, that data is beautiful!