Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BLING "spring cleaning" #778

Open
wants to merge 15 commits into
base: master
Choose a base branch
from
Open

BLING "spring cleaning" #778

wants to merge 15 commits into from

Conversation

averdy
Copy link
Contributor

@averdy averdy commented Oct 3, 2023

What changes does this PR introduce?

  1. Bug fix: add proper initialization of biomass for the 3 phytoplankton groups
  2. New feature: add calculation of chlorophyll concentration as detected by satellite (2D field).

What is the current behaviour?

  1. Unrealistic phytoplankton biomass values at the start of the simulation when using ADVECT_PHYTO flag
  2. Chlorophyll concentration in the top layer is not directly comparable with satellite observations

What is the new behaviour

  1. Realistic phytoplankton biomass values at the start of the simulation
  2. New diagnostic BLGCHLSA is calculated at 13:30 local time based on model chlorophyll concentration in the upper ocean and light attenuation, similarly to how satellites observe ocean color

Does this PR introduce a breaking change?

No

Other information:

Suggested addition to tag-index

(To avoid unnecessary merge conflicts, please don't update tag-index. One of the admins will do that when merging your pull request.)

@jm-c jm-c requested a review from mmazloff October 24, 2023 15:37
@mmazloff
Copy link
Collaborator

This looks great to me!

@jm-c jm-c self-requested a review November 24, 2023 14:55
Copy link
Member

@jm-c jm-c left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I took a quick look, and this looks good. I have 3 comments:

  1. looks like TAF is getting 3 recomputation warnings now, in AD experiment global_oce_biogeo_bling. Now the gradient check did not change at all but seeing some changes to AD-Monitor for adptracer03 and adptracer06 (2 digit match only).
  2. it would be good to have the "local time refinement" to also work when pkg/cal is compiled but not used (useCAL=F) as well as when it's not even compiled, since it's not too difficult to change (see comment below).
  3. the time-window boundaries (12:00 and 13:30) could be changed to parameter, if not run-time params at least fortran PARAMETER (locally in bling_light.F). This could help to follow what the code is doing and also in case this time-window needs to be changed.

C Satellite chlorophyll
#ifdef ALLOW_CAL
c mydate is utc time
diffutc = XC(i,j,bi,bj)/360.*24
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would not be difficult to have something equivalent for the case useCAL = F. I could suggest some changes in this routine if you like.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure what the equivalent of the variable localtime is when pkg/cal is not used. Suggestions appreciated!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In addition, the code only makes sense, if XC holds the longitude, i.e. not for cartesian coordinates. Do you want to catch that case?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ideally yes! And anything else that I've missed, because it's not in the setup I'm used to :)

@@ -1895,6 +1899,8 @@ SUBROUTINE BLING_BIO_NITROGEN(
CALL DIAGNOSTICS_FILL(N_reminp, 'BLGNREM ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(P_reminp, 'BLGPREM ',0,Nr,2,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(no3_adj, 'BLGNONEN',0,Nr,2,bi,bj,myThid)
C 2d global variables
CALL DIAGNOSTICS_FILL(chl_sat, 'BLGCHLSA',0,1,1,bi,bj,myThid)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

An alternative would be to move this DIAGNOSTICS_FILL call directly inside bling_light.F where it's computed.
And I did not check carefully, but if variable chl_sat is just used for diagnostics, then it could be changed to a local variable inside bling_light.F and removed from common block. Now this might not be a good idea if you have in mind of computing some cost function from it.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm moving the diagnostic_fill call inside bling_light.F as suggested.
chl_sat is going to be used in the cost function, so I'm going to keep it as a global variable.

Comment on lines 924 to 927
Phy_lg_local(i,j,k) = max( epsln, Phy_lg_local(i,j,k) )
Phy_sm_local(i,j,k) = max( epsln, Phy_sm_local(i,j,k) )
Phy_diaz_local(i,j,k) = max( epsln, Phy_diaz_local(i,j,k) )

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To avoid the new recomputation warnings (the associated recomputations make the run much longer: ALL increases from 128s to 434s, because ADTHE_MAIN_LOOP increases from 52s to 358s) can be avoided without store directives, see below, but most of this routine is actually recomputed in the AD-code and one may want to fix that.

Suggested change
Phy_lg_local(i,j,k) = max( epsln, Phy_lg_local(i,j,k) )
Phy_sm_local(i,j,k) = max( epsln, Phy_sm_local(i,j,k) )
Phy_diaz_local(i,j,k) = max( epsln, Phy_diaz_local(i,j,k) )
#ifdef ALLOW_AUTODIFF_TAMC
C avoid recomputation warnings
ENDIF
ENDDO
ENDDO
ENDDO
DO k=1,Nr
DO j=jmin,jmax
DO i=imin,imax
IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
#endif
Phy_lg_local(i,j,k) = max( epsln, Phy_lg_local(i,j,k) )
Phy_sm_local(i,j,k) = max( epsln, Phy_sm_local(i,j,k) )
Phy_diaz_local(i,j,k) = max( epsln, Phy_diaz_local(i,j,k) )

@mjlosch
Copy link
Member

mjlosch commented Nov 27, 2023

@averdy

  1. I went through bling_bio_nitrogen.F and adjusted store directives and reorganised the code a little to make it more straightforward for TAF. I have now a version that uses fewer store directives (21 instead of 22) and the recomputations seem to be at bay meaning that I think we avoid the expensive part and still do some recomputations where they do not take more than two lines of code, basically just one statement. Shall I push this version of bling_bio_nitrogen.F it, so that you can see if you like it?

  2. I also noticed this code:

C-----------------------------------------------------------
C  avoid negative nutrient concentrations that can result from
C  advection when low concentrations
#ifdef BLING_NO_NEG
      CALL BLING_MIN_VAL(PTR_O2,  1. _d -11, o2_adj, bi, bj)
[...]
#endif

It's not clear to me why this needs to be "hidden from the adjoint" (as stated the routine BLING_MIN_VAL). It would make the adjoint more accurate, if you don't hide it. I think this can be replaced by just capping the input, and avoiding the extra *_adj variables altogether, but I may not understand the meaning of this entirely.

  1. Also there are some fields that are purely diagnostic, e.g. theta_Fe_inv, I suggest to move them to the diagnostics section, so that they are only computed when diagnostics are used (or even when the particular diagnostic is active).

  2. The long "Remineralisation"-loop could be improved/simplified by turning PONflux_u,POPflux_u, PFEflux_u, CaCO3flux_u, SIflux_u into 2D fields, but that's not a priority here.

@averdy
Copy link
Contributor Author

averdy commented Nov 29, 2023

@jm-c Thank you for the feedback - great suggestions! I’m working on it.

@averdy
Copy link
Contributor Author

averdy commented Nov 29, 2023

@mjlosch Thank you for your help! Yes, if you can push your version of bling_bio_nitrogen.F, that would be great. I'm working on addressing your other suggestions.

This change speeds up ADTHE_MAIN_LOOP in
verification/global_oce_bio_bling by a factor of 7 (!), mainly because
it fixes the extensive recomputations. Relative to the version prior
to this PR the speed up is small.

With this commit (on my laptop):
User/System/Wall clock time:   51.93 sec / 0.30 sec /  52.59 sec
Before:
User/System/Wall clock time:  358.49 sec / 0.35 sec / 360.67 sec
@averdy
Copy link
Contributor Author

averdy commented Mar 26, 2024

@jm-c Regarding your comment "1" above, I think Martin's edits will have solved the recomputations. Should I do something about "changes to AD-Monitor for adptracer03 and adptracer06"?

In response to comment "3" I'm making the time-window boundaries run-time parameters.

@averdy
Copy link
Contributor Author

averdy commented Mar 26, 2024

@mjlosch Regarding your comment "2": the "adj" variables are meant to keep track for the adjustments for budget calculations. If the bling_min_val code doesn't need to be hidden from the adjoint, should I remove it from bling_ad.flow and bling_ad_diff.list?

I'm moving the purely diagnostic fields to the diagnostics section as suggested in comment "3".

@mjlosch
Copy link
Member

mjlosch commented Mar 26, 2024

about comment 2: yes, you could do that, but that will give you recomputation warnings that need to be dealt with via store directives. The question is if you really need that. I can push something, if you like.

Also in bling_ad_diff.list, bling_read_pickup.f (there are flow directives) and bling_write_pickup.f (not called directly) can be removed.

@jm-c
Copy link
Member

jm-c commented Apr 4, 2024

@averdy I saw the 2 new parameters chlsat_tbegin and chlsat_tend, but looks like they are not set in bling_readparms.F (neither default value nor in namelist).

mjlosch and others added 4 commits April 5, 2024 09:47
- remove bling_write_pickup.f from bling_ad_diff.list because the call
is not visible to taf anyway
- remove bling_read_pickup.f from bling_ad_diff.list because there are
flow directives (avoid another warning)
@jm-c
Copy link
Member

jm-c commented Apr 8, 2024

@averdy I am going to revert your last commit, because it "undo" all the changes that have been added since last November (and, as a result, this PR now affects 273 files). I went back 1 commit on my clone and I've checked that I have
all the previous changes, including (a) your latest addition to bling_readparms.F, (b) Martin TAF related adjustment, and all the recently merged PR.
Note that If I do this, it will likely make less easy for you to update your working "bling_cleaning" branch and you may want to consider removing this branch from your local clone (but not from your fork !) and checking-out a new local copy of this branch (after updating your clone with a "git fetch").

@jm-c
Copy link
Member

jm-c commented Apr 15, 2024

I remove the last commit by doing:

git reset --hard HEAD~1
git push -f averdy bling_cleanup

so I used a "forced push", and regarding updates, if you already have a local copy of this branch, using simply "git pull" might not work well. For this case, either (a) removing this branch from your local copy ("git branch -D bling_cleanup") and getting it back after a "git fetch", as mentionned earlier, or (b) with "git fetch [remote] ; git checkout bling_cleanup ; git restore . " might do it too.

@averdy
Copy link
Contributor Author

averdy commented Apr 15, 2024

Thank you @jm-c! I've updated my local copy.

@mmazloff
Copy link
Collaborator

LGTM

@jm-c
Copy link
Member

jm-c commented May 21, 2024

@averdy Sorry for waiting so long.
I would like to rename (+ change units) the 2 new parameters chlsat_tbegin and chlsat_tend and make the little changes to allow to use this code without pkg/cal.
I suggest "chlsat_locTimWindow(2)" (a bit long but not used in too many places in the code), in decimal hours,
to be even more clear about local-time (as almost everywhere else time is UTC).
And since these 2 parameters are new (not even there when this PR was submitted) it should not break any set-up (I hope).
If you don't like this please let us know.

@averdy
Copy link
Contributor Author

averdy commented May 21, 2024

@jm-c That sounds great!

@jm-c
Copy link
Member

jm-c commented May 21, 2024

Also, as I am looking at the details of this local-time window, it seems that the current version is not very accurate:
if it's noon at long=0, I will collect in "chl_sat" the chlorophyll from XC=0 to XC+151.3 (instead of XC+151.5, for a one hour and half longitude band).

@jm-c
Copy link
Member

jm-c commented May 21, 2024

so my little changes are going to fix also this inaccuracy.

jm-c added 4 commits May 21, 2024 16:10
- rename the 2 new params to "chlsat_locTimWindow(2)" (in decimal hour)
- use "decimal hour" consistently (this fix inaccurate time-window criteria).
- allow to be used without pkg/cal
just to avoid un-used variables
Copy link
Member

@jm-c jm-c left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM. I run fwd, adm TAF and Tapenade testreport for exp. global_oce_biogeo_bling and all look good. I see minor changes to AD-monitor for the TAF one but I guess it's expected given the changes in store directives. @mjlosch do you agree ?

@jm-c
Copy link
Member

jm-c commented May 22, 2024

@averdy From my side, this PR is good to be merged. So, when you have time, please check/review my set of recent changes and let us know if you aggree or not.

@averdy
Copy link
Contributor Author

averdy commented May 22, 2024

@jm-c It looks great!! Thank you so much for improving this code!

@mjlosch
Copy link
Member

mjlosch commented May 22, 2024

I run fwd, adm TAF and Tapenade testreport for exp. global_oce_biogeo_bling and all look good. I see minor changes to AD-monitor for the TAF one but I guess it's expected given the changes in store directives. @mjlosch do you agree ?

I agree, this looks like the usual round-off issue with recomputation vs restoring from tape (common block). Nothing to worry about.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants