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.

Sign up for GitHub

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

Jump to bottom

MOM5 updates from ACCESS-ESM1.5 #381

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

MOM5 updates from ACCESS-ESM1.5 #381

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6ffcdac
Add files via upload
matt-csiro May 23, 2023
252ce49
Test upload
matt-csiro May 23, 2023
c9f4ce2
WOMBAT update from ESM
matt-csiro May 25, 2023
abc9060
Update to compile ESM.
matt-csiro May 25, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions bin/environs.nci
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,3 +4,4 @@ module load intel-compiler/2019.5.281
module load netcdf/4.7.4
module load openmpi/4.0.2
setenv mpirunCommand "mpirun --mca orte_base_help_aggregate 0 -np"
setenv OASIS_ROOT /g/data/p66/pbd562/test/t47-hxw/jan20/4.0.2/oasis3-mct
Copy link
Collaborator

@aekiss aekiss May 30, 2023 •
edited

Choose a reason for hiding this comment

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

not sure we want to hardwire a specific path like this, which is only useable for users with membership of p66

Copy link
Author

Choose a reason for hiding this comment

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

I'm needing this setting otherwise the compile for the ESM fails, for instance...

.../src/accesscm_coupler/mom_oasis3_interface.F90(79): error #7002: Error in opening the compiled module file. Check INCLUDE paths. [MOD_PRISM]
use mod_prism
----^

Let me know if there are alternatives to try.

2 changes: 2 additions & 0 deletions exp/compile_esm.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
./MOM_compile.csh --platform=nci --type=ACCESS-ESM

18 changes: 16 additions & 2 deletions src/mom5/ocean_csiro_bgc/bio_v3.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ integer :: index_temp, index_salt

integer :: ts_npzd ! number of time steps within NPZD model
integer :: kmeuph = 20 ! deepest level of euphotic zone
integer :: k100 = 10 ! deepest level less than 100 m
integer :: tn, trn

real :: pi = 3.14159265358979 !yes, this is pi
Expand Down Expand Up @@ -96,8 +97,10 @@ integer :: index_temp, index_salt
! set the maximum index for euphotic depth
do k=1,grid%nk
if (grid%zw(k).le. 400) kmeuph=k
if (grid%zw(k).le. 100) k100=k
enddo
! print*,'rjm euphotic ', kmeuph
! print*,'mac k100 ', k100

