search

NOAA-CEFI-Regional-Ocean-Modeling / ocean_BGC

3 stars 4 forks source link

Divide by zero in iceberg sources #5

Closed andrew-c-ross closed 5 months ago

andrew-c-ross commented 5 months ago

The 4P version of COBALT has the following lines here: 

cobalt%jfe_iceberg(i,j,1) = cobalt%jfe_iceberg_ratio*max(frunoff(i,j),0.0)/rho_dzt(i,j,1)
cobalt%jno3_iceberg(i,j,1) = cobalt%jno3_iceberg_ratio*max(frunoff(i,j),0.0)/rho_dzt(i,j,1)
cobalt%jpo4_iceberg(i,j,1) = cobalt%jpo4_iceberg_ratio*max(frunoff(i,j),0.0)/rho_dzt(i,j,1)

The problem is that rho_dzt(:, :, 1) (density times layer thickness) is 0 over land.

This block probably needs to be wrapped in a check for the index of the bottom layer like line 9788:

if (grid_kmt(i,j) .gt. 0) then

I don't think this affects the solution since it happens over land, but debugging other issues (like I'm currently trying to do) is harder since the model aborts here when it is compiled in debug mode.

yichengt900 commented 5 months ago

Thanks, @andrew-c-ross, for reporting this bug. I can open a pull request to hotfix this issue if you prefer.