From 799031d44d524ec72b1f73a5e3da69c7b611e593 Mon Sep 17 00:00:00 2001 From: Russell Fiedler Date: Fri, 17 Jan 2020 14:32:54 +1100 Subject: [PATCH 1/7] Add CSIRO BGC and allow 10m winds unconditionally in accessom --- exp/MOM_compile.csh | 14 +- exp/ocean_compile.csh | 5 +- src/accessom_coupler/mom_oasis3_interface.F90 | 8 - src/mom5/ocean_core/ocean_model.F90 | 9 + src/mom5/ocean_core/ocean_sbc.F90 | 114 +- src/mom5/ocean_csiro_bgc/bio_v3.inc | 545 ++++ src/mom5/ocean_csiro_bgc/csiro_bgc.F90 | 2523 +++++++++++++++++ src/mom5/ocean_csiro_bgc/ocmip2_co2calc.F90 | 729 +++++ src/mom5/ocean_tracers/ocean_tpm.F90 | 77 +- 9 files changed, 3992 insertions(+), 32 deletions(-) create mode 100755 src/mom5/ocean_csiro_bgc/bio_v3.inc create mode 100755 src/mom5/ocean_csiro_bgc/csiro_bgc.F90 create mode 100755 src/mom5/ocean_csiro_bgc/ocmip2_co2calc.F90 diff --git a/exp/MOM_compile.csh b/exp/MOM_compile.csh index 0f98d6a30e..48aa291178 100755 --- a/exp/MOM_compile.csh +++ b/exp/MOM_compile.csh @@ -54,8 +54,10 @@ if ( $help ) then echo " EBM : ocean-seaice-land-atmosphere coupled model with energy balance atmosphere" echo " ACCESS-CM : ocean component of ACCESS-CM model." echo " ACCESS-OM : ocean component of ACCESS-OM model." + echo " ACCESS-ESM : ocean component of ACCESS-ESM model with CSIRO BGC (Wombat)." + echo " ACCESS-OM-BGC: ocean component of ACCESS-OM model with CSIRO BGC (Wombat)." echo - echo "--platform followed by the platform name that has a corresponding environs file in the ../bin dir, default is gfortran" + echo "--platform followed by the platform name that has a corresponfing environ file in the ../bin dir, default is gfortran" echo echo "--use_netcdf4 use NetCDF4, the default is NetCDF4. Warning: many of the standard experiments don't work with NetCDF4." echo @@ -87,8 +89,12 @@ if ( $type == EBM ) then set cppDefs = ( "-Duse_netCDF -Duse_netCDF3 -Duse_libMPI -DLAND_BND_TRACERS -DOVERLOAD_C8 -DOVERLOAD_C4 -DOVERLOAD_R4" ) else if( $type == ACCESS-OM ) then set cppDefs = ( "-Duse_netCDF -Duse_libMPI -DACCESS -DACCESS_OM" ) +else if( $type == ACCESS-OM-BGC ) then + set cppDefs = ( "-Duse_netCDF -Duse_libMPI -DACCESS -DACCESS_OM -DCSIRO_BGC" ) else if( $type == ACCESS-CM ) then set cppDefs = ( "-Duse_netCDF -Duse_libMPI -DACCESS -DACCESS_CM" ) +else if( $type == ACCESS-ESM ) then + set cppDefs = ( "-Duse_netCDF -Duse_libMPI -DACCESS -DACCESS_CM -DCSIRO_BGC" ) endif if ( $unit_testing ) then @@ -134,7 +140,7 @@ cd $root/exp source ./ocean_compile.csh if ( $status ) exit $status -if( $type != MOM_solo && $type != ACCESS-OM && $type != ACCESS-CM) then +if( $type != MOM_solo && $type != ACCESS-OM && $type != ACCESS-CM && $type != ACCESS-OM-BGC && $type != ACCESS-ESM) then cd $root/exp source ./ice_compile.csh if ( $status ) exit $status @@ -186,12 +192,12 @@ cd $executable:h if( $type == MOM_solo ) then set srcList = ( mom5/drivers ) set libs = "$executable:h:h/lib_ocean/lib_ocean.a $executable:h:h/lib_FMS/lib_FMS.a" -else if( $type == ACCESS-CM ) then +else if( $type == ACCESS-CM || $type == ACCESS-ESM) then set srcList = ( accesscm_coupler ) set includes = "-I$executable:h:h/lib_FMS -I$executable:h:h/$type/lib_ocean" set libs = "$executable:h:h/$type/lib_ocean/lib_ocean.a $executable:h:h/lib_FMS/lib_FMS.a" setenv OASIS true -else if( $type == ACCESS-OM ) then +else if( $type == ACCESS-OM || $type == ACCESS-OM-BGC) then set srcList = ( accessom_coupler ) set includes = "-I$executable:h:h/lib_FMS -I$executable:h:h/$type/lib_ocean" set libs = "$executable:h:h/$type/lib_ocean/lib_ocean.a $executable:h:h/lib_FMS/lib_FMS.a" diff --git a/exp/ocean_compile.csh b/exp/ocean_compile.csh index cd887ea81f..c547360ce8 100644 --- a/exp/ocean_compile.csh +++ b/exp/ocean_compile.csh @@ -4,8 +4,11 @@ set srcList = ( mom5/ocean_core mom5/ocean_diag mom5/ocean_wave mom5/ocean_blobs set lib_name = "lib_ocean" -if( $type == ACCESS-OM || $type == ACCESS-CM ) then +if( $type == ACCESS-OM || $type == ACCESS-CM || $type == ACCESS-OM-BGC || $type == ACCESS-ESM) then set srcList = ( $srcList mom5/ocean_access ) + if( $type == ACCESS-OM-BGC || $type == ACCESS-ESM) then + set srcList = ( $srcList mom5/ocean_csiro_bgc ) + endif mkdir -p $executable:h:h/$type/$lib_name cd $executable:h:h/$type/$lib_name else diff --git a/src/accessom_coupler/mom_oasis3_interface.F90 b/src/accessom_coupler/mom_oasis3_interface.F90 index 46d0a0723f..716284f321 100644 --- a/src/accessom_coupler/mom_oasis3_interface.F90 +++ b/src/accessom_coupler/mom_oasis3_interface.F90 @@ -143,11 +143,7 @@ module mom_oasis3_interface_mod integer iisc,iiec,jjsc,jjec integer iisd,iied,jjsd,jjed -#if defined(ACCESS_WND) integer, parameter :: max_fields_in=21 -#else -integer, parameter :: max_fields_in=20 -#endif integer, parameter :: max_fields_out=8 @@ -358,9 +354,7 @@ subroutine coupler_init(Dom, Time, Time_step_coupled, dt_cpld, Run_len, & mom_name_read(18)='mh_flux' ! Heat flux due to melting mom_name_read(19)='wfimelt' !Water flux due to ice melting mom_name_read(20)='wfiform' !Water flux due to ice forming -#if defined(ACCESS_WND) mom_name_read(21)='wnd_io' ! -#endif !ocn ==> ice mom_name_write(:)='' @@ -741,10 +735,8 @@ subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Time) Ice_ocean_boundary%wfimelt(iisc:iiec,jjsc:jjec) = vwork(iisc:iiec,jjsc:jjec) case('wfiform') Ice_ocean_boundary%wfiform(iisc:iiec,jjsc:jjec) = vwork(iisc:iiec,jjsc:jjec) -#if defined(ACCESS_WND) case('wnd_io') Ice_ocean_boundary%wnd(iisc:iiec,jjsc:jjec) = vwork(iisc:iiec,jjsc:jjec) -#endif case DEFAULT ! Probable error. Leave as warning for the moment. RASF call mpp_error(WARNING,& diff --git a/src/mom5/ocean_core/ocean_model.F90 b/src/mom5/ocean_core/ocean_model.F90 index 04686c07b2..83f8893e3a 100644 --- a/src/mom5/ocean_core/ocean_model.F90 +++ b/src/mom5/ocean_core/ocean_model.F90 @@ -1638,13 +1638,22 @@ subroutine update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, & ! compute ocean tendencies from tracer packages call mpp_clock_begin(id_otpm_source) +#if ! defined (CSIRO_BGC) call ocean_tpm_source(isd, ied, jsd, jed, Domain, Grid, T_prog(:), T_diag(:), & Time, Thickness, Dens, surf_blthick, dtts) +#else + call ocean_tpm_source(isd, ied, jsd, jed, Domain, Grid, T_prog(:), T_diag(:), & + Time, Thickness, Dens, surf_blthick, dtts, swflx, sw_frac_zt) +#endif call mpp_clock_end(id_otpm_source) ! set ocean surface boundary conditions for the tracer packages call mpp_clock_begin(id_otpm_bbc) +#if ! defined (CSIRO_BGC) call ocean_tpm_bbc(Domain, Grid, T_prog(:)) +#else + call ocean_tpm_bbc(Domain, Grid, T_prog(:), Time) +#endif call mpp_clock_end(id_otpm_bbc) ! add sponges to T_prog%th_tendency diff --git a/src/mom5/ocean_core/ocean_sbc.F90 b/src/mom5/ocean_core/ocean_sbc.F90 index c85409b006..dd93e0e808 100644 --- a/src/mom5/ocean_core/ocean_sbc.F90 +++ b/src/mom5/ocean_core/ocean_sbc.F90 @@ -534,6 +534,9 @@ module ocean_sbc_mod use ocean_util_mod, only: diagnose_2d, diagnose_2d_u, diagnose_3d_u, diagnose_sum use ocean_tracer_util_mod, only: diagnose_3d_rho +#if defined(CSIRO_BGC) +use csiro_bgc_mod, only: csiro_bgc_virtual_fluxes, do_csiro_bgc +#endif implicit none private @@ -671,10 +674,13 @@ module ocean_sbc_mod #if defined(ACCESS_CM) || defined(ACCESS_OM) integer :: id_wfimelt =-1 integer :: id_wfiform =-1 +integer :: id_aice =-1 +integer :: id_wnd =-1 #endif #if defined(ACCESS_CM) integer :: id_licefw =-1 integer :: id_liceht =-1 +integer :: id_atm_co2 =-1 #endif @@ -817,6 +823,11 @@ module ocean_sbc_mod real, allocatable, dimension(:,:,:) :: sslope real, allocatable, dimension(:,:) :: aice #endif +#if defined(ACCESS_CM) +real, allocatable, dimension(:,:) :: co2flux +real, allocatable, dimension(:,:) :: ocn_co2 +real, allocatable, dimension(:,:) :: atm_co2 +#endif ! ice-ocean-boundary fields are allocated using absolute @@ -1104,10 +1115,12 @@ subroutine ocean_sbc_init(Grid, Domain, Time, T_prog, T_diag, & #if defined(ACCESS_CM) || defined(ACCESS_OM) allocate ( Ocean_sfc%gradient (isc_bnd:iec_bnd,jsc_bnd:jec_bnd,2)) allocate ( sslope(isc:iec, jsc:jec, 2) ) - allocate ( aice(isc:iec, jsc:jec) ) + allocate ( aice(isd:ied, jsd:jed) ) #if defined(ACCESS_CM) allocate ( Ocean_sfc%co2 (isc_bnd:iec_bnd,jsc_bnd:jec_bnd), & Ocean_sfc%co2flux (isc_bnd:iec_bnd,jsc_bnd:jec_bnd)) + allocate ( co2flux(isd:ied,jsd:jed),ocn_co2(isd:ied,jsd:jed)) + allocate ( atm_co2(isd:ied,jsd:jed)) #endif #endif @@ -1126,6 +1139,9 @@ subroutine ocean_sbc_init(Grid, Domain, Time, T_prog, T_diag, & #if defined(ACCESS_CM) Ocean_sfc%co2 = 0.0 Ocean_sfc%co2flux = 0.0 + co2flux = 0.0 + ocn_co2 = 0.0 + atm_co2 = 0.0 #endif Ocean_sfc%area = Grid%dat(isc:iec, jsc:jec) * Grid%tmask(isc:iec, jsc:jec, 1) !grid cell area @@ -1969,6 +1985,14 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog) '(kg/sec)/1e15', missing_value=missing_value,range=(/-1.e10,1.e10/)) #if defined(ACCESS_CM) || defined(ACCESS_OM) + id_aice = register_diag_field('ocean_model','aice', Grd%tracer_axes(1:2),& + Time%model_time, 'fraction of surface area covered with ice', 'm^2/m^2' , & + missing_value=missing_value,range=(/-1.e1,1.e1/), & + standard_name='areal_ice_concentration' ) + id_wnd = register_diag_field('ocean_model','wnd', Grd%tracer_axes(1:2),& + Time%model_time, 'Wind speed', 'm/s' , & + missing_value=missing_value,range=(/-1.e3,1.e3/), & + standard_name='wind_speed' ) id_total_ocean_wfimelt = register_diag_field('ocean_model','total_ocean_wfimelt', & Time%model_time, 'total icemelt into ocean (>0 enters ocean)', & 'kg/sec/1e15', missing_value=missing_value,range=(/-1.e10,1.e10/)) @@ -1978,6 +2002,10 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog) #endif #if defined(ACCESS_CM) + id_atm_co2 = register_diag_field('ocean_model','atm_co2', Grd%tracer_axes(1:2),& + Time%model_time, 'Atmospheric CO2 content', 'ppm' , & + missing_value=missing_value,range=(/-1.e1,1.e4/), & + standard_name='atmospheric_co2' ) id_total_ocean_licefw = register_diag_field('ocean_model','total_ocean_licefw', & Time%model_time, 'total land icemelt into ocean (>0 enters ocean)', & 'kg/sec/1e15', missing_value=missing_value,range=(/-1.e10,1.e10/)) @@ -2780,6 +2808,10 @@ subroutine initialize_ocean_sfc(Time, Thickness, T_prog, T_diag, Velocity, Ocean Ocean_sfc%v_surf(isc_bnd:iec_bnd,jsc_bnd:jec_bnd) = Velocity%u(isc:iec,jsc:jec,1,2,taup1) Ocean_sfc%sea_lev(isc_bnd:iec_bnd,jsc_bnd:jec_bnd) = Thickness%sea_lev(isc:iec,jsc:jec) Ocean_sfc%frazil(isc_bnd:iec_bnd,jsc_bnd:jec_bnd) = 0.0 +#if defined(ACCESS_CM) + Ocean_sfc%co2flux(isc_bnd:iec_bnd,jsc_bnd:jec_bnd) = co2flux(isc:iec,jsc:jec) !These should really come from a restart RASF + Ocean_sfc%co2(isc_bnd:iec_bnd,jsc_bnd:jec_bnd) = ocn_co2(isc:iec,jsc:jec) +#endif end where ! when enabled, use FAFMIP redistributed heat tracer for sst @@ -2814,6 +2846,11 @@ subroutine initialize_ocean_sfc(Time, Thickness, T_prog, T_diag, Velocity, Ocean id_field = register_restart_field(Sfc_restart, filename, 'v_surf', Ocean_sfc%v_surf,Ocean_sfc%Domain) id_field = register_restart_field(Sfc_restart, filename, 'sea_lev',Ocean_sfc%sea_lev,Ocean_sfc%Domain) id_field = register_restart_field(Sfc_restart, filename, 'frazil', Ocean_sfc%frazil,Ocean_sfc%Domain) +#if defined(ACCESS_CM) +!RASF Make these optional so we don't break existing runs. + id_field = register_restart_field(Sfc_restart, filename, 'co2flux',Ocean_sfc%co2flux,Ocean_sfc%Domain, mandatory=.false.) + id_field = register_restart_field(Sfc_restart, filename, 'ocn_co2', Ocean_sfc%co2,Ocean_sfc%Domain, mandatory=.false.) +#endif if (file_exist('INPUT/ocean_sbc.res.nc')) then call restore_state(Sfc_restart) endif @@ -2898,6 +2935,10 @@ subroutine sum_ocean_sfc(Time, Thickness, T_prog, T_diag, Dens, Velocity, Ocean_ Ocean_sfc%sea_lev(i,j) = Ocean_sfc%sea_lev(i,j) + Thickness%sea_lev(ii,jj) #if defined(ACCESS_CM) || defined(ACCESS_OM) Ocean_sfc%gradient(i,j,:) = Ocean_sfc%gradient(i,j,:) + sslope(ii,jj,:) +#endif +#if defined(ACCESS_CM) + Ocean_sfc%co2flux(i,j) = Ocean_sfc%co2flux(i,j) + co2flux(ii,jj) + Ocean_sfc%co2(i,j) = Ocean_sfc%co2(i,j) + ocn_co2(ii,jj) #endif enddo enddo @@ -2997,6 +3038,10 @@ subroutine zero_ocean_sfc(Ocean_sfc) Ocean_sfc%frazil(i,j) = 0.0 #if defined(ACCESS_CM) || defined(ACCESS_OM) Ocean_sfc%gradient(i,j,:)= 0.0 +#endif +#if defined(ACCESS_CM) + Ocean_sfc%co2flux(i,j) = 0.0 + Ocean_sfc%co2(i,j) = 0.0 #endif enddo enddo @@ -3074,6 +3119,10 @@ subroutine avg_ocean_sfc(Time, Thickness, T_prog, T_diag, Velocity, Ocean_sfc) Ocean_sfc%sea_lev(i,j) = Ocean_sfc%sea_lev(i,j)*divid #if defined(ACCESS_CM) || defined(ACCESS_OM) Ocean_sfc%gradient(i,j,:) = Ocean_sfc%gradient(i,j,:)*divid +#endif +#if defined(ACCESS_CM) + Ocean_sfc%co2flux(i,j) = Ocean_sfc%co2flux(i,j)*divid + Ocean_sfc%co2(i,j) = Ocean_sfc%co2(i,j)*divid #endif endif enddo @@ -3099,6 +3148,10 @@ subroutine avg_ocean_sfc(Time, Thickness, T_prog, T_diag, Velocity, Ocean_sfc) if(Grd%tmask(ii,jj,1) == 1.0) then Ocean_sfc%s_surf(i,j) = T_prog(index_salt)%field(ii,jj,1,taup1) Ocean_sfc%sea_lev(i,j)= Thickness%sea_lev(ii,jj) +#if defined(ACCESS_CM) + Ocean_sfc%co2flux(i,j) = co2flux(ii,jj) + Ocean_sfc%co2(i,j) = ocn_co2(ii,jj) +#endif endif enddo enddo @@ -4282,13 +4335,43 @@ subroutine get_ocean_sbc(Time, Ice_ocean_boundary, Thickness, Dens, Ext_mode, T_ call mpp_update_domains(pme(:,:), Dom%domain2d) endif +#if defined(ACCESS_CM) + do j = jsc_bnd,jec_bnd + do i = isc_bnd,iec_bnd + ii = i + i_shift + jj = j + j_shift + atm_co2(ii,jj) = Ice_ocean_boundary%co2(i,j)*Grd%tmask(ii,jj,1) + enddo + enddo +#endif + +#if defined(ACCESS_CM) || defined(ACCESS_OM) + do j = jsc_bnd,jec_bnd + do i = isc_bnd,iec_bnd + ii = i + i_shift + jj = j + j_shift + aice(ii,jj) = Ice_ocean_boundary%aice(i,j)*Grd%tmask(ii,jj,1) + enddo + enddo +#endif !--------compute surface tracer fluxes from tracer packages------------------- ! +#if defined(ACCESS_CM) && defined(CSIRO_BGC) + call ocean_tpm_sbc(Dom, Grd, T_prog(:), Time, Ice_ocean_boundary%fluxes, runoff, & + isc_bnd, iec_bnd, jsc_bnd, jec_bnd,aice, Velocity%u10, & + use_waterflux, salt_restore_as_salt_flux, atm_co2, co2flux, ocn_co2) +#elif defined(ACCESS_OM) && defined(CSIRO_BGC) +! Do not pass co2flux, ocn_co2 or atm_co2 + call ocean_tpm_sbc(Dom, Grd, T_prog(:), Time, Ice_ocean_boundary%fluxes, runoff, & + isc_bnd, iec_bnd, jsc_bnd, jec_bnd,aice=aice, wnd=Velocity%u10, & + use_waterflux=use_waterflux, salt_restore_as_salt_flux=salt_restore_as_salt_flux) +#else call ocean_tpm_sbc(Dom, Grd, T_prog(:), Time, Ice_ocean_boundary%fluxes, runoff, & isc_bnd, iec_bnd, jsc_bnd, jec_bnd) - +! Leave this case blank at the moment. We want to use the bgc with SIS etc. +#endif !--------send diagnostics------------------- ! @@ -4297,16 +4380,6 @@ subroutine get_ocean_sbc(Time, Ice_ocean_boundary, Thickness, Dens, Ext_mode, T_ melt, liquid_precip, frozen_precip, evaporation, sensible, longwave, & latent, swflx, swflx_vis) -#if defined(ACCESS_CM) || defined(ACCESS_OM) - do j = jsc_bnd,jec_bnd - do i = isc_bnd,iec_bnd - ii = i + i_shift - jj = j + j_shift - aice(ii,jj) = Ice_ocean_boundary%aice(i,j)*Grd%tmask(ii,jj,1) - enddo - enddo -#endif - end subroutine get_ocean_sbc ! NAME="get_ocean_sbc" @@ -4548,6 +4621,14 @@ subroutine flux_adjust(Time, T_diag, Dens, Ext_mode, T_prog, Velocity, river, me enddo enddo +#if defined(CSIRO_BGC) + ! if salt fluxes are used to restore salinity, then virtual fluxes are + ! needed for csiro BGC tracers. mac, dec12. + if (do_csiro_bgc) then + call csiro_bgc_virtual_fluxes(isc, iec, jsc, jec, isd, ied, jsd, jed, flx_restore, T_prog) + endif +#endif + endif ! endif for if (use_waterflux .and. .not. salt_restore_as_salt_flux) then endif ! endif for if(id_restore(index_salt) > 0 .and. salt_damp_factor > 0.0) @@ -5887,6 +5968,12 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou total_stuff = mpp_global_sum(Dom%domain2d,wrk1_2d(:,:), NON_BITWISE_EXACT_SUM) used = send_data (id_total_ocean_wfiform, total_stuff*1e-15, Time%model_time) endif + if (id_aice > 0) used = send_data(id_aice, aice(:,:), & + Time%model_time, rmask=Grd%tmask(:,:,1), & + is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) + if (id_wnd > 0) used = send_data(id_wnd, Velocity%u10(:,:), & + Time%model_time, rmask=Grd%tmask(:,:,1), & + is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) #endif #if defined(ACCESS_CM) ! waterflux associated with land ice melt into ocean (kg/(m2*sec)) @@ -5902,6 +5989,9 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou Time%model_time, rmask=Grd%tmask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif + if (id_atm_co2 > 0) used = send_data(id_atm_co2, atm_co2(:,:), & + Time%model_time, rmask=Grd%tmask(:,:,1), & + is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) #endif diff --git a/src/mom5/ocean_csiro_bgc/bio_v3.inc b/src/mom5/ocean_csiro_bgc/bio_v3.inc new file mode 100755 index 0000000000..59558c5a3c --- /dev/null +++ b/src/mom5/ocean_csiro_bgc/bio_v3.inc @@ -0,0 +1,545 @@ +! + +subroutine bio_v3(isc, iec, jsc, jec, isd, ied, jsd, jed, T_prog, Grid, Time, dtts, Thickness, Dens, swflx, sw_frac_zt) + + +! Based on the pt_npzd.aos version - which is the optimized version + +!----------------------------------------------------------------------- +! arguments +!----------------------------------------------------------------------- +! + +integer, intent(in) :: isc, iec +integer, intent(in) :: jsc, jec +integer, intent(in) :: isd, ied +integer, intent(in) :: jsd, jed +type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog +type(ocean_grid_type), intent(in) :: Grid +type(ocean_time_type), intent(in) :: Time +real, intent(in) :: dtts +type(ocean_thickness_type), intent(in) :: Thickness +type(ocean_density_type), intent(in) :: Dens +real, intent(in), dimension(isd:ied,jsd:jed) :: swflx ! short wave radiation flux (W/m^2) +real, intent(in), dimension(isd:,jsd:,:) :: sw_frac_zt ! short wave fraction on T grid (none) + +!----------------------------------------------------------------------- +! local variables +!----------------------------------------------------------------------- + +integer :: i +integer :: j +integer :: k +integer :: n + +logical :: used +integer :: index_temp, index_salt + +! BGC parameters now read from bgc_param.nc at run time. mac, aug11. +! include "rjm_param_201005.inc" + + integer :: ts_npzd ! number of time steps within NPZD model + integer :: kmeuph = 20 ! deepest level of euphotic zone + integer :: tn, trn + + real :: pi = 3.14159265358979 !yes, this is pi + real :: biotr(isc:iec,grid%nk,ntr_bmax), bioma(isc:iec,grid%nk,ntr_bmax), pprod(isc:iec,jsc:jec,grid%nk) + real :: bion,biop,bioz,biod,bioo,biocaco3,bioi + real :: u_npz,g_npz + real :: fx1,fx2,fx3,fx4,fu1,fu2,radbio, sw_zt, sw_zt1, sw_zt2 + real :: vpbio(isc:iec,grid%nk) + real :: avej(isc:iec,grid%nk) +!chd auxiliary variables to prevent unnecessary computation + real :: fbc + real :: f11,f21,f22,f23,f31,f32,f41,f51 + real :: epsi = 1e-15 + real :: rdtts !1/dtts + real :: dtsb + real :: adv_fb(isc:iec,1:grid%nk+1) + real, dimension(isd:ied,jsd:jed) :: mld + real :: caco3_bgc_change, no3_bgc_change + +! +! ===================================================================== +! begin executable code +! ===================================================================== + +! read biotic parameters + + call time_interp_external(alphabio_id, time%model_time, alphabio) + call time_interp_external(parbio_id, time%model_time, parbio) + call time_interp_external(kwbio_id, time%model_time, kwbio) + call time_interp_external(kcbio_id, time%model_time, kcbio) + call time_interp_external(abio_id, time%model_time, abio) + call time_interp_external(bbio_id, time%model_time, bbio) + call time_interp_external(cbio_id, time%model_time, cbio) + call time_interp_external(k1bio_id, time%model_time, k1bio) + call time_interp_external(muepbio_id, time%model_time, muepbio) + call time_interp_external(muepsbio_id, time%model_time, muepsbio) + call time_interp_external(gam1bio_id, time%model_time, gam1bio) + call time_interp_external(gbio_id, time%model_time, gbio) + call time_interp_external(epsbio_id, time%model_time, epsbio) + call time_interp_external(muezbio_id, time%model_time, muezbio) + call time_interp_external(gam2bio_id, time%model_time, gam2bio) + call time_interp_external(muedbio_id, time%model_time, muedbio) + call time_interp_external(wdetbio_id, time%model_time, wdetbio) + call time_interp_external(muecaco3_id, time%model_time, muecaco3) + call time_interp_external(wcaco3_id, time%model_time, wcaco3) + call time_interp_external(tscav_fe_id, time%model_time, tscav_fe) + call time_interp_external(fe_bkgnd_id, time%model_time, fe_bkgnd) + call time_interp_external(f_inorg_id, time%model_time, f_inorg) + + + + + +! set the maximum index for euphotic depth + do k=1,grid%nk + if (grid%zw(k).le. 400) kmeuph=k + enddo +! print*,'rjm euphotic ', kmeuph + +! +! +!----------------------------------------------------------------------- +! calculate the source terms for BIOTICs +!----------------------------------------------------------------------- +! + + + index_temp = fm_get_index('/ocean_mod/prog_tracers/temp') + index_salt = fm_get_index('/ocean_mod/prog_tracers/salt') + + + ts_npzd= max(1,nint(dtts / 900.)) + rdtts = 1/dtts + +! write (stdout(),*) ' AO-NPZD model will do ',ts_npzd,' time steps' +! write (stdout(),*) ' time step in NPZD model will be ', & +! dtts/ts_npzd,'sec.' + +!chd time step within NPZD model +!chd + dtsb=dtts/float(ts_npzd) + +! Calculate the mixed layer depth. mac, aug11. + call calc_mixed_layer_depth(Thickness, & + T_prog(index_salt)%field(isd:ied,jsd:jed,:,Time%tau), & + T_prog(index_temp)%field(isd:ied,jsd:jed,:,Time%tau), & + Dens%rho(isd:ied,jsd:jed,:,Time%tau), & + Dens%pressure_at_depth(isd:ied,jsd:jed,:), & + mld) + +! +! Loop over multiple instances +! + + + +do n = 1, instances !{ + + pprod_gross(:,:,:) = 0.0 + zprod_gross(:,:,:) = 0.0 + light_limit(:,:) = 0.0 + +! possible tracers in model + 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_o2 = biotic(n)%ind_bgc(id_o2) + + ind_no3 = biotic(n)%ind_bgc(id_no3) + + ind_zoo = biotic(n)%ind_bgc(id_zoo) + ind_phy = biotic(n)%ind_bgc(id_phy) + ind_det = biotic(n)%ind_bgc(id_det) + + ind_caco3= biotic(n)%ind_bgc(id_caco3) + ind_fe = biotic(n)%ind_bgc(id_fe) + +!chd calculate biotic source terms using Euler +!chd forward timestepping + + do j = jsc, jec !{ + + +!chd transfer tracers at tau-1 to temporary arrays + do k=1,grid%nk + do i=isc,iec + biotr(i,k,id_no3) = max(0.0,t_prog(ind_no3)%field(i,j,k,Time%taum1)) + biotr(i,k,id_phy) = max(0.0,t_prog(ind_phy)%field(i,j,k,Time%taum1)) + biotr(i,k,id_zoo) = max(0.0,t_prog(ind_zoo)%field(i,j,k,Time%taum1)) + biotr(i,k,id_det) = max(0.0,t_prog(ind_det)%field(i,j,k,Time%taum1)) + biotr(i,k,id_o2) = max(0.0,t_prog(ind_o2)%field(i,j,k,Time%taum1)) + if (id_caco3.ne.0) & + biotr(i,k,id_caco3) = max(0.0,t_prog(ind_caco3)%field(i,j,k,Time%taum1)) + if (id_fe.ne.0) & + biotr(i,k,id_fe) = max(0.0,t_prog(ind_fe)%field(i,j,k,Time%taum1)) + enddo + enddo + +!chd create mask bioma, which is 0 for biotr < epsi and 1 +!chd otherwise + do trn=1,ntr_bmax + do k=1,grid%nk + do i=isc,iec + bioma(i,k,trn) = 0. + if (biotr(i,k,trn) .gt. epsi) then + bioma(i,k,trn) = 1. + endif + enddo + enddo + enddo + + + avej(:,:)= 0.0 + vpbio(:,:)= 0.0 + + do k=1,kmeuph + do i=isc,iec + +! Shortwave intensity averaged over layer + sw_zt2 = (sw_frac_zt(i,j,k)-sw_frac_zt(i,j,k+1))/(grid%zt(k)-grid%zt(k+1) ) * (grid%zw(k)-grid%zt(k)) + sw_frac_zt(i,j,k) + if(k.eq.1) then + sw_zt1=1. + else + sw_zt1 = (sw_frac_zt(i,j,k-1)-sw_frac_zt(i,j,k))/(grid%zt(k-1)-grid%zt(k) ) * (grid%zw(k-1)-grid%zt(k-1)) + sw_frac_zt(i,j,k-1) + endif + sw_zt = (sw_zt1 + sw_zt2)*.5 * swflx(i,j) +! + +! add ice feedback on growth ...... +! Use sw_thru_ice = .true. to indicate ice term already included in swflx when running with ice model +! otherwise shortwave flux from the ocean model is corrected here. mac, sep11. + if (sw_thru_ice) then + radbio = max(0.0,parbio(i,j)*sw_zt) + else + radbio = max(0.0,parbio(i,j)*(1.0-fice_t(i,j))*sw_zt) + endif + + +! Assume now that shortwave flux from ocean model has been modified for ice. mac, aug12 +! radbio = max(0.0,parbio(i,j)*sw_zt) +! well, now doing an experiment without ice, so assumption no longer true. mac, oct12. + + vpbio(i,k) = abio(i,j)*bbio(i,j)** & + (cbio(i,j)*t_prog(index_temp)%field(i,j,k,time%tau)) + +! vpbio(i,k)= 0.6* (1.006)**T_prog(index_temp)%field(i,j,k,Time%tau) + +! avej(i,k) = vpbio(i,k)*alphabio*radbio/(vpbio(i,k)+(alphabio*radbio)**2) +! Growth term from Evans and Parslow 1985, mac Apr10 +! avej(i,k) = vpbio(i,k)*alphabio*radbio/(vpbio(i,k)**2+(alphabio*radbio)**2)**0.5 + +! Growth term from Brian Griffiths + avej(i,k) = vpbio(i,k)*(1.0-exp(-1.0*(alphabio(i,j)*radbio)/vpbio(i,k))) + + if (Grid%zw(k) .le. mld(i,j)) & + light_limit(i,j)=light_limit(i,j) + & + (1.0-exp(-1.0*(alphabio(i,j)*radbio)/vpbio(i,k))) * & + thickness%dzt(i,j,k) + +! avej(i,k)=0. ! zero phyto growth + enddo + enddo + + +!chd begin time stepping over dtts in NPZD model with Euler forward +!chd This is the NPZD model: +!chd (P: phytoplankton, Z: Zooplankton, N: Nitrate and D: Detritus) +!chd +!chd dP/dt = u(N,Temp.,Light) P - p_P P - g(P) P Z +!chd +!chd dZ/dt = a g(P) P Z - d Z - p_Z Z^2 +!chd +!chd dN/dt = r D + d Z - u(N,Temp.,Light) P [ + r_d DOC ] +!chd +!chd dD/dt = (1-s)[ (1-a) g(P) P Z + p_P P + p_Z Z^2] -r D + w_D dD/dz + + do tn = 1,ts_npzd !{ + do k=1,grid%nk !{ + do i=isc,iec !{ + bion = max(0.0,biotr(i,k,id_no3)) + biop = max(0.0,biotr(i,k,id_phy)) + bioz = max(0.0,biotr(i,k,id_zoo)) + biod = max(0.0,biotr(i,k,id_det)) + bioo = max(0.0,biotr(i,k,id_o2)) + if (id_caco3.ne.0) biocaco3 = max(0.0,biotr(i,k,id_caco3)) + if (id_fe.ne.0) bioi = max(0.0,biotr(i,k,id_fe)) + +!chd -- phytoplankton equation +!chdc +!chdc use Liebigs Law of the Minimum (Liebig, 1845) for growth rate +!chdc (minimum of light-limited and nutrient limited growth rates; +!chdc although chlorophyll is not explicitly considered, this will +!chdc later allow for a diagnostic determination of a Chl:N ratio +!chdc depending on light- or nutrient-limited growth. +!chdc --> Hurtt and Armstrong, 1996) +!chdc saturation growth rate (infinite light, infinite nutrients) +!chd vpbio = abio*bbio**(cbio*t_prog(index_temp)%field(i,j,k,tau)) +!chd growth rate + u_npz = min(avej(i,k),vpbio(i,k)*bion/(k1bio(i,j)+bion)) +! rjm - iron limitation + if(id_fe.ne.0) u_npz = min(u_npz,vpbio(i,k)*bioi/(0.1+bioi)) + +!chd if (k .eq. 1) then +!chd write(stdout(),*) 'ave' +!chd write(stdout(),*) avej(i,k) +!chd endif +!chd grazing function + g_npz = gbio(i,j)*epsbio(i,j)*biop*biop/(gbio(i,j)+epsbio(i,j)*biop*biop) +!chd temperature dependance of growth rates + fbc = bbio(i,j)**(cbio(i,j)*t_prog(index_temp)%field(i,j,k,time%tau)) + + f11 = u_npz*biop * bioma(i,k,id_no3) + f21 = g_npz*bioz * bioma(i,k,id_phy) + f22 = muepbio(i,j)*fbc*biop * bioma(i,k,id_phy) + f23 = muepsbio(i,j)*biop*biop * bioma(i,k,id_phy) + f31 = gam2bio(i,j)*fbc*bioz * bioma(i,k,id_zoo) + f32 = muezbio(i,j)*bioz*bioz * bioma(i,k,id_zoo) + f41 = muedbio(i,j)*fbc*biod * bioma(i,k,id_det) +! if (grid%zw(k) .ge. 180) f41 = f41*.2 ! reduce decay below 300m +! 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 (id_caco3.ne.0) then + f51 = muecaco3(i,j)*biocaco3 *bioma(i,k,id_caco3) + endif + +!chd -- nutrient equation + biotr(i,k,id_no3) = biotr(i,k,id_no3) + dtsb * ( & + f41 + f31 + f22 - f11) + +!chd -- phyto plankton equation + biotr(i,k,id_phy) = biotr(i,k,id_phy) + dtsb * ( & + f11 - f21 - f22 - f23) + +! Estimate primary productivity from phytoplankton growth + pprod_gross(i,j,k) = pprod_gross(i,j,k) + dtsb*f11 + +!chd -- zooplankton equation + biotr(i,k,id_zoo) = biotr(i,k,id_zoo) + dtsb * ( & + gam1bio(i,j)*f21 - f31 - f32) + +! Estimate secondary productivity from zooplankton growth + zprod_gross(i,j,k) = zprod_gross(i,j,k) + dtsb*f21 + +!chd -- detritus equation + biotr(i,k,id_det) = biotr(i,k,id_det) + dtsb * ( & + (1-gam1bio(i,j))*f21 + f23 + f32 - f41) + +!chd -- oxygen equation + biotr(i,k,id_o2) = biotr(i,k,id_o2) - bioma(i,k,id_o2) * 172 / 16 * dtsb * ( & + f41 + f31 + f22 - f11) + +! rjm -- extra equation for caco3 - alkalinity + if (id_caco3.ne.0) & + biotr(i,k,id_caco3) = biotr(i,k,id_caco3) + dtsb * ( ( & + (1-gam1bio(i,j))*f21 + f23 + f32)*f_inorg(i,j)*106/16 - f51) ! 8% of POC 106/16*.08 + +! rjm -- extra equation for iron +! if (id_fe.ne.0) & +! biotr(i,k,id_fe) = biotr(i,k,id_fe) +dtsb * 200./16.* ( & +! f41 + f31 + f22 - f11) + +! mac -- extra equation for iron, molar Fe:N = 1.98e-5:1.0 (Christian et al. 2002), and Fe units are micro mole/m3 cf milli mole/m3 for others. + if (id_fe.ne.0) & + biotr(i,k,id_fe) = biotr(i,k,id_fe) +dtsb * 2.0e-2* ( & + f41 + f31 + f22 - f11) + + + enddo !} i + enddo !} k + enddo !} tn + + +!chd add biotically induced tendency to biotracers, + do k = 1, grid%nk !{ + do i = isc, iec !{ + + no3_bgc_change = grid%tmask(i,j,k)* & + (biotr(i,k,id_no3) - max(0.0,t_prog(ind_no3)%field(i,j,k,time%taum1))) + caco3_bgc_change = grid%tmask(i,j,k)* & + (biotr(i,k,id_caco3) - max(0.0,t_prog(ind_caco3)%field(i,j,k,time%taum1))) + + t_prog(ind_no3)%field(i,j,k,time%taum1) = biotr(i,k,id_no3) + t_prog(ind_phy)%field(i,j,k,time%taum1) = biotr(i,k,id_phy) + t_prog(ind_zoo)%field(i,j,k,time%taum1) = biotr(i,k,id_zoo) + t_prog(ind_det)%field(i,j,k,time%taum1) = biotr(i,k,id_det) + t_prog(ind_o2)%field(i,j,k,time%taum1) = biotr(i,k,id_o2) + if ( id_caco3 .ne. 0 ) & + t_prog(ind_caco3)%field(i,j,k,time%taum1) = biotr(i,k,id_caco3) + if ( id_fe .ne. 0 ) & + t_prog(ind_fe)%field(i,j,k,time%taum1) = biotr(i,k,id_fe) - & + dtts * tscav_fe(i,j) * max(0.0,(biotr(i,k,id_fe) - fe_bkgnd(i,j)) ) + + if (id_dic.ne.0) & + t_prog(ind_dic)%field(i,j,k,time%taum1) = & + t_prog(ind_dic)%field(i,j,k,time%taum1) + & + 106./16. * no3_bgc_change - caco3_bgc_change + + if (id_adic.ne.0) & + t_prog(ind_adic)%field(i,j,k,time%taum1) = & + t_prog(ind_adic)%field(i,j,k,time%taum1) + & + 106./16. * no3_bgc_change - caco3_bgc_change + + if (id_alk.ne.0) & + t_prog(ind_alk)%field(i,j,k,time%taum1) = & + t_prog(ind_alk)%field(i,j,k,time%taum1) + & + ( - 2.0 * caco3_bgc_change - no3_bgc_change) + + pprod_gross(i,j,k)=rdtts*pprod_gross(i,j,k)*grid%tmask(i,j,k) + zprod_gross(i,j,k)=rdtts*zprod_gross(i,j,k)*grid%tmask(i,j,k) + + + enddo !} i + enddo !} k + +!RASF upstream sinking of detritus +! change field(..,k,..) to field(..,k-1,..), mac apr10 + do k=2,grid%nk+1 + do i=isc,iec +! adv_fb(i,k)=wdetbio(i,j)*t_prog(ind_det)%field(i,j,k-1,time%taum1) + adv_fb(i,k)=wdetbio(i,j)*biotr(i,k-1,id_det) + enddo + enddo + +!RASF no flux boundary conditions + do i=isc,iec + adv_fb(i,1) = 0.0 + enddo + +! Deposit tracer to sediment as tracer sinks through base of column. mac, nov12 + do i = isc, iec + k = grid%kmt(i,j) + if (k .gt. 0) then + + biotic(n)%det_sed_depst(i,j) = adv_fb(i,k+1) + + endif ! k .gt. 0 + enddo ! i + + +!!mac remineralise tracer sinking through the base of the column. +! do i = isc, iec +! k = grid%kmt(i,j) +! if (k .gt. 0) then +! +! t_prog(ind_no3)%th_tendency(i,j,k) = t_prog(ind_no3)%th_tendency(i,j,k) + & +! grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,time%tau) & +! * adv_fb(i,k+1)/Thickness%dzt(i,j,k) +! +! t_prog(ind_o2)%th_tendency(i,j,k) = t_prog(ind_o2)%th_tendency(i,j,k) - & +! 172.0 / 16.0 * grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,time%tau) & +! * adv_fb(i,k+1)/Thickness%dzt(i,j,k) +! +! if (id_fe .ne. 0) then +! t_prog(ind_fe)%th_tendency(i,j,k) = t_prog(ind_fe)%th_tendency(i,j,k) + & +! 2.0e-2 * grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,time%tau) & +! * adv_fb(i,k+1)/Thickness%dzt(i,j,k) +! endif +! +! if (id_dic .ne. 0) then +! t_prog(ind_dic)%th_tendency(i,j,k) = t_prog(ind_dic)%th_tendency(i,j,k) + & +! 106.0 / 16.0 * grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,time%tau) & +! * adv_fb(i,k+1)/Thickness%dzt(i,j,k) +! endif +! +! if (id_adic .ne. 0) then +! t_prog(ind_adic)%th_tendency(i,j,k) =t_prog(ind_adic)%th_tendency(i,j,k) + & +! 106.0 / 16.0 * grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,time%tau) & +! * adv_fb(i,k+1)/Thickness%dzt(i,j,k) +! endif +! +! if (id_alk .ne. 0) then +! t_prog(ind_alk)%th_tendency(i,j,k) = t_prog(ind_alk)%th_tendency(i,j,k) - & +! grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,time%tau) & +! *adv_fb(i,k+1)/Thickness%dzt(i,j,k) +! endif +! +! endif ! k .gt. 0 +! +! enddo ! i + + + do k =1,grid%nk !{ + do i =isc,iec !{ + t_prog(ind_det)%field(i,j,k,time%taum1) = t_prog(ind_det)%field(i,j,k,time%taum1) + & + grid%tmask(i,j,k) * dtts * & + (-adv_fb(i,k+1) + adv_fb(i,k))/Thickness%dzt(i,j,k) + enddo !} i + enddo !} k + +!RASF upstream sinking of caco3 +if (id_caco3.ne.0) then + do k=2,grid%nk+1 + do i=isc,iec +! adv_fb(i,k)=wcaco3(i,j)*t_prog(ind_caco3)%field(i,j,k-1,time%taum1) + adv_fb(i,k)=wcaco3(i,j)*biotr(i,k-1,id_caco3) + enddo + enddo +! no flux boundary conditions + do i=isc,iec + adv_fb(i,1) = 0.0 + enddo + + +! Deposit tracer to sediment as tracer sinks through base of column. mac, nov12 + do i = isc, iec + k = grid%kmt(i,j) + if (k .gt. 0) then + + biotic(n)%caco3_sed_depst(i,j) = adv_fb(i,k+1) + + endif ! k .gt. 0 + enddo ! i + + +!!mac remineralise tracer sinking through the base of the column. +! do i = isc, iec +! k = grid%kmt(i,j) +! if (k .gt. 0) then +! +! if (id_dic .ne. 0) then +! t_prog(ind_dic)%th_tendency(i,j,k) = t_prog(ind_dic)%th_tendency(i,j,k) + & +! grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,time%tau) & +! * adv_fb(i,k+1)/Thickness%dzt(i,j,k) +! endif +! +! if (id_adic .ne. 0) then +! t_prog(ind_adic)%th_tendency(i,j,k) =t_prog(ind_adic)%th_tendency(i,j,k) + & +! grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,time%tau) & +! * adv_fb(i,k+1)/Thickness%dzt(i,j,k) +! endif +! +! if (id_alk .ne. 0) then +! t_prog(ind_alk)%th_tendency(i,j,k) = t_prog(ind_alk)%th_tendency(i,j,k) + & +! 2.0 * grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,time%tau) & +! *adv_fb(i,k+1)/Thickness%dzt(i,j,k) +! endif +! endif ! k .gt. 0 +! enddo ! i +! + + + do k =1,grid%nk !{ + do i =isc,iec !{ + t_prog(ind_caco3)%field(i,j,k,time%taum1)=t_prog(ind_caco3)%field(i,j,k,time%taum1) +& + grid%tmask(i,j,k) * dtts * & + (-adv_fb(i,k+1) + adv_fb(i,k))/Thickness%dzt(i,j,k) + enddo !} i + enddo !} k +endif ! end loop for caco3 + + + + + enddo !} j + +enddo !} n + +return +end subroutine bio_v3 diff --git a/src/mom5/ocean_csiro_bgc/csiro_bgc.F90 b/src/mom5/ocean_csiro_bgc/csiro_bgc.F90 new file mode 100755 index 0000000000..af4e7178fa --- /dev/null +++ b/src/mom5/ocean_csiro_bgc/csiro_bgc.F90 @@ -0,0 +1,2523 @@ +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! !! +!! GNU General Public License !! +!! !! +!! This file is part of the Flexible Modeling System (FMS). !! +!! !! +!! FMS is free software; you can redistribute it and/or modify !! +!! it and are expected to follow the terms of the GNU General Public !! +!! License as published by the Free Software Foundation. !! +!! !! +!! FMS is distributed in the hope that it will be useful, !! +!! but WITHOUT ANY WARRANTY; without even the implied warranty of !! +!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! +!! GNU General Public License for more details. !! +!! !! +!! You should have received a copy of the GNU General Public License !! +!! along with FMS; if not, write to: !! +!! Free Software Foundation, Inc. !! +!! 59 Temple Place, Suite 330 !! +!! Boston, MA 02111-1307 USA !! +!! or see: !! +!! http://www.gnu.org/licenses/gpl.txt !! +!! !! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +! +! +! Richard Matear +! +! +! Matt Chamberlain +! +! +! +! +! +! +! CSIRO v3 bgc model +! +! +! +! Generic setup of the CSIRO BGC model; NO3-based NPZD cycle, +! with coupled carbon cycle +! +! +! +! +! +! +! +! +! +! +! +!------------------------------------------------------------------ +! +! Module csiro_bgc_mod +! +! csiro bgc model version 3 +! +!------------------------------------------------------------------ + +module csiro_bgc_mod !{ + +!------------------------------------------------------------------ +! +! Global definitions +! +!------------------------------------------------------------------ + +!---------------------------------------------------------------------- +! +! Modules +! +!---------------------------------------------------------------------- + +use diag_manager_mod, only: send_data +use field_manager_mod, only: fm_field_name_len, fm_path_name_len, fm_string_len +use field_manager_mod, only: fm_get_length, fm_get_value, fm_new_value +use field_manager_mod, only: fm_get_index +use fms_mod, only: write_data, write_version_number +use mpp_mod, only: stdout, stdlog, mpp_error, mpp_sum, FATAL +use mpp_mod, only: mpp_pe, mpp_root_pe, mpp_sync +use mpp_mod, only: mpp_clock_id, mpp_clock_begin, mpp_clock_end, CLOCK_ROUTINE +use fms_io_mod, only: register_restart_field, save_restart, restore_state +use fms_io_mod, only: restart_file_type, reset_field_pointer + +use ocean_types_mod, only: ocean_domain_type, ocean_grid_type +use mpp_domains_mod, only: mpp_global_field, mpp_global_sum, BITWISE_EXACT_SUM + +use time_manager_mod, only: get_date +use time_interp_external_mod, only: time_interp_external + +use ocean_parameters_mod, only: rho0 + +use ocean_tpm_util_mod, only: otpm_set_tracer_package, otpm_set_prog_tracer +use fm_util_mod, only: fm_util_check_for_bad_fields, fm_util_set_value +use fm_util_mod, only: fm_util_get_string, fm_util_get_logical, fm_util_get_integer, fm_util_get_real +use fm_util_mod, only: fm_util_get_logical_array, fm_util_get_real_array, fm_util_get_string_array +use fm_util_mod, only: fm_util_start_namelist, fm_util_end_namelist +! use fm_util_mod, only: domain, grid, time, dtts +! use fm_util_mod, only: isc, iec, jsc, jec, nk, isd, ied, jsd, jed +! use fm_util_mod, only: taum1, tau, taup1 +! use fm_util_mod, only: t_prog, t_diag +! use fm_util_mod, only: indsal, indtemp +! use fm_util_mod, only: end_of_year, end_of_month + +use ocean_types_mod, only: ocean_thickness_type +use ocean_types_mod, only: ocean_density_type +use ocean_types_mod, only: ocean_grid_type, ocean_domain_type +use ocean_types_mod, only: ocean_prog_tracer_type, ocean_diag_tracer_type +use ocean_types_mod, only: ocean_time_type + +use ocean_types_mod, only: ocean_public_type + +use ocean_tracer_diag_mod, only: calc_mixed_layer_depth + +!---------------------------------------------------------------------- +! force all variables to be "typed" +!---------------------------------------------------------------------- + +implicit none + +!---------------------------------------------------------------------- +! +! Make all routines and variables private by default +! +!---------------------------------------------------------------------- +private + +!---------------------------------------------------------------------- +! +! Public routines +! +!---------------------------------------------------------------------- +! order of call init, start, loop(sbc, source, bbc, tracer), end +public :: csiro_bgc_bbc +public :: csiro_bgc_end +public :: csiro_bgc_init +public :: csiro_bgc_sbc +public :: csiro_bgc_source +public :: csiro_bgc_start +public :: csiro_bgc_tracer +public :: csiro_bgc_virtual_fluxes + +!---------------------------------------------------------------------- +! +! Private routines +! +!---------------------------------------------------------------------- + +private :: allocate_arrays + +!---------------------------------------------------------------------- +! +! Private parameters +! +!---------------------------------------------------------------------- + +character(len=fm_field_name_len), parameter :: package_name = 'csiro_bgc' +character(len=48), parameter :: mod_name = 'csiro_bgc_mod' +character(len=fm_string_len), parameter :: default_file_in = 'INPUT/csiro_bgc.res.nc' +character(len=fm_string_len), parameter :: default_file_out = 'RESTART/csiro_bgc.res.nc' + +integer, parameter :: ntr_bmax = 10 + +!------------------------------------------------------- +! private types +!-------------------------------------------------------- + +type biotic_type !{ +! set-up tracer indexes + integer :: ntr_bgc = ntr_bmax + integer, dimension(ntr_bmax) :: ind_bgc + integer, dimension(ntr_bmax) :: id_bgc_stf + integer, dimension(ntr_bmax) :: id_bgc_btf + integer, dimension(ntr_bmax) :: id_bgc_vstf + integer, dimension(ntr_bmax) :: id_bgc_src + + logical :: init + character(len=fm_field_name_len) :: name +! arrays for air-sea fluxes + real :: sal_global + real, allocatable, dimension(:) :: bgc_global + real, allocatable, dimension(:,:) :: htotal + real, allocatable, dimension(:,:) :: alpha + real, allocatable, dimension(:,:) :: csat + real, allocatable, dimension(:,:) :: csat_csurf + real, allocatable, dimension(:,:) :: csat_acsurf + real, allocatable, dimension(:,:) :: csurf + real, allocatable, dimension(:,:) :: acsurf + real, allocatable, dimension(:,:) :: dpco2 + real, allocatable, dimension(:,:) :: pco2surf + real, allocatable, dimension(:,:) :: paco2surf + real, allocatable, dimension(:,:) :: pco2atm + real, allocatable, dimension(:,:) :: paco2atm + real, allocatable, dimension(:,:) :: det_sediment ! mmol(NO3) m-2 in DET sitting at base of column as sediment. + real, allocatable, dimension(:,:) :: caco3_sediment ! mmol(CaCO3) m-2 sitting at base of column as sediment. + real, allocatable, dimension(:,:) :: det_sed_remin ! mmol(NO3) m-2 s-1, rate of remineralisation of DET in sediment. + real, allocatable, dimension(:,:) :: caco3_sed_remin ! mmol m-2 s-1, rate of remineralisation of CaCO3 in sediment. + real, allocatable, dimension(:,:) :: det_sed_depst ! mmol(NO3) m-2 s-1, rate of deposition of DET in sediment. + real, allocatable, dimension(:,:) :: caco3_sed_depst ! mmol m-2 s-1, rate of deposition of CaCO3 in sediment. + real, allocatable, dimension(:,:) :: sio2 + real, allocatable, dimension(:,:) :: po4 + real, allocatable, dimension(:,:,:) :: vstf +end type biotic_type !} + +!---------------------------------------------------------------------- +! Public variables +!---------------------------------------------------------------------- + +logical, public :: do_csiro_bgc + +!---------------------------------------------------------------------- +! Private variables +!---------------------------------------------------------------------- +integer :: package_index + +! set the tracer index for the various tracers +integer :: id_po4, id_dic, id_alk, id_o2, id_no3, id_phy, id_det, id_zoo & + , id_caco3, id_adic, id_fe, id_caco3_sediment, id_det_sediment +! internal pointer to make reading the code easier +integer :: ind_po4, ind_dic, ind_alk, ind_o2, ind_no3, ind_phy, ind_det, ind_zoo & + , ind_caco3, ind_adic, ind_fe +character*6 :: qbio_model +integer :: bio_version ! version of the bgc module to use +logical :: zero_floor ! apply hard floor to bgc tracers +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. +logical :: use_access_co2=.false. ! determine whether to use ACCESS atmospheric CO2 or a CO2 file to drive the ocean's anthropogenic CO2 tracer. mac, may13. + +integer :: id_clock_csiro_obgc + +integer :: atmpress_id +character*128 :: atmpress_file +character*32 :: atmpress_name +real, allocatable, dimension(:,:) :: fice_t +integer :: id_light_limit = -1 +integer :: id_pprod_gross = -1 +integer :: id_pprod_gross_2d = -1 +integer :: id_zprod_gross = -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_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 +integer :: id_total_aco2_flux, id_total_co2_flux +real, allocatable, dimension(:,:) :: kw_co2 +real, allocatable, dimension(:,:) :: kw_o2 +real, allocatable, dimension(:,:) :: patm_t +integer :: pistonveloc_id +integer :: aco2_id +integer :: seaicefract_id +character*128 :: pistonveloc_file +character*32 :: pistonveloc_name +character*128 :: aco2_file +character*32 :: aco2_name +real, allocatable, dimension(:,:) :: sc_o2 +real, allocatable, dimension(:,:) :: sc_co2 +real, allocatable, dimension(:,:) :: o2_saturation +character*128 :: seaicefract_file +character*32 :: seaicefract_name +character*128 :: dust_file +character*32 :: dust_name +integer :: dust_id +real, allocatable, dimension(:,:) :: dust_t + +real, allocatable, dimension(:,:) :: xkw_t +real, allocatable, dimension(:,:) :: aco2 +real, allocatable, dimension(:,:) :: htotalhi +real, allocatable, dimension(:,:) :: htotallo + +type(biotic_type), allocatable, dimension(:) :: biotic +integer :: instances +real, allocatable, dimension(:) :: tk +real, allocatable, dimension(:) :: ts +real, allocatable, dimension(:) :: ts2 +real, allocatable, dimension(:) :: ts3 +real, allocatable, dimension(:) :: ts4 +real, allocatable, dimension(:) :: ts5 +real, allocatable, dimension(:) :: tt + + +! rjm biotic parameters +integer :: n_eudepth=4 +real :: p_k = 0.1 +real :: rain_ratio = 8.48 +real :: s_npp = -1 +real :: ratio_poc(ntr_bmax) +real :: ratio_pic(ntr_bmax) + +real, allocatable, dimension(:,:,:) :: poc +real, allocatable, dimension(:,:,:) :: pic +real, allocatable, dimension(:,:) :: poc_tot +real, allocatable, dimension(:,:) :: pic_tot +real, allocatable, dimension(:,:) :: pmax_growth +real, allocatable, dimension(:,:) :: pp +real, allocatable, dimension(:) :: fmin_poc +real, allocatable, dimension(:) :: fmin_pic +real, allocatable, dimension(:,:,:) :: biotr +real, allocatable, dimension(:,:) :: light_limit +real, allocatable, dimension(:,:,:) :: pprod_gross +real, allocatable, dimension(:,:) :: pprod_gross_2d +real, allocatable, dimension(:,:,:) :: zprod_gross +real, allocatable, dimension(:) :: ray +real, allocatable, dimension(:) :: dummy +real :: dummy1 +real, allocatable, dimension(:) :: tracer_sources +real, allocatable, dimension(:,:) :: area_k +real, allocatable, dimension(:,:) :: tmp + +integer :: global_sum_flag ! flag for mpp_global_sum + + +integer :: alphabio_id +real, allocatable, dimension(:,:) :: alphabio +integer :: parbio_id +real, allocatable, dimension(:,:) :: parbio +integer :: kwbio_id +real, allocatable, dimension(:,:) :: kwbio +integer :: kcbio_id +real, allocatable, dimension(:,:) :: kcbio +integer :: abio_id +real, allocatable, dimension(:,:) :: abio +integer :: bbio_id +real, allocatable, dimension(:,:) :: bbio +integer :: cbio_id +real, allocatable, dimension(:,:) :: cbio +integer :: k1bio_id +real, allocatable, dimension(:,:) :: k1bio +integer :: muepbio_id +real, allocatable, dimension(:,:) :: muepbio +integer :: muepsbio_id +real, allocatable, dimension(:,:) :: muepsbio +integer :: gam1bio_id +real, allocatable, dimension(:,:) :: gam1bio +integer :: gbio_id +real, allocatable, dimension(:,:) :: gbio +integer :: epsbio_id +real, allocatable, dimension(:,:) :: epsbio +integer :: muezbio_id +real, allocatable, dimension(:,:) :: muezbio +integer :: gam2bio_id +real, allocatable, dimension(:,:) :: gam2bio +integer :: muedbio_id +real, allocatable, dimension(:,:) :: muedbio +integer :: muecaco3_id +real, allocatable, dimension(:,:) :: muecaco3 +integer :: muedbio_sed_id +real, allocatable, dimension(:,:) :: muedbio_sed +integer :: muecaco3_sed_id +real, allocatable, dimension(:,:) :: muecaco3_sed +integer :: wdetbio_id +real, allocatable, dimension(:,:) :: wdetbio +integer :: wcaco3_id +real, allocatable, dimension(:,:) :: wcaco3 +integer :: nat_co2_id +real, allocatable, dimension(:,:) :: nat_co2 +integer :: tscav_fe_id +real, allocatable, dimension(:,:) :: tscav_fe +integer :: fe_bkgnd_id +real, allocatable, dimension(:,:) :: fe_bkgnd +integer :: f_inorg_id +real, allocatable, dimension(:,:) :: f_inorg + + +! for extra restart file(s) +integer :: id_restart(2)=0 +type(restart_file_type), save :: sed_restart + + +! Tracer names + +character(5), dimension(10) :: tracer_name + +!----------------------------------------------------------------------- +! +! Subroutine and function definitions +! +!----------------------------------------------------------------------- + +contains + +!####################################################################### +! +! +! +! Dynamically allocate arrays +! +! + +subroutine allocate_arrays (isc, iec, jsc, jec, isd, ied, jsd, jed, nk) !{ + +! arguments + +! local variables + +integer, intent(in) :: isd, isc +integer, intent(in) :: ied, iec +integer, intent(in) :: jsd, jsc +integer, intent(in) :: jed, jec +integer, intent(in) :: nk +integer :: m +integer :: k +integer :: j +integer :: i +integer :: n +integer :: ntr_bgc + +!----------------------------------------------------------------------- +! start executable code +!----------------------------------------------------------------------- + +allocate( xkw_t(isd:ied,jsd:jed) ) +allocate( patm_t(isd:ied,jsd:jed) ) +allocate( fice_t(isd:ied,jsd:jed) ) +allocate( aco2(isd:ied,jsd:jed) ) +allocate( dust_t(isd:ied,jsd:jed) ) + +allocate( sc_o2(isc:iec,jsc:jec) ) +allocate( sc_co2(isc:iec,jsc:jec) ) +allocate( kw_o2(isc:iec,jsc:jec) ) +allocate( kw_co2(isc:iec,jsc:jec) ) +allocate( htotalhi(isd:ied,jsd:jed) ) +allocate( htotallo(isd:ied,jsd:jed) ) + + +allocate( o2_saturation(isc:iec,jsc:jec) ) + +allocate( tt(isc:iec) ) +allocate( tk(isc:iec) ) +allocate( ts(isc:iec) ) +allocate( ts2(isc:iec) ) +allocate( ts3(isc:iec) ) +allocate( ts4(isc:iec) ) +allocate( ts5(isc:iec) ) + +allocate( poc(isc:iec,jsc:jec,nk) ) +allocate( pic(isc:iec,jsc:jec,nk) ) +allocate( poc_tot(isc:iec,jsc:jec) ) +allocate( pic_tot(isc:iec,jsc:jec) ) +allocate( pmax_growth(isc:iec,nk) ) +allocate( pp(isc:iec,nk) ) +allocate( fmin_poc(nk) ) +allocate( fmin_pic(nk) ) + ntr_bgc = biotic(1)%ntr_bgc +allocate( ray(nk) ) +allocate( biotr(isc:iec,nk,ntr_bgc) ) +allocate( light_limit(isc:iec,jsc:jec) ) +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 (tmp(isd:ied,jsd:jed) ) +allocate ( tracer_sources(0:nk) ) +allocate(area_k(isd:ied,jsd:jed) ) +allocate( dummy(isc:iec) ) + +! initialize some arrays + +xkw_t(:,:) = 0.0 +patm_t(:,:) = 0.0 +fice_t(:,:) = 0.0 +dust_t(:,:) = 0.0 +aco2(:,:) = 0.0 +sc_co2(:,:) = 0.0 +kw_co2(:,:) = 0.0 +sc_o2(:,:) = 0.0 +kw_o2(:,:) = 0.0 +o2_saturation(:,:) = 0.0 + +tt(:) = 0.0 +tk(:) = 0.0 +ts(:) = 0.0 +ts2(:) = 0.0 +ts3(:) = 0.0 +ts4(:) = 0.0 +ts5(:) = 0.0 + +! allocate biotic array elements + +do n = 1, instances !{ +! rjm problem - the allocation of the arrays is done +! after the initialization routine which cause problems +! solution - move the call to array allocation to the init subroutine +! ntr_bgc = biotic(n)%ntr_bgc +! allocate( biotic(n)%ind_bgc(1:ntr_bgc) ) + + allocate( biotic(n)%htotal(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) ) + allocate( biotic(n)%csat_acsurf(isd:ied,jsd:jed) ) + allocate( biotic(n)%csurf(isd:ied,jsd:jed) ) + allocate( biotic(n)%acsurf(isd:ied,jsd:jed) ) + allocate( biotic(n)%dpco2(isd:ied,jsd:jed) ) + allocate( biotic(n)%pco2atm(isd:ied,jsd:jed) ) + allocate( biotic(n)%paco2atm(isd:ied,jsd:jed) ) + allocate( biotic(n)%pco2surf(isd:ied,jsd:jed) ) + allocate( biotic(n)%paco2surf(isd:ied,jsd:jed) ) + allocate( biotic(n)%caco3_sediment(isd:ied,jsd:jed) ) + allocate( biotic(n)%det_sediment(isd:ied,jsd:jed) ) + allocate( biotic(n)%caco3_sed_remin(isd:ied,jsd:jed) ) + allocate( biotic(n)%det_sed_remin(isd:ied,jsd:jed) ) + allocate( biotic(n)%caco3_sed_depst(isd:ied,jsd:jed) ) + allocate( biotic(n)%det_sed_depst(isd:ied,jsd:jed) ) + allocate( biotic(n)%sio2(isd:ied,jsd:jed) ) + allocate( biotic(n)%po4(isd:ied,jsd:jed) ) + + ntr_bgc = biotic(n)%ntr_bgc + allocate(biotic(n)%bgc_global(1:ntr_bgc) ) + allocate(biotic(n)%vstf(isc:iec,jsc:jec,1:ntr_bgc) ) + +enddo !} + +! allocate biotic parameters + +allocate( alphabio(isd:ied,jsd:jed) ) +allocate( parbio(isd:ied,jsd:jed) ) +allocate( kwbio(isd:ied,jsd:jed) ) +allocate( kcbio(isd:ied,jsd:jed) ) +allocate( abio(isd:ied,jsd:jed) ) +allocate( bbio(isd:ied,jsd:jed) ) +allocate( cbio(isd:ied,jsd:jed) ) +allocate( k1bio(isd:ied,jsd:jed) ) +allocate( muepbio(isd:ied,jsd:jed) ) +allocate( muepsbio(isd:ied,jsd:jed) ) +allocate( gam1bio(isd:ied,jsd:jed) ) +allocate( gbio(isd:ied,jsd:jed) ) +allocate( epsbio(isd:ied,jsd:jed) ) +allocate( muezbio(isd:ied,jsd:jed) ) +allocate( gam2bio(isd:ied,jsd:jed) ) +allocate( muedbio(isd:ied,jsd:jed) ) +allocate( muecaco3(isd:ied,jsd:jed) ) +allocate( muedbio_sed(isd:ied,jsd:jed) ) +allocate( muecaco3_sed(isd:ied,jsd:jed) ) +allocate( wdetbio(isd:ied,jsd:jed) ) +allocate( wcaco3(isd:ied,jsd:jed) ) +allocate( nat_co2(isd:ied,jsd:jed) ) +allocate( tscav_fe(isd:ied,jsd:jed) ) +allocate( fe_bkgnd(isd:ied,jsd:jed) ) +allocate( f_inorg(isd:ied,jsd:jed) ) + + +! initialize some arrays + +do n = 1, instances !{ + biotic(n)%htotal(:,:) = 1.e-8 + biotic(n)%sio2(:,:) = 35. *1e-3 + biotic(n)%bgc_global(:) = 0. ! this will make vstf zero + biotic(n)%sal_global = 35. + biotic(n)%vstf(:,:,:) = 0. + biotic(n)%caco3_sediment(:,:) = 0.0 + biotic(n)%det_sediment(:,:) = 0.0 +enddo !} n + +return +end subroutine allocate_arrays !} +! NAME="allocate_arrays" + + +!####################################################################### +! +! +! +! compute bottom boundary conditions e.g. fluxes due to +! interaction with the sediment +! +! + +subroutine csiro_bgc_bbc (isc, iec, jsc, jec, T_prog, grid, Time) !{ + +integer, intent(in) :: isc, iec +integer, intent(in) :: jsc, jec +type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog +type(ocean_grid_type), intent(in) :: Grid +type(ocean_time_type), intent(in) :: Time + + +!----------------------------------------------------------------------- +! local definitions +!----------------------------------------------------------------------- + +character(len=64), parameter :: sub_name = 'csiro_bgc_bbc' +character(len=256), parameter :: error_header = & + '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: warn_header = & + '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: note_header = & + '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + +integer :: i, j, n, k, nn +integer :: ind_temp +real :: fbc +logical :: used + +! ===================================================================== +! begin executable code +! ===================================================================== + +! +!! rjm -- relax Fe to 1 at ocean bottom with 1 day time-constant +!! Huh? This is setting it to zero. RASF +!! Should we delete commented out code? RASF +!! Yes, fe%btf() values are left as zero... +!! bottom values for fe are assigned in csiro_bgc_tracer subroutine... +!! so commenting out this block. mac, nov10. +! +!do n = 1, instances !{ +! if (id_fe.ne.0) then +! ind_fe = biotic(n)%ind_bgc(id_fe) +! do j = jsc, jec !{ +! do i = isc, iec !{ +! k = grid%kmt(i,j) +! t_prog(ind_fe)%btf(i,j) = rho0 * -1./86400. *0 ! negative into the ocean +!!grid%tmask(i,j,k) & +!! / (-1.) * & +!! (1 - max(0.0,t_prog(ind_fe)%field(i,j,k,taum1))) +! enddo !} i +! enddo !} j +! endif +!enddo !} n +! + +!! Trial fe source within csiro_bgc_bbc. mac, nov12. +!do n = 1, instances !{ +! if (id_fe.ne.0) then +! do j = jsc, jec !{ +! do i = isc, iec !{ +! k = grid%kmt(i,j) +! if (grid%zw(k) .le. 200) & +!t_prog(ind_fe)%btf(i,j) = -1.0 * rho0 * 1000. / 86400. +!! t_prog(ind_fe)%btf(i,j) = -1.0 * rho0 * (0.999 - t_prog(ind_fe)%field(i,j,k,time%taum1)) * & +!! grid%dzt(k) / (0.1 * 86400) ! flux ~ density * diff(concentration) * dz / time_scale +! ! setting time scale initial to 0.1 of a day, since at present, the concentration of the bottom is fixed in csiro_bgc_tracer; +! ! and want the time scale to be shorter (stronger) than the scavenging to 0.6 within csiro_bgc_source. +! ! mac, nov12. +! enddo !} i +! enddo !} j +! endif +!enddo !} n + + + + + + ! find remineralisation rate of sediment tracers. mac, nov12. + call time_interp_external(muedbio_sed_id, time%model_time, muedbio_sed) + call time_interp_external(muecaco3_sed_id, time%model_time, muecaco3_sed) + ind_temp = fm_get_index('/ocean_mod/prog_tracers/temp') + + do n = 1, instances !{ + do j = jsc, jec !{ + do i = isc, iec !{ + k = grid%kmt(i,j) + if (k .gt. 0) then + fbc= bbio(i,j)**(cbio(i,j)*T_prog(ind_temp)%field(i,j,k,time%taum1)) + biotic(n)%det_sed_remin(i,j) = muedbio_sed(i,j)*fbc*biotic(n)%det_sediment(i,j) + biotic(n)%caco3_sed_remin(i,j) = muecaco3_sed(i,j)*fbc*biotic(n)%caco3_sediment(i,j) + +! remineralisation of sediments to supply nutrient fields. +! NB, btf values are positive from the water column into the sediment. mac, nov12. + T_prog(ind_no3)%btf(i,j) = -1.0 * rho0 * biotic(n)%det_sed_remin(i,j) + if (id_o2 .ne. 0) & + T_prog(ind_o2)%btf(i,j) = -172./16. * T_prog(ind_no3)%btf(i,j) + if (id_dic .ne. 0) & + T_prog(ind_dic)%btf(i,j) = 106./16. * T_prog(ind_no3)%btf(i,j) - & + rho0 * biotic(n)%caco3_sed_remin(i,j) + if (id_adic .ne. 0) & + T_prog(ind_adic)%btf(i,j) = T_prog(ind_dic)%btf(i,j) + if (id_fe .ne. 0) & + T_prog(ind_fe)%btf(i,j) = 2.0e-2 * T_prog(ind_no3)%btf(i,j) + if (id_alk .ne. 0) & + T_prog(ind_alk)%btf(i,j) = -2.0 * rho0 * biotic(n)%caco3_sed_remin(i,j) - & + T_prog(ind_no3)%btf(i,j) + + endif !} k > 0 + enddo !} i + enddo !} j + enddo !} n + + + + ! send rate of remineralisation of sediment tracers to output. mac, nov12. + do n = 1, instances !{ + if (id_caco3_sed_remin .gt. 0) then + used = send_data(id_caco3_sed_remin, biotic(n)%caco3_sed_remin(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) + endif + if (id_det_sed_remin .gt. 0) then + used = send_data(id_det_sed_remin, biotic(n)%det_sed_remin(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) + endif + enddo + +do n = 1, instances !{ + do nn=1,biotic(n)%ntr_bgc +! save the tracer bottom fluxes + if (biotic(n)%id_bgc_btf(nn) .gt. 0) then + used = send_data(biotic(n)%id_bgc_btf(nn), & + t_prog(biotic(n)%ind_bgc(nn))%btf(isc:iec,jsc:jec)/rho0, & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) + endif + enddo !} nn +enddo !} n + + +return +end subroutine csiro_bgc_bbc !} +! NAME="csiro_bgc_bbc" + +!####################################################################### +! +! +! +! Clean up various BIOTIC quantities for this run. +! +! + +subroutine csiro_bgc_end(isc, iec, jsc, jec, taup1, Thickness, T_prog, grid) !{ + +type(ocean_thickness_type), intent(in) :: Thickness +integer, intent(in) :: isc +integer, intent(in) :: iec +integer, intent(in) :: jsc +integer, intent(in) :: jec +integer, intent(in) :: taup1 +type(ocean_prog_tracer_type), dimension(:), intent(in) :: T_prog +type(ocean_grid_type), intent(in) :: Grid + +!----------------------------------------------------------------------- +! local definitions +!----------------------------------------------------------------------- + +character(len=64), parameter :: sub_name = 'csiro_bgc_end' +character(len=256), parameter :: error_header = & + '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: warn_header = & + '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: note_header = & + '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + +integer :: i +integer :: j +integer :: k +integer :: lun +integer :: n +real :: total_det +real :: total_zoo +real :: total_phy +real :: total_o2 +real :: total_no3 + +! ===================================================================== +! begin executable code +! ===================================================================== + +! integrate the total concentrations of some tracers +! for the end of the run + +total_no3 = 0.0 +total_phy = 0.0 +total_o2 = 0.0 +total_zoo = 0.0 +total_det = 0.0 + +do n = 1, instances !{ + do k = 1,grid%nk !{ + do j = jsc, jec !{ + do i = isc, iec !{ + total_no3 = total_no3 + & + t_prog(biotic(n)%ind_bgc(1))%field(i,j,k,taup1) * & + grid%dat(i,j) * grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,taup1) + total_phy = total_phy + & + t_prog(biotic(n)%ind_bgc(2))%field(i,j,k,taup1) * & + grid%dat(i,j) * grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,taup1) + total_o2 = total_o2 + & + t_prog(biotic(n)%ind_bgc(3))%field(i,j,k,taup1) * & + grid%dat(i,j) * grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,taup1) + total_zoo = total_zoo + & + t_prog(biotic(n)%ind_bgc(4))%field(i,j,k,taup1) * & + grid%dat(i,j) * grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,taup1) + total_det = total_det + & + t_prog(biotic(n)%ind_bgc(5))%field(i,j,k,taup1) * & + grid%dat(i,j) * grid%tmask(i,j,k) * Thickness%rho_dzt(i,j,k,taup1) + enddo !} i + enddo !} j + enddo !} k + + call mpp_sum(total_no3) + call mpp_sum(total_phy) + call mpp_sum(total_o2) + call mpp_sum(total_zoo) + call mpp_sum(total_det) + + write (stdout(),*) ' Instance ', trim(biotic(n)%name) + write (stdout(), & + '(/'' Total NO3 = '',es19.12,'' Gmol N'')') & + total_no3 * 1.0e-12 + write (stdout(), & + '(/'' Total PHY = '',es19.12,'' Gmol N'')') & + total_phy * 1.0e-12 + write (stdout(), & + '(/'' Total O2 = '',es19.12,'' Gmol O2'')') & + total_o2 * 1.0e-12 + write (stdout(), & + '(/'' Total ZOO = '',es19.12,'' Gmol N'')') & + total_zoo * 1.0e-12 + write (stdout(), & + '(/'' Total DET = '',es19.12,'' Gmol N'')') & + total_det * 1.0e-12 + write (stdout(), & + '(/'' Total N = '',es19.12,'' Gmol N'')') & + (total_no3+total_phy+total_zoo+total_det) * 1.0e-12 + +! +! Write out extra restart file(s) +! + + call reset_field_pointer(sed_restart, id_restart(1), biotic(n)%caco3_sediment(:,:)) + call reset_field_pointer(sed_restart, id_restart(2), biotic(n)%det_sediment(:,:)) + + call save_restart(sed_restart) + +enddo !} n + +return +end subroutine csiro_bgc_end !} +! NAME="csiro_bgc_end" + + +!####################################################################### +! +! +! +! Calculate the surface boundary conditions +! +! + +subroutine csiro_bgc_sbc(isc, iec, jsc, jec, isd, ied, jsd, jed, & + T_prog, aice, wnd, grid, time, use_waterflux, salt_restore_as_salt_flux, atm_co2, co2flux, sfc_co2) + +use ocmip2_co2calc_mod +use mpp_mod, only : mpp_sum +use time_interp_external_mod, only: time_interp_external +use time_manager_mod, only: days_in_year, days_in_month, & + get_date, set_date +use field_manager_mod, only: fm_get_index + +integer, intent(in) :: isc, iec +integer, intent(in) :: jsc, jec +integer, intent(in) :: isd, ied +integer, intent(in) :: jsd, jed +type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog +type(ocean_grid_type), intent(in) :: Grid +type(ocean_time_type), intent(in) :: Time +real, intent(in), dimension(isd:ied,jsd:jed) :: aice, wnd +logical, intent(in) :: use_waterflux, salt_restore_as_salt_flux + +real, intent(in), dimension(isd:ied,jsd:jed), optional :: atm_co2 +real, intent(out), dimension(isd:ied,jsd:jed), optional :: co2flux +real, intent(out), dimension(isd:ied,jsd:jed), optional :: sfc_co2 + +real :: total_co2_flux, total_aco2_flux +logical :: used + +!----------------------------------------------------------------------- +! local definitions +!----------------------------------------------------------------------- + +character(len=64), parameter :: sub_name = 'csiro_bgc_sbc' +character(len=256), parameter :: error_header = & + '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: warn_header = & + '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: note_header = & + '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + +! coefficients for O2 saturation +! + +real, parameter :: a_0 = 2.00907 +real, parameter :: a_1 = 3.22014 +real, parameter :: a_2 = 4.05010 +real, parameter :: a_3 = 4.94457 +real, parameter :: a_4 = -2.56847e-01 +real, parameter :: a_5 = 3.88767 +real, parameter :: b_0 = -6.24523e-03 +real, parameter :: b_1 = -7.37614e-03 +real, parameter :: b_2 = -1.03410e-02 +real, parameter :: b_3 = -8.17083e-03 +real, parameter :: c_0 = -4.88682e-07 + +integer :: i +integer :: j +integer :: k +integer :: n +integer :: nn +integer :: ntr_bgc +integer :: ind_trc +integer :: m +integer :: kz +integer :: hour +integer :: day +real :: days_in_this_year +integer :: minute +integer :: month +integer :: num_days +integer :: second +integer :: year +integer :: indtemp, indsal + +! ===================================================================== +! begin executable code +! ===================================================================== + + call get_date(time%model_time, & + year, month, day, hour, minute, second) + num_days = days_in_year(time%model_time) + days_in_this_year = 0.0 + do m = 1, month - 1 + days_in_this_year = days_in_this_year + & + days_in_month(set_date(year, m, 1)) + enddo + days_in_this_year = days_in_this_year + day - 1 + hour/24.0 + & + minute/1440.0 + second/86400.0 + +!--------------------------------------------------------------------- +! calculate interpolated xkw, seaice fraction and atmospheric +! pressure & solar radiation +!--------------------------------------------------------------------- +if (gasx_from_file) then + call time_interp_external(pistonveloc_id, time%model_time, xkw_t) +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 3.6e5 is a conversion of units; cm/hr to m/s. +! mac, may13. + xkw_t(i,j)=(0.31 * wnd(i,j)**2.0 ) /3.6e5 + enddo ! i + enddo ! j +endif ! if (gasx_from_file) + +call time_interp_external(nat_co2_id, time%model_time, nat_co2) +call time_interp_external(atmpress_id, time%model_time, patm_t) +call time_interp_external(dust_id, time%model_time, dust_t) +if (id_adic .ne. 0) then +! The atmospheric co2 value for the anthropogenic+natural carbon tracer +! is either read from a file or a value from the access atmospheric model, +! as determined by the flag use_access_co2. mac, may13. + if (use_access_co2) then + aco2(isc:iec,jsc:jec) = atm_co2(isc:iec,jsc:jec) + else + call time_interp_external(aco2_id, time%model_time, aco2) + endif +endif +if (ice_file4gasx) then + call time_interp_external(seaicefract_id, time%model_time, fice_t) +else + fice_t(isc:iec,jsc:jec) = aice(isc:iec,jsc:jec) +endif + +indtemp = fm_get_index('/ocean_mod/prog_tracers/temp') +indsal = fm_get_index('/ocean_mod/prog_tracers/salt') + +! ------------------------------------------------------------------- +! rjm: start by zeroing the fluxes +if ((.not. use_waterflux) .or. (salt_restore_as_salt_flux)) then !{ + do n = 1, instances !{ + ntr_bgc = biotic(n)%ntr_bgc + do nn=1,ntr_bgc + ind_trc = biotic(n)%ind_bgc(nn) + t_prog(ind_trc)%stf(:,:) = 0 + enddo + enddo +endif + +!--------------------------------------------------------------------- +! Compute the Schmidt number of CO2 in seawater using the +! formulation presented by Wanninkhof (1992, J. Geophys. Res., 97, +! 7373-7382). +!--------------------------------------------------------------------- +! +! no point doing the following calculations if no dic and o2 tracer +if (id_dic .eq. 0 .and. id_o2.eq. 0) return + +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) * & + (-0.043219))) * grid%tmask(i,j,1) + enddo !} i +enddo !} j + +!--------------------------------------------------------------------- +! Compute the Schmidt number of O2 in seawater using the +! formulation proposed by Keeling et al. (1998, Global Biogeochem. +! Cycles, 12, 141-163). +!--------------------------------------------------------------------- + +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) * & + (-0.008004))) * grid%tmask(i,j,1) + enddo !} i +enddo !} j + + +!--------------------------------------------------------------------- +! calculate csurf, csat and csat - csurf via the routine co2calc +! input and output units are in mol/m^3 +!--------------------------------------------------------------------- + +! determine the atmospheric pCO2 concetration for this time-step +do n = 1, instances !{ + do j = jsc, jec + do i = isc, iec + biotic(n)%pco2atm(i,j) = nat_co2(i,j) + enddo + enddo +enddo !} n + +if (id_adic .ne. 0) then + do n = 1, instances !{ + do j = jsc, jec + do i = isc, iec + biotic(n)%paco2atm(i,j) = aco2(i,j) + enddo + enddo + enddo !} n +endif + +!call init_ocmip2_co2calc( & +! time%model_time, isc, iec, jsc, jec, 1, & +! t_prog(indtemp)%field(isc:iec,jsc:jec,1,time%taum1), & +! t_prog(indsal)%field(isc:iec,jsc:jec,1,time%taum1)) + +do n = 1, instances !{ + do j = jsc, jec !{ + do i = isc, iec !{ + htotallo(i,j) = .1 + htotalhi(i,j) = 100.0 + enddo !} i + enddo !} j + +! ocmip2+co2cal expects tracers in mol/m3 !!!! + if (id_dic.ne.0) then + biotic(n)%po4(:,:) = t_prog(ind_dic)%field(isd:ied,jsd:jed,1,time%taum1)*0 + if (id_no3.ne.0) biotic(n)%po4(:,:) = t_prog(ind_no3)%field(isd:ied,jsd:jed,1,time%taum1)/16.*1e-3 + if (id_po4.ne.0) biotic(n)%po4(:,:) = t_prog(ind_po4)%field(isd:ied,jsd:jed,1,time%taum1)*1e-3 + + call ocmip2_co2calc(isd, jsd, & + isc, iec, jsc, jec, & + grid%tmask(isd:ied,jsd:jed,1), & + t_prog(indtemp)%field(isd:ied,jsd:jed,1,time%taum1), & + t_prog(indsal)%field(isd:ied,jsd:jed,1,time%taum1), & + t_prog(ind_dic)%field(isd:ied,jsd:jed,1,time%taum1)*1e-3,& + t_prog(ind_alk)%field(isd:ied,jsd:jed,1,time%taum1)*1e-3,& + 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)%csurf(isc:iec,jsc:jec), & + alpha=biotic(n)%alpha(isc:iec,jsc:jec) , & + pco2surf = biotic(n)%pco2surf(isc:iec,jsc:jec), & + scale= 1.0/1024.5 ) + + if (id_adic .ne. 0) then ! calculate CO2 flux including anthropogenic CO2 + call ocmip2_co2calc(isd, jsd, & + isc, iec, jsc, jec, & + grid%tmask(isd:ied,jsd:jed,1), & + t_prog(indtemp)%field(isd:ied,jsd:jed,1,time%taum1), & + t_prog(indsal)%field(isd:ied,jsd:jed,1,time%taum1), & + t_prog(ind_adic)%field(isd:ied,jsd:jed,1,time%taum1)*1e-3,& + t_prog(ind_alk)%field(isd:ied,jsd:jed,1,time%taum1)*1e-3,& + 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)%acsurf(isc:iec,jsc:jec), & + alpha=biotic(n)%alpha(isc:iec,jsc:jec) , & + pco2surf = biotic(n)%paco2surf(isc:iec,jsc:jec), & + scale = 1.0/1024.5 ) + endif ! no adic + endif ! no dic +enddo !} n + + +!--------------------------------------------------------------------- +! Compute the oxygen saturation concentration at 1 atm total +! pressure in mol/m^3 given the temperature (t, in deg C) and +! the salinity (s, in permil) +! +! From Garcia and Gordon (1992), Limnology and Oceonography. +! The formula used is from page 1310, eq (8). +! +! *** Note: the "a3*ts^2" term (in the paper) is incorrect. *** +! *** It shouldn't be there. *** +! +! o2_saturation is defiend between T(freezing) <= T <= 40 deg C and +! 0 permil <= S <= 42 permil +! +! check value: T = 10 deg C, S = 35 permil, +! o2_saturation = 0.282015 mol/m^3 +!--------------------------------------------------------------------- +! + +do j = jsc, jec !{ + do i = isc, iec !{ + tt(i) = 298.15 - t_prog(indtemp)%field(i,j,1,time%taum1) + tk(i) = 273.15 + t_prog(indtemp)%field(i,j,1,time%taum1) + ts(i) = log(tt(i) / tk(i)) + ts2(i) = ts(i) * ts(i) + ts3(i) = ts2(i) * ts(i) + ts4(i) = ts3(i) * ts(i) + ts5(i) = ts4(i) * ts(i) + o2_saturation(i,j) = & + exp(a_0 + a_1*ts(i) + a_2*ts2(i) + & + a_3*ts3(i) + a_4*ts4(i) + a_5*ts5(i) + & + t_prog(indsal)%field(i,j,1,time%taum1) * & + (b_0 + b_1*ts(i) + b_2*ts2(i) + b_3*ts3(i) + & + c_0*t_prog(indsal)%field(i,j,1,time%taum1))) + enddo !} i +enddo !} j + +!chd +!chd convert from ml/l to mmol/m^3 +!chd + +do j = jsc, jec !{ + do i = isc, iec !{ + o2_saturation(i,j) = o2_saturation(i,j) * & + (1000.*1000.0/22391.6) + enddo !} i +enddo !} j + +! +!--------------------------------------------------------------------- +! calculate piston-velocities +! including effect of sea-ice +! xkw is given in cm/s (converted in read_biotic_bc), therefore +! kw is also in cm/s +!--------------------------------------------------------------------- + +do j = jsc, jec !{ + do i = isc, iec !{ + kw_co2(i,j) = (1.0 - fice_t(i,j)) * xkw_t(i,j) * & + sqrt(660.0/sc_co2(i,j)) * grid%tmask(i,j,1) + + kw_o2(i,j) = (1.0 - fice_t(i,j)) * xkw_t(i,j) * & + sqrt(660.0/sc_o2(i,j)) * grid%tmask(i,j,1) + enddo !} i +enddo !} j + +! +!chd--------------------------------------------------------------------- +!chd calculate surface fluxes for BIOTICs +!chd kw is in m/s, o2_sat is in mmol/m^3, +!chd hence stf should be in mmol/m^2/s +!chd--------------------------------------------------------------------- +! + +total_co2_flux = 0.0 + +if (id_dic.ne.0) then + do n = 1, instances !{ + do j = jsc, jec !{ + do i = isc, iec !{ + ! These calculations used to be done in ocmip2_co2calc in mom4p1_2007, + ! but not with mom4p1_2009. + ! Some extra calculations now required in csiro_bgc_sbc for co2 air-sea + ! alpha is in units of mol/m3/atm (atm -> partial pressure of CO2), and 1e6 is to convert micro-atm to atm. + ! flux. mac, apr11. + biotic(n)%csat(i,j) = biotic(n)%pco2atm(i,j) / 1e6 * biotic(n)%alpha(i,j) * patm_t(i,j) + biotic(n)%csat_csurf(i,j) = biotic(n)%csat(i,j) - biotic(n)%csurf(i,j) + + t_prog(ind_dic)%stf(i,j) = rho0 * kw_co2(i,j) * & + biotic(n)%csat_csurf(i,j)*1e3 !convert from mol/m^3 to mmol/m^3 + + total_co2_flux = total_co2_flux + kw_co2(i,j) * & + biotic(n)%csat_csurf(i,j) * 3.7843e-7 * grid%dat(i,j) * grid%tmask(i,j,1) ! convert from mol/s to Pg/year (12.0*1e-15*86400*365=3.78e-7) + enddo !} i + enddo !} j + enddo !} n +endif + +if (id_total_co2_flux .gt. 0) then + call mpp_sum(total_co2_flux) + used = send_data(id_total_co2_flux,total_co2_flux,Time%model_time) +endif + + +total_aco2_flux = 0.0 + +if (id_adic.ne.0) then + do n = 1, instances !{ + do j = jsc, jec !{ + do i = isc, iec !{ + ! These calculations used to be done in ocmip2_co2calc in mom4p1_2007, + ! but not with mom4p1_2009. + ! Some extra calculations now required in csiro_bgc_sbc for co2 air-sea + ! flux. mac, apr11. + ! Note, csat has been used as temporary variable, alpha is only f(T,S), + ! csurf had been used as a temporary variable, but now needs to be saved + ! for preindustrial vs anthropogenic. + biotic(n)%csat(i,j) = biotic(n)%paco2atm(i,j) / 1e6 * biotic(n)%alpha(i,j) * patm_t(i,j) + biotic(n)%csat_acsurf(i,j) = biotic(n)%csat(i,j) - biotic(n)%acsurf(i,j) + + t_prog(ind_adic)%stf(i,j) = rho0 * kw_co2(i,j) * & + biotic(n)%csat_acsurf(i,j)*1e3 !convert from mol/m^3 to mmol/m^3 + + total_aco2_flux = total_aco2_flux + kw_co2(i,j) * & + biotic(n)%csat_acsurf(i,j) * 3.7843e-7 * grid%dat(i,j) * grid%tmask(i,j,1) ! convert from mol/s to Pg/year (12.0*1e-15*86400*365=3.78e-7) +! send the anthropogenic Pco2 and co2 flux into the ocean back to the atmospheric model. mac, may13. +!RASF avoid using Ocean_sfc. +! if(present(Ocean_sfc)) then +! Ocean_sfc%co2flux(i,j) = kw_co2(i,j) * & +! biotic(n)%csat_acsurf(i,j)*0.04401 !convert from mol/m^2/s to kg(CO2)/m^2/s, 0.04401 kg(CO2)/mole +! Ocean_sfc%co2(i,j) = biotic(n)%paco2surf(i,j) +! endif + enddo !} i + enddo !} j + if(present(co2flux) .and. present(sfc_co2)) then + do j = jsc, jec !{ + do i = isc, iec !{ + co2flux(i,j) = kw_co2(i,j) * & + biotic(n)%csat_acsurf(i,j)*0.04401 !convert from mol/m^2/s to kg(CO2)/m^2/s, 0.04401 kg(CO2)/mole + sfc_co2(i,j) = biotic(n)%paco2surf(i,j) + enddo !} i + enddo !} j + endif + + enddo !} n +endif + +if (id_total_aco2_flux .gt. 0) then + call mpp_sum(total_aco2_flux) + used = send_data(id_total_aco2_flux,total_aco2_flux,Time%model_time) +endif + + +if (id_o2.ne.0) then + do n = 1, instances !{ + do j = jsc, jec !{ + do i = isc, iec !{ + t_prog(ind_o2)%stf(i,j) = rho0 * kw_o2(i,j) * & + (o2_saturation(i,j) * patm_t(i,j) - & + t_prog(ind_o2)%field(i,j,1,time%taum1)) + enddo !} i + enddo !} j + enddo !} n +endif + +!rjm: Do atmospheric iron input +! for Fe units of nmol/m3, Fe input must be in nmol/m2/s +if (id_fe.ne.0) then + do n = 1, instances !{ + do j = jsc, jec !{ + do i = isc, iec !{ + t_prog(ind_fe)%stf(i,j) = rho0 * dust_t(i,j) + enddo !} i + enddo !} j + enddo !} n +endif + +if (.not. use_waterflux) then !{ +! rjm - One only needs to compute virtual fluxes if waterflux is not used +! NB, when this routine is called, t_prog(ind_sal)%stf() only has applied fluxes +! and not any restoring fluxes. +! csiro_bgc_virtual_fluxes() now calculated BGC virtual fluxes if any virtual flux +! is use to restore salinity. mac, dec12. + do n = 1, instances !{ + ntr_bgc = biotic(n)%ntr_bgc + + do nn=1,ntr_bgc + ind_trc = biotic(n)%ind_bgc(nn) + + do j = jsc, jec !{ + do i = isc, iec !{ + biotic(n)%vstf(i,j,nn) = & + t_prog(indsal)%stf(i,j) * & + biotic(n)%bgc_global(nn) / biotic(n)%sal_global + t_prog(ind_trc)%stf(i,j) = & + t_prog(ind_trc)%stf(i,j) + & + biotic(n)%vstf(i,j,nn) + enddo !} i + enddo !} j + enddo !} nn + enddo !} n +endif !} + +return + +end subroutine csiro_bgc_sbc !} +! NAME="csiro_bgc_sbc" + +!####################################################################### +! +! +! +! Apply virtual fluxes to BGC tracers when sea surface salinity +! is restored as a salt flux within ocean_sbc//flux_adjust() +! This subroutine only adds the virtual flux associated with +! any virtual salt fluxes from salinity restoring; any virtual fluxes +! from applied forcing (i.e. when use_waterflux = .false.) are applied to +! bgc virtual fluxes with the csiro_bgc_sbc() subroutine. +! mac, dec12. +! + +subroutine csiro_bgc_virtual_fluxes(isc, iec, jsc, jec, isd, ied, jsd, jed, salt_vstf, T_prog) !{ + +!----------------------------------------------------------------------- +! arguments +!----------------------------------------------------------------------- + +integer, intent(in) :: isc, iec, jsc, jec, isd, ied, jsd, jed +real, intent(in), dimension(isd:ied,jsd:jed) :: salt_vstf ! virtual salt flux. +type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog + + +!----------------------------------------------------------------------- +! local definitions +!----------------------------------------------------------------------- + +integer :: i, j, n, nn, ntr_bgc, ind_trc + +!----------------------------------------------------------------------- +! start executable code +!----------------------------------------------------------------------- + + do n = 1, instances !{ + ntr_bgc = biotic(n)%ntr_bgc + + do nn=1,ntr_bgc + ind_trc = biotic(n)%ind_bgc(nn) + + do j = jsc, jec !{ + do i = isc, iec !{ + biotic(n)%vstf(i,j,nn) = & + salt_vstf(i,j) * & + biotic(n)%bgc_global(nn) / biotic(n)%sal_global + t_prog(ind_trc)%stf(i,j) = & + t_prog(ind_trc)%stf(i,j) + & + biotic(n)%vstf(i,j,nn) + enddo !} i + enddo !} j + enddo !} nn + enddo !} n + + +end subroutine csiro_bgc_virtual_fluxes !} +! NAME="csiro_bgc_virtual_fluxes" + + +!####################################################################### +! +! +! +! Set up any extra fields needed by the tracer packages +! +! Save pointers to various "types", such as Grid and Domains. +! + +subroutine csiro_bgc_init !{ + +!----------------------------------------------------------------------- +! local definitions +!----------------------------------------------------------------------- + +character(len=64), parameter :: sub_name = 'csiro_bgc_init' +character(len=256), parameter :: error_header = & + '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: warn_header = & + '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: note_header = & + '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=128) :: version = "git:/mom_mac/src/mom5csiro_bgc.F90" +character(len=128) :: tagname = "mac, sep16." + +integer :: n +character(len=fm_field_name_len) :: name +character(len=fm_path_name_len) :: path_to_names +character(len=fm_field_name_len+1) :: suffix +character(len=fm_string_len) :: string +character(len=fm_field_name_len+3) :: long_suffix +character(len=256) :: caller_str +character(len=fm_string_len), pointer, dimension(:) :: good_list + +integer :: nn, ntr_bgc +real :: min_range=0.0, max_range=1.e4 +real :: sum_ntr = 0.0 +character(len=8) :: bgc_trc='tracer00' + +! Initialize the csiro_bgc package + +!package_index = otpm_set_tracer_package(package_name, & +! file_in = default_file_in, file_out = default_file_out, & +! caller = trim(mod_name) // '(' // trim(sub_name) // ')') +package_index = otpm_set_tracer_package(package_name, & + restart_file = default_file_in, & + caller = trim(mod_name) // '(' // trim(sub_name) // ')') + +! Check whether to use this package + +path_to_names = '/ocean_mod/tracer_packages/' // trim(package_name) // '/names' +instances = fm_get_length(path_to_names) +if (instances .lt. 0) then !{ + call mpp_error(FATAL, trim(error_header) // ' Could not get number of instances') +endif !} + +! Check some things + +if (instances .eq. 0) then !{ + write (stdout(),*) trim(note_header), ' No instances' + do_csiro_bgc = .false. +else !}{ + if (instances .eq. 1) then !{ + write (stdout(),*) trim(note_header), ' ', instances, ' instance' + else !}{ + write (stdout(),*) trim(note_header), ' ', instances, ' instances' + call mpp_error(FATAL, trim(error_header) // ' CSIRO_bgc code not yet ready for instances > 1, mac, nov12') + endif !} + do_csiro_bgc = .true. +endif !} + +! Return if we don't want to use this package, +! after changing the list back + +if (.not. do_csiro_bgc) then !{ + return +endif !} + +! allocate storage for biotic array +allocate ( biotic(instances) ) + +! loop over the names, saving them into the biotic array + +do n = 1, instances !{ + if (fm_get_value(path_to_names, name, index = n)) then !{ + biotic(n)%name = name + else !}{ + write (name,*) n + call mpp_error(FATAL, trim(error_header) // & + 'Bad field name for index ' // trim(name)) + endif !} +enddo !} + + +!----------------------------------------------------------------------- +! Process the namelists +!----------------------------------------------------------------------- + +! Add the package name to the list of good namelists, to be used +! later for a consistency check + +if (fm_new_value('/ocean_mod/GOOD/good_namelists', package_name, append = .true.) .le. 0) then !{ + call mpp_error(FATAL, trim(error_header) // & + ' Could not add ' // trim(package_name) // ' to "good_namelists" list') +endif !} + +!----------------------------------------------------------------------- +! Set up the *global* namelist +!----------------------------------------------------------------------- + +caller_str = trim(mod_name) // '(' // trim(sub_name) // ')' + +call fm_util_start_namelist(package_name, '*global*', caller = caller_str, no_overwrite = .true., & + check = .true.) + + call fm_util_set_value('atmpress_file', 'INPUT/atmpress_ocmip2.nc') + call fm_util_set_value('atmpress_name', 'atmpress') + call fm_util_set_value('pistonveloc_file', 'INPUT/pistonveloc_ocmip2.nc') + call fm_util_set_value('pistonveloc_name', 'pistonveloc') + call fm_util_set_value('aco2_file', 'INPUT/co2_b35_2D.nc') + call fm_util_set_value('aco2_name', 'co2') + call fm_util_set_value('seaicefract_file', 'INPUT/f_ice_ocmip2.nc') + call fm_util_set_value('seaicefract_name', 'f_ice') + call fm_util_set_value('dust_file', 'INPUT/dust.nc') + call fm_util_set_value('dust_name', 'DUST') + +! additional information + call fm_util_set_value('s_npp', -.1) ! scale factor for NP + 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('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? +! Set defaults to use the ACCESS atmospheric CO2 for the anthropogenic CO2 in the ocean if this is compiled for ACCESS-CM/ESM, otherwise to read CO2 from files provided. mac, may13. +#if defined(ACCESS_CM) + call fm_util_set_value('use_access_co2', .true.)! use access model co2 values. mac, may13. +#else + call fm_util_set_value('use_access_co2', .false.)! use file for atmospheric co2 values. +#endif + call fm_util_set_value('id_po4',0) + call fm_util_set_value('id_dic',0) + call fm_util_set_value('id_adic',0) + call fm_util_set_value('id_alk',0) + call fm_util_set_value('id_o2',0) + call fm_util_set_value('id_no3',0) + call fm_util_set_value('id_phy',0) + call fm_util_set_value('id_zoo',0) + call fm_util_set_value('id_det',0) + call fm_util_set_value('id_caco3',0) + call fm_util_set_value('id_fe',0) + + atmpress_file = fm_util_get_string ('atmpress_file', scalar = .true.) + atmpress_name = fm_util_get_string ('atmpress_name', scalar = .true.) + pistonveloc_file = fm_util_get_string ('pistonveloc_file', scalar = .true.) + pistonveloc_name = fm_util_get_string ('pistonveloc_name', scalar = .true.) + aco2_file = fm_util_get_string ('aco2_file', scalar = .true.) + aco2_name = fm_util_get_string ('aco2_name', scalar = .true.) + seaicefract_file = fm_util_get_string ('seaicefract_file', scalar = .true.) + seaicefract_name = fm_util_get_string ('seaicefract_name', scalar = .true.) + dust_file = fm_util_get_string ('dust_file', scalar = .true.) + dust_name = fm_util_get_string ('dust_name', scalar = .true.) + + qbio_model = fm_util_get_string ('qbio_model', scalar = .true.) + s_npp = fm_util_get_real ('s_npp', scalar = .true.) + bio_version = fm_util_get_integer ('bio_version', scalar = .true.) + ! rjm: Tracer to use + zero_floor = fm_util_get_logical ('zero_floor', 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.) + use_access_co2 = fm_util_get_logical ('use_access_co2', scalar = .true.) + + id_dic = fm_util_get_integer ('id_dic', scalar = .true.) + id_adic = fm_util_get_integer ('id_adic', scalar = .true.) + id_alk = fm_util_get_integer ('id_alk', scalar = .true.) + id_po4 = fm_util_get_integer ('id_po4', scalar = .true.) + id_o2 = fm_util_get_integer ('id_o2', scalar = .true.) + id_no3 = fm_util_get_integer ('id_no3', scalar = .true.) + id_zoo = fm_util_get_integer ('id_zoo', scalar = .true.) + id_phy = fm_util_get_integer ('id_phy', scalar = .true.) + id_det = fm_util_get_integer ('id_det', scalar = .true.) + id_caco3 = fm_util_get_integer ('id_caco3', scalar = .true.) + id_fe = fm_util_get_integer ('id_fe', scalar = .true.) + + +call fm_util_end_namelist(package_name, '*global*', caller = caller_str, check = .true.) + + +sum_ntr = min(1,id_dic)+min(1,id_adic)+min(1,id_po4)+min(1,id_alk) & + +min(1,id_o2)+min(1,id_no3)+min(1,id_phy)+min(1,id_zoo)+min(1,id_det) & + +min(1,id_caco3)+min(1,id_fe) +if (mpp_pe() == mpp_root_pe() ) print*,'csiro_bgc_init: Number bgc tracers = ',sum_ntr + + +!----------------------------------------------------------------------- +! Set up the prognostic tracers. +!----------------------------------------------------------------------- + +! RASF: Use sensible names for tracers +! For default case names will be simply adic, o2 etc +! If we have multiple instances the names will take the form +! adic_instancename, o2_instancename +! Different input files etc can be set in the field_table like for temp, salt + +do n = 1, instances !{ + + biotic(n)%ntr_bgc = min(1,id_dic)+min(1,id_adic)+min(1,id_po4)+min(1,id_alk) & + +min(1,id_o2)+min(1,id_no3)+min(1,id_phy)+min(1,id_zoo)+min(1,id_det) & + +min(1,id_caco3)+min(1,id_fe) + if (mpp_pe() == mpp_root_pe() ) print*,'Number bgc tracers = ',biotic(n)%ntr_bgc + + + name = biotic(n)%name + if (name(1:1) .eq. '_') then !{ + suffix = '' + long_suffix = '' + else !}{ + suffix = '_' // name + long_suffix = ' (' // trim(name) // ')' + endif !} + + ntr_bgc=biotic(n)%ntr_bgc + + do nn=1,ntr_bgc + if ( nn == id_no3 ) then + bgc_trc='no3' + min_range=-1e-5 + max_range=100.0 + else if ( nn == id_phy ) then + bgc_trc='phy' + min_range=-1e-7 + max_range=10.0 + else if (nn == id_o2 ) then + bgc_trc='o2' + min_range=-1e-7 + max_range=10.0 + else if (nn == id_zoo ) then + bgc_trc='zoo' + min_range=-1e-7 + max_range=10.0 + else if (nn == id_det ) then + bgc_trc='det' + min_range=0.0 + max_range=10.0 + else if (nn == id_caco3 ) then + bgc_trc='caco3' + min_range=-1e-5 + max_range=10.0 + else if (nn == id_dic ) then + bgc_trc='dic' + min_range=0.0 + max_range=3000.0 + else if (nn == id_alk ) then + bgc_trc='alk' + min_range=0.0 + max_range=3000.0 + else if (nn == id_adic ) then + bgc_trc='adic' + min_range=-1e-5 + max_range=3000.0 + else if (nn == id_fe ) then + bgc_trc='fe' + min_range=-0.0001 + max_range=10.0 + else if (nn == id_po4 ) then + bgc_trc='po4' + min_range=0.0 + max_range=100.0 + else + bgc_trc(1:6)='dummy' + write(bgc_trc(7:8),'(i2.2)') nn + min_range=0.0 + max_range=100.0 + endif + + + biotic(n)%ind_bgc(nn) = otpm_set_prog_tracer(trim(bgc_trc) // trim(suffix), & + package_name, & + longname = trim(bgc_trc) // trim(long_suffix), & + units = 'mmol/m^3', flux_units = 'mmol/m^2/s', & + caller = trim(mod_name)//'('//trim(sub_name)//')', & + min_range=min_range,max_range=max_range) + enddo !} nn +enddo !} n + +!----------------------------------------------------------------------- +! Set up the instance namelists +!----------------------------------------------------------------------- + + +do n = 1, instances !{ + +! create the instance namelist + + call fm_util_start_namelist(package_name, trim(biotic(n)%name) , caller = caller_str, no_overwrite = .true., & + check = .true.) + + + call fm_util_set_value('file_in', default_file_in) + call fm_util_set_value('file_out', default_file_out) + call fm_util_set_value('init', .false.) + + call fm_util_end_namelist(package_name, biotic(n)%name, check = .true., caller = caller_str) + +enddo !} n + +! Check for any errors in the number of fields in the namelists for this package + +good_list => fm_util_get_string_array('/ocean_mod/GOOD/namelists/' // trim(package_name) // '/good_values', & + caller = trim(mod_name) // '(' // trim(sub_name) // ')') +if (associated(good_list)) then !{ + call fm_util_check_for_bad_fields('/ocean_mod/namelists/' // trim(package_name), good_list, & + caller = trim(mod_name) // '(' // trim(sub_name) // ')') + deallocate(good_list) +else !}{ + call mpp_error(FATAL,trim(error_header) // ' Empty "' // trim(package_name) // '" list') +endif !} + + +! possible tracers in model +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) + 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) + + ind_zoo = biotic(n)%ind_bgc(id_zoo) + ind_phy = biotic(n)%ind_bgc(id_phy) + ind_det = biotic(n)%ind_bgc(id_det) + + ind_caco3= biotic(n)%ind_bgc(id_caco3) + if(id_fe .ne. 0 ) ind_fe= biotic(n)%ind_bgc(id_fe) + +enddo !} n + +! Write SVN version numbers to both the logfile and the standard output +! mac, mar12. +call write_version_number(version, tagname) + +if (mpp_pe() == mpp_root_pe() ) then + print* + print*,"Repo. details of csiro_bgc compiled into executable" + print*,"git:/mom_mac/src/mom5csiro_bgc.F90" + print*,"mac, sep16." + print* +endif + +id_clock_csiro_obgc = mpp_clock_id('(CSIRO ocean biogeochemistry)',grain=CLOCK_ROUTINE) + +return + +end subroutine csiro_bgc_init !} +! NAME="csiro_bgc_init" + + +!####################################################################### +! +! +! +! compute the source terms for the BIOTICs, including boundary +! conditions (not done in setvbc, to minimize number +! of hooks required in MOM base code) +! +! + +subroutine csiro_bgc_source(isc, iec, jsc, jec, isd, ied, jsd, jed, T_prog, grid, time, dtts, Thickness, Dens, swflx, sw_frac_zt) !{ + +use field_manager_mod, only: fm_get_index + +!----------------------------------------------------------------------- +! arguments +!----------------------------------------------------------------------- +! + +integer, intent(in) :: isc, iec +integer, intent(in) :: jsc, jec +integer, intent(in) :: isd, ied +integer, intent(in) :: jsd, jed +type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog +type(ocean_grid_type), intent(in) :: Grid +type(ocean_time_type), intent(in) :: Time +real, intent(in) :: dtts +type(ocean_thickness_type), intent(in) :: Thickness +type(ocean_density_type), intent(in) :: Dens +real, intent(in), dimension(isd:ied,jsd:jed) :: swflx ! short wave radiation flux (W/m^2) +real, intent(in), dimension(isd:,jsd:,:) :: sw_frac_zt ! fraction of short wave radiation flux (none) + + +!----------------------------------------------------------------------- +! local definitions +!----------------------------------------------------------------------- + +character(len=64), parameter :: sub_name = 'csiro_bgc_source' +character(len=256), parameter :: error_header = & + '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: warn_header = & + '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: note_header = & + '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + +integer :: i +integer :: j +integer :: k +integer :: n +logical :: used +integer :: day +integer :: month +integer :: year +integer :: hour +integer :: minute +integer :: second +integer :: indtemp + +integer :: nn, ntr_bgc + +! include "bio_v1.par" + +real :: total_trc +! +! ===================================================================== +! begin executable code +! ===================================================================== + +! get the model month + +call mpp_clock_begin(id_clock_csiro_obgc) + +call get_date(time%model_time, year, month, day, & + hour, minute, second) + +!----------------------------------------------------------------------- +! calculate the source terms for BIOTICs +!----------------------------------------------------------------------- + + + select case(bio_version) +case(0) +! call bio_v0a(Thickness) +!include "bio_v0.f90" + +case(1) +! call bio_v1a(Thickness) +!include "bio_v1.f90" + +case(2) +!include "bio_v2.f90" + +case(3) +call bio_v3(isc, iec, jsc, jec, isd, ied, jsd, jed, T_prog, grid, time, dtts, Thickness, Dens, swflx, sw_frac_zt) + +case default +if (mpp_pe() == mpp_root_pe() )print*,'rjm do no bio version' + +end select + +!----------------------------------------------------------------------- +! Save variables for diagnostics +!----------------------------------------------------------------------- +if (id_pco2 .gt. 0) then + used = send_data(id_pco2, biotic(1)%pco2surf(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) +endif +if (id_paco2 .gt. 0) then + 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_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)) +endif +if (id_aco2_sat .gt. 0) then + used = send_data(id_aco2_sat, biotic(1)%csat_acsurf(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) +endif + +if (id_sc_o2 .gt. 0) then + used = send_data(id_sc_o2, sc_o2(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) +endif +if (id_o2_sat .gt. 0) then + used = send_data(id_o2_sat, o2_saturation(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) +endif +if (id_kw_o2 .gt. 0) then + used = send_data(id_kw_o2, kw_o2(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) +endif + +! Gross production of phytoplankton + +if (id_pprod_gross .gt. 0) then + used = send_data(id_pprod_gross, pprod_gross(isc:iec,jsc:jec,:), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,:)) +endif + +if (id_pprod_gross_2d .gt. 0) then + pprod_gross_2d(:,:)=0.0 + do k=1,grid%nk + do j=jsc,jec + do i=isc,iec + pprod_gross_2d(i,j)=pprod_gross_2d(i,j) + pprod_gross(i,j,k)*Thickness%dzt(i,j,k) + enddo + enddo + enddo + used = send_data(id_pprod_gross_2d, pprod_gross_2d(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 + used = send_data(id_zprod_gross, zprod_gross(isc:iec,jsc:jec,:), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,:)) +endif + +! growth limit of phytoplankton for light. + +if (id_light_limit .gt. 0) then + used = send_data(id_light_limit, light_limit(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) +endif + +do n = 1, instances !{ + + ntr_bgc=biotic(n)%ntr_bgc + do nn=1,ntr_bgc + +! save the tracer surface fluxes + if (biotic(n)%id_bgc_stf(nn) .gt. 0) then + used = send_data(biotic(n)%id_bgc_stf(nn), & + t_prog(biotic(n)%ind_bgc(nn))%stf(isc:iec,jsc:jec)/rho0, & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) + endif +! save the tracer sources + if (biotic(n)%id_bgc_src(nn) .gt. 0) then + used = send_data(biotic(n)%id_bgc_src(nn), & + t_prog(biotic(n)%ind_bgc(nn))%th_tendency(isc:iec,jsc:jec,:)/ & + Thickness%rho_dzt(isc:iec,jsc:jec,:,time%tau), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,:)) + endif + enddo !} nn + + ! rate of deposition of sinking tracers to sediment + if (id_caco3_sed_depst .gt. 0) then + used = send_data(id_caco3_sed_depst, biotic(n)%caco3_sed_depst(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) + endif + if (id_det_sed_depst .gt. 0) then + used = send_data(id_det_sed_depst, biotic(n)%det_sed_depst(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) + endif + + +enddo !} n + + +call mpp_clock_end(id_clock_csiro_obgc) + +return +end subroutine csiro_bgc_source !} +! NAME="csiro_bgc_source" + + +!####################################################################### +! +! +! +! Initialize variables, read in namelists, calculate constants +! for a given run and allocate diagnostic arrays +! +! + +subroutine csiro_bgc_start (time, domain, grid) !{ + +use time_manager_mod, only: days_in_year, days_in_month, & + get_date, set_date +use time_interp_external_mod, only: init_external_field +use diag_manager_mod, only: register_diag_field, diag_axis_init +use fms_mod, only : read_data + +! +!----------------------------------------------------------------------- +! arguments +!----------------------------------------------------------------------- +! + +type(ocean_time_type), intent(in) :: Time +type(ocean_domain_type), intent(in) :: Domain +type(ocean_grid_type), intent(in) :: Grid + +! +!----------------------------------------------------------------------- +! local definitions +!----------------------------------------------------------------------- + +character(len=64), parameter :: sub_name = 'csiro_bgc_start' +character(len=256), parameter :: error_header = & + '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: warn_header = & + '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: note_header = & + '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + +character(len=fm_string_len) :: file_in +character(len=fm_string_len) :: file_out +integer :: index +logical :: init +integer :: check +integer :: m, n +integer :: second +integer :: minute +integer :: hour +integer :: day +integer :: month +integer :: year +character(len=fm_field_name_len) :: name +real :: days_in_this_year +integer :: num_days +character(len=fm_field_name_len+1) :: str +integer, dimension(:) :: tracer_axes_1(3) +integer :: id_zt_1 +character(len=256) :: caller_str + +integer :: nn, ntr_bgc +character(len=5) :: bgc_stf='stf00' +character(len=5) :: bgc_btf='btf00' +character(len=5) :: bgc_src='src00' +character(len=64) :: name1, name2, name3, name4 + + +! ===================================================================== +! begin of executable code +! ===================================================================== + +write(stdout(),*) +write(stdout(),*) trim(note_header), & + 'Starting ', trim(package_name), ' module' + +!----------------------------------------------------------------------- +! dynamically allocate the global BIOTIC arrays +!----------------------------------------------------------------------- + +call allocate_arrays(Domain%isc, Domain%iec, Domain%jsc, Domain%jec, Domain%isd, Domain%ied, Domain%jsd, Domain%jed, grid%nk) + + +!----------------------------------------------------------------------- +! Open up the files for boundary conditions +!----------------------------------------------------------------------- + +atmpress_id = init_external_field(atmpress_file, & + atmpress_name, & + domain = Domain%domain2d) +if (atmpress_id .eq. 0) then !{ + call mpp_error(FATAL, trim(error_header) // & + 'Could not open atmpress file: ' // & + trim(atmpress_file)) +endif !} + +pistonveloc_id = init_external_field(pistonveloc_file, & + pistonveloc_name, & + domain = Domain%domain2d) +if (pistonveloc_id .eq. 0) then !{ + call mpp_error(FATAL, trim(error_header) // & + 'Could not open pistonveloc file: ' // & + trim(pistonveloc_file)) +endif !} + +!RASF I think the ifdafs are redundant +if (id_adic .ne. 0 .and. .not. use_access_co2) then + aco2_id = init_external_field(aco2_file, & + aco2_name, & + domain = Domain%domain2d) + if (aco2_id .eq. 0) then !{ + call mpp_error(FATAL, trim(error_header) // & + 'Could not open aco2 file: ' // & + trim(aco2_file)) + endif !} +endif + +seaicefract_id = init_external_field(seaicefract_file, & + seaicefract_name, & + domain = Domain%domain2d) +if (seaicefract_id .eq. 0) then !{ + call mpp_error(FATAL, trim(error_header) // & + 'Could not open seaicefract file: ' // & + trim(seaicefract_file)) +endif !} + +dust_id = init_external_field(dust_file, & + dust_name, & + domain = Domain%domain2d) +if (dust_id .eq. 0) then !{ + call mpp_error(FATAL, trim(error_header) // & + 'Could not open dust file: ' // & + trim(dust_file)) +endif !} + + +alphabio_id = init_external_field("INPUT/bgc_param.nc", & + "alphabio", domain = Domain%domain2d) + +if (alphabio_id .eq. 0) then !{ + call mpp_error(FATAL, trim(error_header) // & + 'Could not open bgc_param.nc file' ) +endif !} + +parbio_id = init_external_field("INPUT/bgc_param.nc", & + "parbio", domain = Domain%domain2d) +kwbio_id = init_external_field("INPUT/bgc_param.nc", & + "kwbio", domain = Domain%domain2d) +kcbio_id = init_external_field("INPUT/bgc_param.nc", & + "kcbio", domain = Domain%domain2d) +abio_id = init_external_field("INPUT/bgc_param.nc", & + "abio", domain = Domain%domain2d) +bbio_id = init_external_field("INPUT/bgc_param.nc", & + "bbio", domain = Domain%domain2d) +cbio_id = init_external_field("INPUT/bgc_param.nc", & + "cbio", domain = Domain%domain2d) +k1bio_id = init_external_field("INPUT/bgc_param.nc", & + "k1bio", domain = Domain%domain2d) +muepbio_id = init_external_field("INPUT/bgc_param.nc", & + "muepbio", domain = Domain%domain2d) +muepsbio_id = init_external_field("INPUT/bgc_param.nc", & + "muepsbio", domain = Domain%domain2d) +gam1bio_id = init_external_field("INPUT/bgc_param.nc", & + "gam1bio", domain = Domain%domain2d) +gbio_id = init_external_field("INPUT/bgc_param.nc", & + "gbio", domain = Domain%domain2d) +epsbio_id = init_external_field("INPUT/bgc_param.nc", & + "epsbio", domain = Domain%domain2d) +muezbio_id = init_external_field("INPUT/bgc_param.nc", & + "muezbio", domain = Domain%domain2d) +gam2bio_id = init_external_field("INPUT/bgc_param.nc", & + "gam2bio", domain = Domain%domain2d) +muedbio_id = init_external_field("INPUT/bgc_param.nc", & + "muedbio", domain = Domain%domain2d) +muecaco3_id = init_external_field("INPUT/bgc_param.nc", & + "muecaco3", domain = Domain%domain2d) +muedbio_sed_id = init_external_field("INPUT/bgc_param.nc", & + "muedbio_sed", domain = Domain%domain2d) +muecaco3_sed_id = init_external_field("INPUT/bgc_param.nc", & + "muecaco3_sed", domain = Domain%domain2d) +wdetbio_id = init_external_field("INPUT/bgc_param.nc", & + "wdetbio", domain = Domain%domain2d) +wcaco3_id = init_external_field("INPUT/bgc_param.nc", & + "wcaco3", domain = Domain%domain2d) +nat_co2_id = init_external_field("INPUT/bgc_param.nc", & + "nat_co2", domain = Domain%domain2d) +tscav_fe_id = init_external_field("INPUT/bgc_param.nc", & + "tscav_fe", domain = Domain%domain2d) +fe_bkgnd_id = init_external_field("INPUT/bgc_param.nc", & + "fe_bkgnd", domain = Domain%domain2d) +f_inorg_id = init_external_field("INPUT/bgc_param.nc", & + "f_inorg", domain = Domain%domain2d) + +! --------------------------------- +! +! Read extra restart field(s) +! +! --------------------------------- + +do n = 1, instances !{ + id_restart(1) = register_restart_field(sed_restart, "csiro_bgc_sediment.res.nc", "caco3_sediment", biotic(n)%caco3_sediment(:,:), domain=Domain%domain2d) + id_restart(2) = register_restart_field(sed_restart, "csiro_bgc_sediment.res.nc", "det_sediment", biotic(n)%det_sediment(:,:), domain=Domain%domain2d) +enddo !} n + +call restore_state(sed_restart) + + +!----------------------------------------------------------------------- +! do some calendar calculation for the next section, +! calculate the number of days in this year as well as +! those that have already passed +!----------------------------------------------------------------------- + +call get_date(time%model_time, & + year, month, day, hour, minute, second) +num_days = days_in_year(time%model_time) +days_in_this_year = 0.0 +do m = 1, month - 1 + days_in_this_year = days_in_this_year + & + days_in_month(set_date(year, m, 1)) +enddo +days_in_this_year = days_in_this_year + day - 1 + hour/24.0 + & + minute/1440.0 + second/86400.0 + +!----------------------------------------------------------------------- +! read in additional information for a restart +!----------------------------------------------------------------------- + +write(stdout(),*) + +!----------------------------------------------------------------------- +! +! initialize some arrays which are held constant for this +! simulation +! +!----------------------------------------------------------------------- + +!----------------------------------------------------------------------- +! register the fields +!----------------------------------------------------------------------- +id_pco2 = register_diag_field('ocean_model', & + 'pco2', grid%tracer_axes(1:2), & + Time%model_time, 'pCO2', ' ', & + missing_value = -1.0e+10) +id_paco2 = register_diag_field('ocean_model', & + 'paco2', grid%tracer_axes(1:2), & + Time%model_time, 'pCO2 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', & + missing_value = -1.0e+10) +id_aco2_sat = register_diag_field('ocean_model', & + 'aco2_saturation', grid%tracer_axes(1:2), & + Time%model_time, 'CO2 saturation inc. anthropogenic', 'mmol/m^3', & + missing_value = -1.0e+10) + +id_sc_o2 = register_diag_field('ocean_model', & + 'sc_o2', grid%tracer_axes(1:2), & + Time%model_time, 'Schmidt number - O2', ' ', & + missing_value = -1.0e+10) + +id_o2_sat = register_diag_field('ocean_model', & + 'o2_saturation', grid%tracer_axes(1:2), & + Time%model_time, 'O2 saturation', 'mmolO2/m^3', & + missing_value = -1.0e+10) + +id_kw_o2 = register_diag_field('ocean_model', & + 'kw_o2', grid%tracer_axes(1:2), & + Time%model_time, 'Piston velocity - O2', ' ', & + missing_value = -1.0e+10) + +id_light_limit = register_diag_field('ocean_model','light_limit', & + grid%tracer_axes(1:2),Time%model_time, 'Integrated light limitation of phytoplankton growth', & + ' ',missing_value = -1.0e+10) + +id_pprod_gross = register_diag_field('ocean_model','pprod_gross', & + grid%tracer_axes(1:3),Time%model_time, 'Gross PHY production', & + 'mmolN/m^3/s',missing_value = -1.0e+10) + +id_pprod_gross_2d = register_diag_field('ocean_model','pprod_gross_2d', & + grid%tracer_axes(1:2),Time%model_time, 'Vertically integrated Gross PHY production', & + 'mmolN/m^2/s',missing_value = -1.0e+10) + +id_zprod_gross = register_diag_field('ocean_model','zprod_gross', & + grid%tracer_axes(1:3),Time%model_time, 'Gross ZOO production', & + 'mmolN/m^3/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) + +id_det_sediment = register_diag_field('ocean_model','det_sediment', & + grid%tracer_axes(1:2),Time%model_time, 'Accumulated DET in sediment at base of water column', & + 'mmolN/m^2',missing_value = -1.0e+10) + +id_caco3_sed_remin = register_diag_field('ocean_model','caco3_sed_remin', & + grid%tracer_axes(1:2),Time%model_time, 'Rate of remineralisation of CaCO3 in accumulated sediment', & + 'mmolN/m^2',missing_value = -1.0e+10) + +id_det_sed_remin = register_diag_field('ocean_model','det_sed_remin', & + grid%tracer_axes(1:2),Time%model_time, 'Rate of remineralisation of DET in accumulated sediment', & + 'mmolN/m^2',missing_value = -1.0e+10) + +id_caco3_sed_depst = register_diag_field('ocean_model','caco3_sed_depst', & + grid%tracer_axes(1:2),Time%model_time, 'Rate of deposition of CaCO3 to sediment at base of water column', & + 'mmolN/m^2',missing_value = -1.0e+10) + +id_det_sed_depst = register_diag_field('ocean_model','det_sed_depst', & + grid%tracer_axes(1:2),Time%model_time, 'Rate of deposition of DET to sediment at base of water column', & + 'mmolN/m^2',missing_value = -1.0e+10) + +id_total_co2_flux = register_diag_field('ocean_model','total_co2_flux', & + Time%model_time, 'Total surface flux of inorganic C (natural) into ocean', & + 'Pg/yr',missing_value = -1.0e+30) + +id_total_aco2_flux = register_diag_field('ocean_model','total_aco2_flux', & + Time%model_time, 'Total surface flux of inorganic C (natural + anthropogenic) into ocean', & + 'Pg/yr',missing_value = -1.0e+30) + +do n = 1, instances !{ + + if (instances .eq. 1) then !{ + str = ' ' + else !}{ + str = '_' // biotic(n)%name + endif !} + +! use generic definition for +! surface tracer fluxes -stf1, stf2 +! virtual tracer flux - vstf1 +! sources - src1 + ntr_bgc=biotic(n)%ntr_bgc + + do nn=1,ntr_bgc + if (nn .ge. 10) then + write(bgc_stf(4:5),'(i2)') nn + write(bgc_src(4:5),'(i2)') nn + write(bgc_btf(4:5),'(i2)') nn + else + write(bgc_stf(4:4),'(i1)') 0 + write(bgc_src(4:4),'(i1)') 0 + write(bgc_btf(4:4),'(i1)') 0 + write(bgc_stf(5:5),'(i1)') nn + write(bgc_src(5:5),'(i1)') nn + write(bgc_btf(5:5),'(i1)') nn + endif + +! default field names... + name1 = bgc_stf // ' flux into ocean' + name2 = bgc_stf // ' virtual flux into ocean' + name3 = bgc_src // ' source term' + name4 = bgc_btf // ' flux into sediment' +! specified field names + if (nn .eq. id_no3) then + name1 = 'Flux into ocean - nitrate' + name2 = 'Virtual flux into ocean - nitrate' + name3 = 'Source term - nitrate' + name4 = 'Flux into sediment - nitrate' + endif + if (nn .eq. id_phy) then + name1 = 'Flux into ocean - phytoplankton' + name2 = 'Virtual flux into ocean - phytoplankton' + name3 = 'Source term - phytoplankton' + name4 = 'Flux into sediment - phytoplankton' + endif + if (nn .eq. id_zoo) then + name1 = 'Flux into ocean - zooplankton' + name2 = 'Virtual flux into ocean - zooplankton' + name3 = 'Source term - zooplankton' + name4 = 'Flux into sediment - zooplankton' + endif + if (nn .eq. id_det) then + name1 = 'Flux into ocean - detritus' + name2 = 'Virtual flux into ocean - detritus' + name3 = 'Source term - detritus' + name4 = 'Flux into sediment - detritus' + endif + if (nn .eq. id_o2) then + name1 = 'Flux into ocean - oxygen' + name2 = 'Virtual flux into ocean - oxygen' + name3 = 'Source term - oxygen' + name4 = 'Flux into sediment - oxygen' + endif + if (nn .eq. id_po4) then + name1 = 'Flux into ocean - phosphate' + name2 = 'Virtual flux into ocean - phosphate' + name3 = 'Source term - phosphate' + name4 = 'Flux into sediment - phosphate' + endif + if (nn .eq. id_dic) then + name1 = 'Flux into ocean - DIC, PI' + name2 = 'Virtual flux into ocean - DIC, PI' + name3 = 'Source term - DIC, PI' + name4 = 'Flux into sediment - DIC, PI' + endif + if (nn .eq. id_adic) then + name1 = 'Flux into ocean - DIC, inc. anth.' + name2 = 'Virtual flux into ocean - DIC, inc. anth.' + name3 = 'Source term - DIC, inc. anth.' + name4 = 'Flux into sediment - DIC, inc. anth.' + endif + if (nn .eq. id_alk) then + name1 = 'Flux into ocean - alkalinity' + name2 = 'Virtual flux into ocean - alkalinity' + name3 = 'Source term - alkalinity' + name4 = 'Flux into sediment - alkalinity' + endif + if (nn .eq. id_caco3) then + name1 = 'Flux into ocean - calcium carbonate' + name2 = 'Virtual flux into ocean - calcium carbonate' + name3 = 'Source term - calcium carbonate' + name4 = 'Flux into sediment - calcium carbonate' + endif + if (nn .eq. id_fe) then + name1 = 'Flux into ocean - iron' + name2 = 'Virtual flux into ocean - iron' + name3 = 'Source term - iron' + name4 = 'Flux into sediment - iron' + endif + if (mpp_pe() == mpp_root_pe() )print*,'rjm bio',bgc_stf,'v'//bgc_stf + + biotic(n)%id_bgc_stf(nn) = register_diag_field('ocean_model', & + bgc_stf//str, grid%tracer_axes(1:2), & + Time%model_time, name1, 'mmol/m^2/s', & +! Time%model_time, bgc_stf//'flux into ocean', 'mmol/m^2/s', & + missing_value = -1.0e+10) + + biotic(n)%id_bgc_vstf(nn) = register_diag_field('ocean_model', & + 'v'//bgc_stf//str, grid%tracer_axes(1:2), & + Time%model_time, name2, 'mmol/m^3/s', & +! Time%model_time, bgc_stf//'virtual flux into ocean', 'mmol/m^3/s', & + missing_value = -1.0e+10) + + biotic(n)%id_bgc_src(nn) = register_diag_field('ocean_model', & + bgc_src//str, grid%tracer_axes(1:3), & + Time%model_time, name3, 'mmolN/m^3/s', & +! Time%model_time, bgc_src, 'mmolN/m^3/s', & + missing_value = -1.0e+10) + + biotic(n)%id_bgc_btf(nn) = register_diag_field('ocean_model', & + bgc_btf//str, grid%tracer_axes(1:2), & + Time%model_time, name4, 'mmol/m^2/s', & +! Time%model_time, bgc_btf//'flux into sediment', 'mmol/m^2/s', & + missing_value = -1.0e+10) + + enddo !} nn + +enddo !} n + +!----------------------------------------------------------------------- +! give info +!----------------------------------------------------------------------- + +do n = 1, instances !{ + biotic(n)%sal_global=34.6 + if (id_dic.ne.0) biotic(n)%bgc_global(id_dic)=1966. + if (id_adic.ne.0) biotic(n)%bgc_global(id_adic)=1966. + if (id_alk.ne.0) biotic(n)%bgc_global(id_alk)=2303. + if (id_no3.ne.0) biotic(n)%bgc_global(id_no3)=4.67 +enddo + +write(stdout(),*) +write(stdout(),*) trim(note_header), 'Tracer runs initialized' +write(stdout(),*) + + +return +end subroutine csiro_bgc_start !} +! NAME="csiro_bgc_start" + + +!####################################################################### +! +! +! +! Perform things that should be done in tracer, but are done here +! in order to minimize the number of hooks necessary in the MOM4 basecode +! +! + +subroutine csiro_bgc_tracer (isc, iec, jsc, jec, t_prog, grid, time, dtts) !{ + +use mpp_mod, only : mpp_sum + + +type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog +type(ocean_grid_type), intent(in) :: Grid +type(ocean_time_type), intent(in) :: Time +integer, intent(in) :: isc, iec +integer, intent(in) :: jsc, jec +real, intent(in) :: dtts + +!----------------------------------------------------------------------- +! local definitions +!----------------------------------------------------------------------- + +character(len=64), parameter :: sub_name = 'csiro_bgc_tracer' +character(len=256), parameter :: error_header = & + '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: warn_header = & + '==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):' +character(len=256), parameter :: note_header = & + '==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + +integer :: i +integer :: j +integer :: n +integer :: k +integer :: nn, ntr_bgc, ind_trc +real :: total_trc +integer :: indsal +logical :: used + + +indsal = fm_get_index('/ocean_mod/prog_tracers/salt') + +!----------------------------------------------------------------------- +! accumulate global annual means +!----------------------------------------------------------------------- + + +do n = 1, instances !{ +! use generic definition for bgc tracers -tracer1, tracer2 + ntr_bgc=biotic(n)%ntr_bgc + + ! Use zero_floor = .false. to skip this floor to BGC tracers and allow the time_tendency calculation to non-negative values bring any small negatives up. + ! This then improves the 'mismatch' in diagnostic output. + ! mac, aug11. + if (zero_floor) then + do nn=1,ntr_bgc + + + do k = 1,grid%nk !{ + do j = jsc, jec !{ + do i = isc, iec !{ + t_prog(biotic(n)%ind_bgc(nn))%field(i,j,k,time%taup1)= & + max(0.,t_prog(biotic(n)%ind_bgc(nn))%field(i,j,k,time%taup1) ) + enddo !} i + enddo !} j + enddo !} k + + enddo !} nn + endif ! zero_floor + +! comment out this sediment source of fe while testing equivalent code in csiro_bgc_bbc. mac, nov12. + +! rjm bottom Fe fix +! mac aug10, only apply this fix when the water is <= 200 m deep. + if (id_fe.ne.0) then + do j = jsc, jec !{ + do i = isc, iec !{ + if (grid%kmt(i,j) .gt. 0) then + k = grid%kmt(i,j) + if (grid%zw(k) .le. 200) & + t_prog(biotic(n)%ind_bgc(id_fe))%field(i,j,k,time%taup1)= 0.999 + endif + enddo !} i + enddo !} j + endif + +! apply deposition and remineralisation rates to sediment tracers. mac, nov12. + do j = jsc, jec !{ + do i = isc, iec !{ + if (grid%kmt(i,j) .gt. 0) then + biotic(n)%det_sediment(i,j) = biotic(n)%det_sediment(i,j) + dtts* & + ( biotic(n)%det_sed_depst(i,j) - & + biotic(n)%det_sed_remin(i,j) ) + biotic(n)%caco3_sediment(i,j) = biotic(n)%caco3_sediment(i,j) + dtts* & + ( biotic(n)%caco3_sed_depst(i,j) - & + biotic(n)%caco3_sed_remin(i,j) ) + endif + enddo !} i + enddo !} j + + ! send sediment tracers to output. mac, nov12. + if (id_caco3_sediment .gt. 0) then + used = send_data(id_caco3_sediment, biotic(n)%caco3_sediment(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) + endif + if (id_det_sediment .gt. 0) then + used = send_data(id_det_sediment, biotic(n)%det_sediment(isc:iec,jsc:jec), & + time%model_time, rmask = grid%tmask(isc:iec,jsc:jec,1)) + endif + +enddo !} n + +return +end subroutine csiro_bgc_tracer !} +! NAME="csiro_bgc_tracer" + +include "bio_v3.inc" + +end module csiro_bgc_mod !} diff --git a/src/mom5/ocean_csiro_bgc/ocmip2_co2calc.F90 b/src/mom5/ocean_csiro_bgc/ocmip2_co2calc.F90 new file mode 100755 index 0000000000..fa3ac5ee7f --- /dev/null +++ b/src/mom5/ocean_csiro_bgc/ocmip2_co2calc.F90 @@ -0,0 +1,729 @@ +! ---------------------------------------------------------------- +! GNU General Public License +! This file is a part of MOM. +! +! MOM is free software; you can redistribute it and/or modify it and +! are expected to follow the terms of the GNU General Public License +! as published by the Free Software Foundation; either version 2 of +! the License, or (at your option) any later version. +! +! MOM is distributed in the hope that it will be useful, but WITHOUT +! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +! License for more details. +! +! For the full text of the GNU General Public License, +! write to: Free Software Foundation, Inc., +! 675 Mass Ave, Cambridge, MA 02139, USA. +! or see: http://www.gnu.org/licenses/gpl.html +!----------------------------------------------------------------------- +! +! +! Richard D. Slater +! +! +! John P. Dunne +! +! +! +! Surface fCO2 calculation +! +! +! +! Calculate the fugacity of CO2 at the surface in thermodynamic +! equilibrium with the current alkalinity (Alk) and total dissolved +! inorganic carbon (DIC) at a particular temperature and salinity +! using an initial guess for the total hydrogen +! ion concentration (htotal) +! +! + +module ocmip2_co2calc_mod + +implicit none + +private + +public :: ocmip2_co2calc +public :: ocmip2_co2_alpha + +character(len=128) :: version = '$Id: ocmip2_co2calc.F90,v 20.0 2013/12/14 00:09:44 fms Exp $' +character(len=128) :: tagname = '$Name: tikal $' + +contains + + +!####################################################################### +! +! +! +! Calculate CO2 solubility, alpha, from +! temperature (t) and salinity (s). +! +! INPUT +! +! isd = first i-limit of the arrays with halo +! ied = last i-limit of the arrays with halo +! jsd = first j-limit of the arrays with halo +! jed = last j-limit of the arrays with halo +! isc = first i-limit of the arrays for computation +! iec = last i-limit of the arrays for computation +! jsc = first j-limit of the arrays for computation +! jec = last j-limit of the arrays for computation +! +! t = temperature (degrees C) +! +! s = salinity (PSU) +! +! mask = land mask array (0.0 = land) +! +! OUTPUT +! alpha = Solubility of CO2 for air? (mol/kg/atm unless scaled) +! +! IMPORTANT: Some words about units - (JCO, 4/4/1999) +! +! - Models may carry tracers in mol/m^3 (on a per volume basis) +! +! - Conversely, this routine, which was written by observationalists +! (C. Sabine and R. Key), passes input arguments in umol/kg +! (i.e., on a per mass basis) +! +! - Thus, if the input or output units need to be changed from/to mol/m^3 +! then set scale to an appropriate value. For example, if the model +! uses mol/m^3, then scale should be set to something like 1.0/1024.5 +! to convert from mol/m^3 to mol/kg. +! +! +subroutine ocmip2_co2_alpha(isd, jsd, isc, iec, jsc, jec, t, s, mask, alpha, scale) + +integer, intent(in) :: isd +integer, intent(in) :: jsd +integer, intent(in) :: isc +integer, intent(in) :: iec +integer, intent(in) :: jsc +integer, intent(in) :: jec +real, dimension(isd:,jsd:), intent(in) :: t +real, dimension(isd:,jsd:), intent(in) :: s +real, dimension(isd:,jsd:), intent(in) :: mask +real, dimension(isc:,jsc:), intent(out) :: alpha +real, intent(in), optional :: scale + +integer :: i, j +real :: log100 +real :: tk +real :: tk100 +real :: tk1002 +real :: logtk +real :: ff +real :: scale_factor + +! set the scale factor for unit conversion +if (present(scale)) then + scale_factor = scale +else + scale_factor = 1.0 +endif + +! Calculate all constants needed to convert between various measured +! carbon species. References for each equation are noted in the code. +! Once calculated, the constants are stored and passed in the common +! block "const". The original version of this code was based on +! the code by Dickson in Version 2 of "Handbook of Methods for the +! Analysis of the Various Parameters of the Carbon Dioxide System +! in Seawater", DOE, 1994 (SOP No. 3, p25-26). +! +! Derive simple terms used more than once +log100 = log(100.0) +do j = jsc, jec + do i = isc, iec + + if (mask(i,j) .ne. 0.0) then + + tk = 273.15 + t(i,j) + tk100 = tk / 100.0 + tk1002 = tk100 * tk100 + logtk = log(tk) + + ! f = k0(1-pH2O)*correction term for non-ideality + + ! Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 + ! values) + ff = exp(-162.8301 + 218.2968 / tk100 + & + 90.9241 * (logtk - log100) - 1.47696 * tk1002 + & + s(i,j) * (0.025695 - 0.025225 * tk100 + 0.0049867 * tk1002)) + + ! Should we be using K0 or ff for the solubility here? + ! convert to output units + alpha(i,j) = ff / scale_factor + else + alpha(i,j) = 0.0 + endif + enddo +enddo + +return + +end subroutine ocmip2_co2_alpha +! NAME="ocmip2_co2_alpha" + + +!####################################################################### +! +! +! +! Calculate co2* from total alkalinity and total CO2 at +! temperature (t) and salinity (s). +! It is assumed that init_ocmip2_co2calc has already been called with +! the T and S to calculate the various coefficients. +! +! INPUT +! +! isd = first i-limit of the arrays with halo +! ied = last i-limit of the arrays with halo +! jsd = first j-limit of the arrays with halo +! jed = last j-limit of the arrays with halo +! isc = first i-limit of the arrays for computation +! iec = last i-limit of the arrays for computation +! jsc = first j-limit of the arrays for computation +! jec = last j-limit of the arrays for computation +! +! mask = land mask array (0.0 = land) +! +! t = temperature (degrees C) +! +! s = salinity (PSU) +! +! dic_in = total inorganic carbon (mol/kg unless scaled) +! where 1 T = 1 metric ton = 1000 kg +! +! ta_in = total alkalinity (eq/kg unless scaled) +! +! pt_in = inorganic phosphate (mol/kg unless scaled) +! +! sit_in = inorganic silicate (mol/kg unless scaled) +! +! htotallo = factor to set lower limit of htotal range +! +! htotalhi = factor to set upper limit of htotal range +! +! htotal = H+ concentraion +! +! OUTPUT +! co2star = CO2*water (kg/kg unless scaled) +! alpha = Solubility of CO2 for air? (kg/kg/atm unless scaled) +! pco2surf = oceanic pCO2 (ppmv) +! +! k1 = activity factors for carbonate species +! +! k2 (see below) +! +! invtk = 1/(t+273.15) +! +! IMPORTANT: Some words about units - (JCO, 4/4/1999) +! +! - Models may carry tracers in mol/m^3 (on a per volume basis) +! +! - Conversely, this routine, which was written by observationalists +! (C. Sabine and R. Key), passes input arguments in umol/kg +! (i.e., on a per mass basis) +! +! - Thus, if the input or output units need to be changed from/to mol/m^3 +! then set scale to an appropriate value. For example, if the model +! uses mol/m^3, then scale should be set to something like 1.0/1024.5 +! to convert from mol/m^3 to mol/kg. +! +! +subroutine ocmip2_co2calc(isd, jsd, isc, iec, jsc, jec, mask, t, s, & + dic_in, ta_in, pt_in, sit_in, htotallo, htotalhi, htotal, & + co2star, co3_ion, alpha, pCO2surf, k1_out, k2_out, invtk_out, scale) + +real, parameter :: permeg = 1.e-6 +real, parameter :: xacc = 1.0e-10 + +integer, intent(in) :: isd +integer, intent(in) :: jsd +integer, intent(in) :: isc +integer, intent(in) :: iec +integer, intent(in) :: jsc +integer, intent(in) :: jec +real, dimension(isd:,jsd:), intent(in) :: mask +real, dimension(isd:,jsd:), intent(in) :: t +real, dimension(isd:,jsd:), intent(in) :: s +real, dimension(isd:,jsd:), intent(in) :: dic_in +real, dimension(isd:,jsd:), intent(in) :: ta_in +real, dimension(isd:,jsd:), intent(in) :: pt_in +real, dimension(isd:,jsd:), intent(in) :: sit_in +real, dimension(isc:,jsc:), intent(in) :: htotallo +real, dimension(isc:,jsc:), intent(in) :: htotalhi +real, dimension(isc:,jsc:), intent(inout) :: htotal +real, dimension(isc:,jsc:), intent(out), optional :: co2star +real, dimension(isc:,jsc:), intent(out), optional :: co3_ion +real, dimension(isc:,jsc:), intent(out), optional :: alpha +real, dimension(isc:,jsc:), intent(out), optional :: pCO2surf +real, dimension(isc:,jsc:), intent(out), optional :: invtk_out +real, dimension(isc:,jsc:), intent(out), optional :: k1_out +real, dimension(isc:,jsc:), intent(out), optional :: k2_out +real, intent(in), optional :: scale + +integer :: i +integer :: j +real :: log100 +real :: pt +real :: sit +real :: ta +real :: dic +real :: tk +real :: invtk +real :: tk100 +real :: tk1002 +real :: logtk +real :: is +real :: is2 +real :: sqrtis +real :: s2 +real :: sqrts +real :: s15 +real :: scl +real :: logf_of_s +real :: ff +!real :: k0 +real :: k1 +real :: k2 +real :: kb +real :: k1p +real :: k2p +real :: k3p +real :: ksi +real :: kw +real :: ks +real :: kf +real :: bt +real :: st +real :: ft +real :: htotal2 +real :: co2star_internal +real :: scale_factor + +! set the scale factor for unit conversion +if (present(scale)) then + scale_factor = scale +else + scale_factor = 1.0 +endif + +! Calculate all constants needed to convert between various measured +! carbon species. References for each equation are noted in the code. +! Once calculated, the constants are stored and passed in the common +! block "const". The original version of this code was based on +! the code by Dickson in Version 2 of "Handbook of Methods for the +! Analysis of the Various Parameters of the Carbon Dioxide System +! in Seawater", DOE, 1994 (SOP No. 3, p25-26). +! +! Derive simple terms used more than once +log100 = log(100.0) +do j = jsc, jec + do i = isc, iec + + if (mask(i,j) .ne. 0.0) then + + tk = 273.15 + t(i,j) + tk100 = tk / 100.0 + tk1002 = tk100 * tk100 + invtk = 1.0 / tk + logtk = log(tk) + is = 19.924 * s(i,j) / (1000.0 - 1.005 * s(i,j)) + is2 = is * is + sqrtis = sqrt(is) + s2 = s(i,j) * s(i,j) + sqrts = sqrt(s(i,j)) + s15 = sqrts ** 3 + scl = s(i,j) / 1.80655 + logf_of_s = log(1.0 - 0.001005 * s(i,j)) + + ! k0 from Weiss 1974 + + !k0 = exp(93.4517/tk100 - 60.2409 + 23.3585 * log(tk100) + & + !s(i,j) * (0.023517 - 0.023656 * tk100 + 0.0047036 * tk1002)) + + + ! k1 = [H][HCO3]/[H2CO3] + ! k2 = [H][CO3]/[HCO3] + + ! Millero p.664 (1995) using Mehrbach et al. data on seawater scale + k1 = 10.0 ** (-(3670.7 * invtk - 62.008 + 9.7944 * logtk - & + 0.0118 * s(i,j) + 0.000116 * s2)) + + k2 = 10.0 ** (-(1394.7 * invtk + 4.777 - 0.0184 * s(i,j) + 0.000118 * s2)) + + ! kb = [H][BO2]/[HBO2] + + ! Millero p.669 (1995) using data from Dickson (1990) + kb = exp((-8966.90 - 2890.53 * sqrts - 77.942 * s(i,j) + & + 1.728 * s15 - 0.0996 * s2) * invtk + & + (148.0248 + 137.1942 * sqrts + 1.62142 * s(i,j)) + & + (-24.4344 - 25.085 * sqrts - 0.2474 * s(i,j)) * logtk + & + 0.053105 * sqrts * tk) + + ! k1p = [H][H2PO4]/[H3PO4] + + ! DOE(1994) eq 7.2.20 with footnote using data from Millero (1974) + k1p = exp(-4576.752 * invtk + 115.525 - 18.453 * logtk + & + (-106.736 * invtk + 0.69171) * sqrts + & + (-0.65643 * invtk - 0.01844) * s(i,j)) + + ! k2p = [H][HPO4]/[H2PO4] + + ! DOE(1994) eq 7.2.23 with footnote using data from Millero (1974)) + k2p = exp(-8814.715 * invtk + 172.0883 - 27.927 * logtk + & + (-160.340 * invtk + 1.3566) * sqrts + & + (0.37335 * invtk - 0.05778) * s(i,j)) + + ! k3p = [H][PO4]/[HPO4] + + ! DOE(1994) eq 7.2.26 with footnote using data from Millero (1974) + k3p = exp(-3070.75 * invtk - 18.141 + & + (17.27039 * invtk + 2.81197) * sqrts + & + (-44.99486 * invtk - 0.09984) * s(i,j)) + + ! ksi = [H][SiO(OH)3]/[Si(OH)4] + + ! Millero p.671 (1995) using data from Yao and Millero (1995) + ksi = exp(-8904.2 * invtk + 117.385 - 19.334 * logtk + & + (-458.79 * invtk + 3.5913) * sqrtis + & + (188.74 * invtk - 1.5998) * is + & + (-12.1652 * invtk + 0.07871) * is2 + logf_of_s) + + ! kw = [H][OH] + + ! Millero p.670 (1995) using composite data + kw = exp(-13847.26 * invtk + 148.9652 - 23.6521 * logtk + & + (118.67 * invtk - 5.977 + 1.0495 * logtk) * sqrts - 0.01615 * s(i,j)) + + ! ks = [H][SO4]/[HSO4] + + ! Dickson (1990, J. chem. Thermodynamics 22, 113) + ks = exp(-4276.1 * invtk + 141.328 - 23.093 * logtk + & + (-13856.0 * invtk + 324.57 - 47.986 * logtk) * sqrtis + & + (35474.0 * invtk - 771.54 + 114.723 * logtk) * is - & + 2698.0 * invtk * sqrtis ** 3 + 1776.0 * invtk * is2 + logf_of_s) + + ! kf = [H][F]/[HF] + + ! Dickson and Riley (1979) -- change pH scale to total + kf = exp(1590.2 * invtk - 12.641 + 1.525 * sqrtis + logf_of_s + & + log(1.0 + (0.1400 / 96.062) * scl / ks)) + + ! Calculate concentrations for borate, sulfate, and fluoride + + ! Uppstrom (1974) + bt = 0.000232 / 10.811 * scl + + ! Morris & Riley (1966) + st = 0.14 / 96.062 * scl + + ! Riley (1965) + ft = 0.000067 / 18.9984 * scl + + ! set some stuff to pass back, if requested + if (present(k1_out)) then + k1_out(i,j) = k1 + endif + + if (present(k2_out)) then + k2_out(i,j) = k2 + endif + + if (present(invtk_out)) then + invtk_out(i,j) = invtk + endif + + ! Possibly convert input in mol/m^3 -> mol/kg + sit = sit_in(i,j) * scale_factor + ta = ta_in(i,j) * scale_factor + dic = dic_in(i,j) * scale_factor + pt = pt_in(i,j) * scale_factor + +! +!*********************************************************************** +! +! Calculate [H+] total when DIC and TA are known at T, S and 1 atm. +! The solution converges to err of xacc. The solution must be within +! the range x1 to x2. +! +! If DIC and TA are known then either a root finding or iterative method +! must be used to calculate htotal. In this case we use the +! Newton-Raphson "safe" method taken from "Numerical Recipes" +! (function "rtsafe.f" with error trapping removed). +! +! As currently set, this procedure iterates about 12 times. The x1 +! and x2 values set below will accomodate ANY oceanographic values. +! If an initial guess of the pH is known, then the number of +! iterations can be reduced to about 5 by narrowing the gap between +! x1 and x2. It is recommended that the first few time steps be run +! with x1 and x2 set as below. After that, set x1 and x2 to the +! previous value of the pH +/- ~0.5. The current setting of xacc will +! result in co2star accurate to 3 significant figures (xx.y). Making +! xacc bigger will result in faster convergence also, but this is not +! recommended (xacc of 10**-9 drops precision to 2 significant +! figures). +! + + htotal(i,j) = drtsafe(htotalhi(i,j)*htotal(i,j), & + htotallo(i,j)*htotal(i,j), & + dic, ta, pt, sit, k1, k2, k1p, k2p, k3p, & + bt, ft, st, kb, kw, kf, ks, ksi, xacc) + + ! Calculate [CO2*] as defined in DOE Methods Handbook 1994 Ver.2, + ! ORNL/CDIAC-74, Dickson and Goyet, eds. (Ch 2 p 10, Eq A.49) + ! Convert units of output arguments + ! Note: co2star is calculated in + ! mol/kg within this routine + ! + htotal2 = htotal(i,j) * htotal(i,j) + co2star_internal = dic * htotal2 / (htotal2 + k1 * (htotal(i,j) + k2)) / scale_factor + if (present(co2star)) then + co2star(i,j) = co2star_internal + endif + if (present(co3_ion)) then + co3_ion(i,j) = co2star_internal * k1 * k2 / htotal2 + endif + !ph = -log10(htotal(i,j)) + + ! f = k0(1-pH2O)*correction term for non-ideality + + ! Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 + ! values) + if (present(alpha) .or. present(pco2surf)) then + + ff = exp(-162.8301 + 218.2968 / tk100 + & + 90.9241 * (logtk - log100) - 1.47696 * tk1002 + & + s(i,j) * (0.025695 - 0.025225 * tk100 + 0.0049867 * tk1002)) + + ! Add two output arguments for storing pCO2surf + ! Should we be using K0 or ff for the solubility here? + ! Convert pCO2surf from atm to uatm + ! possibly convert output units + if (present(alpha)) then + alpha(i,j) = ff / scale_factor + endif + if (present(pco2surf)) then + pCO2surf(i,j) = co2star_internal / (ff / scale_factor) / permeg + endif + + endif + + else + + if (present(co2star)) then + co2star(i,j) = 0.0 + endif + if (present(k1_out)) then + k1_out(i,j) = 0.0 + endif + if (present(k2_out)) then + k2_out(i,j) =0.0 + endif + if (present(invtk_out)) then + invtk_out(i,j) = 0.0 + endif + if (present(alpha)) then + alpha(i,j) = 0.0 + endif + if (present(pco2surf)) then + pCO2surf(i,j) = 0.0 + endif + + endif + + enddo +enddo + +return + +end subroutine ocmip2_co2calc +! NAME="ocmip2_co2calc" + + +!####################################################################### +! +! +! +! File taken from Numerical Recipes. Modified R. M. Key 4/94 +! + +function drtsafe(x1, x2, dic, ta, pt, sit, k1, k2, k1p, k2p, k3p, & + bt, ft, st, kb, kw, kf, ks, ksi, xacc) + +real :: drtsafe +real, intent(in) :: x1 +real, intent(in) :: x2 +real, intent(in) :: dic +real, intent(in) :: ta +real, intent(in) :: pt +real, intent(in) :: sit +real, intent(in) :: k1 +real, intent(in) :: k2 +real, intent(in) :: k1p +real, intent(in) :: k2p +real, intent(in) :: k3p +real, intent(in) :: bt +real, intent(in) :: ft +real, intent(in) :: st +real, intent(in) :: kb +real, intent(in) :: kw +real, intent(in) :: kf +real, intent(in) :: ks +real, intent(in) :: ksi +real, intent(in) :: xacc + +integer, parameter :: maxit = 100 + +integer :: j +real :: fl +real :: df +real :: fh +real :: swap +real :: xl +real :: xh +real :: dxold +real :: dx +real :: f +real :: temp + +call ta_iter_1(x1, fl, df, dic, ta, pt, sit, k1, k2, & + k1p, k2p, k3p, bt, ft, st, kb, kw, kf, ks, ksi) +call ta_iter_1(x2, fh, df, dic, ta, pt, sit, k1, k2, & + k1p, k2p, k3p, bt, ft, st, kb, kw, kf, ks, ksi) +if(fl .lt. 0.0) then + xl=x1 + xh=x2 +else + xh=x1 + xl=x2 + swap=fl + fl=fh + fh=swap +end if +drtsafe=0.5*(x1+x2) +dxold=abs(x2-x1) +dx=dxold +call ta_iter_1(drtsafe, f, df, dic, ta, pt, sit, k1, k2, & + k1p, k2p, k3p, bt, ft, st, kb, kw, kf, ks, ksi) +do j=1,maxit + if (((drtsafe-xh)*df-f)*((drtsafe-xl)*df-f) .ge. 0.0 .or. & + abs(2.0*f) .gt. abs(dxold*df)) then + dxold=dx + dx=0.5*(xh-xl) + drtsafe=xl+dx + if (xl .eq. drtsafe) then + return + endif + else + dxold=dx + dx=f/df + temp=drtsafe + drtsafe=drtsafe-dx + if (temp .eq. drtsafe) then + return + endif + end if + if (abs(dx) .lt. xacc) then + return + endif + call ta_iter_1(drtsafe, f, df, dic, ta, pt, sit, k1, k2, & + k1p, k2p, k3p, bt, ft, st, kb, kw, kf, ks, ksi) + if(f .lt. 0.0) then + xl=drtsafe + fl=f + else + xh=drtsafe + fh=f + end if +enddo + +return + +end function drtsafe +! NAME="drtsafe" + + +!####################################################################### +! +! +! +! This routine expresses TA as a function of DIC, htotal and constants. +! It also calculates the derivative of this function with respect to +! htotal. It is used in the iterative solution for htotal. In the call +! "x" is the input value for htotal, "fn" is the calculated value for TA +! and "df" is the value for dTA/dhtotal +! + +subroutine ta_iter_1(x, fn, df, dic, ta, pt, sit, k1, k2, & + k1p, k2p, k3p, bt, ft, st, kb, kw, kf, ks, ksi) + +real, intent(in) :: x +real, intent(out) :: fn +real, intent(out) :: df +real, intent(in) :: dic +real, intent(in) :: ta +real, intent(in) :: pt +real, intent(in) :: sit +real, intent(in) :: k1 +real, intent(in) :: k2 +real, intent(in) :: k1p +real, intent(in) :: k2p +real, intent(in) :: k3p +real, intent(in) :: bt +real, intent(in) :: ft +real, intent(in) :: st +real, intent(in) :: kb +real, intent(in) :: kw +real, intent(in) :: kf +real, intent(in) :: ks +real, intent(in) :: ksi + +real :: x2 +real :: x3 +real :: k12 +real :: k12p +real :: k123p +real :: c +real :: a +real :: a2 +real :: da +real :: b +real :: b2 +real :: db + +x2 = x*x +x3 = x2*x +k12 = k1*k2 +k12p = k1p*k2p +k123p = k12p*k3p +c = 1.0 + st/ks +a = x3 + k1p*x2 + k12p*x + k123p +a2 = a*a +da = 3.0*x2 + 2.0*k1p*x + k12p +b = x2 + k1*x + k12 +b2 = b*b +db = 2.0*x + k1 + +! fn = hco3+co3+borate+oh+hpo4+2*po4+silicate+hfree+hso4+hf+h3po4-ta +fn = k1*x*dic/b + 2.0*dic*k12/b + bt/(1.0 + x/kb) + kw/x + & + pt*k12p*x/a + 2.0*pt*k123p/a + sit/(1.0 + x/ksi) - & + x/c - st/ (1.0 + ks/x/c) - ft/(1.0 + kf/x) - & + pt*x3/a - ta + +! df = dfn/dx +df = ((k1*dic*b) - k1*x*dic*db)/b2 - 2.0*dic*k12*db/b2 - & + bt/kb/(1.0+x/kb)**2 - kw/x2 + (pt*k12p*(a - x*da))/a2 - & + 2.0*pt*k123p*da/a2 - sit/ksi/(1.0+x/ksi)**2 - 1.0/c + & + st*(1.0 + ks/x/c)**(-2)*(ks/c/x2) + ft*(1.0 + kf/x)**(-2)*kf/x2 - & + pt*x2*(3.0*a-x*da)/a2 + +return + +end subroutine ta_iter_1 +! NAME="ta_iter_1" + +end module ocmip2_co2calc_mod diff --git a/src/mom5/ocean_tracers/ocean_tpm.F90 b/src/mom5/ocean_tracers/ocean_tpm.F90 index e9837d2d4d..d836419038 100644 --- a/src/mom5/ocean_tracers/ocean_tpm.F90 +++ b/src/mom5/ocean_tracers/ocean_tpm.F90 @@ -225,7 +225,7 @@ module ocean_tpm_mod !{ use ocean_generic_mod, only: ocean_generic_column_physics use ocean_generic_mod, only: ocean_generic_end use ocean_generic_mod, only: ocean_generic_flux_init -#endif +#endif use ocean_frazil_mod, only: ocean_frazil_init use ocean_tempsalt_mod, only: ocean_tempsalt_init @@ -237,6 +237,10 @@ module ocean_tpm_mod !{ use transport_matrix_mod, only: transport_matrix_start use transport_matrix_mod, only: transport_matrix_store_implicit +#if defined(CSIRO_BGC) +use csiro_bgc_mod +#endif + ! ! force all variables to be "typed" ! @@ -380,7 +384,7 @@ subroutine do_time_calc(time, dtts) !{ #ifdef USE_OCEAN_OCMIP2 call mpp_error(FATAL, & '==>Error in ocean_tmp_mod: cpp option USE_OCEAN_OCMIP2 is now called USE_OCEAN_BGC. Please recompile...sorry.') -#endif +#endif time_tau = time%tau @@ -570,7 +574,7 @@ end subroutine do_time_calc !} ! calculations ! -subroutine ocean_tpm_bbc(Domain, Grid, T_prog) !{ +subroutine ocean_tpm_bbc(Domain, Grid, T_prog, Time) !{ ! !----------------------------------------------------------------------- @@ -581,6 +585,7 @@ subroutine ocean_tpm_bbc(Domain, Grid, T_prog) !{ type(ocean_domain_type), intent(in) :: Domain type(ocean_grid_type), intent(in) :: Grid type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog +type(ocean_time_type), intent(in), optional :: Time ! !----------------------------------------------------------------------- @@ -605,7 +610,13 @@ subroutine ocean_tpm_bbc(Domain, Grid, T_prog) !{ Domain%isd, Domain%ied, Domain%jsd, Domain%jed, T_prog, Grid%kmt) endif !} -#endif +#endif + +#if defined(CSIRO_BGC) +if (do_csiro_bgc) then !{ + call csiro_bgc_bbc(Domain%isc, Domain%iec, Domain%jsc, Domain%jec, T_prog, Grid, Time) +endif !} +#endif return @@ -729,6 +740,13 @@ subroutine ocean_tpm_end(Domain, Grid, T_prog, T_diag, Time, Thickness) !{ #endif +#if defined(CSIRO_BGC) +if (do_csiro_bgc) then !{ + call csiro_bgc_end(Domain%isc, Domain%iec, Domain%jsc, Domain%jec, Time%taup1, & + Thickness, T_prog, Grid) +endif !} +#endif + return end subroutine ocean_tpm_end !} @@ -1200,7 +1218,9 @@ end subroutine ocean_tpm_sfc_end !} ! subroutine ocean_tpm_sbc(Domain, Grid, T_prog, Time, Ice_ocean_boundary_fluxes, & - runoff, isc_bnd, iec_bnd, jsc_bnd, jec_bnd) !{ + runoff, isc_bnd, iec_bnd, jsc_bnd, jec_bnd, aice, wnd, & + use_waterflux, salt_restore_as_salt_flux, atm_co2, co2flux, ocn_co2) + use coupler_types_mod, only: coupler_2d_bc_type @@ -1223,6 +1243,14 @@ subroutine ocean_tpm_sbc(Domain, Grid, T_prog, Time, Ice_ocean_boundary_fluxes, integer, intent(in) :: jec_bnd real, dimension(Domain%isd:,Domain%jsd:), intent(in) :: runoff +! Optional arguments. Passed through for csiro BGC. + +real, intent(in), dimension(Domain%isd:,Domain%jsd:), optional :: aice +real, intent(in), dimension(Domain%isd:,Domain%jsd:), optional :: atm_co2 +real, intent(in), dimension(Domain%isd:,Domain%jsd:), optional :: wnd +logical, intent(in), optional :: use_waterflux, salt_restore_as_salt_flux + +real, intent(out), dimension(Domain%isd:,Domain%jsd:), optional :: co2flux, ocn_co2 ! !----------------------------------------------------------------------- ! local parameters @@ -1289,9 +1317,15 @@ subroutine ocean_tpm_sbc(Domain, Grid, T_prog, Time, Ice_ocean_boundary_fluxes, Grid, Time, Ice_ocean_boundary_fluxes) endif !} - #endif +#if defined(CSIRO_BGC) +if (do_csiro_bgc) then !{ + call csiro_bgc_sbc(Domain%isc, Domain%iec, Domain%jsc, Domain%jec, Domain%isd, Domain%ied, Domain%jsd, Domain%jed, & + T_prog, aice, wnd, Grid, Time, use_waterflux, salt_restore_as_salt_flux, atm_co2, co2flux, ocn_co2) +endif !} +#endif + return end subroutine ocean_tpm_sbc !} @@ -1382,6 +1416,10 @@ subroutine ocean_tpm_init(Domain, Grid, Time, Time_steps, & #endif +#if defined(CSIRO_BGC) +call csiro_bgc_init +#endif + call transport_matrix_init @@ -1474,7 +1512,7 @@ end subroutine ocean_tpm_flux_init !} ! subroutine ocean_tpm_source(isd, ied, jsd, jed, Domain, Grid, T_prog, T_diag, & - Time, Thickness, Dens, hblt_depth, dtts) + Time, Thickness, Dens, hblt_depth, dtts, swflx, sw_frac_zt) implicit none @@ -1497,6 +1535,11 @@ subroutine ocean_tpm_source(isd, ied, jsd, jed, Domain, Grid, T_prog, T_diag, type(ocean_density_type), intent(in) :: Dens real, intent(in), dimension(isd:,jsd:) :: hblt_depth real, intent(in) :: dtts + +! Optional input for csiro bgc +real, intent(in), dimension(isd:,jsd:), optional :: swflx ! short wave radiation flux (W/m^2) +real, intent(in), dimension(isd:,jsd:,:), optional :: sw_frac_zt ! short wave radiation fraction + ! !----------------------------------------------------------------------- ! local parameters @@ -1555,6 +1598,14 @@ subroutine ocean_tpm_source(isd, ied, jsd, jed, Domain, Grid, T_prog, T_diag, #endif +#if defined(CSIRO_BGC) +if (do_csiro_bgc) then + call csiro_bgc_source(Domain%isc, Domain%iec, Domain%jsc, Domain%jec, & + Domain%isd, Domain%ied, Domain%jsd, Domain%jed, & + T_prog, grid, Time, dtts, Thickness, Dens, swflx, sw_frac_zt) +endif +#endif + if (do_ocean_residency) then !{ ! must come last call ocean_residency_source(Domain%isc, Domain%iec, Domain%jsc, Domain%jec, & isd, ied, jsd, jed, Grid%nk, T_prog, T_diag, Time, Thickness, Dens, & @@ -1689,6 +1740,12 @@ subroutine ocean_tpm_start(Domain, Grid, T_prog, T_diag, Time, Thickness) !{ #endif +#if defined(CSIRO_BGC) +if (do_csiro_bgc) then !{ + call csiro_bgc_start(Time, Domain, Grid) +endif !} +#endif + if (do_transport_matrix) then !{ call transport_matrix_start(Time, T_prog, Domain%isd, Domain%ied, Domain%jsd, & Domain%jed, Grid%nk, Grid%tracer_axes) @@ -1792,6 +1849,12 @@ subroutine ocean_tpm_tracer(Domain, T_prog, T_diag, Grid, Time, Thickness, Dens, #endif +#if defined(CSIRO_BGC) +if (do_csiro_bgc) then !{ + call csiro_bgc_tracer(Domain%isc, Domain%iec, Domain%jsc, Domain%jec, T_prog, grid, Time, dtts) +endif !} +#endif + if (do_transport_matrix) then !{ call transport_matrix_store_implicit(Time, T_prog, Domain%isd, Domain%ied, Domain%jsd, Domain%jed, & Grid%nk, Domain%isc, Domain%iec, Domain%jsc, Domain%jec, Grid%tmask) From c456d943f00e7d5e6259de202820873ddfc4ef77 Mon Sep 17 00:00:00 2001 From: Russell Fiedler Date: Fri, 24 Jan 2020 16:22:42 +1100 Subject: [PATCH 2/7] update bin templates --- bin/environs.nci | 9 ++++----- bin/mkmf.template.nci | 2 +- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/bin/environs.nci b/bin/environs.nci index be2c1138b5..f1f4e84d93 100644 --- a/bin/environs.nci +++ b/bin/environs.nci @@ -1,7 +1,6 @@ -source /etc/profile.d/nf_csh_modules +source /etc/profile.d/modules.sh module purge -module load intel-fc/17.0.1.132 -module load intel-cc/17.0.1.132 -module load netcdf/4.3.2 -module load openmpi/1.10.2 +module load intel-compiler/2020.0.166 +module load netcdf/4.7.1 +module load openmpi/4.0.2 setenv mpirunCommand "mpirun --mca orte_base_help_aggregate 0 -np" diff --git a/bin/mkmf.template.nci b/bin/mkmf.template.nci index a9a1e7a614..031a41fcf5 100644 --- a/bin/mkmf.template.nci +++ b/bin/mkmf.template.nci @@ -38,7 +38,7 @@ endif FPPFLAGS := -fpp -Wp,-w $(INCLUDE) FFLAGS := -fno-alias -safe-cray-ptr -fpe0 -ftz -assume byterecl -i4 -r8 -traceback -nowarn -check noarg_temp_created -assume nobuffered_io -convert big_endian -grecord-gcc-switches -align all -FFLAGS_OPT = -O2 -debug minimal -xHost +FFLAGS_OPT = g -O2 -debug minimal -xHost FFLAGS_DEBUG = -g -O0 -debug all -check -check noarg_temp_created -check nopointer -warn -warn noerrors -ftrapuv FFLAGS_REPRO = -fp-model precise FFLAGS_VERBOSE = -v -V -what From ed589ae40f4717ab8032444f15511219e4d67ae1 Mon Sep 17 00:00:00 2001 From: Russell Fiedler Date: Mon, 27 Apr 2020 14:02:08 +1000 Subject: [PATCH 3/7] co2 only in accesscm config --- src/mom5/ocean_core/ocean_sbc.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mom5/ocean_core/ocean_sbc.F90 b/src/mom5/ocean_core/ocean_sbc.F90 index fb68c473f2..d6f200df8b 100644 --- a/src/mom5/ocean_core/ocean_sbc.F90 +++ b/src/mom5/ocean_core/ocean_sbc.F90 @@ -1995,10 +1995,12 @@ subroutine ocean_sbc_diag_init(Time, Dens, T_prog) id_total_ocean_wfiform = register_diag_field('ocean_model','total_ocean_wfiform', & Time%model_time, 'total iceform outof ocean (>0 enters ocean)', & 'kg/sec/1e15', missing_value=missing_value,range=(/-1.e10,1.e10/)) +#if defined(ACCESS_CM) id_atm_co2 = register_diag_field('ocean_model','atm_co2', Grd%tracer_axes(1:2),& Time%model_time, 'Atmospheric CO2 content', 'ppm' , & missing_value=missing_value,range=(/-1.e1,1.e4/), & standard_name='atmospheric_co2' ) +#endif id_total_ocean_licefw = register_diag_field('ocean_model','total_ocean_licefw', & Time%model_time, 'total land icemelt into ocean (>0 enters ocean)', & 'kg/sec/1e15', missing_value=missing_value,range=(/-1.e10,1.e10/)) From 3e5e3bfd4a2b6d72133404f032a5ca88c83bebad Mon Sep 17 00:00:00 2001 From: Russell Fiedler Date: Thu, 30 Apr 2020 12:52:21 +1000 Subject: [PATCH 4/7] atm_co2 only used in ACCESS_CM --- src/mom5/ocean_core/ocean_sbc.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mom5/ocean_core/ocean_sbc.F90 b/src/mom5/ocean_core/ocean_sbc.F90 index d6f200df8b..28021783d1 100644 --- a/src/mom5/ocean_core/ocean_sbc.F90 +++ b/src/mom5/ocean_core/ocean_sbc.F90 @@ -5977,10 +5977,12 @@ subroutine ocean_sbc_diag(Time, Velocity, Thickness, Dens, T_prog, Ice_ocean_bou Time%model_time, rmask=Grd%tmask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) endif +#if defined(ACCESS_CM) if (id_atm_co2 > 0) used = send_data(id_atm_co2, atm_co2(:,:), & Time%model_time, rmask=Grd%tmask(:,:,1), & is_in=isc, js_in=jsc, ie_in=iec, je_in=jec) #endif +#endif ! output saline contraction coeff (1/rho drho/dS) at the ocean surface (1/psu) From 2c877109be97828b37ef7bf6aaa28e42ef6c2512 Mon Sep 17 00:00:00 2001 From: Russell Fiedler Date: Mon, 25 May 2020 11:25:42 +1000 Subject: [PATCH 5/7] typo fix in compilation script --- exp/MOM_compile.csh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exp/MOM_compile.csh b/exp/MOM_compile.csh index 48aa291178..b63526e668 100755 --- a/exp/MOM_compile.csh +++ b/exp/MOM_compile.csh @@ -57,7 +57,7 @@ if ( $help ) then echo " ACCESS-ESM : ocean component of ACCESS-ESM model with CSIRO BGC (Wombat)." echo " ACCESS-OM-BGC: ocean component of ACCESS-OM model with CSIRO BGC (Wombat)." echo - echo "--platform followed by the platform name that has a corresponfing environ file in the ../bin dir, default is gfortran" + echo "--platform followed by the platform name that has a corresponding environ file in the ../bin dir, default is gfortran" echo echo "--use_netcdf4 use NetCDF4, the default is NetCDF4. Warning: many of the standard experiments don't work with NetCDF4." echo From 82321af302eb5084c2b091e794f20712bfde645f Mon Sep 17 00:00:00 2001 From: Russell Fiedler Date: Mon, 25 May 2020 12:53:53 +1000 Subject: [PATCH 6/7] add travis testing for wombat --- .travis.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.travis.yml b/.travis.yml index 51136f4d80..f0348fb5f2 100644 --- a/.travis.yml +++ b/.travis.yml @@ -19,6 +19,8 @@ env: - TYPE=EBM - TYPE=ACCESS-OM - TYPE=ACCESS-CM + - TYPE=ACCESS-OM-BGC + - TYPE=ACCESS-ESM global: - FC=gfortran-4.8 - OMPI_FC=${FC} From 1aea40c429af6fc75ff3ceae5012d8e39eee0200 Mon Sep 17 00:00:00 2001 From: Aidan Heerdegen Date: Tue, 26 May 2020 15:21:58 +1000 Subject: [PATCH 7/7] Added ACCESS-CM-BGC to allow_failures --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index f0348fb5f2..aeefe48761 100644 --- a/.travis.yml +++ b/.travis.yml @@ -29,7 +29,7 @@ env: matrix: allow_failures: - env: TYPE=ACCESS-OM - # - env: TYPE=ACCESS-CM + - env: TYPE=ACCESS-CM-BGC install: - sudo apt-add-repository --yes ppa:ubuntu-toolchain-r/test