!
!
Expand Down Expand Up @@ -140,6 +143,8 @@ do n = 1, instances !{

pprod_gross(:,:,:) = 0.0
zprod_gross(:,:,:) = 0.0
export_prod(:,:) = 0.0
export_inorg(:,:) = 0.0
light_limit(:,:) = 0.0
radbio3d(:,:,:) = 0.0
npp3d(:,:,:) = 0.0
Expand Down Expand Up @@ -169,7 +174,8 @@ do n = 1, instances !{
ind_dic = biotic(n)%ind_bgc(id_dic)
ind_adic = biotic(n)%ind_bgc(id_adic)
ind_alk = biotic(n)%ind_bgc(id_alk)
ind_po4 = biotic(n)%ind_bgc(id_po4)
! ind_po4 = biotic(n)%ind_bgc(id_po4)
if (id_po4 .ne. 0) ind_po4 = biotic(n)%ind_bgc(id_po4)
ind_o2 = biotic(n)%ind_bgc(id_o2)

ind_no3 = biotic(n)%ind_bgc(id_no3)
Expand Down Expand Up @@ -328,7 +334,7 @@ do n = 1, instances !{
! change the ratio of remineralisation of det in upper to lower ocean from 5 to 2,
! but keep the remin rate the same at depth, so change the value in rjm_param_2010.... as well
! mac, jun10.
if (grid%zw(k) .ge. 180) f41 = f41*.5 ! reduce decay below 300m
if (grid%zw(k) .ge. 180) f41 = f41*.5 ! reduce decay at depth

if (id_caco3.ne.0) then
f51 = muecaco3(i,j)*biocaco3 *bioma(i,k,id_caco3)
Expand Down Expand Up @@ -383,6 +389,14 @@ do n = 1, instances !{
enddo !} tn


! Calculate export production, organic and inorganic. mac, sep18.
k=k100
do i = isc, iec !{
export_prod(i,j)=t_prog(ind_det)%field(i,j,k,time%taum1)*wdetbio(i,j)*grid%tmask(i,j,k)
export_inorg(i,j)=t_prog(ind_caco3)%field(i,j,k,time%taum1)*wcaco3(i,j)*grid%tmask(i,j,k)
enddo !} i


!chd add biotically induced tendency to biotracers,
do k = 1, grid%nk !{
do i = isc, iec !{
Expand Down
112 changes: 81 additions & 31 deletions src/mom5/ocean_csiro_bgc/csiro_bgc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,6 +184,7 @@ module csiro_bgc_mod !{
real :: sal_global
real, allocatable, dimension(:) :: bgc_global
real, allocatable, dimension(:,:) :: htotal
real, allocatable, dimension(:,:) :: ahtotal
real, allocatable, dimension(:,:) :: alpha
real, allocatable, dimension(:,:) :: csat
real, allocatable, dimension(:,:) :: csat_csurf
Expand Down Expand Up @@ -235,6 +236,7 @@ module csiro_bgc_mod !{
character*6 :: qbio_model
integer :: bio_version ! version of the bgc module to use
logical :: zero_floor ! apply hard floor to bgc tracers
logical :: sea_ice_bgc ! exchange algae/phytoplankton and nutrients from sea ice BGC model
logical :: sw_thru_ice ! make this true in a coupled model, so bgc knows swflx is already modified for the presense of ice.
logical :: gasx_from_file ! use gasx exchange coefficients from provided files. mac, may13.
logical :: ice_file4gasx ! make this true in a ocean-only model, option to control whether to use ice file or the model ice when determining masking effect of ice on co2 gas exchange.
Expand Down Expand Up @@ -276,10 +278,13 @@ module csiro_bgc_mod !{
integer :: id_pprod_gross = -1
integer :: id_pprod_gross_2d = -1
integer :: id_zprod_gross = -1
integer :: id_export_prod = -1
integer :: id_export_inorg = -1
integer :: id_kw_o2 = -1
integer :: id_o2_sat = -1
integer :: id_sc_o2 = -1
integer :: id_pco2 = -1, id_paco2 = -1
integer :: id_htotal = -1, id_ahtotal = -1
integer :: id_co2_sat = -1, id_aco2_sat = -1
integer :: id_caco3_sed_remin, id_det_sed_remin
integer :: id_caco3_sed_depst, id_det_sed_depst
Expand Down Expand Up @@ -347,7 +352,7 @@ module csiro_bgc_mod !{
real, allocatable, dimension(:,:) :: npp2d
real, allocatable, dimension(:,:,:) :: npp3d
real, allocatable, dimension(:,:,:) :: pprod_gross
real, allocatable, dimension(:,:) :: pprod_gross_2d
real, allocatable, dimension(:,:) :: pprod_gross_2d, export_prod, export_inorg
real, allocatable, dimension(:,:,:) :: zprod_gross
real, allocatable, dimension(:) :: ray
real, allocatable, dimension(:) :: dummy
Expand Down Expand Up @@ -521,6 +526,8 @@ subroutine allocate_arrays (isc, iec, jsc, jec, isd, ied, jsd, jed, nk) !{
allocate( pprod_gross(isc:iec,jsc:jec,nk) )
allocate( pprod_gross_2d(isc:iec,jsc:jec) )
allocate( zprod_gross(isc:iec,jsc:jec,nk) )
allocate( export_prod(isc:iec,jsc:jec) )
allocate( export_inorg(isc:iec,jsc:jec) )

allocate (tmp(isd:ied,jsd:jed) )
allocate ( tracer_sources(0:nk) )
Expand Down Expand Up @@ -558,6 +565,7 @@ subroutine allocate_arrays (isc, iec, jsc, jec, isd, ied, jsd, jed, nk) !{
! allocate( biotic(n)%ind_bgc(1:ntr_bgc) )

allocate( biotic(n)%htotal(isd:ied,jsd:jed) )
allocate( biotic(n)%ahtotal(isd:ied,jsd:jed) )
allocate( biotic(n)%alpha(isd:ied,jsd:jed) )
allocate( biotic(n)%csat(isd:ied,jsd:jed) )
allocate( biotic(n)%csat_csurf(isd:ied,jsd:jed) )
Expand Down Expand Up @@ -617,6 +625,7 @@ subroutine allocate_arrays (isc, iec, jsc, jec, isd, ied, jsd, jed, nk) !{

do n = 1, instances !{
biotic(n)%htotal(:,:) = 1.e-8
biotic(n)%ahtotal(:,:) = 1.e-8
biotic(n)%sio2(:,:) = 35. *1e-3
biotic(n)%bgc_global(:) = 0. ! this will make vstf zero
biotic(n)%sal_global = 35.
Expand Down Expand Up @@ -981,6 +990,7 @@ subroutine csiro_bgc_sbc(isc, iec, jsc, jec, isd, ied, jsd, jed, &
integer :: second
integer :: year
integer :: indtemp, indsal
real :: schmidt_temp

! =====================================================================
! begin executable code
Expand All @@ -1006,9 +1016,10 @@ subroutine csiro_bgc_sbc(isc, iec, jsc, jec, isd, ied, jsd, jed, &
else
do j=jsc, jec
do i=isc, iec
! the relations 0.31 * u^2 comes from Wanninkhof 1992, and is the coefficient for steady wind speed; the equation is 0.39 * u^2 for instaneous wind speeds. I don't know what exactly the timescale is between steady and instaneous...
! the relations 0.31 * u^2 comes from Wanninkhof 1992, and is the coefficient for steady wind speeds, suitable for shipboard spot
! measurements or derived from scatterometers; the equation is 0.39 * u^2 for long-term wind speeds (e.g. climatology).
! the 3.6e5 is a conversion of units; cm/hr to m/s.
! mac, may13.
! mac, may13 (comments updated feb22).
xkw_t(i,j)=(0.31 * wnd(i,j)**2.0 ) /3.6e5
enddo ! i
enddo ! j
Expand Down Expand Up @@ -1064,9 +1075,10 @@ subroutine csiro_bgc_sbc(isc, iec, jsc, jec, isd, ied, jsd, jed, &

do j = jsc, jec !{
do i = isc, iec !{
sc_co2(i,j) = 2073.1 + t_prog(indtemp)%field(i,j,1,time%taum1) * &
(-125.62 + t_prog(indtemp)%field(i,j,1,time%taum1) * &
(3.6276 + t_prog(indtemp)%field(i,j,1,time%taum1) * &
schmidt_temp = min(35.0, t_prog(indtemp)%field(i,j,1,time%taum1))! add max limit to temp used for Schmidt number calculation. mac, jan21.
sc_co2(i,j) = 2073.1 + schmidt_temp * &
Copy link
Collaborator

Choose a reason for hiding this comment

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

will this give numerically identical results to the old sc_co2 calculation?

Copy link
Author

Choose a reason for hiding this comment

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

Good point, the Schmidt Number and gas exchange will not be identical wherever there is a surface temperature over 35 degC.

The change was made for the ESM after some points in extended warm climate CMIP6 experiments reached about 42 degC, where the Schmidt number went negative and blew up the model, as found by Tilo Z. and Martin D. in late 2020.
Checking the reference to the Schmidt formula (Wanninkhof 1992, https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/92JC00188), the temperature range of the calculation should only be up to 30 degC, which we have probably been abusing forever.
Martin found that the limit of 35 degC was already being used in the equivalent component of NEMO, so we adopted the same in the ESM1.5.
This way the model is stable for these hot climate experiments.

(-125.62 + schmidt_temp * &
(3.6276 + schmidt_temp * &
(-0.043219))) * grid%tmask(i,j,1)
enddo !} i
enddo !} j
Expand All @@ -1079,9 +1091,10 @@ subroutine csiro_bgc_sbc(isc, iec, jsc, jec, isd, ied, jsd, jed, &

do j = jsc, jec !{
do i = isc, iec !{
sc_o2(i,j) = 1638.0 + t_prog(indtemp)%field(i,j,1,time%taum1) * &
(-81.83 + t_prog(indtemp)%field(i,j,1,time%taum1) * &
(1.483 + t_prog(indtemp)%field(i,j,1,time%taum1) * &
schmidt_temp = min(35.0, t_prog(indtemp)%field(i,j,1,time%taum1))! add max limit to temp used for Schmidt number calculation. mac, jan21.
sc_o2(i,j) = 1638.0 + schmidt_temp * &
Copy link
Collaborator

Choose a reason for hiding this comment

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

as above, will this give numerically identical results to the old sc_co2 calculation?

Copy link
Author

Choose a reason for hiding this comment

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

Good point, the Schmidt Number and gas exchange will not be identical wherever there is a surface temperature over 35 degC.

The change was made for the ESM after some points in extended warm climate CMIP6 experiments reached about 42 degC, where the Schmidt number went negative and blew up the model, as found by Tilo Z. and Martin D. in late 2020.
Checking the reference to the Schmidt formula (Wanninkhof 1992, https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/92JC00188), the temperature range of the calculation should only be up to 30 degC, which we have probably been abusing forever.
Martin found that the limit of 35 degC was already being used in the equivalent component of NEMO, so we adopted the same in the ESM1.5.
This way the model is stable for these hot climate experiments.

(-81.83 + schmidt_temp * &
(1.483 + schmidt_temp * &
(-0.008004))) * grid%tmask(i,j,1)
enddo !} i
enddo !} j
Expand Down Expand Up @@ -1157,7 +1170,7 @@ subroutine csiro_bgc_sbc(isc, iec, jsc, jec, isd, ied, jsd, jed, &
biotic(n)%po4(isd:ied,jsd:jed),&
biotic(n)%sio2(isd:ied,jsd:jed), &
htotallo(isc:iec,jsc:jec), htotalhi(isc:iec,jsc:jec), &
biotic(n)%htotal(isc:iec,jsc:jec), &
biotic(n)%ahtotal(isc:iec,jsc:jec), &
biotic(n)%acsurf(isc:iec,jsc:jec), &
alpha=biotic(n)%alpha(isc:iec,jsc:jec) , &
pco2surf = biotic(n)%paco2surf(isc:iec,jsc:jec), &
Expand Down Expand Up @@ -1344,27 +1357,28 @@ subroutine csiro_bgc_sbc(isc, iec, jsc, jec, isd, ied, jsd, jed, &
enddo !} n
endif

!ice-to-ocean flux of algae
if (id_phy.ne.0) then
do n = 1, instances !{
do j = jsc, jec !{
do i = isc, iec !{
t_prog(ind_phy)%stf(i,j) = iof_alg(i,j)
enddo !} i
enddo !} j
enddo !} n
endif
!ice-to-ocean flux of nitrate
if (id_no3.ne.0) then
do n = 1, instances !{
do j = jsc, jec !{
do i = isc, iec !{
t_prog(ind_no3)%stf(i,j) = iof_nit(i,j)
enddo !} i
enddo !} j
enddo !} n
endif

if (sea_ice_bgc) then
!ice-to-ocean flux of algae
if (id_phy.ne.0) then
do n = 1, instances !{
do j = jsc, jec !{
do i = isc, iec !{
t_prog(ind_phy)%stf(i,j) = iof_alg(i,j)
enddo !} i
enddo !} j
enddo !} n
endif
!ice-to-ocean flux of nitrate
if (id_no3.ne.0) then
do n = 1, instances !{
do j = jsc, jec !{
do i = isc, iec !{
t_prog(ind_no3)%stf(i,j) = iof_nit(i,j)
enddo !} i
enddo !} j
enddo !} n
endif
endif ! sea_ice_bgc

if (.not. use_waterflux) then !{
! rjm - One only needs to compute virtual fluxes if waterflux is not used
Expand Down Expand Up @@ -1587,6 +1601,7 @@ subroutine csiro_bgc_init !{
call fm_util_set_value('qbio_model', 'bio_vx') ! version name
call fm_util_set_value('bio_version',1) ! version of the bgc module
call fm_util_set_value('zero_floor', .false.) ! apply hard floor to bgc tracers
call fm_util_set_value('sea_ice_bgc', .true.) ! exchange algae/phytoplankton and nutrients from sea ice BGC model
call fm_util_set_value('sw_thru_ice', .true.) ! is shortwave flux modified by an ice model?
call fm_util_set_value('gasx_from_file', .true.)! use file with gas exchange coefficients?
call fm_util_set_value('ice_file4gasx', .true.)! use file with ice cover for gas exchange?
Expand Down Expand Up @@ -1624,6 +1639,7 @@ subroutine csiro_bgc_init !{
bio_version = fm_util_get_integer ('bio_version', scalar = .true.)
! rjm: Tracer to use
zero_floor = fm_util_get_logical ('zero_floor', scalar = .true.)
sea_ice_bgc = fm_util_get_logical ('sea_ice_bgc', scalar = .true.)
sw_thru_ice = fm_util_get_logical ('sw_thru_ice', scalar = .true.)
gasx_from_file = fm_util_get_logical ('gasx_from_file', scalar = .true.)
ice_file4gasx = fm_util_get_logical ('ice_file4gasx', scalar = .true.)
Expand Down Expand Up @@ -1934,6 +1950,14 @@ subroutine csiro_bgc_source(isc, iec, jsc, jec, isd, ied, jsd, jed, T_prog, grid
used = send_data(id_paco2, biotic(1)%paco2surf(isc:iec,jsc:jec), &
time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
endif
if (id_htotal .gt. 0) then
used = send_data(id_htotal, biotic(1)%htotal(isc:iec,jsc:jec), &
time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
endif
if (id_ahtotal .gt. 0) then
used = send_data(id_ahtotal, biotic(1)%ahtotal(isc:iec,jsc:jec), &
time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
endif
if (id_co2_sat .gt. 0) then
used = send_data(id_co2_sat, biotic(1)%csat_csurf(isc:iec,jsc:jec), &
time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
Expand Down Expand Up @@ -2005,6 +2029,16 @@ subroutine csiro_bgc_source(isc, iec, jsc, jec, isd, ied, jsd, jed, T_prog, grid
time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
endif

if (id_export_prod .gt. 0) then
used = send_data(id_export_prod, export_prod(isc:iec,jsc:jec), &
time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
endif

if (id_export_inorg .gt. 0) then
used = send_data(id_export_inorg, export_inorg(isc:iec,jsc:jec), &
time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1))
endif

! Gross production of zooplankton

if (id_zprod_gross .gt. 0) then
Expand Down Expand Up @@ -2412,6 +2446,14 @@ subroutine csiro_bgc_start (time, domain, grid) !{
'paco2', grid%tracer_axes(1:2), &
Time%model_time, 'pCO2 inc. anthropogenic', ' ', &
missing_value = -1.0e+10)
id_htotal = register_diag_field('ocean_model', &
'htotal', grid%tracer_axes(1:2), &
Time%model_time, 'Htotal', ' ', &
missing_value = -1.0e+10)
id_ahtotal = register_diag_field('ocean_model', &
'ahtotal', grid%tracer_axes(1:2), &
Time%model_time, 'Htotal inc. anthropogenic', ' ', &
missing_value = -1.0e+10)
id_co2_sat = register_diag_field('ocean_model', &
'co2_saturation', grid%tracer_axes(1:2), &
Time%model_time, 'CO2 saturation', 'mmol/m^3', &
Expand Down Expand Up @@ -2558,6 +2600,14 @@ subroutine csiro_bgc_start (time, domain, grid) !{
grid%tracer_axes(1:3),Time%model_time, 'Gross ZOO production', &
'mmolN/m^3/s',missing_value = -1.0e+10)

id_export_prod = register_diag_field('ocean_model','export_prod', &
grid%tracer_axes(1:2),Time%model_time, 'Organic export through 100m', &
'mmolN/m^2/s',missing_value = -1.0e+10)

id_export_inorg = register_diag_field('ocean_model','export_inorg', &
grid%tracer_axes(1:2),Time%model_time, 'Inorganic export through 100m', &
'mmolN/m^2/s',missing_value = -1.0e+10)

id_caco3_sediment = register_diag_field('ocean_model','caco3_sediment', &
grid%tracer_axes(1:2),Time%model_time, 'Accumulated CaCO3 in sediment at base of water column', &
'mmolN/m^2',missing_value = -1.0e+10)
Expand Down