!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module ecosys_mod 2,31
!BOP
! !MODULE: ecosys_mod
!
! !DESCRIPTION:
!
! Multispecies ecosystem based on Doney et al. 1996, Moore et al., 2002
! Based on POP Global NCAR Nitrogen Ecosystem Model
! version 0.0 (June 15th, 1998) from S.C. Doney.
! Based on Doney et al., 1996 model.
! Climate and Global Dynamics, NCAR
! (doney@whoi.edu)
!
! Version 1.0
! Multispecies, multiple limiting nutrient version of ecosystem
! based on mixed layer model of Moore et al.(2002). Implemented here with
! fixed elemental ratios and including only the diatoms and small
! phytoplankton, with a parameterization of calcification,
! by Keith Lindsay and Keith Moore, Fall 2001 - Spring 2002.
! Calcification parameterization based on Moore et al. 2002.
!
! Version 2.0, January 2003
! Adds diazotrophs as a phytoplankton group, (based on Moore et al., 2002a)
! Allows for variable fe/C for all phytoplankton groups
! Allows for variable si/C for the diatoms
! Adds explicit tracers for DON, DOP, DOFe
! variable remin length scale for detrital soft POM and bSi f(temperature)
! Extensive modifications to iron scavenging parameterization
! Addition of a sedimentary dissolved iron source,
! (implemented in ballast code as excess remin in bottom cell)
! coded by J.K. Moore, (jkmoore@uci.edu)
!
! Version 2.01. March 2003
! corrected O2 bug
! corrected grazing parameter z_grz bug at depth
! dust dissolution at depth releases iron,
! increased length scale for dust diss., increased hard fraction dust
! no deep ocean reduction in scavenging rates,
! increase bSi OC/ballast ratio 0.3 -> 0.35,
! corrected bug in diazotroph photoadaptation, and diat and sp adapatation
!
! Version 2.02.
! corrected bug in Fe_scavenge (units for dust), May 2003
! changed C/N/P ratios to 117/16/1 (Anderson & Sarmiento, 1994)
!
! Version 2.03., July 2003
! Remin of DOM no longer temperature dependent,
! new iron scavenging parameterization added,
! some dissolution of hard fraction of ballast materials added
!
! Version 2.1, September 2003
! modfied iron scavenging and dust dissolution at depth
!
! Version 2.11, March 2004
! fixed bug in iron scavenging code, replace dust and POC flux_in w/ flux_out
!
! Version 2.12, April 2004 - Final version for GBC paper revision,
! (Questions/comments, Keith Moore - jkmoore@uci.edu
!
! References
! Doney, S.C., Glover, D.M., Najjar, R.G., 1996. A new coupled, one-dimensional
! biological-physical model for the upper ocean: applications to the JGOFS
! Bermuda Time-Series Study (BATS) site. Deep-Sea Res. II, 43: 591-624.
!
! Moore, JK, Doney, SC, Kleypas, JA, Glover, DM, Fung, IY, 2002. An intermediate
! complexity marine ecosystem model for the global domain. Deep-Sea Res. II, 49:
! 403-462.
!
! Moore, JK, Doney, SC, Glover, DM, Fung, IY, 2002. Iron cycling and nutrient
! limitation patterns in surface waters of the world ocean. Deep-Sea Res. II,
! 49: 463-507.
! !REVISION HISTORY:
! SVN:$Id: $
!-----------------------------------------------------------------------
! variables/subroutines/function used from other modules
! The following are used extensively in this ecosys, so are used at
! the module level. The use statements for variables that are only needed
! locally are located at the module subprogram level.
!-----------------------------------------------------------------------
! !USES:
use POP_KindsMod
use POP_ErrorMod
use POP_CommMod
use POP_GridHorzMod
use POP_FieldMod
use POP_HaloMod
use kinds_mod
use constants
use communicate
use broadcast
use global_reductions
use blocks
use domain_size
use domain
use exit_mod
use prognostic
use grid
use io
use io_types
use io_tools
use tavg
use timers
use passive_tracer_tools
use named_field_mod
use forcing_tools
use time_management
use ecosys_parms
use registry
use named_field_mod
use co2calc
#ifdef CCSMCOUPLED
use POP_MCT_vars_mod
use shr_strdata_mod
#endif
! !INPUT PARAMETERS:
!-----------------------------------------------------------------------
! include ecosystem parameters
! all variables from this modules have a parm_ prefix
!-----------------------------------------------------------------------
implicit none
save
private
!-----------------------------------------------------------------------
! public/private declarations
!-----------------------------------------------------------------------
public :: &
ecosys_tracer_cnt, &
ecosys_init, &
ecosys_tracer_ref_val, &
ecosys_set_sflux, &
ecosys_tavg_forcing, &
ecosys_set_interior, &
ecosys_write_restart, &
ecosys_qsw_distrb_const
!-----------------------------------------------------------------------
! module variables required by forcing_passive_tracer
!-----------------------------------------------------------------------
integer (int_kind), parameter :: &
ecosys_tracer_cnt = 24
!-----------------------------------------------------------------------
! flags controlling which portion of code are executed
! usefull for debugging
!-----------------------------------------------------------------------
logical (log_kind) :: &
lsource_sink, &
lflux_gas_o2, &
lflux_gas_co2,&
locmip_k1_k2_bug_fix
logical (log_kind), dimension(:,:,:), allocatable :: &
LAND_MASK
!-----------------------------------------------------------------------
! relative tracer indices
!-----------------------------------------------------------------------
integer (int_kind), parameter :: &
po4_ind = 1, & ! dissolved inorganic phosphate
no3_ind = 2, & ! dissolved inorganic nitrate
sio3_ind = 3, & ! dissolved inorganic silicate
nh4_ind = 4, & ! dissolved ammonia
fe_ind = 5, & ! dissolved inorganic iron
o2_ind = 6, & ! dissolved oxygen
dic_ind = 7, & ! dissolved inorganic carbon
alk_ind = 8, & ! alkalinity
doc_ind = 9, & ! dissolved organic carbon
spC_ind = 10, & ! small phytoplankton carbon
spChl_ind = 11, & ! small phytoplankton chlorophyll
spCaCO3_ind = 12, & ! small phytoplankton caco3
diatC_ind = 13, & ! diatom carbon
diatChl_ind = 14, & ! diatom chlorophyll
zooC_ind = 15, & ! zooplankton carbon
spFe_ind = 16, & ! small phytoplankton iron
diatSi_ind = 17, & ! diatom silicon
diatFe_ind = 18, & ! diatom iron
diazC_ind = 19, & ! diazotroph carbon
diazChl_ind = 20, & ! diazotroph Chlorophyll
diazFe_ind = 21, & ! diazotroph iron
don_ind = 22, & ! dissolved organic nitrogen
dofe_ind = 23, & ! dissolved organic iron
dop_ind = 24 ! dissolved organic phosphorus
!-----------------------------------------------------------------------
! derived type & parameter for tracer index lookup
!-----------------------------------------------------------------------
type(ind_name_pair), dimension(ecosys_tracer_cnt) :: &
ind_name_table
! type(ind_name_pair), dimension(ecosys_tracer_cnt) :: &
! ind_name_table = (/ &
! ind_name_pair(po4_ind, 'PO4'), &
! ind_name_pair(no3_ind, 'NO3'), &
! ind_name_pair(sio3_ind, 'SiO3'), &
! ind_name_pair(nh4_ind, 'NH4'), &
! ind_name_pair(fe_ind, 'Fe'), &
! ind_name_pair(o2_ind, 'O2'), &
! ind_name_pair(dic_ind, 'DIC'), &
! ind_name_pair(alk_ind, 'ALK'), &
! ind_name_pair(doc_ind, 'DOC'), &
! ind_name_pair(spC_ind, 'spC'), &
! ind_name_pair(spChl_ind, 'spChl'), &
! ind_name_pair(spCaCO3_ind, 'spCaCO3'), &
! ind_name_pair(diatC_ind, 'diatC'), &
! ind_name_pair(diatChl_ind, 'diatChl'), &
! ind_name_pair(zooC_ind, 'zooC'), &
! ind_name_pair(spFe_ind, 'spFe'), &
! ind_name_pair(diatSi_ind, 'diatSi'), &
! ind_name_pair(diatFe_ind, 'diatFe'), &
! ind_name_pair(diazC_ind, 'diazC'), &
! ind_name_pair(diazChl_ind, 'diazChl'), &
! ind_name_pair(diazFe_ind, 'diazFe'), &
! ind_name_pair(don_ind, 'DON'), &
! ind_name_pair(dofe_ind, 'DOFe'), &
! ind_name_pair(dop_ind, 'DOP') /)
!
!-----------------------------------------------------------------------
! options for forcing of gas fluxes
!-----------------------------------------------------------------------
integer (int_kind), parameter :: &
gas_flux_forcing_iopt_drv = 1, &
gas_flux_forcing_iopt_file = 2, &
atm_co2_iopt_const = 1, &
atm_co2_iopt_drv_prog = 2, &
atm_co2_iopt_drv_diag = 3
integer (int_kind) :: &
gas_flux_forcing_iopt, &
atm_co2_iopt
real (r8) :: &
atm_co2_const ! value of atmospheric co2 (ppm, dry-air, 1 atm)
character(char_len) :: &
gas_flux_forcing_file ! file containing gas flux forcing fields
!-----------------------------------------------------------------------
real (r8) :: &
rest_time_inv_surf, & ! inverse restoring timescale at surface
rest_time_inv_deep, & ! inverse restoring timescale at depth
rest_z0, & ! shallow end of transition regime
rest_z1 ! deep end of transition regime
type(tracer_read) :: &
po4_rest, & ! restoring data for PO4
no3_rest, & ! restoring data for NO3
sio3_rest, & ! restoring data for SiO3
gas_flux_fice, & ! ice fraction for gas fluxes
gas_flux_ws, & ! wind speed for gas fluxes
gas_flux_ap, & ! atmospheric pressure for gas fluxes
fesedflux_input ! namelist input for iron_flux
!-----------------------------------------------------------------------
! module variables related to ph computations
!-----------------------------------------------------------------------
real (r8), dimension(:,:,:), allocatable, target :: &
PH_PREV, & ! computed ph from previous time step
IRON_PATCH_FLUX ! localized iron patch flux
real (r8), dimension(:,:,:), allocatable :: &
dust_FLUX_IN ! dust flux not stored in STF since dust is not prognostic
real (r8), dimension(:,:,:,:), allocatable :: &
PH_PREV_3D ! computed pH_3D from previous time step
!-----------------------------------------------------------------------
! restoring climatologies for nutrients
!-----------------------------------------------------------------------
logical (log_kind) :: &
lrest_po4, & ! restoring on po4 ?
lrest_no3, & ! restoring on no3 ?
lrest_sio3 ! restoring on sio3 ?
real (r8), dimension(km) :: &
nutr_rest_time_inv ! inverse restoring time scale for nutrients (1/secs)
real (r8), dimension(:,:,:,:), allocatable, target :: &
PO4_CLIM, NO3_CLIM, SiO3_CLIM
real (r8), dimension(:,:,:,:), allocatable, target :: &
FESEDFLUX
character(char_len) :: &
nutr_rest_file ! file containing nutrient fields
!maltrud variable restoring
logical (log_kind) :: &
lnutr_variable_restore ! geographically varying nutrient restoring
character(char_len) :: &
nutr_variable_rest_file, & ! file containing variable restoring info
nutr_variable_rest_file_fmt ! format of file containing variable restoring info
real (r8), dimension(:,:,:), allocatable, target :: &
NUTR_RESTORE_RTAU ! inverse restoring timescale for variable
! interior restoring
integer (int_kind), dimension(:,:,:), allocatable :: &
NUTR_RESTORE_MAX_LEVEL ! maximum level for applying variable
! interior restoring
real (r8), dimension(:,:,:,:), allocatable :: &
INTERP_WORK ! temp array for interpolate_forcing output
type(forcing_monthly_every_ts) :: &
dust_flux, & ! surface dust flux
iron_flux, & ! iron component of surface dust flux
fice_file, & ! ice fraction, if read from file
xkw_file, & ! a * wind-speed ** 2, if read from file
ap_file ! atmoshperic pressure, if read from file
character(char_len) :: &
ndep_data_type ! type of ndep forcing
type(forcing_monthly_every_ts) :: &
nox_flux_monthly, & ! surface NOx species flux, added to nitrate pool
nhy_flux_monthly ! surface NHy species flux, added to ammonium pool
integer (int_kind) :: &
ndep_shr_stream_year_first, & ! first year in stream to use
ndep_shr_stream_year_last, & ! last year in stream to use
ndep_shr_stream_year_align ! align ndep_shr_stream_year_first with this model year
integer (int_kind), parameter :: &
ndep_shr_stream_var_cnt = 2, & ! number of variables in ndep shr_stream
ndep_shr_stream_no_ind = 1, & ! index for NO forcing
ndep_shr_stream_nh_ind = 2 ! index for NH forcing
character(char_len) :: &
ndep_shr_stream_file ! file containing domain and input data
real (r8) :: &
ndep_shr_stream_scale_factor ! unit conversion factor
#ifdef CCSMCOUPLED
type(shr_strdata_type) :: ndep_sdat ! input data stream for ndep
#endif
!-----------------------------------------------------------------------
! derived type for implicit handling of sinking particulate matter
!-----------------------------------------------------------------------
type sinking_particle
real (r8) :: &
diss, & ! dissolution length for soft subclass
gamma, & ! fraction of production -> hard subclass
mass, & ! mass of 1e9 base units in g
rho ! QA mass ratio of POC to this particle class
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
sflux_in, & ! incoming flux of soft subclass (base units/cm^2/sec)
hflux_in, & ! incoming flux of hard subclass (base units/cm^2/sec)
prod, & ! production term (base units/cm^3/sec)
sflux_out, & ! outgoing flux of soft subclass (base units/cm^2/sec)
hflux_out, & ! outgoing flux of hard subclass (base units/cm^2/sec)
remin ! remineralization term (base units/cm^3/sec)
end type sinking_particle
!-----------------------------------------------------------------------
! define tavg id for 2d fields related to surface fluxes
!-----------------------------------------------------------------------
integer (int_kind) :: &
tavg_ECOSYS_IFRAC, &! tavg id for ice fraction
tavg_ECOSYS_IFRAC_2,&! tavg id for duplicate of ice fraction
tavg_ECOSYS_XKW, &! tavg id for xkw
tavg_ECOSYS_XKW_2, &! tavg id for duplicate of xkw
tavg_ECOSYS_ATM_PRESS, &! tavg id for atmospheric pressure
tavg_PV_O2, &! tavg id for o2 piston velocity
tavg_SCHMIDT_O2, &! tavg id for O2 schmidt number
tavg_O2SAT, &! tavg id for O2 saturation
tavg_CO2STAR, &! tavg id for co2star
tavg_DCO2STAR, &! tavg id for dco2star
tavg_pCO2SURF, &! tavg id for surface pco2
tavg_DpCO2, &! tavg id for delta pco2
tavg_DpCO2_2, &! tavg id for duplicate of delta pco2
tavg_PV_CO2, &! tavg id for co2 piston velocity
tavg_SCHMIDT_CO2, &! tavg id for co2 schmidt number
tavg_DIC_GAS_FLUX, &! tavg id for dic flux
tavg_DIC_GAS_FLUX_2,&! tavg id for duplicate of dic flux
tavg_O2_GAS_FLUX_2,&! tavg id for duplicate of O2 flux
tavg_PH, &! tavg id for surface pH
tavg_ATM_CO2, &! tavg id for atmospheric CO2
tavg_IRON_FLUX, &! tavg id for dust flux
tavg_DUST_FLUX, &! tavg id for dust flux
tavg_NOx_FLUX, &! tavg id for nox flux
tavg_NHy_FLUX ! tavg id for nhy flux
!-----------------------------------------------------------------------
! define tavg id for nonstandard 2d fields
!-----------------------------------------------------------------------
integer (int_kind) :: &
tavg_O2_ZMIN, &! tavg id for vertical minimum of O2
tavg_O2_ZMIN_DEPTH ! tavg id for depth of vertical minimum of O2
!-----------------------------------------------------------------------
! define tavg id for nonstandard 3d fields
!-----------------------------------------------------------------------
integer (int_kind) :: &
tavg_O2_PRODUCTION,&! tavg id for o2 production
tavg_O2_CONSUMPTION,&! tavg id for o2 consumption
tavg_AOU, &! tavg id for AOU
tavg_PO4_RESTORE, &! tavg id for po4 restoring
tavg_NO3_RESTORE, &! tavg id for no3 restoring
tavg_SiO3_RESTORE, &! tavg id for sio3 restoring
tavg_PAR_avg, &! tavg id for available radiation avg over mixed layer
tavg_POC_FLUX_IN, &! tavg id for poc flux into cell
tavg_POC_PROD, &! tavg id for poc production
tavg_POC_REMIN, &! tavg id for poc remineralization
tavg_POC_ACCUM, &! tavg id for poc accumulation
tavg_CaCO3_FLUX_IN,&! tavg id for caco3 flux into cell
tavg_CaCO3_PROD, &! tavg id for caco3 production
tavg_CaCO3_REMIN, &! tavg id for caco3 remineralization
tavg_SiO2_FLUX_IN, &! tavg id for sio2 flux into cell
tavg_SiO2_PROD, &! tavg id for sio2 production
tavg_SiO2_REMIN, &! tavg id for sio2 remineralization
tavg_dust_FLUX_IN, &! tavg id for dust flux into cell
tavg_dust_REMIN, &! tavg id for dust remineralization
tavg_P_iron_FLUX_IN, &! tavg id for p_iron flux into cell
tavg_P_iron_PROD, &! tavg id for p_iron production
tavg_P_iron_REMIN, &! tavg id for p_iron remineralization
tavg_graze_sp, &! tavg id for small phyto grazing
tavg_graze_diat, &! tavg id for diatom grazing
tavg_graze_diaz, &! tavg id for diaz grazing
tavg_graze_TOT, &! tavg id for total grazing
tavg_sp_loss, &! tavg id for small phyto loss
tavg_diat_loss, &! tavg id for diatom loss
tavg_diaz_loss, &! tavg id for diaz loss
tavg_zoo_loss, &! tavg id for zooplankton loss
tavg_sp_agg, &! tavg id for small phyto aggregate
tavg_diat_agg ! tavg id for diatom aggregate
!-----------------------------------------------------------------------
! define tavg id for MORE nonstandard 3d fields
!-----------------------------------------------------------------------
integer (int_kind) :: &
tavg_photoC_sp, &! tavg id for small phyto C fixation
tavg_CaCO3_form, &! tavg id for CaCO3 formation
tavg_photoC_diat, &! tavg id for diatom C fixation
tavg_photoC_diaz, &! tavg id for diaz C fixation
tavg_photoC_TOT, &! tavg id for total C fixation
tavg_photoC_sp_zint, &! tavg id for small phyto C fixation vertical integral
tavg_CaCO3_form_zint, &! tavg id for CaCO3 formation vertical integral
tavg_photoC_diat_zint, &! tavg id for diatom C fixation vertical integral
tavg_photoC_diaz_zint, &! tavg id for diaz C fixation vertical integral
tavg_photoC_TOT_zint, &! tavg id for total C fixation vertical integral
tavg_photoC_NO3_sp, &! tavg id for small phyto C fixation from NO3
tavg_photoC_NO3_sp_zint, &! tavg id for small phyto C fixation from NO3 vertical integral
tavg_photoC_NO3_diat, &! tavg id for diatom C fixation from NO3
tavg_photoC_NO3_diat_zint, &! tavg id for diatom C fixation from NO3 vertical integral
tavg_photoC_NO3_diaz, &! tavg id for diaz C fixation from NO3
tavg_photoC_NO3_diaz_zint, &! tavg id for diaz C fixation from NO3 vertical integral
tavg_photoC_NO3_TOT, &! tavg id for total C fixation from NO3
tavg_photoC_NO3_TOT_zint ! tavg id for total C fixation from NO3 vertical integral
!-----------------------------------------------------------------------
! define tavg id for MORE nonstandard 3d fields
!-----------------------------------------------------------------------
integer (int_kind) :: &
tavg_DOC_prod, &! tavg id for doc production
tavg_DOC_remin, &! tavg id for doc remineralization
tavg_DON_prod, &! tavg id for don production
tavg_DON_remin, &! tavg id for don remineralization
tavg_DOFe_prod, &! tavg id for dofe production
tavg_DOFe_remin, &! tavg id for dofe remineralization
tavg_DOP_prod, &! tavg id for dop production
tavg_DOP_remin, &! tavg id for dop remineralization
tavg_Fe_scavenge, &! tavg id for iron scavenging
tavg_Fe_scavenge_rate, &! tavg id for iron scavenging rate
tavg_sp_N_lim, &! tavg id for small phyto N limitation
tavg_sp_Fe_lim, &! tavg id for small phyto Fe limitation
tavg_sp_PO4_lim, &! tavg id for small phyto PO4 limitation
tavg_sp_light_lim, &! tavg id for small phyto light limitation
tavg_diat_N_lim, &! tavg id for diatom N limitation
tavg_diat_Fe_lim, &! tavg id for diatom Fe limitation
tavg_diat_PO4_lim, &! tavg id for diatom PO4 limitation
tavg_diat_SiO3_lim, &! tavg id for diatom SiO3 limitation
tavg_diat_light_lim, &! tavg id for diatom light limitation
tavg_diaz_Fe_lim, &! tavg id for diaz Fe limitation
tavg_diaz_P_lim, &! tavg id for diaz PO4 limitation
tavg_diaz_light_lim, &! tavg id for diaz light limitation
tavg_diaz_Nfix, &! tavg id for diaz N fixation
tavg_bSi_form, &! tavg id for diatom Si uptake
tavg_NITRIF, &! tavg id for nitrification
tavg_DENITRIF, &! tavg id for denitrification
tavg_photoFe_sp, &! tavg id for small phyto Fe uptake
tavg_photoFe_diat, &! tavg id for diatom Fe uptake
tavg_photoFe_diaz ! tavg id for diaz Fe uptake
integer (int_kind) :: &
tavg_CO3, &! tavg id for 3D carbonate ion
tavg_HCO3, &! tavg id for 3D bicarbonate ion
tavg_H2CO3, &! tavg id for 3D carbonic acid
tavg_pH_3D, &! tavg id for 3D pH
tavg_co3_sat_calc, &! tavg id for co3 concentration at calcite saturation
tavg_zsatcalc, &! tavg id for calcite saturation depth
tavg_co3_sat_arag, &! tavg id for co3 concentration at aragonite saturation
tavg_zsatarag ! tavg id for aragonite saturation depth
!-----------------------------------------------------------------------
! define array for holding flux-related quantities that need to be time-averaged
! this is necessary since the forcing routines are called before tavg flags
!-----------------------------------------------------------------------
real (r8), dimension(:,:,:,:), allocatable :: &
ECO_SFLUX_TAVG
!-----------------------------------------------------------------------
! average surface tracer value related variables
! used as reference value for virtual flux computations
!-----------------------------------------------------------------------
logical (log_kind), dimension(ecosys_tracer_cnt) :: &
vflux_flag ! which tracers get virtual fluxes applied
integer (int_kind) :: &
comp_surf_avg_flag ! time flag id for computing average
! surface tracer values
real (r8), dimension(ecosys_tracer_cnt) :: &
surf_avg ! average surface tracer values
logical (log_kind) :: &
ecosys_qsw_distrb_const
!-----------------------------------------------------------------------
! iron patch fertilization
!-----------------------------------------------------------------------
logical (log_kind) :: &
liron_patch ! flag for iron patch fertilization
character(char_len) :: &
iron_patch_flux_filename ! file containing name of iron patch file
integer (int_kind) :: &
iron_patch_month ! integer month to add patch flux
!-----------------------------------------------------------------------
! timers
!-----------------------------------------------------------------------
integer (int_kind) :: &
ecosys_shr_strdata_advance_timer, &
ecosys_comp_CO3terms_timer, &
ecosys_interior_timer, &
ecosys_sflux_timer
!-----------------------------------------------------------------------
! named field indices
!-----------------------------------------------------------------------
integer (int_kind) :: &
totChl_surf_nf_ind = 0, & ! total chlorophyll in surface layer
sflux_co2_nf_ind = 0, & ! air-sea co2 gas flux
atm_co2_nf_ind = 0 ! atmospheric co2
!-----------------------------------------------------------------------
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
PAR_out ! photosynthetically available radiation (W/m^2)
!-----------------------------------------------------------------------
real (r8), parameter :: &
phlo_init = 7.0_r8, & ! low bound for ph for no prev soln
phhi_init = 9.0_r8, & ! high bound for ph for no prev soln
del_ph = 0.20_r8 ! delta-ph for prev soln
!*****************************************************************************
contains
!*****************************************************************************
!BOP
! !IROUTINE: ecosys_init
! !INTERFACE:
subroutine ecosys_init(init_ts_file_fmt, read_restart_filename, & 1,148
tracer_d_module, TRACER_MODULE, tadvect_ctype, &
errorCode)
! !DESCRIPTION:
! Initialize ecosys tracer module. This involves setting metadata, reading
! the module namelist, setting initial conditions, setting up forcing,
! and defining additional tavg variables.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
character (*), intent(in) :: &
init_ts_file_fmt, & ! format (bin or nc) for input file
read_restart_filename ! file name for restart file
! !INPUT/OUTPUT PARAMETERS:
type (tracer_field), dimension(ecosys_tracer_cnt), intent(inout) :: &
tracer_d_module ! descriptors for each tracer
real (r8), dimension(nx_block,ny_block,km,ecosys_tracer_cnt,3,max_blocks_clinic), &
intent(inout) :: TRACER_MODULE
! !OUTPUT PARAMETERS:
character (char_len), dimension(ecosys_tracer_cnt), intent(out) :: &
tadvect_ctype ! advection method for ecosys tracers
integer (POP_i4), intent(out) :: &
errorCode
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
character(*), parameter :: subname = 'ecosys_mod:ecosys_init'
character(char_len) :: &
init_ecosys_option, & ! option for initialization of bgc
init_ecosys_init_file, & ! filename for option 'file'
init_ecosys_init_file_fmt, & ! file format for option 'file'
comp_surf_avg_freq_opt, & ! choice for freq of comp_surf_avg
gas_flux_forcing_opt, & ! option for forcing gas fluxes
atm_co2_opt, & ! option for atmospheric co2 concentration
ecosys_tadvect_ctype ! advection method for ecosys tracers
type(tracer_read), dimension(ecosys_tracer_cnt) :: &
tracer_init_ext ! namelist variable for initializing tracers
type(tracer_read) :: &
dust_flux_input, & ! namelist input for dust_flux
iron_flux_input, & ! namelist input for iron_flux
nox_flux_monthly_input, & ! namelist input for nox_flux_monthly
nhy_flux_monthly_input ! namelist input for nhy_flux_monthly
logical (log_kind) :: &
default, & ! arg to init_time_flag
lnml_found, & ! Was ecosys_nml found ?
lmarginal_seas ! Is ecosystem active in marginal seas ?
integer (int_kind) :: &
n, & ! index for looping over tracers
k, & ! index for looping over depth levels
l, & ! index for looping over time levels
ind, & ! tracer index for tracer name from namelist
iblock, & ! index for looping over blocks
nml_error ! namelist i/o error flag
integer (int_kind) :: &
freq_opt, freq, & ! args for init_time_flag
comp_surf_avg_freq_iopt, & ! choice for freq of comp_surf_avg
comp_surf_avg_freq ! choice for freq of comp_surf_avg
logical (log_kind) :: &
use_nml_surf_vals ! do namelist surf values override values from restart file
logical (log_kind) :: &
lecovars_full_depth_tavg ! should ecosystem vars be written full depth
!-----------------------------------------------------------------------
! values to be used when comp_surf_avg_freq_opt==never
!-----------------------------------------------------------------------
real (r8) :: &
surf_avg_dic_const, surf_avg_alk_const
namelist /ecosys_nml/ &
init_ecosys_option, init_ecosys_init_file, tracer_init_ext, &
init_ecosys_init_file_fmt, &
dust_flux_input, iron_flux_input, fesedflux_input, &
ndep_data_type, nox_flux_monthly_input, nhy_flux_monthly_input, &
ndep_shr_stream_year_first, ndep_shr_stream_year_last, &
ndep_shr_stream_year_align, ndep_shr_stream_file, &
ndep_shr_stream_scale_factor, &
gas_flux_forcing_opt, gas_flux_forcing_file, &
gas_flux_fice, gas_flux_ws, gas_flux_ap, &
lrest_po4, lrest_no3, lrest_sio3, &
rest_time_inv_surf, rest_time_inv_deep, rest_z0, rest_z1, &
nutr_rest_file, po4_rest, no3_rest, sio3_rest, &
comp_surf_avg_freq_opt, comp_surf_avg_freq, &
use_nml_surf_vals, surf_avg_dic_const, surf_avg_alk_const, &
ecosys_qsw_distrb_const, lmarginal_seas, &
lsource_sink, lflux_gas_o2, lflux_gas_co2, locmip_k1_k2_bug_fix, &
lnutr_variable_restore, nutr_variable_rest_file, &
nutr_variable_rest_file_fmt,atm_co2_opt,atm_co2_const, &
ecosys_tadvect_ctype, &
liron_patch,iron_patch_flux_filename,iron_patch_month, &
lecovars_full_depth_tavg
character (char_len) :: &
ecosys_restart_filename ! modified file name for restart file
real (r8), dimension (nx_block,ny_block) :: WORK
!-----------------------------------------------------------------------
! initialize name table
!-----------------------------------------------------------------------
errorCode = POP_Success
ind_name_table( 1) = ind_name_pair(po4_ind, 'PO4')
ind_name_table( 2) = ind_name_pair(no3_ind, 'NO3')
ind_name_table( 3) = ind_name_pair(sio3_ind, 'SiO3')
ind_name_table( 4) = ind_name_pair(nh4_ind, 'NH4')
ind_name_table( 5) = ind_name_pair(fe_ind, 'Fe')
ind_name_table( 6) = ind_name_pair(o2_ind, 'O2')
ind_name_table( 7) = ind_name_pair(dic_ind, 'DIC')
ind_name_table( 8) = ind_name_pair(alk_ind, 'ALK')
ind_name_table( 9) = ind_name_pair(doc_ind, 'DOC')
ind_name_table(10) = ind_name_pair(spC_ind, 'spC')
ind_name_table(11) = ind_name_pair(spChl_ind, 'spChl')
ind_name_table(12) = ind_name_pair(spCaCO3_ind, 'spCaCO3')
ind_name_table(13) = ind_name_pair(diatC_ind, 'diatC')
ind_name_table(14) = ind_name_pair(diatChl_ind, 'diatChl')
ind_name_table(15) = ind_name_pair(zooC_ind, 'zooC')
ind_name_table(16) = ind_name_pair(spFe_ind, 'spFe')
ind_name_table(17) = ind_name_pair(diatSi_ind, 'diatSi')
ind_name_table(18) = ind_name_pair(diatFe_ind, 'diatFe')
ind_name_table(19) = ind_name_pair(diazC_ind, 'diazC')
ind_name_table(20) = ind_name_pair(diazChl_ind, 'diazChl')
ind_name_table(21) = ind_name_pair(diazFe_ind, 'diazFe')
ind_name_table(22) = ind_name_pair(don_ind, 'DON')
ind_name_table(23) = ind_name_pair(dofe_ind, 'DOFe')
ind_name_table(24) = ind_name_pair(dop_ind, 'DOP')
!-----------------------------------------------------------------------
! initialize forcing_monthly_every_ts variables
!-----------------------------------------------------------------------
call init_forcing_monthly_every_ts(dust_flux)
call init_forcing_monthly_every_ts(iron_flux)
call init_forcing_monthly_every_ts(fice_file)
call init_forcing_monthly_every_ts(xkw_file)
call init_forcing_monthly_every_ts(ap_file)
call init_forcing_monthly_every_ts(nox_flux_monthly)
call init_forcing_monthly_every_ts(nhy_flux_monthly)
!-----------------------------------------------------------------------
! initialize ecosystem parameters
!-----------------------------------------------------------------------
call ecosys_parms_init
!-----------------------------------------------------------------------
! initialize tracer_d values
!-----------------------------------------------------------------------
do n = 1, ecosys_tracer_cnt
tracer_d_module(n)%short_name = ind_name_table(n)%name
if (any(n == (/ spChl_ind, diatChl_ind, diazChl_ind /))) then
tracer_d_module(n)%units = 'mg/m^3'
tracer_d_module(n)%tend_units = 'mg/m^3/s'
tracer_d_module(n)%flux_units = 'mg/m^3 cm/s'
else if (n == alk_ind) then
tracer_d_module(n)%units = 'meq/m^3'
tracer_d_module(n)%tend_units = 'meq/m^3/s'
tracer_d_module(n)%flux_units = 'meq/m^3 cm/s'
else
tracer_d_module(n)%units = 'mmol/m^3'
tracer_d_module(n)%tend_units = 'mmol/m^3/s'
tracer_d_module(n)%flux_units = 'mmol/m^3 cm/s'
endif
end do
tracer_d_module(po4_ind )%long_name='Dissolved Inorganic Phosphate'
tracer_d_module(no3_ind )%long_name='Dissolved Inorganic Nitrate'
tracer_d_module(sio3_ind )%long_name='Dissolved Inorganic Silicate'
tracer_d_module(nh4_ind )%long_name='Dissolved Ammonia'
tracer_d_module(fe_ind )%long_name='Dissolved Inorganic Iron'
tracer_d_module(o2_ind )%long_name='Dissolved Oxygen'
tracer_d_module(dic_ind )%long_name='Dissolved Inorganic Carbon'
tracer_d_module(alk_ind )%long_name='Alkalinity'
tracer_d_module(doc_ind )%long_name='Dissolved Organic Carbon'
tracer_d_module(spC_ind )%long_name='Small Phytoplankton Carbon'
tracer_d_module(spChl_ind )%long_name='Small phytoplankton Chlorophyll'
tracer_d_module(spCaCO3_ind)%long_name='Small Phytoplankton CaCO3'
tracer_d_module(diatC_ind )%long_name='Diatom Carbon'
tracer_d_module(diatChl_ind)%long_name='Diatom Chlorophyll'
tracer_d_module(zooC_ind )%long_name='Zooplankton Carbon'
tracer_d_module(spFe_ind )%long_name='Small Phytoplankton Iron'
tracer_d_module(diatSi_ind )%long_name='Diatom Silicon'
tracer_d_module(diatFe_ind )%long_name='Diatom Iron'
tracer_d_module(diazC_ind )%long_name='Diazotroph Carbon'
tracer_d_module(diazChl_ind)%long_name='Diazotroph Chlorophyll'
tracer_d_module(diazFe_ind )%long_name='Diazotroph Iron'
tracer_d_module(don_ind )%long_name='Dissolved Organic Nitrogen'
tracer_d_module(dofe_ind )%long_name='Dissolved Organic Iron'
tracer_d_module(dop_ind )%long_name='Dissolved Organic Phosphorus'
!-----------------------------------------------------------------------
! default namelist settings
!-----------------------------------------------------------------------
init_ecosys_option = 'unknown'
init_ecosys_init_file = 'unknown'
init_ecosys_init_file_fmt = 'bin'
gas_flux_forcing_opt = 'drv'
gas_flux_forcing_file = 'unknown'
gas_flux_fice%filename = 'unknown'
gas_flux_fice%file_varname = 'FICE'
gas_flux_fice%scale_factor = c1
gas_flux_fice%default_val = c0
gas_flux_fice%file_fmt = 'bin'
gas_flux_ws%filename = 'unknown'
gas_flux_ws%file_varname = 'XKW'
gas_flux_ws%scale_factor = c1
gas_flux_ws%default_val = c0
gas_flux_ws%file_fmt = 'bin'
gas_flux_ap%filename = 'unknown'
gas_flux_ap%file_varname = 'P'
gas_flux_ap%scale_factor = c1
gas_flux_ap%default_val = c0
gas_flux_ap%file_fmt = 'bin'
lrest_po4 = .false.
lrest_no3 = .false.
lrest_sio3 = .false.
nutr_rest_file = 'unknown'
rest_time_inv_surf = c0
rest_time_inv_deep = c0
rest_z0 = c1000
rest_z1 = c2 * c1000
po4_rest%filename = 'unknown'
po4_rest%file_varname = tracer_d_module(po4_ind)%short_name
po4_rest%scale_factor = c1
po4_rest%default_val = c0
po4_rest%file_fmt = 'bin'
no3_rest%filename = 'unknown'
no3_rest%file_varname = tracer_d_module(no3_ind)%short_name
no3_rest%scale_factor = c1
no3_rest%default_val = c0
no3_rest%file_fmt = 'bin'
sio3_rest%filename = 'unknown'
sio3_rest%file_varname = tracer_d_module(sio3_ind)%short_name
sio3_rest%scale_factor = c1
sio3_rest%default_val = c0
sio3_rest%file_fmt = 'bin'
!maltrud variable restoring
lnutr_variable_restore = .false.
nutr_variable_rest_file = 'unknown'
nutr_variable_rest_file_fmt = 'bin'
dust_flux_input%filename = 'unknown'
dust_flux_input%file_varname = 'dust_flux'
dust_flux_input%scale_factor = c1
dust_flux_input%default_val = c0
dust_flux_input%file_fmt = 'bin'
iron_flux_input%filename = 'unknown'
iron_flux_input%file_varname = 'iron_flux'
iron_flux_input%scale_factor = c1
iron_flux_input%default_val = c0
iron_flux_input%file_fmt = 'bin'
fesedflux_input%filename = 'unknown'
fesedflux_input%file_varname = 'FESEDFLUXIN'
fesedflux_input%scale_factor = c1
fesedflux_input%default_val = c0
fesedflux_input%file_fmt = 'bin'
ndep_data_type = 'monthly-calendar'
nox_flux_monthly_input%filename = 'unknown'
nox_flux_monthly_input%file_varname = 'nox_flux'
nox_flux_monthly_input%scale_factor = c1
nox_flux_monthly_input%default_val = c0
nox_flux_monthly_input%file_fmt = 'bin'
nhy_flux_monthly_input%filename = 'unknown'
nhy_flux_monthly_input%file_varname = 'nhy_flux'
nhy_flux_monthly_input%scale_factor = c1
nhy_flux_monthly_input%default_val = c0
nhy_flux_monthly_input%file_fmt = 'bin'
ndep_shr_stream_year_first = 1
ndep_shr_stream_year_last = 1
ndep_shr_stream_year_align = 1
ndep_shr_stream_file = 'unknown'
ndep_shr_stream_scale_factor = c1
do n = 1,ecosys_tracer_cnt
tracer_init_ext(n)%mod_varname = 'unknown'
tracer_init_ext(n)%filename = 'unknown'
tracer_init_ext(n)%file_varname = 'unknown'
tracer_init_ext(n)%scale_factor = c1
tracer_init_ext(n)%default_val = c0
tracer_init_ext(n)%file_fmt = 'bin'
end do
lmarginal_seas = .true.
lsource_sink = .true.
lflux_gas_o2 = .true.
lflux_gas_co2 = .true.
locmip_k1_k2_bug_fix = .true.
comp_surf_avg_freq_opt = 'never'
comp_surf_avg_freq = 1
use_nml_surf_vals = .false.
surf_avg_dic_const = 1944.0_r8
surf_avg_alk_const = 2225.0_r8
ecosys_qsw_distrb_const = .true.
liron_patch = .false.
iron_patch_flux_filename = 'unknown_iron_patch_filename'
iron_patch_month = 1
atm_co2_opt = 'const'
atm_co2_const = 280.0_r8
ecosys_tadvect_ctype = 'base_model'
lecovars_full_depth_tavg = .false.
if (my_task == master_task) then
open (nml_in, file=nml_filename, status='old',iostat=nml_error)
if (nml_error /= 0) then
nml_error = -1
else
nml_error = 1
endif
do while (nml_error > 0)
read(nml_in, nml=ecosys_nml,iostat=nml_error)
end do
if (nml_error == 0) close(nml_in)
endif
call broadcast_scalar(nml_error, master_task)
if (nml_error /= 0) then
call document(subname, 'ecosys_nml not found')
call exit_POP(sigAbort, 'ERROR : stopping in '/&
&/ subname)
endif
if (my_task == master_task) then
write(stdout,blank_fmt)
write(stdout,ndelim_fmt)
write(stdout,blank_fmt)
write(stdout,*) ' ecosys:'
write(stdout,blank_fmt)
write(stdout,*) ' ecosys_nml namelist settings:'
write(stdout,blank_fmt)
write(stdout,ecosys_nml)
write(stdout,blank_fmt)
write(stdout,delim_fmt)
endif
!-----------------------------------------------------------------------
! broadcast all namelist variables
!-----------------------------------------------------------------------
call broadcast_scalar(init_ecosys_option, master_task)
call broadcast_scalar(init_ecosys_init_file, master_task)
call broadcast_scalar(init_ecosys_init_file_fmt, master_task)
call broadcast_scalar(gas_flux_forcing_opt, master_task)
if (trim(gas_flux_forcing_opt) == 'drv') then
gas_flux_forcing_iopt = gas_flux_forcing_iopt_drv
else if (trim(gas_flux_forcing_opt) == 'file') then
gas_flux_forcing_iopt = gas_flux_forcing_iopt_file
else
call document(subname, 'gas_flux_forcing_opt', gas_flux_forcing_opt)
call exit_POP(sigAbort, 'unknown gas_flux_forcing_opt')
endif
call broadcast_scalar(gas_flux_forcing_file, master_task)
call broadcast_scalar(gas_flux_fice%filename, master_task)
call broadcast_scalar(gas_flux_fice%file_varname, master_task)
call broadcast_scalar(gas_flux_fice%scale_factor, master_task)
call broadcast_scalar(gas_flux_fice%default_val, master_task)
call broadcast_scalar(gas_flux_fice%file_fmt, master_task)
fice_file%input = gas_flux_fice
call broadcast_scalar(gas_flux_ws%filename, master_task)
call broadcast_scalar(gas_flux_ws%file_varname, master_task)
call broadcast_scalar(gas_flux_ws%scale_factor, master_task)
call broadcast_scalar(gas_flux_ws%default_val, master_task)
call broadcast_scalar(gas_flux_ws%file_fmt, master_task)
xkw_file%input = gas_flux_ws
call broadcast_scalar(gas_flux_ap%filename, master_task)
call broadcast_scalar(gas_flux_ap%file_varname, master_task)
call broadcast_scalar(gas_flux_ap%scale_factor, master_task)
call broadcast_scalar(gas_flux_ap%default_val, master_task)
call broadcast_scalar(gas_flux_ap%file_fmt, master_task)
ap_file%input = gas_flux_ap
call broadcast_scalar(lrest_po4, master_task)
call broadcast_scalar(lrest_no3, master_task)
call broadcast_scalar(lrest_sio3, master_task)
call broadcast_scalar(nutr_rest_file, master_task)
call broadcast_scalar(rest_time_inv_surf, master_task)
call broadcast_scalar(rest_time_inv_deep, master_task)
call broadcast_scalar(rest_z0, master_task)
call broadcast_scalar(rest_z1, master_task)
call broadcast_scalar(po4_rest%filename, master_task)
call broadcast_scalar(po4_rest%file_varname, master_task)
call broadcast_scalar(po4_rest%scale_factor, master_task)
call broadcast_scalar(po4_rest%default_val, master_task)
call broadcast_scalar(po4_rest%file_fmt, master_task)
call broadcast_scalar(no3_rest%filename, master_task)
call broadcast_scalar(no3_rest%file_varname, master_task)
call broadcast_scalar(no3_rest%scale_factor, master_task)
call broadcast_scalar(no3_rest%default_val, master_task)
call broadcast_scalar(no3_rest%file_fmt, master_task)
call broadcast_scalar(sio3_rest%filename, master_task)
call broadcast_scalar(sio3_rest%file_varname, master_task)
call broadcast_scalar(sio3_rest%scale_factor, master_task)
call broadcast_scalar(sio3_rest%default_val, master_task)
call broadcast_scalar(sio3_rest%file_fmt, master_task)
!maltrud variable restoring
call broadcast_scalar(lnutr_variable_restore, master_task)
call broadcast_scalar(nutr_variable_rest_file, master_task)
call broadcast_scalar(nutr_variable_rest_file_fmt, master_task)
call broadcast_scalar(dust_flux_input%filename, master_task)
call broadcast_scalar(dust_flux_input%file_varname, master_task)
call broadcast_scalar(dust_flux_input%scale_factor, master_task)
call broadcast_scalar(dust_flux_input%default_val, master_task)
call broadcast_scalar(dust_flux_input%file_fmt, master_task)
dust_flux%input = dust_flux_input
call broadcast_scalar(iron_flux_input%filename, master_task)
call broadcast_scalar(iron_flux_input%file_varname, master_task)
call broadcast_scalar(iron_flux_input%scale_factor, master_task)
call broadcast_scalar(iron_flux_input%default_val, master_task)
call broadcast_scalar(iron_flux_input%file_fmt, master_task)
iron_flux%input = iron_flux_input
call broadcast_scalar(fesedflux_input%filename, master_task)
call broadcast_scalar(fesedflux_input%file_varname, master_task)
call broadcast_scalar(fesedflux_input%scale_factor, master_task)
call broadcast_scalar(fesedflux_input%default_val, master_task)
call broadcast_scalar(fesedflux_input%file_fmt, master_task)
call broadcast_scalar(ndep_data_type, master_task)
call broadcast_scalar(nox_flux_monthly_input%filename, master_task)
call broadcast_scalar(nox_flux_monthly_input%file_varname, master_task)
call broadcast_scalar(nox_flux_monthly_input%scale_factor, master_task)
call broadcast_scalar(nox_flux_monthly_input%default_val, master_task)
call broadcast_scalar(nox_flux_monthly_input%file_fmt, master_task)
nox_flux_monthly%input = nox_flux_monthly_input
call broadcast_scalar(nhy_flux_monthly_input%filename, master_task)
call broadcast_scalar(nhy_flux_monthly_input%file_varname, master_task)
call broadcast_scalar(nhy_flux_monthly_input%scale_factor, master_task)
call broadcast_scalar(nhy_flux_monthly_input%default_val, master_task)
call broadcast_scalar(nhy_flux_monthly_input%file_fmt, master_task)
nhy_flux_monthly%input = nhy_flux_monthly_input
call broadcast_scalar(ndep_shr_stream_year_first, master_task)
call broadcast_scalar(ndep_shr_stream_year_last, master_task)
call broadcast_scalar(ndep_shr_stream_year_align, master_task)
call broadcast_scalar(ndep_shr_stream_file, master_task)
call broadcast_scalar(ndep_shr_stream_scale_factor, master_task)
do n = 1,ecosys_tracer_cnt
call broadcast_scalar(tracer_init_ext(n)%mod_varname, master_task)
call broadcast_scalar(tracer_init_ext(n)%filename, master_task)
call broadcast_scalar(tracer_init_ext(n)%file_varname, master_task)
call broadcast_scalar(tracer_init_ext(n)%scale_factor, master_task)
call broadcast_scalar(tracer_init_ext(n)%default_val, master_task)
call broadcast_scalar(tracer_init_ext(n)%file_fmt, master_task)
end do
call broadcast_scalar(comp_surf_avg_freq_opt, master_task)
call broadcast_scalar(comp_surf_avg_freq, master_task)
call broadcast_scalar(use_nml_surf_vals, master_task)
call broadcast_scalar(surf_avg_dic_const, master_task)
call broadcast_scalar(surf_avg_alk_const, master_task)
call broadcast_scalar(ecosys_qsw_distrb_const, master_task)
call broadcast_scalar(lmarginal_seas, master_task)
call broadcast_scalar(lsource_sink, master_task)
call broadcast_scalar(lflux_gas_o2, master_task)
call broadcast_scalar(lflux_gas_co2, master_task)
call broadcast_scalar(locmip_k1_k2_bug_fix, master_task)
call broadcast_scalar(liron_patch, master_task)
call broadcast_scalar(iron_patch_flux_filename, master_task)
call broadcast_scalar(iron_patch_month, master_task)
call broadcast_scalar(atm_co2_opt, master_task)
call broadcast_scalar(atm_co2_const, master_task)
call broadcast_scalar(ecosys_tadvect_ctype, master_task)
tadvect_ctype = ecosys_tadvect_ctype
call broadcast_scalar(lecovars_full_depth_tavg, master_task)
!-----------------------------------------------------------------------
! set variables immediately dependent on namelist variables
!-----------------------------------------------------------------------
select case (comp_surf_avg_freq_opt)
case ('never')
comp_surf_avg_freq_iopt = freq_opt_never
case ('nyear')
comp_surf_avg_freq_iopt = freq_opt_nyear
case ('nmonth')
comp_surf_avg_freq_iopt = freq_opt_nmonth
case default
call document(subname, 'comp_surf_avg_freq_opt', comp_surf_avg_freq_opt)
call exit_POP(sigAbort, 'unknown comp_surf_avg_freq_opt')
end select
call init_time_flag('ecosys_comp_surf_avg', comp_surf_avg_flag, &
default=.false., freq_opt=comp_surf_avg_freq_iopt, &
freq=comp_surf_avg_freq, owner='ecosys_init')
select case (atm_co2_opt)
case ('const')
atm_co2_iopt = atm_co2_iopt_const
case ('drv_prog')
atm_co2_iopt = atm_co2_iopt_drv_prog
case ('drv_diag')
atm_co2_iopt = atm_co2_iopt_drv_diag
case default
call document(subname, 'atm_co2_opt', atm_co2_opt)
call exit_POP(sigAbort, 'unknown atm_co2_opt')
end select
!-----------------------------------------------------------------------
! namelist consistency checking
!-----------------------------------------------------------------------
if (use_nml_surf_vals .and. comp_surf_avg_freq_iopt /= freq_opt_never) then
call document(subname, 'use_nml_surf_vals', use_nml_surf_vals)
call document(subname, 'comp_surf_avg_freq_opt', comp_surf_avg_freq_opt)
call exit_POP(sigAbort, 'use_nml_surf_vals can only be .true. if ' /&
&/ ' comp_surf_avg_freq_opt is never')
endif
!-----------------------------------------------------------------------
! initialize virtual flux flag array
!-----------------------------------------------------------------------
vflux_flag = .false.
vflux_flag(dic_ind) = .true.
vflux_flag(alk_ind) = .true.
!-----------------------------------------------------------------------
! allocate various ecosys allocatable module variables
!-----------------------------------------------------------------------
allocate( PH_PREV(nx_block,ny_block,max_blocks_clinic) )
allocate( dust_FLUX_IN(nx_block,ny_block,max_blocks_clinic) )
allocate( PH_PREV_3D(nx_block,ny_block,km,max_blocks_clinic) )
PH_PREV_3D = c0
!-----------------------------------------------------------------------
! allocate and initialize LAND_MASK
!-----------------------------------------------------------------------
allocate( LAND_MASK(nx_block,ny_block,nblocks_clinic) )
if (lmarginal_seas) then
LAND_MASK = REGION_MASK /= 0
else
LAND_MASK = REGION_MASK > 0
endif
!-----------------------------------------------------------------------
! initialize tracers
!-----------------------------------------------------------------------
select case (init_ecosys_option)
case ('restart', 'ccsm_continue', 'ccsm_branch', 'ccsm_hybrid')
ecosys_restart_filename = char_blank
if (init_ecosys_init_file == 'same_as_TS') then
if (read_restart_filename == 'undefined') then
call document(subname, 'no restart file to read ecosys from')
call exit_POP(sigAbort, 'stopping in ' /&
&/ subname)
endif
ecosys_restart_filename = read_restart_filename
init_ecosys_init_file_fmt = init_ts_file_fmt
else ! do not read from TS restart file
ecosys_restart_filename = trim(init_ecosys_init_file)
endif
call rest_read_tracer_block(init_ecosys_init_file_fmt, &
ecosys_restart_filename, &
tracer_d_module, &
TRACER_MODULE)
call read_field(init_ecosys_init_file_fmt, &
ecosys_restart_filename, &
'PH_SURF', PH_PREV)
if (use_nml_surf_vals) then
surf_avg = c0
surf_avg(dic_ind) = surf_avg_dic_const
surf_avg(alk_ind) = surf_avg_alk_const
else
call extract_surf_avg(init_ecosys_init_file_fmt, &
ecosys_restart_filename)
endif
call eval_time_flag(comp_surf_avg_flag) ! evaluates time_flag(comp_surf_avg_flag)%value via time_to_do
if (check_time_flag(comp_surf_avg_flag)) &
call comp_surf_avg(TRACER_MODULE(:,:,1,:,oldtime,:), &
TRACER_MODULE(:,:,1,:,curtime,:))
case ('file', 'ccsm_startup')
call document(subname, 'ecosystem vars being read from separate files')
call file_read_tracer_block(init_ecosys_init_file_fmt, &
init_ecosys_init_file, &
tracer_d_module, &
ind_name_table, &
tracer_init_ext, &
TRACER_MODULE)
if (n_topo_smooth > 0) then
do n = 1, ecosys_tracer_cnt
do k=1,km
call fill_points(k,TRACER_MODULE(:,:,k,n,oldtime,:), &
errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'ecosys_init: error in fill points for tracers(oldtime)')
return
endif
call fill_points(k,TRACER_MODULE(:,:,k,n,curtime,:), &
errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'ecosys_init: error in fill points for tracers(newtime)')
return
endif
enddo
enddo
endif
PH_PREV = c0
if (use_nml_surf_vals) then
surf_avg = c0
surf_avg(dic_ind) = surf_avg_dic_const
surf_avg(alk_ind) = surf_avg_alk_const
else
call comp_surf_avg(TRACER_MODULE(:,:,1,:,oldtime,:), &
TRACER_MODULE(:,:,1,:,curtime,:))
endif
case default
call document(subname, 'init_ecosys_option', init_ecosys_option)
call exit_POP(sigAbort, 'unknown init_ecosys_option')
end select
!-----------------------------------------------------------------------
! register Chl field for short-wave absorption
! apply land mask to tracers
! set Chl field for short-wave absorption
!-----------------------------------------------------------------------
call named_field_register('model_chlorophyll', totChl_surf_nf_ind)
!$OMP PARALLEL DO PRIVATE(iblock,n,k,WORK)
do iblock=1,nblocks_clinic
do n = 1,ecosys_tracer_cnt
do k = 1,km
where (.not. LAND_MASK(:,:,iblock) .or. k > KMT(:,:,iblock))
TRACER_MODULE(:,:,k,n,curtime,iblock) = c0
TRACER_MODULE(:,:,k,n,oldtime,iblock) = c0
end where
end do
end do
WORK = max(c0,p5*(TRACER_MODULE(:,:,1,spChl_ind,oldtime,iblock) + &
TRACER_MODULE(:,:,1,spChl_ind,curtime,iblock))) + &
max(c0,p5*(TRACER_MODULE(:,:,1,diatChl_ind,oldtime,iblock) + &
TRACER_MODULE(:,:,1,diatChl_ind,curtime,iblock))) + &
max(c0,p5*(TRACER_MODULE(:,:,1,diazChl_ind,oldtime,iblock) + &
TRACER_MODULE(:,:,1,diazChl_ind,curtime,iblock)))
call named_field_set(totChl_surf_nf_ind, iblock, WORK)
enddo
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
! timer init
!-----------------------------------------------------------------------
call get_timer(ecosys_comp_CO3terms_timer, 'comp_CO3terms', &
nblocks_clinic, distrb_clinic%nprocs)
call get_timer(ecosys_interior_timer, 'ECOSYS_INTERIOR', &
nblocks_clinic, distrb_clinic%nprocs)
call get_timer(ecosys_sflux_timer, 'ECOSYS_SFLUX',1, &
distrb_clinic%nprocs)
if (ndep_data_type == 'shr_stream') then
call get_timer(ecosys_shr_strdata_advance_timer, &
'ecosys_shr_strdata_advance',1, distrb_clinic%nprocs)
endif
!-----------------------------------------------------------------------
! call other initialization subroutines
!-----------------------------------------------------------------------
call ecosys_init_tavg
call ecosys_init_sflux
call ecosys_init_interior_restore
!-----------------------------------------------------------------------
! set lfull_depth_tavg flag for short-lived ecosystem tracers
!-----------------------------------------------------------------------
tracer_d_module(spC_ind )%lfull_depth_tavg = lecovars_full_depth_tavg
tracer_d_module(spChl_ind )%lfull_depth_tavg = lecovars_full_depth_tavg
tracer_d_module(spCaCO3_ind)%lfull_depth_tavg = lecovars_full_depth_tavg
tracer_d_module(diatC_ind )%lfull_depth_tavg = lecovars_full_depth_tavg
tracer_d_module(diatChl_ind)%lfull_depth_tavg = lecovars_full_depth_tavg
tracer_d_module(zooC_ind )%lfull_depth_tavg = lecovars_full_depth_tavg
tracer_d_module(spFe_ind )%lfull_depth_tavg = lecovars_full_depth_tavg
tracer_d_module(diatSi_ind )%lfull_depth_tavg = lecovars_full_depth_tavg
tracer_d_module(diatFe_ind )%lfull_depth_tavg = lecovars_full_depth_tavg
tracer_d_module(diazC_ind )%lfull_depth_tavg = lecovars_full_depth_tavg
tracer_d_module(diazChl_ind)%lfull_depth_tavg = lecovars_full_depth_tavg
tracer_d_module(diazFe_ind )%lfull_depth_tavg = lecovars_full_depth_tavg
!-----------------------------------------------------------------------
!EOC
end subroutine ecosys_init
!***********************************************************************
!BOP
! !IROUTINE: extract_surf_avg
! !INTERFACE:
subroutine extract_surf_avg(init_ecosys_init_file_fmt, & 1,6
ecosys_restart_filename)
! !DESCRIPTION:
! Extract average surface values from restart file.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
character (*), intent(in) :: &
init_ecosys_init_file_fmt, & ! file format (bin or nc)
ecosys_restart_filename ! file name for restart file
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
type (datafile) ::&
restart_file ! io file descriptor
integer (int_kind) :: &
n ! tracer index
character (char_len) :: &
short_name ! tracer name temporaries
!-----------------------------------------------------------------------
surf_avg = c0
restart_file = construct_file(init_ecosys_init_file_fmt, &
full_name=trim(ecosys_restart_filename), &
record_length=rec_type_dbl, &
recl_words=nx_global*ny_global)
do n = 1, ecosys_tracer_cnt
if (vflux_flag(n)) then
short_name = 'surf_avg_' /&
&/ ind_name_table(n)%name
call add_attrib_file(restart_file, trim(short_name), surf_avg(n))
endif
end do
call data_set(restart_file, 'open_read')
do n = 1, ecosys_tracer_cnt
if (vflux_flag(n)) then
short_name = 'surf_avg_' /&
&/ ind_name_table(n)%name
call extract_attrib_file(restart_file, trim(short_name), surf_avg(n))
endif
end do
call data_set (restart_file, 'close')
call destroy_file (restart_file)
!-----------------------------------------------------------------------
!EOC
end subroutine extract_surf_avg
!***********************************************************************
!BOP
! !IROUTINE: ecosys_init_tavg
! !INTERFACE:
subroutine ecosys_init_tavg 1,113
! !DESCRIPTION:
! call define_tavg_field for nonstandard tavg fields
!
! !REVISION HISTORY:
! same as module
!
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
integer (int_kind) :: &
var_cnt ! how many tavg variables are defined
!-----------------------------------------------------------------------
! 2D fields related to surface fluxes
!-----------------------------------------------------------------------
var_cnt = 0
call define_tavg_field(tavg_ECOSYS_IFRAC,'ECOSYS_IFRAC',2, &
long_name='Ice Fraction for ecosys fluxes', &
units='fraction', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_ECOSYS_IFRAC_2,'ECOSYS_IFRAC_2',2, &
long_name='Ice Fraction for ecosys fluxes', &
units='fraction', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_ECOSYS_XKW,'ECOSYS_XKW',2, &
long_name='XKW for ecosys fluxes', &
units='cm/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_ECOSYS_XKW_2,'ECOSYS_XKW_2',2, &
long_name='XKW for ecosys fluxes', &
units='cm/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_ECOSYS_ATM_PRESS,'ECOSYS_ATM_PRESS',2, &
long_name='Atmospheric Pressure for ecosys fluxes', &
units='atmospheres', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_PV_O2,'PV_O2',2, &
long_name='PV_O2', &
units='cm/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_SCHMIDT_O2,'SCHMIDT_O2',2, &
long_name='O2 Schmidt Number', &
units='none', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_O2SAT,'O2SAT',2, &
long_name='O2 Saturation', &
units='mmol/m^3', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_CO2STAR,'CO2STAR',2, &
long_name='CO2 Star', &
units='mmol/m^3', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_DCO2STAR,'DCO2STAR',2, &
long_name='D CO2 Star', &
units='mmol/m^3', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_pCO2SURF,'pCO2SURF',2, &
long_name='surface pCO2', &
units='ppmv', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_DpCO2,'DpCO2',2, &
long_name='D pCO2', &
units='ppmv', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_DpCO2_2,'DpCO2_2',2, &
long_name='D pCO2', &
units='ppmv', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_PV_CO2,'PV_CO2',2, &
long_name='CO2 Piston Velocity', &
units='cm/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_SCHMIDT_CO2,'SCHMIDT_CO2',2, &
long_name='CO2 Schmidt Number', &
units='none', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_DIC_GAS_FLUX,'FG_CO2',2, &
long_name='DIC Surface Gas Flux', &
units='mmol/m^3 cm/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_DIC_GAS_FLUX_2,'FG_CO2_2',2, &
long_name='DIC Surface Gas Flux', &
units='mmol/m^3 cm/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_PH,'PH',2, &
long_name='Surface pH', &
units='none', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_ATM_CO2,'ATM_CO2',2, &
long_name='Atmospheric CO2', &
units='ppmv', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
call define_tavg_field(tavg_O2_GAS_FLUX_2,'STF_O2_2',2, &
long_name='Dissolved Oxygen Surface Flux', &
units='mmol/m^3 cm/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
var_cnt = var_cnt+1
!-----------------------------------------------------------------------
allocate(ECO_SFLUX_TAVG(nx_block,ny_block,var_cnt,max_blocks_clinic))
ECO_SFLUX_TAVG = c0
!-----------------------------------------------------------------------
! The follow surface flux related fields are not with the above because
! their values are directly available from STF in ecosys_tavg_forcing
! and do not need to be stored in ECO_SFLUX_TAVG.
!-----------------------------------------------------------------------
call define_tavg_field(tavg_IRON_FLUX,'IRON_FLUX',2, &
long_name='Iron Flux', &
units='nmol/cm^2/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_DUST_FLUX,'DUST_FLUX',2, &
long_name='Dust Flux', &
units='g/cm^2/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_NOx_FLUX,'NOx_FLUX',2, &
long_name='Flux of NOx from Atmosphere', &
units='nmol/cm^2/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_NHy_FLUX,'NHy_FLUX',2, &
long_name='Flux of NHy from Atmosphere', &
units='nmol/cm^2/s', grid_loc='2110', &
coordinates='TLONG TLAT time')
!-----------------------------------------------------------------------
! nonstandard 2D fields
!-----------------------------------------------------------------------
call define_tavg_field(tavg_O2_ZMIN,'O2_ZMIN',2, &
long_name='Vertical Minimum of O2', &
units='mmol/m^3', grid_loc='2111', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_O2_ZMIN_DEPTH,'O2_ZMIN_DEPTH',2, &
long_name='Depth of Vertical Minimum of O2', &
units='cm', grid_loc='2111', &
coordinates='TLONG TLAT time')
!-----------------------------------------------------------------------
! nonstandard 3D fields
!-----------------------------------------------------------------------
call define_tavg_field(tavg_O2_PRODUCTION,'O2_PRODUCTION',3, &
long_name='O2 Production', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_O2_CONSUMPTION,'O2_CONSUMPTION',3, &
long_name='O2 Consumption', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_AOU,'AOU',3, &
long_name='Apparent O2 Utilization ', &
units='mmol/m^3', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_PO4_RESTORE,'PO4_RESTORE',3, &
long_name='PO4 Restoring', &
units='mmol/m^3', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_NO3_RESTORE,'NO3_RESTORE',3, &
long_name='NO3 Restoring', &
units='mmol/m^3', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_SiO3_RESTORE,'SiO3_RESTORE',3, &
long_name='SiO3 Restoring', &
units='mmol/m^3', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_PAR_avg,'PAR_avg',3, &
long_name='PAR Average over Model Cell', &
units='w/m^2', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_POC_FLUX_IN,'POC_FLUX_IN',3, &
long_name='POC Flux into Cell', &
units='mmol/m^3 cm/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_POC_PROD,'POC_PROD',3, &
long_name='POC Production', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_POC_REMIN,'POC_REMIN',3, &
long_name='POC Remineralization', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_POC_ACCUM,'POC_ACCUM',3, &
long_name='POC Accumulation', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_CaCO3_FLUX_IN,'CaCO3_FLUX_IN',3, &
long_name='CaCO3 flux into cell', &
units='mmol/m^3 cm/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_CaCO3_PROD,'CaCO3_PROD',3, &
long_name='CaCO3 Production', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_CaCO3_REMIN,'CaCO3_REMIN',3, &
long_name='CaCO3 Remineralization', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_SiO2_FLUX_IN,'SiO2_FLUX_IN',3, &
long_name='SiO2 Flux into Cell', &
units='mmol/m^3 cm/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_SiO2_PROD,'SiO2_PROD',3, &
long_name='SiO2 Production', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_SiO2_REMIN,'SiO2_REMIN',3, &
long_name='SiO2 Remineralization', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_dust_FLUX_IN,'dust_FLUX_IN',3, &
long_name='Dust Flux into Cell', &
units='ng/s/m^2', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_dust_REMIN,'dust_REMIN',3, &
long_name='Dust Remineralization', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_P_iron_FLUX_IN,'P_iron_FLUX_IN',3, &
long_name='P_iron Flux into Cell', &
units='mmol/m^3 cm/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_P_iron_PROD,'P_iron_PROD',3, &
long_name='P_iron Production', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_P_iron_REMIN,'P_iron_REMIN',3, &
long_name='P_iron Remineralization', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_graze_sp,'graze_sp',3, &
long_name='Small Phyto Grazing', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_graze_diat,'graze_diat',3, &
long_name='Diatom Grazing', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_graze_diaz,'graze_diaz',3, &
long_name='Diazotroph Grazing', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_graze_TOT,'graze_TOT',3, &
long_name='Total Grazing', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_sp_loss,'sp_loss',3, &
long_name='Small Phyto Loss', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diat_loss,'diat_loss',3, &
long_name='Diatom Loss', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diaz_loss,'diaz_loss',3, &
long_name='Diaz Loss', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_zoo_loss,'zoo_loss',3, &
long_name='Zooplankton Loss', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_sp_agg,'sp_agg',3, &
long_name='Small Phyto Aggregate', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diat_agg,'diat_agg',3, &
long_name='Diatom Aggregate', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_photoC_sp,'photoC_sp',3, &
long_name='Small Phyto C Fixation', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_CaCO3_form,'CaCO3_form',3, &
long_name='CaCO3 Formation', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_photoC_diat,'photoC_diat',3, &
long_name='Diatom C Fixation', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_photoC_diaz,'photoC_diaz',3, &
long_name='Diaz C Fixation', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_photoC_TOT,'photoC_TOT',3, &
long_name='Total C Fixation', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_photoC_sp_zint,'photoC_sp_zint',2, &
long_name='Small Phyto C Fixation Vertical Integral',&
units='mmol/m^3 cm/s', grid_loc='2111', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_CaCO3_form_zint,'CaCO3_form_zint',2, &
long_name='CaCO3 Formation Vertical Integral',&
units='mmol/m^3 cm/s', grid_loc='2111', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_photoC_diat_zint,'photoC_diat_zint',2, &
long_name='Diatom C Fixation Vertical Integral',&
units='mmol/m^3 cm/s', grid_loc='2111', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_photoC_diaz_zint,'photoC_diaz_zint',2, &
long_name='Diaz C Fixation Vertical Integral',&
units='mmol/m^3 cm/s', grid_loc='2111', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_photoC_TOT_zint,'photoC_TOT_zint',2, &
long_name='Total C Fixation Vertical Integral',&
units='mmol/m^3 cm/s', grid_loc='2111', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_photoC_NO3_sp,'photoC_NO3_sp',3, &
long_name='Small Phyto C Fixation from NO3', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_photoC_NO3_sp_zint,'photoC_NO3_sp_zint',2,&
long_name='Small Phyto C Fixation from NO3 Vertical Integral',&
units='mmol/m^3 cm/s', grid_loc='2111', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_photoC_NO3_diat,'photoC_NO3_diat',3, &
long_name='Diatom C Fixation from NO3', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_photoC_NO3_diat_zint,'photoC_NO3_diat_zint',2,&
long_name='Diatom C Fixation from NO3 Vertical Integral',&
units='mmol/m^3 cm/s', grid_loc='2111', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_photoC_NO3_diaz,'photoC_NO3_diaz',3, &
long_name='Diaz C Fixation from NO3', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_photoC_NO3_diaz_zint,'photoC_NO3_diaz_zint',2,&
long_name='Diaz C Fixation from NO3 Vertical Integral',&
units='mmol/m^3 cm/s', grid_loc='2111', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_photoC_NO3_TOT,'photoC_NO3_TOT',3, &
long_name='Total C Fixation from NO3', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_photoC_NO3_TOT_zint,'photoC_NO3_TOT_zint',2,&
long_name='Total C Fixation from NO3 Vertical Integral',&
units='mmol/m^3 cm/s', grid_loc='2111', &
coordinates='TLONG TLAT time')
!-----------------------------------------------------------------------
! MORE nonstandard 3D fields
!-----------------------------------------------------------------------
call define_tavg_field(tavg_DOC_prod,'DOC_prod',3, &
long_name='DOC Production', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_DOC_remin,'DOC_remin',3, &
long_name='DOC Remineralization', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_DON_prod,'DON_prod',3, &
long_name='DON Production', &
units='mmols/m^3', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_DON_remin,'DON_remin',3, &
long_name='DON Remineralization', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_DOFe_prod,'DOFe_prod',3, &
long_name='DOFe Production', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_DOFe_remin,'DOFe_remin',3, &
long_name='DOFe Remineralization', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_DOP_prod,'DOP_prod',3, &
long_name='DOP Production', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_DOP_remin,'DOP_remin',3, &
long_name='DOP Remineralization', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_Fe_scavenge,'Fe_scavenge',3, &
long_name='Iron Scavenging', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_Fe_scavenge_rate,'Fe_scavenge_rate',3, &
long_name='Iron Scavenging Rate', &
units='1/y', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_sp_N_lim,'sp_N_lim',3, &
long_name='Small Phyto N Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_sp_Fe_lim,'sp_Fe_lim',3, &
long_name='Small Phyto Fe Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_sp_PO4_lim,'sp_PO4_lim',3, &
long_name='Small Phyto PO4 Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_sp_light_lim,'sp_light_lim',3, &
long_name='Small Phyto Light Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diat_N_lim,'diat_N_lim',3, &
long_name='Diatom N Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diat_Fe_lim,'diat_Fe_lim',3, &
long_name='Diatom Fe Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diat_PO4_lim,'diat_PO4_lim',3, &
long_name='Diatom PO4 Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diat_SiO3_lim,'diat_SiO3_lim',3, &
long_name='Diatom SiO3 Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diat_light_lim,'diat_light_lim',3, &
long_name='Diatom Light Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diaz_Fe_lim,'diaz_Fe_lim',3, &
long_name='Diaz Fe Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diaz_P_lim,'diaz_P_lim',3, &
long_name='Diaz PO4 Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diaz_light_lim,'diaz_light_lim',3, &
long_name='Diaz LIGHT Limitation', &
units='none', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_diaz_Nfix,'diaz_Nfix',3, &
long_name='Diaz N Fixation', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_bSi_form,'bSi_form',3, &
long_name='Diatom Si Uptake', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_NITRIF,'NITRIF',3, &
long_name='Nitrification', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_DENITRIF,'DENITRIF',3, &
long_name='Denitrification', &
units='mmol/m^3/s', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_photoFe_sp,'photoFe_sp',3, &
long_name='Small Phyto Fe Uptake', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_photoFe_diat,'photoFe_diat',3, &
long_name='Diatom Fe Uptake', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_photoFe_diaz,'photoFe_diaz',3, &
long_name='Diaz Fe Uptake', &
units='mmol/m^3/s', grid_loc='3114', &
coordinates='TLONG TLAT z_t_150m time')
call define_tavg_field(tavg_CO3,'CO3',3, &
long_name='Carbonate Ion Concentration', &
units='mmol/m^3', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_HCO3,'HCO3',3, &
long_name='Bicarbonate Ion Concentration', &
units='mmol/m^3', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_H2CO3,'H2CO3',3, &
long_name='Carbonic Acid Concentration', &
units='mmol/m^3', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_pH_3D,'pH_3D',3, &
long_name='pH', &
units='none', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_co3_sat_calc,'co3_sat_calc',3, &
long_name='CO3 concentration at calcite saturation', &
units='mmol/m^3', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_zsatcalc,'zsatcalc',2, &
long_name='Calcite Saturation Depth', &
units='cm', grid_loc='2110', &
coordinates='TLONG TLAT time')
call define_tavg_field(tavg_co3_sat_arag,'co3_sat_arag',3, &
long_name='CO3 concentration at aragonite saturation', &
units='mmol/m^3', grid_loc='3111', &
coordinates='TLONG TLAT z_t time')
call define_tavg_field(tavg_zsatarag,'zsatarag',2, &
long_name='Aragonite Saturation Depth', &
units='cm', grid_loc='2110', &
coordinates='TLONG TLAT time')
!-----------------------------------------------------------------------
!EOC
end subroutine ecosys_init_tavg
!***********************************************************************
!BOP
! !IROUTINE: ecosys_set_interior
! !INTERFACE:
subroutine ecosys_set_interior(k, TEMP_OLD, TEMP_CUR, SALT_OLD, SALT_CUR, & 1,178
TRACER_MODULE_OLD, TRACER_MODULE_CUR, DTRACER_MODULE, this_block)
! !DESCRIPTION:
! Compute time derivatives for ecosystem state variables
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: &
k ! vertical level index
real (r8), dimension(nx_block,ny_block), intent(in) :: &
TEMP_OLD, &! old potential temperature (C)
TEMP_CUR, &! current potential temperature (C)
SALT_OLD, &! old salinity (msu)
SALT_CUR ! current salinity (msu)
real (r8), dimension(nx_block,ny_block,km,ecosys_tracer_cnt), intent(in) :: &
TRACER_MODULE_OLD, &! old tracer values
TRACER_MODULE_CUR ! current tracer values
type (block), intent(in) :: &
this_block ! block info for the current block
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,ecosys_tracer_cnt), intent(out) :: &
DTRACER_MODULE ! computed source/sink terms
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
character(*), parameter :: &
subname = 'ecosys_mod:ecosys_set_interior'
real (r8), parameter :: &
epsC = 1.00e-8, & ! small C concentration (mmol C/m^3)
epsTinv = 3.17e-8, & ! small inverse time scale (1/year) (1/sec)
epsnondim = 1.00e-6 ! small non-dimensional number (non-dim)
type(sinking_particle), save :: &
POC, & ! base units = nmol C
P_CaCO3, & ! base units = nmol CaCO3
P_SiO2, & ! base units = nmol SiO2
dust, & ! base units = g
P_iron ! base units = nmol Fe
real (r8), dimension(nx_block,ny_block,max_blocks_clinic), save :: &
QA_dust_def, & ! incoming deficit in the QA(dust) POC flux
ZSATCALC, & ! Calcite Saturation Depth
ZSATARAG, & ! Aragonite Saturation Depth
CO3_CALC_ANOM_km1,&! CO3 concentration above calcite saturation at k-1
CO3_ARAG_ANOM_km1 ! CO3 concentration above aragonite saturation at k-1
real (r8), dimension(nx_block,ny_block) :: &
TEMP, & ! local copy of model TEMP
SALT, & ! local copy of model SALT
DIC_loc, & ! local copy of model DIC
ALK_loc, & ! local copy of model ALK
PO4_loc, & ! local copy of model PO4
NO3_loc, & ! local copy of model NO3
SiO3_loc, & ! local copy of model SiO3
NH4_loc, & ! local copy of model NH4
Fe_loc, & ! local copy of model Fe
O2_loc, & ! local copy of model O2
DOC_loc, & ! local copy of model DOC
spC_loc, & ! local copy of model spC
spChl_loc, & ! local copy of model spChl
spCaCO3_loc, & ! local copy of model spCaCO3
diatC_loc, & ! local copy of model diatC
diatChl_loc, & ! local copy of model diatChl
zooC_loc, & ! local copy of model zooC
spFe_loc, & ! local copy of model spFe
diatSi_loc, & ! local copy of model diatSi
diatFe_loc, & ! local copy of model diatFe
diazC_loc, & ! local copy of model diazC
diazChl_loc, & ! local copy of model diazChl
diazFe_loc, & ! local copy of model diazFe
DON_loc, & ! local copy of model DON
DOFe_loc, & ! local copy of model DOFe
DOP_loc ! local copy of model DOP
real (r8), dimension(nx_block,ny_block) :: &
WORK1,WORK2,WORK3 ! temporaries
real (r8) :: &
C_loss_thres, & ! bio-C threshold at which losses go to zero (mmol C/m^3)
f_loss_thres ! fraction of grazing loss reduction at depth
real (r8), dimension(nx_block,ny_block) :: &
PAR_in, & ! photosynthetically available radiation (W/m^2)
KPARdz, & ! PAR adsorption coefficient (non-dim)
PAR_avg, & ! average PAR over mixed layer depth (W/m^2)
DOC_prod, & ! production of DOC (mmol C/m^3/sec)
DOC_remin, & ! remineralization of DOC (mmol C/m^3/sec)
NITRIF, & ! nitrification (NH4 -> NO3) (mmol N/m^3/sec)
DENITRIF, & ! nitrification (NO3 -> N2) (mmol N/m^3/sec)
RESTORE ! restoring terms for nutrients (mmol ./m^3/sec)
real (r8), dimension(nx_block,ny_block) :: &
z_umax, & ! max. zoo growth rate on sp at local T (1/sec)
diat_umax, & ! max. zoo growth rate on diatoms at local T (1/sec)
z_mort, & ! zoo respiration loss, (1/sec/((mmol C/m3))
C_loss_diaz, & ! bio-C threshold at which losses go to zero (mmol C/m^3)
z_mort2, & ! zoo quad mort rate, tohigherlevels (1/sec/((mmol C/m3))
diaz_umax ! max. zoo growth rate on diazotrophs at local T (1/sec)
real (r8), dimension(nx_block,ny_block) :: &
thetaC_sp, & ! local Chl/C ratio in small phyto (mg Chl/mmol C)
thetaC_diat, & ! local Chl/C ratio in diatoms (mg Chl/mmol C)
QCaCO3, & ! small phyto CaCO3/C ratio (mmol CaCO3/mmol C)
Tfunc, & ! temp response function GD98 (non-dim)
VNO3_sp, & ! small phyto NO3 uptake rate (non-dim)
VNH4_sp, & ! small phyto NH4 uptake rate (non-dim)
VNtot_sp, & ! small phyto total N uptake rate (non-dim)
VFeC_sp, & ! ??? small phyto C-specific iron uptake (non-dim)
VPO4_sp, & ! ??? small phyto C-specific po4 uptake (non-dim)
f_nut, & ! nut limitation factor, modifies C fixation (non-dim)
PCmax, & ! max value of PCphoto at temperature TEMP (1/sec)
PCphoto_sp, & ! small C-specific rate of photosynth. (1/sec)
photoC_sp, & ! small phyto C-fixation (mmol C/m^3/sec)
NO3_V_sp, & ! nitrate uptake by small phyto (mmol NO3/m^3/sec)
NH4_V_sp, & ! ammonium uptake by small phyto (mmol NH4/m^3/sec)
VNC_sp, & ! small phyto C-specific N uptake rate (mmol N/mmol C/sec)
pChl, & ! Chl synth. regulation term (mg Chl/mmol N)
photoacc_sp, & ! Chl synth. term in photoadapt. (GD98) (mg Chl/m^3/sec)
CaCO3_PROD, & ! prod. of CaCO3 by small phyto (mmol CaCO3/m^3/sec)
VNO3_diat, & ! diatom nitrate uptake rate (non-dim)
VNH4_diat, & ! diatom ammonium uptake rate (non-dim)
VNtot_diat, & ! diatom total N uptake rate (non-dim)
VFeC_diat, & ! diatom C-specific iron uptake (non-dim)
VPO4_diat, & ! diatom C-specific PO4 uptake (non-dim)
VSiO3_diat, & ! C-specific SiO3 uptake for diatoms (non-dim)
PCphoto_diat, & ! diatom C-specific rate of photosynth. (1/sec)
photoC_diat, & ! diatom C-fixation (mmol C/m^3/sec)
NO3_V_diat, & ! nitrate uptake by diatoms (mmol NO3/m^3/sec)
NH4_V_diat, & ! ammonium uptake by diatoms (mmol NH4/m^3/sec)
VNC_diat, & ! diatom C-specific N uptake rate (mmol N/mmol C/sec)
photoacc_diat, & ! Chl synth. term in photoadapt. (GD98) (mg Chl/m^3/sec)
reduceV, & ! factor in nutrient uptake (mmol C/m^3)^2
graze_sp, & ! grazing rate on small phyto (mmol C/m^3/sec)
graze_sp_zoo, & ! graze_sp routed to zoo (mmol C/m^3/sec)
graze_sp_poc, & ! graze_sp routed to poc (mmol C/m^3/sec)
graze_sp_doc, & ! graze_sp routed to doc (mmol C/m^3/sec)
graze_sp_dic ! graze_sp routed to dic (mmol C/m^3/sec)
real (r8), dimension(nx_block,ny_block) :: & ! max of 39 continuation lines
graze_diat, & ! grazing rate on diatoms (mmol C/m^3/sec)
graze_diat_zoo, & ! graze_diat routed to zoo (mmol C/m^3/sec)
graze_diat_poc, & ! graze_diat routed to poc (mmol C/m^3/sec)
graze_diat_doc, & ! graze_diat routed to doc (mmol C/m^3/sec)
graze_diat_dic, & ! graze_diat routed to dic (mmol C/m^3/sec)
Pprime, & ! used to limit phyto mort at low biomass (mmol C/m^3)
sp_loss, & ! small phyto non-grazing mort (mmol C/m^3/sec)
sp_loss_poc, & ! sp_loss routed to poc (mmol C/m^3/sec)
sp_loss_doc, & ! sp_loss routed to doc (mmol C/m^3/sec)
sp_loss_dic, & ! sp_loss routed to dic (mmol C/m^3/sec)
sp_agg, & ! small phyto agg loss (mmol C/m^3/sec)
diat_loss, & ! diatom non-grazing mort (mmol C/m^3/sec)
diat_loss_poc, & ! diat_loss routed to poc (mmol C/m^3/sec)
diat_loss_doc, & ! diat_loss routed to doc (mmol C/m^3/sec)
diat_loss_dic, & ! diat_loss routed to dic (mmol C/m^3/sec)
diat_agg, & ! diatom agg (mmol C/m^3/sec)
f_zoo_detr, & ! frac of zoo losses into large detrital pool (non-dim)
Fe_scavenge, & ! loss of dissolved iron, scavenging (mmol Fe/m^3/sec)
Zprime, & ! used to limit zoo mort at low biomass (mmol C/m^3)
zoo_loss, & ! mortality & higher trophic grazing on zooplankton (mmol C/m^3/sec)
zoo_loss_doc, & ! zoo_loss routed to doc (mmol C/m^3/sec)
zoo_loss_dic, & ! zoo_loss routed to dic (mmol C/m^3/sec)
light_lim, & ! light limitation factor
Qsi, & ! Diatom initial Si/C ratio (mmol Si/mmol C)
gQsi, & ! diatom Si/C ratio for growth (new biomass)
Qfe_sp, & ! small phyto init fe/C ratio (mmolFe/mmolC)
gQfe_sp, & ! small phyto fe/C for growth
Qfe_diat, & ! diatom init fe/C ratio
gQfe_diat, & ! diatom fe/C ratio for growth
Qfe_diaz, & ! diazotrophs init fe/C ratio
gQfe_diaz ! diazotroph fe/C ratio for new growth
real (r8), dimension(nx_block,ny_block) :: & ! max of 39 continuation lines
PCphoto_diaz, & ! diazotroph C-specific rate of photosynth. (1/sec)
photoC_diaz, & ! diazotroph C-fixation (mmol C/m^3/sec)
Vfec_diaz, & ! diazotroph C-specific iron uptake (non-dim)
Vpo4_diaz, & ! diazotroph C-specific po4 uptake (non-dim)
photoacc_diaz, & ! Chl synth. term in photoadapt. (GD98) (mg Chl/m^3/sec)
Vnc_diaz, & ! diazotroph C-specific N uptake rate (mmol N/mmol C/sec)
diaz_Nexcrete, & ! diazotroph fixed N excretion
diaz_Nfix, & ! diazotroph total Nitrogen fixation (mmol N/m^3/sec)
thetaC_diaz, & ! local Chl/C ratio in diazotrophs (mg Chl/mmol C)
photoFe_diaz, & ! iron uptake by diazotrophs (mmolFe/m^3/sec)
photoFe_diat, & ! iron uptake by diatoms
photoFe_sp, & ! iron uptake by small phyto
photoSi_diat, & ! silicon uptake by diatoms (mmol Si/m^3/sec)
remaining_diazP,& ! used in routing P from diazotroph losses
diaz_loss, & ! diazotroph non-grazing mort (mmol C/m^3/sec)
diaz_loss_doc, & ! mortality routed to DOM pool
diaz_loss_poc, & ! mortalitiy routed to POC pool
diaz_loss_dic, & ! mortality routed to remin
diaz_loss_dop, & ! P from mort routed to DOP pool
diaz_loss_dip, & ! P from mort routed to remin
! P from diaz losses must be routed differently than
! other elements to ensure that sinking detritus and
! zooplankton pools get their fixed P/C ratios, the
! remaining P is split evenly between DOP and PO4
graze_diaz, & ! grazing rate on diazotrophs (mmol C/m^3/sec)
graze_diaz_poc, & ! grazing routed to sinking detr (mmol C/m^3/sec)
graze_diaz_doc, & ! grazing routed to DOC (mmol C/m^3/sec)
graze_diaz_dic, & ! grazing routed to remin (mmol C/m^3/sec)
graze_diaz_zoo, & ! grazing routed to new zoo biomass
DON_remin, & ! portion of DON remineralized
DOFe_remin, & ! portion of DOFe remineralized
DOP_remin, & ! portion of DOP remineralized
DOM_remin, & ! fraction of DOM remineralized at current TEMP
Fe_scavenge_rate ! annual scavenging rate of iron as % of ambient
real (r8), dimension(nx_block,ny_block) :: & ! max of 39 continuation lines
DON_prod, & ! production of dissolved organic N
DOFe_prod, & ! produciton of dissolved organic Fe
DOP_prod, & ! production of dissolved organic P
po4_V_sp, & ! sp uptake of po4 (mmol po4/m3/s)
po4_V_diaz, & ! diaz uptake of po4 (mmol po4/m3/s)
dop_V_sp, & ! sp uptake of dop (mmol dop/m3/s)
dop_V_diaz, & ! diaz uptake of dop (mmol dop/m3/s)
cks, & ! constant used in Fe quota modification
cksi, & ! constant used in Si quota modification
VNO3_diaz, & ! relative NO3 uptake by diazotrophs
VNH4_diaz, & ! relative NH4 uptake by diazotrophs
VNtot_diaz, & ! total relative N uptake by diazotrophs
photoNO3_diaz, & ! total nitrate uptake by diazotrophs
photoNH4_diaz, & ! total NH4 uptake by diazotrophs
O2_PRODUCTION, & ! O2 production
O2_CONSUMPTION ! O2 consumption
real (r8), dimension(nx_block,ny_block) :: &
CO3, &! carbonate ion
HCO3, &! bicarbonate ion
H2CO3, &! carbonic acid
OMEGA_CALC, &! solubility ratio for aragonite
OMEGA_ARAG ! solubility ratio for calcite
integer (int_kind) :: &
bid, & ! local_block id
kk ! index for looping over k levels
!-----------------------------------------------------------------------
bid = this_block%local_id
call timer_start(ecosys_interior_timer, block_id=bid)
DTRACER_MODULE = c0
!-----------------------------------------------------------------------
! exit immediately if computations are not to be performed
!-----------------------------------------------------------------------
if (.not. lsource_sink) then
call timer_stop(ecosys_interior_timer, block_id=bid)
return
endif
!-----------------------------------------------------------------------
! create local copies of model tracers
! treat negative values as zero
! apply mask to local copies
!-----------------------------------------------------------------------
TEMP = p5*(TEMP_OLD + TEMP_CUR)
SALT = p5*(SALT_OLD + SALT_CUR)*salt_to_ppt
DIC_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,dic_ind) + &
TRACER_MODULE_CUR(:,:,k,dic_ind)))
ALK_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,alk_ind) + &
TRACER_MODULE_CUR(:,:,k,alk_ind)))
PO4_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,po4_ind) + &
TRACER_MODULE_CUR(:,:,k,po4_ind)))
NO3_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,no3_ind) + &
TRACER_MODULE_CUR(:,:,k,no3_ind)))
SiO3_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,sio3_ind) + &
TRACER_MODULE_CUR(:,:,k,sio3_ind)))
NH4_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,nh4_ind) + &
TRACER_MODULE_CUR(:,:,k,nh4_ind)))
Fe_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,fe_ind) + &
TRACER_MODULE_CUR(:,:,k,fe_ind)))
O2_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,o2_ind) + &
TRACER_MODULE_CUR(:,:,k,o2_ind)))
DOC_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,doc_ind) + &
TRACER_MODULE_CUR(:,:,k,doc_ind)))
spC_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,spC_ind) + &
TRACER_MODULE_CUR(:,:,k,spC_ind)))
spChl_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,spChl_ind) + &
TRACER_MODULE_CUR(:,:,k,spChl_ind)))
spCaCO3_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,spCaCO3_ind) + &
TRACER_MODULE_CUR(:,:,k,spCaCO3_ind)))
diatC_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diatC_ind) + &
TRACER_MODULE_CUR(:,:,k,diatC_ind)))
diatChl_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diatChl_ind) + &
TRACER_MODULE_CUR(:,:,k,diatChl_ind)))
zooC_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,zooC_ind) + &
TRACER_MODULE_CUR(:,:,k,zooC_ind)))
spFe_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,spFe_ind) + &
TRACER_MODULE_CUR(:,:,k,spFe_ind)))
diatSi_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diatSi_ind) + &
TRACER_MODULE_CUR(:,:,k,diatSi_ind)))
diatFe_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diatFe_ind) + &
TRACER_MODULE_CUR(:,:,k,diatFe_ind)))
diazC_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diazC_ind) + &
TRACER_MODULE_CUR(:,:,k,diazC_ind)))
diazChl_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diazChl_ind) + &
TRACER_MODULE_CUR(:,:,k,diazChl_ind)))
diazFe_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,diazFe_ind) + &
TRACER_MODULE_CUR(:,:,k,diazFe_ind)))
DON_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,don_ind) + &
TRACER_MODULE_CUR(:,:,k,don_ind)))
DOFe_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,dofe_ind) + &
TRACER_MODULE_CUR(:,:,k,dofe_ind)))
DOP_loc = max(c0, p5*(TRACER_MODULE_OLD(:,:,k,dop_ind) + &
TRACER_MODULE_CUR(:,:,k,dop_ind)))
where (.not. LAND_MASK(:,:,bid) .or. k > KMT(:,:,bid))
PO4_loc = c0
NO3_loc = c0
SiO3_loc = c0
NH4_loc = c0
Fe_loc = c0
O2_loc = c0
DOC_loc = c0
spC_loc = c0
spChl_loc = c0
spCaCO3_loc = c0
diatC_loc = c0
diatChl_loc = c0
zooC_loc = c0
spFe_loc = c0
diatSi_loc = c0
diatFe_loc = c0
diazC_loc = c0
diazChl_loc = c0
diazFe_loc = c0
DON_loc = c0
DOFe_loc = c0
DOP_loc = c0
end where
!-----------------------------------------------------------------------
! If any phyto box are zero, set others to zeros.
!-----------------------------------------------------------------------
where (spC_loc == c0 .or. spChl_loc == c0 .or. spFe_loc == c0)
spC_loc = c0
spChl_loc = c0
spCaCO3_loc = c0
spFe_loc = c0
end where
where (diatC_loc == c0 .or. diatChl_loc == c0 .or. diatFe_loc == c0 &
.or. diatSi_loc == c0)
diatC_loc = c0
diatChl_loc = c0
diatFe_loc = c0
diatSi_loc = c0
end where
where (diazC_loc == c0 .or. diazChl_loc == c0 .or. diazFe_loc == c0)
diazC_loc = c0
diazChl_loc = c0
diazFe_loc = c0
end where
!-----------------------------------------------------------------------
! set local variables, with incoming ratios
!-----------------------------------------------------------------------
thetaC_sp = spChl_loc / (spC_loc + epsC)
thetaC_diat = diatChl_loc / (diatC_loc + epsC)
thetaC_diaz = diazChl_loc / (diazC_loc + epsC)
Qsi = diatSi_loc / (diatC_loc + epsC)
Qfe_diat = diatFe_loc / (diatC_loc + epsC)
Qfe_sp = spFe_loc / (spC_loc + epsC)
Qfe_diaz = diazFe_loc / (diazC_loc + epsC)
where (Qsi > gQsi_max) Qsi = gQsi_max
!-----------------------------------------------------------------------
! DETERMINE NEW ELEMENTAL RATIOS FOR GROWTH (NEW BIOMASS)
!-----------------------------------------------------------------------
gQsi = gQsi_0
gQfe_diat = gQfe_diat_0
gQfe_sp = gQfe_sp_0
gQfe_diaz = gQfe_diaz_0
cks = 10._r8
cksi = 10._r8
!-----------------------------------------------------------------------
! Modify these initial ratios under low ambient iron conditions
!-----------------------------------------------------------------------
where (Fe_loc < cks * parm_diat_kfe)
gQfe_diat = max((gQfe_diat * Fe_loc /(cks * parm_diat_kfe)), &
gQfe_diat_min)
end where
where ((Fe_loc < cksi * parm_diat_kfe) .and. (Fe_loc > c0) .and. &
(SiO3_loc > (cksi * parm_diat_kSiO3)))
gQsi = min((gQsi*cksi*parm_diat_kfe/Fe_loc), gQsi_max)
end where
where (Fe_loc == c0)
gQsi = gQsi_max
end where
where (Fe_loc < cks * parm_sp_kfe)
gQfe_sp = max((gQfe_sp*Fe_loc/(cks * parm_sp_kfe)), &
gQfe_sp_min)
end where
where (Fe_loc < cks * diaz_kFe)
gQfe_diaz = max((gQfe_diaz*Fe_loc/(cks * diaz_kFe)), &
gQfe_diaz_min)
end where
!-----------------------------------------------------------------------
! Modify the initial si/C ratio under low ambient Si conditions
!-----------------------------------------------------------------------
where (SiO3_loc < (cksi * parm_diat_kSiO3))
gQsi = max((gQsi*SiO3_loc/(cksi * parm_diat_kSiO3)), &
gQsi_min)
end where
!-----------------------------------------------------------------------
! various k==1 initializations
!
! 0.45 fraction of incoming SW -> PAR (non-dim)
!-----------------------------------------------------------------------
if (k == 1) then
call init_particulate_terms(POC, P_CaCO3, P_SiO2, dust, P_iron, &
QA_dust_def, dust_FLUX_IN, this_block)
endif
!-----------------------------------------------------------------------
! QCaCO3 is the percentage of sp organic matter which is associated
! with coccolithophores
!-----------------------------------------------------------------------
QCaCO3 = spCaCO3_loc / (spC_loc + epsC)
where (QCaCO3 > QCaCO3_max) QCaCO3 = QCaCO3_max
!-----------------------------------------------------------------------
! compute PAR related quantities
!
! 0.03e-2 atten. coeff. per unit chlorophyll (1/cm/(mg Chl/m^3))
! 0.04e-2 atten. coeff. for water (1/cm)
!-----------------------------------------------------------------------
PAR_in = PAR_out(:,:,bid)
where (.not. LAND_MASK(:,:,bid) .or. k > KMT(:,:,bid)) PAR_in = c0
if (partial_bottom_cells) then
KPARdz = (k_chl * (spChl_loc + diatChl_loc + &
diazChl_loc) + k_h2o) * DZT(:,:,k,bid)
else
KPARdz = (k_chl * (spChl_loc + diatChl_loc + &
diazChl_loc) + k_h2o) * dz(k)
endif
PAR_out(:,:,bid) = PAR_in * exp(-KPARdz)
PAR_avg = PAR_in * (c1 - exp(-KPARdz)) / KPARdz
!-----------------------------------------------------------------------
! Tref = 30.0 reference temperature (deg. C)
!
! Using q10 formulation with Q10 value of 2.0 (Doney et al., 1996).
!-----------------------------------------------------------------------
Tfunc = Q_10**(((TEMP + T0_Kelvin)-(Tref + T0_Kelvin))/c10)
!-----------------------------------------------------------------------
! modify growth mort rates by Tfunc
!-----------------------------------------------------------------------
z_umax = parm_z_umax_0 * Tfunc
diat_umax = parm_diat_umax_0 * Tfunc
z_mort2 = parm_z_mort2_0 * Tfunc
z_mort = parm_z_mort_0 * Tfunc
diaz_umax = parm_diaz_umax_0 * Tfunc
DOM_remin = parm_sd_remin_0
!-----------------------------------------------------------------------
! Get relative nutrient uptake rates for phytoplankton,
! min. relative uptake rate modifies C fixation in the manner
! that the min. cell quota does in GD98.
!-----------------------------------------------------------------------
VNO3_sp = (NO3_loc / parm_sp_kNO3) / &
(c1 + (NO3_loc / parm_sp_kNO3) + (NH4_loc / parm_sp_kNH4))
VNH4_sp = (NH4_loc / parm_sp_kNH4) / &
(c1 + (NO3_loc / parm_sp_kNO3) + (NH4_loc / parm_sp_kNH4))
VNtot_sp = VNO3_sp + VNH4_sp
if (tavg_requested(tavg_sp_N_lim)) then
call accumulate_tavg_field(VNtot_sp, tavg_sp_N_lim,bid,k)
endif
!-----------------------------------------------------------------------
! get relative Fe uptake by phytoplankton
! get relative P uptake rates for phytoplankton
!-----------------------------------------------------------------------
VFeC_sp = Fe_loc / (Fe_loc + parm_sp_kFe)
if (tavg_requested(tavg_sp_Fe_lim)) then
call accumulate_tavg_field(VFeC_sp, tavg_sp_Fe_lim,bid,k)
endif
!-----------------------------------------------------------------------
! Partition diaz P uptake in same manner to no3/nh4 partition
!-----------------------------------------------------------------------
VPO4_sp = PO4_loc / (PO4_loc + parm_sp_kPO4)
if (tavg_requested(tavg_sp_PO4_lim)) then
call accumulate_tavg_field(VPO4_sp, tavg_sp_PO4_lim,bid,k)
endif
!-----------------------------------------------------------------------
! Small Phytoplankton C-fixation - given light and relative uptake rates
! determine small phyto nutrient limitation factor for carbon fixation
!-----------------------------------------------------------------------
f_nut = min(VNtot_sp, VFeC_sp)
f_nut = min(f_nut, VPO4_sp)
!-----------------------------------------------------------------------
! get small phyto photosynth. rate, phyto C biomass change, photoadapt
!-----------------------------------------------------------------------
PCmax = PCref * f_nut * Tfunc
light_lim = (c1 - exp((-c1 * parm_alphaChlsp * thetaC_sp * PAR_avg) / &
(PCmax + epsTinv)))
PCphoto_sp = PCmax * light_lim
if (tavg_requested(tavg_sp_light_lim)) then
call accumulate_tavg_field(light_lim, tavg_sp_light_lim,bid,k)
endif
photoC_sp = PCphoto_sp * spC_loc
!-----------------------------------------------------------------------
! Get nutrient uptakes by small phyto based on calculated C fixation
! total N uptake (VNC_sp) is used in photoadaption
!-----------------------------------------------------------------------
where (VNtot_sp > c0)
NO3_V_sp = (VNO3_sp / VNtot_sp) * photoC_sp * Q
NH4_V_sp = (VNH4_sp / VNtot_sp) * photoC_sp * Q
VNC_sp = PCphoto_sp * Q
elsewhere
NO3_V_sp = c0
NH4_V_sp = c0
VNC_sp = c0
photoC_sp = c0
end where
po4_V_sp = photoC_sp * Qp
photoFe_sp = photoC_sp * gQfe_sp
if (tavg_requested(tavg_photoFe_sp)) then
call accumulate_tavg_field(photoFe_sp, tavg_photoFe_sp,bid,k)
endif
!-----------------------------------------------------------------------
! calculate pChl, (used in photoadapt., GD98)
! 2.3 max value of thetaN (Chl/N ratio) (mg Chl/mmol N)
! GD 98 Chl. synth. term
!-----------------------------------------------------------------------
WORK1 = parm_alphaChlsp * thetaC_sp * PAR_avg
where (WORK1 > c0)
pChl = thetaN_max_sp * PCphoto_sp / WORK1
photoacc_sp = (pChl * VNC_sp / thetaC_sp) * spChl_loc
elsewhere
photoacc_sp = c0
end where
!-----------------------------------------------------------------------
! CaCO3 Production, parameterized as function of small phyto production
! decrease CaCO3 as function of nutrient limitation decrease CaCO3 prod
! at low temperatures increase CaCO3 prod under bloom conditions
! maximum calcification rate is 40% of primary production
!-----------------------------------------------------------------------
CaCO3_PROD = f_prod_sp_CaCO3 * photoC_sp
CaCO3_PROD = CaCO3_PROD * f_nut
where (TEMP < CaCO3_temp_thres1) &
CaCO3_PROD = CaCO3_PROD * max((TEMP-CaCO3_temp_thres2), c0) / &
(CaCO3_temp_thres1-CaCO3_temp_thres2)
where (spC_loc > CaCO3_sp_thres) &
CaCO3_PROD = min((CaCO3_PROD*spC_loc/CaCO3_sp_thres), &
(f_photosp_CaCO3*photoC_sp))
if (tavg_requested(tavg_CaCO3_form)) then
call accumulate_tavg_field(CaCO3_PROD, tavg_CaCO3_form,bid,k)
endif
if (tavg_requested(tavg_CaCO3_form_zint)) then
if (partial_bottom_cells) then
WORK1 = DZT(:,:,k,bid) * CaCO3_PROD
else
WORK1 = dz(k) * CaCO3_PROD
endif
call accumulate_tavg_field(WORK1, tavg_CaCO3_form_zint,bid,k)
endif
!-----------------------------------------------------------------------
! Relative uptake rates for diatoms nitrate is VNO3, ammonium is VNH4
!-----------------------------------------------------------------------
VNO3_diat = (NO3_loc / parm_diat_kNO3) / &
(c1 + (NO3_loc / parm_diat_kNO3) + (NH4_loc / parm_diat_kNH4))
VNH4_diat = (NH4_loc / parm_diat_kNH4) / &
(c1 + (NO3_loc / parm_diat_kNO3) + (NH4_loc / parm_diat_kNH4))
VNtot_diat = VNO3_diat + VNH4_diat
if (tavg_requested(tavg_diat_N_lim)) then
call accumulate_tavg_field(VNtot_diat, tavg_diat_N_lim,bid,k)
endif
!-----------------------------------------------------------------------
! get relative Fe uptake by diatoms
! get relative P uptake rates for diatoms
! get relative SiO3 uptake rate for diatoms
!-----------------------------------------------------------------------
VFeC_diat = Fe_loc / (Fe_loc + parm_diat_kFe)
if (tavg_requested(tavg_diat_Fe_lim)) then
call accumulate_tavg_field(VFeC_diat, tavg_diat_Fe_lim,bid,k)
endif
VPO4_diat = PO4_loc / (PO4_loc + parm_diat_kPO4)
if (tavg_requested(tavg_diat_PO4_lim)) then
call accumulate_tavg_field(VPO4_diat, tavg_diat_PO4_lim,bid,k)
endif
VSiO3_diat = SiO3_loc / (SiO3_loc + parm_diat_kSiO3)
if (tavg_requested(tavg_diat_SiO3_lim)) then
call accumulate_tavg_field(VSiO3_diat, tavg_diat_SiO3_lim,bid,k)
endif
!-----------------------------------------------------------------------
! Diatom carbon fixation and photoadapt.
! determine diatom nutrient limitation factor for carbon fixation
!-----------------------------------------------------------------------
f_nut = min(VNtot_diat, VFeC_diat)
f_nut = min(f_nut, VSiO3_diat)
f_nut = min(f_nut, VPO4_diat)
!-----------------------------------------------------------------------
! get diatom photosynth. rate, phyto C biomass change, photoadapt
!-----------------------------------------------------------------------
PCmax = PCref * f_nut * Tfunc
light_lim = (c1 - exp((-c1 * parm_alphaChl * thetaC_diat * PAR_avg) / &
(PCmax + epsTinv)))
PCphoto_diat = PCmax * light_lim
if (tavg_requested(tavg_diat_light_lim)) then
call accumulate_tavg_field(light_lim, tavg_diat_light_lim,bid,k)
endif
photoC_diat = PCphoto_diat * diatC_loc
!-----------------------------------------------------------------------
! Get nutrient uptake by diatoms based on C fixation
!-----------------------------------------------------------------------
where (VNtot_diat > c0)
NO3_V_diat = (VNO3_diat / VNtot_diat) * photoC_diat * Q
NH4_V_diat = (VNH4_diat / VNtot_diat) * photoC_diat * Q
VNC_diat = PCphoto_diat * Q
elsewhere
NO3_V_diat = c0
NH4_V_diat = c0
VNC_diat = c0
photoC_diat = c0
end where
photoFe_diat = photoC_diat * gQfe_diat
photoSi_diat = photoC_diat * gQsi
if (tavg_requested(tavg_photoFe_diat)) then
call accumulate_tavg_field(photoFe_diat, tavg_photoFe_diat,bid,k)
endif
if (tavg_requested(tavg_bSi_form)) then
call accumulate_tavg_field(photoSi_diat, tavg_bSi_form,bid,k)
endif
!-----------------------------------------------------------------------
! calculate pChl, (used in photoadapt., GD98)
! 3.0 max value of thetaN (Chl/N ratio) (mg Chl/mmol N)
! GD 98 Chl. synth. term
!-----------------------------------------------------------------------
WORK1 = parm_alphaChl * thetaC_diat * PAR_avg
where (WORK1 > c0)
pChl = thetaN_max_diat * PCphoto_diat / WORK1
photoacc_diat = (pChl * VNC_diat / thetaC_diat) * diatChl_loc
elsewhere
photoacc_diat = c0
end where
!-----------------------------------------------------------------------
! get relative Fe uptake by diazotrophs
! get relative P uptake rates for diazotrophs
!-----------------------------------------------------------------------
Vfec_diaz = Fe_loc/(Fe_loc + diaz_kFe)
if (tavg_requested(tavg_diaz_Fe_lim)) then
call accumulate_tavg_field(Vfec_diaz, tavg_diaz_Fe_lim,bid,k)
endif
Vpo4_diaz = PO4_loc / (PO4_loc + diaz_kPO4)
if (tavg_requested(tavg_diaz_P_lim)) then
call accumulate_tavg_field(Vpo4_diaz, tavg_diaz_P_lim,bid,k)
endif
f_nut = min(Vpo4_diaz, Vfec_diaz)
!-----------------------------------------------------------------------
! get diazotroph photosynth. rate, phyto C biomass change
!-----------------------------------------------------------------------
PCmax = PCrefDiaz * f_nut * Tfunc
where (TEMP < diaz_temp_thres) PCmax = c0
light_lim = (c1 - exp((-c1 * parm_alphaDiaz * thetaC_diaz * PAR_avg) / &
(PCmax + epsTinv)))
PCphoto_diaz = PCmax * light_lim
if (tavg_requested(tavg_diaz_light_lim)) then
call accumulate_tavg_field(light_lim, tavg_diaz_light_lim,bid,k)
endif
photoC_diaz = PCphoto_diaz * diazC_loc
!-----------------------------------------------------------------------
! Get Fe and po4 uptake by diazotrophs based on C fixation
!-----------------------------------------------------------------------
po4_V_diaz = photoC_diaz * Qp_diaz
photoFe_diaz = photoC_diaz * gQfe_diaz
if (tavg_requested(tavg_photoFe_diaz)) then
call accumulate_tavg_field(photoFe_diaz, tavg_photoFe_diaz,bid,k)
endif
!-----------------------------------------------------------------------
! Relative uptake rates for diazotrophs nitrate is VNO3, ammonium is VNH4
!-----------------------------------------------------------------------
VNO3_diaz = (NO3_loc / diaz_kNO3) / &
(c1 + (NO3_loc / diaz_kNO3) + (NH4_loc / diaz_kNH4))
VNH4_diaz = (NH4_loc / diaz_kNH4) / &
(c1 + (NO3_loc / diaz_kNO3) + (NH4_loc / diaz_kNH4))
VNtot_diaz = c1
!------------------------------------------------------------------------
! Compute inoranic N uptake by diazotrophs.
!------------------------------------------------------------------------
WORK1 = photoC_diaz * Q
photoNO3_diaz = VNO3_diaz / VNtot_diaz *WORK1
photoNH4_diaz = VNH4_diaz / VNtot_diaz *WORK1
!-----------------------------------------------------------------------
! Get N fixation by diazotrophs based on C fixation,
! Diazotrophs fix more than they need then 30% is excreted
!-----------------------------------------------------------------------
diaz_Nfix = (WORK1 * r_Nfix_photo) - photoNO3_diaz - photoNH4_diaz
diaz_Nexcrete = diaz_Nfix + photoNO3_diaz + photoNH4_diaz - WORK1
Vnc_diaz = PCphoto_diaz * Q
if (tavg_requested(tavg_diaz_Nfix)) then
call accumulate_tavg_field(diaz_Nfix, tavg_diaz_Nfix,bid,k)
endif
!-----------------------------------------------------------------------
! calculate pChl, (used in photoadapt., GD98)
! 3.4 max value of thetaN (Chl/N ratio) (mg Chl/mmol N)
! GD 98 Chl. synth. term
!-----------------------------------------------------------------------
WORK1 = parm_alphaDiaz * thetaC_diaz * PAR_avg
where (WORK1 > c0)
pChl = thetaN_max_diaz * PCphoto_diaz / WORK1
photoacc_diaz = (pChl * Vnc_diaz / thetaC_diaz) * diazChl_loc
elsewhere
photoacc_diaz = c0
end where
!-----------------------------------------------------------------------
! CALCULATE GRAZING AND OTHER MORT
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! calculate the loss threshold interpolation factor
!-----------------------------------------------------------------------
if (zt(k) > thres_z1) then
if (zt(k) < thres_z2) then
f_loss_thres = (thres_z2 - zt(k))/(thres_z2 - thres_z1)
else
f_loss_thres = c0
endif
else
f_loss_thres = c1
endif
!-----------------------------------------------------------------------
! 0.001 small phytoplankton threshold C concentration (mmol C/m^3)
! get small phyto loss(in C units)
! small phyto agg loss
! get grazing rate (graze_sp) on small phyto (in C units)
!-----------------------------------------------------------------------
C_loss_thres = f_loss_thres * loss_thres_sp
Pprime = max(spC_loc - C_loss_thres, c0)
sp_loss = sp_mort * Pprime * Tfunc
sp_agg = min((sp_agg_rate_max * dps) * Pprime, sp_mort2 * Pprime * Pprime)
reduceV = Pprime * Pprime
graze_sp = z_umax * zooC_loc * (reduceV / (reduceV + parm_z_grz))
!-----------------------------------------------------------------------
! routing of graze_sp & sp_loss
! sp_agg all goes to POC
! currently assumes that 33% of grazed caco3 is remineralized
! if z_ingest ever changes, coefficients on routing grazed sp must change!
! min.%C routed to POC from grazing for ballast requirements = 0.4 * Qcaco3
! min.%C routed from sp_loss = 0.59 * QCaCO3, or P_CaCO3%rho
!-----------------------------------------------------------------------
graze_sp_zoo = z_ingest * graze_sp
graze_sp_poc = graze_sp * max((caco3_poc_min * QCaCO3), &
min(max(0.05_r8,(spc_poc_fac * Pprime)),&
f_graze_sp_poc_lim))
graze_sp_doc = f_graze_sp_doc * graze_sp
graze_sp_dic = f_graze_sp_dic * graze_sp - graze_sp_poc
sp_loss_poc = QCaCO3 * sp_loss
sp_loss_doc = (c1 - parm_labile_ratio) * (sp_loss - sp_loss_poc)
sp_loss_dic = parm_labile_ratio * (sp_loss - sp_loss_poc)
!-----------------------------------------------------------------------
! 0.01 small diatom threshold C concentration (mmol C/m^3)
! get diatom loss(in C units)
! Diatom agg loss, min. 1%/day
! get grazing rate (graze_diat) on diatoms (in C units)
!-----------------------------------------------------------------------
C_loss_thres = f_loss_thres * loss_thres_diat
Pprime = max(diatC_loc - C_loss_thres, c0)
diat_loss = diat_mort * Pprime * Tfunc
diat_agg = min((diat_agg_rate_max * dps) * Pprime, diat_mort2 * Pprime * Pprime)
diat_agg = max((diat_agg_rate_min * dps) * Pprime, diat_agg)
!-----------------------------------------------------------------------
! Lower z_grz term for diatoms and diazotrophs, larger, more mobile predators
!-----------------------------------------------------------------------
reduceV = Pprime * Pprime
graze_diat = diat_umax * zooC_loc * &
(reduceV / (reduceV + 0.7_r8))
!-----------------------------------------------------------------------
! routing of graze_diat & diat_loss
! diat_agg all goes to POC
! NOTE: if z_ingest is changed, coeff.s for poc,doc and dic must change!
!-----------------------------------------------------------------------
graze_diat_zoo = z_ingest * graze_diat
graze_diat_poc = f_graze_diat_poc * graze_diat
graze_diat_doc = f_graze_diat_doc * graze_diat
graze_diat_dic = f_graze_diat_dic * graze_diat
diat_loss_poc = f_diat_loss_poc * diat_loss
diat_loss_doc = (c1-parm_labile_ratio) * f_diat_loss_dc * diat_loss
diat_loss_dic = parm_labile_ratio * f_diat_loss_dc * diat_loss
!-----------------------------------------------------------------------
! 0.03 small diazotroph threshold C concentration (mmol C/m^3)
! Lower value used at temperatures < 16 deg. C, negligible biomass
! get diazotroph loss(in C units)
! get grazing rate (graze_diaz) on diazotrophs (in C units)
! no aggregation loss for diazotrophs
!-----------------------------------------------------------------------
C_loss_diaz = f_loss_thres * loss_thres_diaz
where (TEMP .lt. diaz_temp_thres) &
C_loss_diaz = f_loss_thres * loss_thres_diaz2
Pprime = max(diazC_loc - C_loss_diaz, c0)
diaz_loss = diaz_mort * Pprime * Tfunc
reduceV = Pprime * Pprime
graze_diaz = diaz_umax * zooC_loc * (reduceV / (reduceV + parm_z_grz))
!-----------------------------------------------------------------------
! routing of graze_diaz & diaz_loss
!-----------------------------------------------------------------------
graze_diaz_zoo = f_graze_diaz_zoo * graze_diaz
graze_diaz_poc = f_graze_diaz_poc * graze_diaz
graze_diaz_doc = f_graze_diaz_doc * graze_diaz
graze_diaz_dic = f_graze_diaz_dic * graze_diaz
diaz_loss_poc = f_diaz_loss_poc * diaz_loss
diaz_loss_doc = (c1 - parm_labile_ratio) * (diaz_loss - diaz_loss_poc)
diaz_loss_dic = parm_labile_ratio * (diaz_loss - diaz_loss_poc)
!-----------------------------------------------------------------------
! Note as diazotrophs have different Qp, we must route enough P into zoopl
! and sinking detritus pools to fill their fixed p/C ratios. The left over
! P (remaining_diazP) is split between DOP and DIP pools
!-----------------------------------------------------------------------
remaining_diazP = ((graze_diaz + diaz_loss) * Qp_diaz) - &
((graze_diaz_poc+graze_diaz_zoo) * Qp)
diaz_loss_dop = (c1 - parm_labile_ratio) * remaining_diazP
diaz_loss_dip = parm_labile_ratio * remaining_diazP
!-----------------------------------------------------------------------
! get fractional factor for routing of zoo losses, based on food supply
! more material is routed to large detrital pool when diatoms eaten
!-----------------------------------------------------------------------
f_zoo_detr = (f_diat_zoo_detr * (graze_diat + epsC * epsTinv) + &
f_sp_zoo_detr * (graze_sp + epsC * epsTinv) + &
f_diaz_zoo_detr * (graze_diaz + epsC * epsTinv)) / &
(graze_diat + graze_sp + graze_diaz + c3 * epsC * epsTinv)
!-----------------------------------------------------------------------
! 0.01 small zoo threshold C concentration (mmol C/m^3)
! zoo losses
!-----------------------------------------------------------------------
C_loss_thres = f_loss_thres * loss_thres_zoo
Zprime = max(zooC_loc - C_loss_thres, c0)
zoo_loss = z_mort2 * Zprime**1.4_r8 + z_mort * Zprime
zoo_loss_doc = (c1 - parm_labile_ratio) * (c1 - f_zoo_detr) * zoo_loss
zoo_loss_dic = parm_labile_ratio * (c1 - f_zoo_detr) * zoo_loss
!-----------------------------------------------------------------------
! compute terms for DOM
!-----------------------------------------------------------------------
DOC_prod = sp_loss_doc + graze_sp_doc + zoo_loss_doc + diat_loss_doc &
+ graze_diat_doc + diaz_loss_doc + graze_diaz_doc
DON_prod = DOC_prod * Q
DOP_prod = (sp_loss_doc + graze_sp_doc + zoo_loss_doc + diat_loss_doc &
+ graze_diat_doc) * Qp + diaz_loss_dop
DOFe_prod = (zoo_loss_doc * Qfe_zoo) &
+ (Qfe_sp * (graze_sp_doc + sp_loss_doc)) &
+ (Qfe_diat * (graze_diat_doc + diat_loss_doc)) &
+ (Qfe_diaz * (graze_diaz_doc + diaz_loss_doc))
DOC_remin = DOC_loc * DOM_remin
DON_remin = DON_loc * DOM_remin
DOFe_remin = DOFe_loc * DOM_remin
DOP_remin = DOP_loc * DOM_remin
!-----------------------------------------------------------------------
! large detritus C
!-----------------------------------------------------------------------
POC%prod(:,:,bid) = sp_agg + graze_sp_poc + sp_loss_poc + &
f_zoo_detr * zoo_loss + diat_loss_poc + &
diat_agg + graze_diat_poc + graze_diaz_poc
!-----------------------------------------------------------------------
! large detrital CaCO3
! 33% of CaCO3 is remin when phyto are grazed
!-----------------------------------------------------------------------
P_CaCO3%prod(:,:,bid) = ((c1 - f_graze_CaCO3_REMIN) * graze_sp + &
sp_loss + sp_agg) * QCaCO3
!-----------------------------------------------------------------------
! large detritus SiO2
! grazed diatom SiO2, 60% is remineralized
!-----------------------------------------------------------------------
P_SiO2%prod(:,:,bid) = ((c1 - f_graze_si_remin) * graze_diat + &
diat_agg + f_diat_loss_poc * diat_loss) * Qsi
dust%prod(:,:,bid) = c0
!-----------------------------------------------------------------------
! Compute iron scavenging :
! 1) compute in terms of loss per year per unit iron (%/year/fe)
! 2) scale by sinking POMx6 + Dust + bSi + CaCO3 flux
! 3) increase scavenging at higher iron (>0.6nM)
! 4) convert to net loss per second
!-----------------------------------------------------------------------
Fe_scavenge_rate = Fe_scavenge_rate0
Fe_scavenge_rate = Fe_scavenge_rate * &
((POC%sflux_out(:,:,bid) + POC%hflux_out(:,:,bid)) * 72.06_r8 + &
(P_CaCO3%sflux_out(:,:,bid) + P_CaCO3%hflux_out(:,:,bid)) * P_CaCO3%mass + &
(P_SiO2%sflux_out(:,:,bid)+P_SiO2%hflux_out(:,:,bid))*P_SiO2%mass + &
(dust%sflux_out(:,:,bid) + dust%hflux_out(:,:,bid)) * dust_fescav_scale)
where (Fe_loc > Fe_scavenge_thres1) &
Fe_scavenge_rate = Fe_scavenge_rate + &
(Fe_loc - Fe_scavenge_thres1) * fe_max_scale2
Fe_scavenge = yps * Fe_loc * Fe_scavenge_rate
P_iron%prod(:,:,bid) = &
((sp_agg + graze_sp_poc + sp_loss_poc) * Qfe_sp) &
+ (zoo_loss * f_zoo_detr * Qfe_zoo) &
+ ((diat_agg + graze_diat_poc + diat_loss_poc) * Qfe_diat) &
+ (graze_diaz_poc * Qfe_diaz) + (f_fescav_P_iron * Fe_scavenge)
call compute_particulate_terms(k, POC, P_CaCO3, P_SiO2, dust, P_iron, &
QA_dust_def, TEMP, O2_loc, this_block)
!-----------------------------------------------------------------------
! nitrate & ammonium
! nitrification in low light
! use exponential decay of PAR across model level to compute taper factor
!-----------------------------------------------------------------------
if (lrest_no3) then
if (lnutr_variable_restore) then
RESTORE = NUTR_RESTORE_RTAU(:,:,bid) * &
merge((NO3_CLIM(:,:,k,bid) - NO3_loc), &
c0, k <= NUTR_RESTORE_MAX_LEVEL(:,:,bid))
else
RESTORE = (NO3_CLIM(:,:,k,bid) - NO3_loc) * nutr_rest_time_inv(k)
endif
else
RESTORE = c0
endif
if (tavg_requested(tavg_NO3_RESTORE)) then
call accumulate_tavg_field(RESTORE, tavg_NO3_RESTORE,bid,k)
endif
where (PAR_out(:,:,bid) < parm_nitrif_par_lim)
NITRIF = parm_kappa_nitrif * NH4_loc
where (PAR_in > parm_nitrif_par_lim)
NITRIF = NITRIF * log(PAR_out(:,:,bid) / parm_nitrif_par_lim) / (-KPARdz)
end where
elsewhere
NITRIF = c0
end where
if (tavg_requested(tavg_NITRIF)) then
call accumulate_tavg_field(NITRIF, tavg_NITRIF,bid,k)
endif
!-----------------------------------------------------------------------
! Compute denitrification under low O2 conditions
!-----------------------------------------------------------------------
WORK1 = ((parm_o2_min+parm_o2_min_delta) - O2_loc) / parm_o2_min_delta
WORK1 = min(max(WORK1,c0),c1)
where (NO3_loc .lt. parm_no3_min)
WORK1 = WORK1 * (NO3_loc / parm_no3_min)
end where
DENITRIF = WORK1 * (DOC_remin + POC%remin(:,:,bid)) / denitrif_C_N
if (tavg_requested(tavg_DENITRIF)) then
call accumulate_tavg_field(DENITRIF, tavg_DENITRIF,bid,k)
endif
!-----------------------------------------------------------------------
! nitrate & ammonium
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,no3_ind) = RESTORE + NITRIF - NO3_V_diat - &
NO3_V_sp - DENITRIF - photoNO3_diaz
DTRACER_MODULE(:,:,nh4_ind) = -NH4_V_diat - NH4_V_sp - NITRIF + &
Q * (zoo_loss_dic + sp_loss_dic + graze_sp_dic + diat_loss_dic + &
graze_diat_dic + POC%remin(:,:,bid) + diaz_loss_dic + &
graze_diaz_dic) + DON_remin + diaz_Nexcrete - photoNH4_diaz
!-----------------------------------------------------------------------
! dissolved iron
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,fe_ind) = P_iron%remin(:,:,bid) &
+ (Qfe_zoo * zoo_loss_dic) + DOFe_remin &
+ (Qfe_sp * (sp_loss_dic + graze_sp_dic)) &
+ (Qfe_diat * (diat_loss_dic + graze_diat_dic)) &
+ (Qfe_diaz * (diaz_loss_dic + graze_diaz_dic)) &
+ graze_diaz_zoo * (Qfe_diaz-Qfe_zoo) &
+ graze_diat_zoo * (Qfe_diat-Qfe_zoo) &
+ graze_sp_zoo * (Qfe_sp-Qfe_zoo) &
- photoFe_sp - photoFe_diat - photoFe_diaz - Fe_scavenge
!-----------------------------------------------------------------------
! dissolved SiO3
!-----------------------------------------------------------------------
if (lrest_sio3) then
if (lnutr_variable_restore) then
RESTORE = NUTR_RESTORE_RTAU(:,:,bid) * &
merge((SiO3_CLIM(:,:,k,bid) - SiO3_loc), &
c0, k <= NUTR_RESTORE_MAX_LEVEL(:,:,bid))
else
RESTORE = (SiO3_CLIM(:,:,k,bid) - SiO3_loc) * nutr_rest_time_inv(k)
endif
else
RESTORE = c0
endif
if (tavg_requested(tavg_SiO3_RESTORE)) then
call accumulate_tavg_field(RESTORE, tavg_SiO3_RESTORE,bid,k)
endif
DTRACER_MODULE(:,:,sio3_ind) = RESTORE + P_SiO2%remin(:,:,bid) + &
Qsi * (f_graze_si_remin * graze_diat +f_diat_loss_dc * diat_loss) &
- photoSi_diat
!-----------------------------------------------------------------------
! phosphate
!-----------------------------------------------------------------------
if (lrest_po4) then
if (lnutr_variable_restore) then
RESTORE = NUTR_RESTORE_RTAU(:,:,bid) * &
merge((PO4_CLIM(:,:,k,bid) - PO4_loc), &
c0, k <= NUTR_RESTORE_MAX_LEVEL(:,:,bid))
else
RESTORE = (PO4_CLIM(:,:,k,bid) - PO4_loc) * nutr_rest_time_inv(k)
endif
else
RESTORE = c0
endif
if (tavg_requested(tavg_PO4_RESTORE)) then
call accumulate_tavg_field(RESTORE, tavg_PO4_RESTORE,bid,k)
endif
DTRACER_MODULE(:,:,po4_ind) = RESTORE + (Qp * (POC%remin(:,:,bid) + &
zoo_loss_dic + sp_loss_dic + graze_sp_dic + diat_loss_dic + &
graze_diat_dic - photoC_diat)) &
+ DOP_remin + diaz_loss_dip - po4_V_sp - po4_V_diaz
!-----------------------------------------------------------------------
! small phyto Carbon
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,spC_ind) = photoC_sp - graze_sp - sp_loss - sp_agg
!-----------------------------------------------------------------------
! small phyto Chlorophyll
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,spChl_ind) = photoacc_sp - &
thetaC_sp * (graze_sp + sp_loss + sp_agg)
!-----------------------------------------------------------------------
! small phytoplankton CaCO3
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,spCaCO3_ind) = CaCO3_PROD - &
(graze_sp + sp_loss + sp_agg) * QCaCO3
!-----------------------------------------------------------------------
! diatom Carbon
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,diatC_ind) = photoC_diat - graze_diat - &
diat_loss - diat_agg
!-----------------------------------------------------------------------
! diatom Chlorophyll
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,diatChl_ind) = photoacc_diat - &
thetaC_diat * (graze_diat + diat_loss + diat_agg)
!-----------------------------------------------------------------------
! zoo Carbon
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,zooC_ind) = graze_sp_zoo + graze_diat_zoo &
+ graze_diaz_zoo - zoo_loss
!-----------------------------------------------------------------------
! dissolved organic Matter
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,doc_ind) = DOC_prod - DOC_remin
DTRACER_MODULE(:,:,don_ind) = DON_prod - DON_remin
DTRACER_MODULE(:,:,dop_ind) = DOP_prod - DOP_remin
DTRACER_MODULE(:,:,dofe_ind) = DOFe_prod - DOFe_remin
!-----------------------------------------------------------------------
! small phyto Fe
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,spFe_ind) = photoFe_sp &
- (Qfe_sp * (graze_sp+sp_loss+sp_agg))
!-----------------------------------------------------------------------
! Diatom Fe
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,diatFe_ind) = photoFe_diat &
- (Qfe_diat * (graze_diat+diat_loss+diat_agg))
!-----------------------------------------------------------------------
! Diatom Si
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,diatSi_ind) = photoSi_diat &
- (Qsi * (graze_diat+diat_loss+diat_agg))
!-----------------------------------------------------------------------
! Diazotroph C
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,diazC_ind) = photoC_diaz - graze_diaz - diaz_loss
!-----------------------------------------------------------------------
! diazotroph Chlorophyll
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,diazChl_ind) = photoacc_diaz - &
thetaC_diaz * (graze_diaz + diaz_loss)
!-----------------------------------------------------------------------
! Diazotroph Fe
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,diazFe_ind) = photoFe_diaz &
- (Qfe_diaz * (graze_diaz + diaz_loss))
!-----------------------------------------------------------------------
! dissolved inorganic Carbon
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,dic_ind) = &
DOC_remin + POC%remin(:,:,bid) + P_CaCO3%remin(:,:,bid) + &
f_graze_CaCO3_REMIN * graze_sp * QCaCO3 + zoo_loss_dic + sp_loss_dic + &
graze_sp_dic + diat_loss_dic + graze_diat_dic - photoC_sp - photoC_diat &
- CaCO3_PROD + graze_diaz_dic + diaz_loss_dic - photoC_diaz
!-----------------------------------------------------------------------
! alkalinity
!-----------------------------------------------------------------------
DTRACER_MODULE(:,:,alk_ind) = -DTRACER_MODULE(:,:,no3_ind) + &
DTRACER_MODULE(:,:,nh4_ind) + &
c2 * (P_CaCO3%remin(:,:,bid) + f_graze_CaCO3_REMIN * graze_sp * QCaCO3 &
- CaCO3_PROD)
!-----------------------------------------------------------------------
! oxygen
!-----------------------------------------------------------------------
! DTRACER_MODULE(:,:,o2_ind) = (photoC_sp + photoC_diat) / &
! parm_Red_D_C_O2 + photoC_diaz/parm_Red_D_C_O2_diaz
! WORK1 = (O2_loc - parm_o2_min) / parm_o2_min_delta
! WORK1 = min(max(WORK1,c0),c1)
! DTRACER_MODULE(:,:,o2_ind) = DTRACER_MODULE(:,:,o2_ind) + WORK1 &
! * ((- POC%remin(:,:,bid) - DOC_remin - zoo_loss_dic - sp_loss_dic &
! - graze_sp_dic - diat_loss_dic - graze_diat_dic &
! - graze_diaz_dic - diaz_loss_dic) / parm_Red_P_C_O2)
DTRACER_MODULE(:,:,o2_ind) = c0
O2_PRODUCTION = c0
where (photoC_diat > c0)
O2_PRODUCTION = O2_PRODUCTION + photoC_diat * &
((NO3_V_diat/(NO3_V_diat+NH4_V_diat)) / parm_Red_D_C_O2 + &
(NH4_V_diat/(NO3_V_diat+NH4_V_diat)) / parm_Remin_D_C_O2)
end where
where (photoC_sp > c0)
O2_PRODUCTION = O2_PRODUCTION + photoC_sp * &
((NO3_V_sp/(NO3_V_sp+NH4_V_sp)) / parm_Red_D_C_O2 + &
(NH4_V_sp/(NO3_V_sp+NH4_V_sp)) / parm_Remin_D_C_O2)
end where
where (photoC_diaz > c0)
O2_PRODUCTION = O2_PRODUCTION + photoC_diaz * &
((photoNO3_diaz/(photoNO3_diaz+photoNH4_diaz+diaz_Nfix)) / parm_Red_D_C_O2 + &
(photoNH4_diaz/(photoNO3_diaz+photoNH4_diaz+diaz_Nfix)) / parm_Remin_D_C_O2 + &
(diaz_Nfix/(photoNO3_diaz+photoNH4_diaz+diaz_Nfix)) / parm_Red_D_C_O2_diaz)
end where
WORK1 = (O2_loc - parm_o2_min) / parm_o2_min_delta
WORK1 = min(max(WORK1,c0),c1)
O2_CONSUMPTION = WORK1 * &
((POC%remin(:,:,bid) + DOC_remin + zoo_loss_dic + &
sp_loss_dic + graze_sp_dic + diat_loss_dic + graze_diat_dic + &
graze_diaz_dic + diaz_loss_dic)/ parm_Remin_D_C_O2 + (c2*NITRIF))
DTRACER_MODULE(:,:,o2_ind) = O2_PRODUCTION - O2_CONSUMPTION
!-----------------------------------------------------------------------
! various tavg/history variables
!-----------------------------------------------------------------------
if (k == 1) then
if (tavg_requested(tavg_O2_ZMIN) .or. tavg_requested(tavg_O2_ZMIN_DEPTH)) then
! WORK1 = O2 at this level
! WORK2 = vertical min of O2
! WORK3 = depth of min
kk = 1
WORK1 = p5*(TRACER_MODULE_OLD(:,:,kk,o2_ind) + &
TRACER_MODULE_CUR(:,:,kk,o2_ind))
WORK2 = WORK1
WORK3 = zt(kk)
do kk = 2,km
WORK1 = p5*(TRACER_MODULE_OLD(:,:,kk,o2_ind) + &
TRACER_MODULE_CUR(:,:,kk,o2_ind))
where (kk <= KMT(:,:,bid) .and. (WORK1 < WORK2))
WORK2 = WORK1
WORK3 = zt(kk)
endwhere
end do
if (tavg_requested(tavg_O2_ZMIN)) then
call accumulate_tavg_field(WORK2, tavg_O2_ZMIN,bid,k)
endif
if (tavg_requested(tavg_O2_ZMIN_DEPTH)) then
call accumulate_tavg_field(WORK3, tavg_O2_ZMIN_DEPTH,bid,k)
endif
endif
endif
if (tavg_requested(tavg_O2_PRODUCTION)) then
call accumulate_tavg_field(O2_PRODUCTION, tavg_O2_PRODUCTION,bid,k)
endif
if (tavg_requested(tavg_O2_CONSUMPTION)) then
call accumulate_tavg_field(O2_CONSUMPTION, tavg_O2_CONSUMPTION,bid,k)
endif
if (tavg_requested(tavg_AOU)) then
WORK1 = O2SAT(TEMP, SALT, &
LAND_MASK(:,:,bid) .and. (k .le. KMT(:,:,bid)))
WORK1 = WORK1 - p5*(TRACER_MODULE_OLD(:,:,k,o2_ind) + &
TRACER_MODULE_CUR(:,:,k,o2_ind))
call accumulate_tavg_field(WORK1, tavg_AOU,bid,k)
endif
if (tavg_requested(tavg_PAR_avg)) then
call accumulate_tavg_field(PAR_avg, tavg_PAR_avg,bid,k)
endif
if (tavg_requested(tavg_graze_sp)) then
call accumulate_tavg_field(graze_sp, tavg_graze_sp,bid,k)
endif
if (tavg_requested(tavg_graze_diat)) then
call accumulate_tavg_field(graze_diat, tavg_graze_diat,bid,k)
endif
if (tavg_requested(tavg_graze_diaz)) then
call accumulate_tavg_field(graze_diaz, tavg_graze_diaz,bid,k)
endif
if (tavg_requested(tavg_graze_TOT)) then
WORK1 = graze_sp + graze_diat + graze_diaz
call accumulate_tavg_field(WORK1, tavg_graze_TOT,bid,k)
endif
if (tavg_requested(tavg_sp_loss)) then
call accumulate_tavg_field(sp_loss, tavg_sp_loss,bid,k)
endif
if (tavg_requested(tavg_diat_loss)) then
call accumulate_tavg_field(diat_loss, tavg_diat_loss,bid,k)
endif
if (tavg_requested(tavg_diaz_loss)) then
call accumulate_tavg_field(diaz_loss, tavg_diaz_loss,bid,k)
endif
if (tavg_requested(tavg_zoo_loss)) then
call accumulate_tavg_field(zoo_loss, tavg_zoo_loss,bid,k)
endif
if (tavg_requested(tavg_sp_agg)) then
call accumulate_tavg_field(sp_agg, tavg_sp_agg,bid,k)
endif
if (tavg_requested(tavg_diat_agg)) then
call accumulate_tavg_field(diat_agg, tavg_diat_agg,bid,k)
endif
if (tavg_requested(tavg_photoC_sp)) then
call accumulate_tavg_field(photoC_sp, tavg_photoC_sp,bid,k)
endif
if (tavg_requested(tavg_photoC_diat)) then
call accumulate_tavg_field(photoC_diat, tavg_photoC_diat,bid,k)
endif
if (tavg_requested(tavg_photoC_diaz)) then
call accumulate_tavg_field(photoC_diaz, tavg_photoC_diaz,bid,k)
endif
if (tavg_requested(tavg_photoC_TOT)) then
WORK1 = photoC_sp + photoC_diat + photoC_diaz
call accumulate_tavg_field(WORK1, tavg_photoC_TOT,bid,k)
endif
if (tavg_requested(tavg_photoC_sp_zint)) then
if (partial_bottom_cells) then
WORK1 = DZT(:,:,k,bid) * photoC_sp
else
WORK1 = dz(k) * photoC_sp
endif
call accumulate_tavg_field(WORK1, tavg_photoC_sp_zint,bid,k)
endif
if (tavg_requested(tavg_photoC_diat_zint)) then
if (partial_bottom_cells) then
WORK1 = DZT(:,:,k,bid) * photoC_diat
else
WORK1 = dz(k) * photoC_diat
endif
call accumulate_tavg_field(WORK1, tavg_photoC_diat_zint,bid,k)
endif
if (tavg_requested(tavg_photoC_diaz_zint)) then
if (partial_bottom_cells) then
WORK1 = DZT(:,:,k,bid) * photoC_diaz
else
WORK1 = dz(k) * photoC_diaz
endif
call accumulate_tavg_field(WORK1, tavg_photoC_diaz_zint,bid,k)
endif
if (tavg_requested(tavg_photoC_TOT_zint)) then
if (partial_bottom_cells) then
WORK1 = DZT(:,:,k,bid) * (photoC_sp + photoC_diat + photoC_diaz)
else
WORK1 = dz(k) * (photoC_sp + photoC_diat + photoC_diaz)
endif
call accumulate_tavg_field(WORK1, tavg_photoC_TOT_zint,bid,k)
endif
if (tavg_requested(tavg_photoC_NO3_sp) .or. &
tavg_requested(tavg_photoC_NO3_sp_zint)) then
where (VNtot_sp > c0)
WORK1 = (VNO3_sp / VNtot_sp) * photoC_sp
elsewhere
WORK1 = c0
end where
if (tavg_requested(tavg_photoC_NO3_sp)) &
call accumulate_tavg_field(WORK1, tavg_photoC_NO3_sp,bid,k)
if (tavg_requested(tavg_photoC_NO3_sp_zint)) then
if (partial_bottom_cells) then
WORK1 = DZT(:,:,k,bid) * WORK1
else
WORK1 = dz(k) * WORK1
endif
call accumulate_tavg_field(WORK1, tavg_photoC_NO3_sp_zint,bid,k)
endif
endif
if (tavg_requested(tavg_photoC_NO3_diat) .or. &
tavg_requested(tavg_photoC_NO3_diat_zint)) then
where (VNtot_diat > c0)
WORK1 = (VNO3_diat / VNtot_diat) * photoC_diat
elsewhere
WORK1 = c0
end where
if (tavg_requested(tavg_photoC_NO3_diat)) &
call accumulate_tavg_field(WORK1, tavg_photoC_NO3_diat,bid,k)
if (tavg_requested(tavg_photoC_NO3_diat_zint)) then
if (partial_bottom_cells) then
WORK1 = DZT(:,:,k,bid) * WORK1
else
WORK1 = dz(k) * WORK1
endif
call accumulate_tavg_field(WORK1, tavg_photoC_NO3_diat_zint,bid,k)
endif
endif
if (tavg_requested(tavg_photoC_NO3_diaz) .or. &
tavg_requested(tavg_photoC_NO3_diaz_zint)) then
where (VNtot_diaz > c0)
WORK1 = (VNO3_diaz / VNtot_diaz) * photoC_diaz
elsewhere
WORK1 = c0
end where
if (tavg_requested(tavg_photoC_NO3_diaz)) &
call accumulate_tavg_field(WORK1, tavg_photoC_NO3_diaz,bid,k)
if (tavg_requested(tavg_photoC_NO3_diaz_zint)) then
if (partial_bottom_cells) then
WORK1 = DZT(:,:,k,bid) * WORK1
else
WORK1 = dz(k) * WORK1
endif
call accumulate_tavg_field(WORK1, tavg_photoC_NO3_diaz_zint,bid,k)
endif
endif
if (tavg_requested(tavg_photoC_NO3_TOT) .or. &
tavg_requested(tavg_photoC_NO3_TOT_zint)) then
where (VNtot_sp > c0)
WORK1 = WORK1 + (VNO3_sp / VNtot_sp) * photoC_sp
elsewhere
WORK1 = c0
end where
where (VNtot_diat > c0)
WORK1 = WORK1 + (VNO3_diat / VNtot_diat) * photoC_diat
end where
where (VNtot_diaz > c0)
WORK1 = WORK1 + (VNO3_diaz / VNtot_diaz) * photoC_diaz
end where
if (tavg_requested(tavg_photoC_NO3_TOT)) &
call accumulate_tavg_field(WORK1, tavg_photoC_NO3_TOT,bid,k)
if (tavg_requested(tavg_photoC_NO3_TOT_zint)) then
if (partial_bottom_cells) then
WORK1 = DZT(:,:,k,bid) * WORK1
else
WORK1 = dz(k) * WORK1
endif
call accumulate_tavg_field(WORK1, tavg_photoC_NO3_TOT_zint,bid,k)
endif
endif
if (tavg_requested(tavg_DOC_prod)) then
call accumulate_tavg_field(DOC_prod, tavg_DOC_prod,bid,k)
endif
if (tavg_requested(tavg_DOC_remin)) then
call accumulate_tavg_field(DOC_remin, tavg_DOC_remin,bid,k)
endif
if (tavg_requested(tavg_DON_prod)) then
call accumulate_tavg_field(DON_prod, tavg_DON_prod,bid,k)
endif
if (tavg_requested(tavg_DON_remin)) then
call accumulate_tavg_field(DON_remin, tavg_DON_remin,bid,k)
endif
if (tavg_requested(tavg_DOP_prod)) then
call accumulate_tavg_field(DOP_prod, tavg_DOP_prod,bid,k)
endif
if (tavg_requested(tavg_DOP_remin)) then
call accumulate_tavg_field(DOP_remin, tavg_DOP_remin,bid,k)
endif
if (tavg_requested(tavg_DOFe_prod)) then
call accumulate_tavg_field(DOFe_prod, tavg_DOFe_prod,bid,k)
endif
if (tavg_requested(tavg_DOFe_remin)) then
call accumulate_tavg_field(DOFe_remin, tavg_DOFe_remin,bid,k)
endif
if (tavg_requested(tavg_Fe_scavenge)) then
call accumulate_tavg_field(Fe_scavenge, tavg_Fe_scavenge,bid,k)
endif
if (tavg_requested(tavg_Fe_scavenge_rate)) then
call accumulate_tavg_field(Fe_scavenge_rate, tavg_Fe_scavenge_rate,bid,k)
endif
if (tavg_requested(tavg_CO3) .or. &
tavg_requested(tavg_HCO3) .or. &
tavg_requested(tavg_H2CO3) .or. &
tavg_requested(tavg_pH_3D) .or. &
tavg_requested(tavg_zsatcalc) .or. &
tavg_requested(tavg_zsatarag)) then
where (PH_PREV_3D(:,:,k,bid) /= c0)
WORK1 = PH_PREV_3D(:,:,k,bid) - del_ph
WORK2 = PH_PREV_3D(:,:,k,bid) + del_ph
elsewhere
WORK1 = phlo_init
WORK2 = phhi_init
end where
call timer_start(ecosys_comp_CO3terms_timer, block_id=bid)
call comp_CO3terms(bid, k, LAND_MASK(:,:,bid) .and. k <= KMT(:,:,bid), &
TEMP, SALT, DIC_loc, ALK_loc, PO4_loc, SiO3_loc, &
WORK1, WORK2, WORK3, H2CO3, HCO3, CO3)
call timer_stop(ecosys_comp_CO3terms_timer, block_id=bid)
if (tavg_requested(tavg_CO3)) then
call accumulate_tavg_field(CO3, tavg_CO3,bid,k)
endif
if (tavg_requested(tavg_HCO3)) then
call accumulate_tavg_field(HCO3, tavg_HCO3,bid,k)
endif
if (tavg_requested(tavg_H2CO3)) then
call accumulate_tavg_field(H2CO3, tavg_H2CO3,bid,k)
endif
if (tavg_requested(tavg_pH_3D)) then
call accumulate_tavg_field(WORK3, tavg_pH_3D,bid,k)
endif
PH_PREV_3D(:,:,k,bid) = WORK3
endif
if (tavg_requested(tavg_co3_sat_calc) .or. &
tavg_requested(tavg_zsatcalc) .or. &
tavg_requested(tavg_co3_sat_arag) .or. &
tavg_requested(tavg_zsatarag)) then
call comp_co3_sat_vals(k, LAND_MASK(:,:,bid) .and. k <= KMT(:,:,bid), &
TEMP, SALT, WORK1, WORK2)
if (tavg_requested(tavg_co3_sat_calc)) then
call accumulate_tavg_field(WORK1, tavg_co3_sat_calc,bid,k)
endif
if (tavg_requested(tavg_zsatcalc)) then
if (k == 1) then
! set to -1, i.e. depth not found yet,
! if mask == .true. and surface supersaturated to -1
ZSATCALC(:,:,bid) = merge(-c1, c0, LAND_MASK(:,:,bid) .and. CO3 > WORK1)
else
where (ZSATCALC(:,:,bid) == -c1 .and. CO3 <= WORK1)
ZSATCALC(:,:,bid) = zt(k-1) + (zt(k)-zt(k-1)) * &
CO3_CALC_ANOM_km1(:,:,bid) / (CO3_CALC_ANOM_km1(:,:,bid) - (CO3 - WORK1))
endwhere
where (ZSATCALC(:,:,bid) == -c1 .and. KMT(:,:,bid) == k)
ZSATCALC(:,:,bid) = zw(k)
endwhere
endif
CO3_CALC_ANOM_km1(:,:,bid) = CO3 - WORK1
if (k == km) then
call accumulate_tavg_field(ZSATCALC(:,:,bid), tavg_zsatcalc,bid,k)
endif
endif
if (tavg_requested(tavg_co3_sat_arag)) then
call accumulate_tavg_field(WORK2, tavg_co3_sat_arag,bid,k)
endif
if (tavg_requested(tavg_zsatarag)) then
if (k == 1) then
! set to -1, i.e. depth not found yet,
! if mask == .true. and surface supersaturated to -1
ZSATARAG(:,:,bid) = merge(-c1, c0, LAND_MASK(:,:,bid) .and. CO3 > WORK2)
else
where (ZSATARAG(:,:,bid) == -c1 .and. CO3 <= WORK2)
ZSATARAG(:,:,bid) = zt(k-1) + (zt(k)-zt(k-1)) * &
CO3_ARAG_ANOM_km1(:,:,bid) / (CO3_ARAG_ANOM_km1(:,:,bid) - (CO3 - WORK2))
endwhere
where (ZSATARAG(:,:,bid) == -c1 .and. KMT(:,:,bid) == k)
ZSATARAG(:,:,bid) = zw(k)
endwhere
endif
CO3_ARAG_ANOM_km1(:,:,bid) = CO3 - WORK2
if (k == km) then
call accumulate_tavg_field(ZSATARAG(:,:,bid), tavg_zsatarag,bid,k)
endif
endif
endif
call timer_stop(ecosys_interior_timer, block_id=bid)
!-----------------------------------------------------------------------
!EOC
end subroutine ecosys_set_interior
!***********************************************************************
!BOP
! !IROUTINE: init_particulate_terms
! !INTERFACE:
subroutine init_particulate_terms(POC, P_CaCO3, P_SiO2, dust, P_iron, & 1
QA_dust_def, NET_DUST_IN, this_block)
! !DESCRIPTION:
! Set incoming fluxes (put into outgoing flux for first level usage).
! Set dissolution length, production fraction and mass terms.
!
! The first 6 arguments are intent(inout) in
! order to preserve contents on other blocks.
! !INPUT/OUTPUT PARAMETERS:
type(sinking_particle), intent(inout) :: &
POC, & ! base units = nmol C
P_CaCO3, & ! base units = nmol CaCO3
P_SiO2, & ! base units = nmol SiO2
dust, & ! base units = g
P_iron ! base units = nmol Fe
real (r8), dimension(nx_block,ny_block,max_blocks_clinic), intent(inout) :: &
QA_dust_def ! incoming deficit in the QA(dust) POC flux
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,max_blocks_clinic), intent(in) :: &
NET_DUST_IN ! dust flux
type (block), intent(in) :: &
this_block ! block info for the current block
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
integer (int_kind) :: &
bid ! local_block id
character (char_len) :: &
tracer_data_label ! label for what is being updated
character (char_len), dimension(:), allocatable :: &
tracer_data_names ! short names for input data fields
integer (int_kind), dimension(:), allocatable :: &
tracer_bndy_loc, & ! location and field type for ghost
tracer_bndy_type ! cell updates
!-----------------------------------------------------------------------
! parameters, from Armstrong et al. 2000
!
! July 2002, length scale for excess POC and bSI modified by temperature
! Value given here is at Tref of 30 deg. C, JKM
!-----------------------------------------------------------------------
POC%diss = 21000.0_r8 ! diss. length (cm), modified by TEMP
POC%gamma = c0 ! not used
POC%mass = 12.01_r8 ! molecular weight of POC
POC%rho = c0 ! not used
P_CaCO3%diss = 80000.0_r8 ! diss. length (cm)
P_CaCO3%gamma = 0.55_r8 ! prod frac -> hard subclass
P_CaCO3%mass = 100.09_r8 ! molecular weight of CaCO
P_CaCO3%rho = 0.05_r8 * P_CaCO3%mass / POC%mass
! QA mass ratio for CaCO3
! This ratio is used in ecos_set_interior
P_SiO2%diss = 21000.0_r8 ! diss. length (cm), modified by TEMP
P_SiO2%gamma = 0.25_r8 ! prod frac -> hard subclass
P_SiO2%mass = 60.08_r8 ! molecular weight of SiO2
P_SiO2%rho = 0.05_r8 * P_SiO2%mass / POC%mass
! QA mass ratio for SiO2
dust%diss = 60000.0_r8 ! diss. length (cm)
dust%gamma = 0.97_r8 ! prod frac -> hard subclass
dust%mass = 1.0e9_r8 ! base units are already grams
dust%rho = 0.05_r8 * dust%mass / POC%mass
! QA mass ratio for dust
P_iron%diss = 60000.0_r8 ! diss. length (cm) - not used
P_iron%gamma = c0 ! prod frac -> hard subclass - not used
P_iron%mass = c0 ! not used
P_iron%rho = c0 ! not used
!-----------------------------------------------------------------------
! Set incoming fluxes
!-----------------------------------------------------------------------
bid = this_block%local_id
P_CaCO3%sflux_out(:,:,bid) = c0
P_CaCO3%hflux_out(:,:,bid) = c0
P_SiO2%sflux_out(:,:,bid) = c0
P_SiO2%hflux_out(:,:,bid) = c0
if (dust_flux%has_data) then
dust%sflux_out(:,:,bid) = (c1 - dust%gamma) * NET_DUST_IN(:,:,bid)
dust%hflux_out(:,:,bid) = dust%gamma * NET_DUST_IN(:,:,bid)
else
dust%sflux_out(:,:,bid) = c0
dust%hflux_out(:,:,bid) = c0
endif
P_iron%sflux_out(:,:,bid) = c0
P_iron%hflux_out(:,:,bid) = c0
!-----------------------------------------------------------------------
! Hard POC is QA flux and soft POC is excess POC.
!-----------------------------------------------------------------------
POC%sflux_out(:,:,bid) = c0
POC%hflux_out(:,:,bid) = c0
!-----------------------------------------------------------------------
! Compute initial QA(dust) POC flux deficit.
!-----------------------------------------------------------------------
QA_dust_def(:,:,bid) = dust%rho * &
(dust%sflux_out(:,:,bid) + dust%hflux_out(:,:,bid))
!-----------------------------------------------------------------------
!EOC
end subroutine init_particulate_terms
!***********************************************************************
!BOP
! !IROUTINE: compute_particulate_terms
! !INTERFACE:
subroutine compute_particulate_terms(k, POC, P_CaCO3, P_SiO2, dust, P_iron, & 1,30
QA_dust_def, TEMP, O2_loc, this_block)
! !DESCRIPTION:
! Compute outgoing fluxes and remineralization terms. Assumes that
! production terms have been set. Incoming fluxes are assumed to be the
! outgoing fluxes from the previous level.
!
! It is assumed that there is no production of dust.
!
! Instantaneous remineralization in the bottom cell is implemented by
! setting the outgoing flux to zero.
!
! For POC, the hard subclass is the POC flux qualitatively associated
! with the ballast flux. The soft subclass is the excess POC flux.
!
! Remineralization for the non-iron particulate pools is computing
! by first computing the outgoing flux and then computing the
! remineralization from conservation, i.e.
! flux_in - flux_out + prod * dz - remin * dz == 0.
!
! For iron, remineralization is first computed from POC remineralization
! and then flux_out is computed from conservation. If the resulting
! flux_out is negative or should be zero because of the sea floor, the
! remineralization is adjusted.
! Note: all the sinking iron is in the P_iron%sflux pool, hflux Fe not
! explicitly tracked, it is assumed that total iron remin is
! proportional to total POC remin.
!
! Based upon Armstrong et al. 2000
!
! July 2002, added temperature effect on remin length scale of
! excess POC (all soft POM& Iron) and on SiO2.
! new variable passed into ballast, Tfunc, main Temperature function
! computed in ecosystem routine. scaling factor for dissolution
! of excess POC, Fe, and Bsi now varies with location (f(temperature)).
!
! Added diffusive iron flux from sediments at depths < 1100m,
! based on Johnson et al., 1999, value of 5 umolFe/m2/day,
! this value too high, using 2 umolFe/m2/day here
!
! Allow hard fraction of ballast to remin with long length scale 40,000m
! thus ~ 10% of hard ballast remins over 4000m water column.
!
! Sinking dust flux is decreased by assumed instant solubility/dissolution
! at ocean surface from the parm_Fe_bioavail.
!
! Modified to allow different Q10 factors for soft POM and bSI remin,
! water TEMP is now passed in instead of Tfunc (1/2005, JKM)
! !USES:
#ifdef CCSMCOUPLED
use shr_sys_mod, only: shr_sys_abort
#endif
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: k ! vertical model level
real (r8), dimension(nx_block,ny_block), intent(in) :: &
TEMP, & ! temperature for scaling functions
O2_loc ! dissolved oxygen used to change/increase POC%diss
type (block), intent(in) :: &
this_block ! block info for the current block
! !INPUT/OUTPUT PARAMETERS:
type(sinking_particle), intent(inout) :: &
POC, & ! base units = nmol C
P_CaCO3, & ! base units = nmol CaCO3
P_SiO2, & ! base units = nmol SiO2
dust, & ! base units = g
P_iron ! base units = nmol Fe
real (r8), dimension(nx_block,ny_block,max_blocks_clinic), intent(inout) :: &
QA_dust_def ! incoming deficit in the QA(dust) POC flux
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
real (r8) :: poc_diss ! diss. length used (cm)
character(*), parameter :: &
subname = 'ecosys_mod:compute_particulate_terms'
real (r8), dimension(nx_block,ny_block) :: &
WORK, & ! temporary for summed quantities to be averaged
TfuncP, & ! temperature scaling from soft POM remin
TfuncS, & ! temperature scaling from soft POM remin
DECAY_CaCO3, & ! scaling factor for dissolution of CaCO3
DECAY_dust, & ! scaling factor for dissolution of dust
DECAY_Hard, & ! scaling factor for dissolution of Hard Ballast
DECAY_HardDust ! scaling factor for dissolution of Hard dust
real (r8) :: &
decay_POC_E, & ! scaling factor for dissolution of excess POC
decay_SiO2, & ! scaling factor for dissolution of SiO2
POC_PROD_avail, & ! POC production available for excess POC flux
new_QA_dust_def, & ! outgoing deficit in the QA(dust) POC flux
dz_loc, dzr_loc ! dz, dzr at a particular i,j location
integer (int_kind) :: &
i, j, & ! loop indices
bid ! local_block id
logical (log_kind) :: &
poc_error ! POC error flag
!-----------------------------------------------------------------------
! incoming fluxes are outgoing fluxes from previous level
!-----------------------------------------------------------------------
bid = this_block%local_id
P_CaCO3%sflux_in(:,:,bid) = P_CaCO3%sflux_out(:,:,bid)
P_CaCO3%hflux_in(:,:,bid) = P_CaCO3%hflux_out(:,:,bid)
P_SiO2%sflux_in(:,:,bid) = P_SiO2%sflux_out(:,:,bid)
P_SiO2%hflux_in(:,:,bid) = P_SiO2%hflux_out(:,:,bid)
dust%sflux_in(:,:,bid) = dust%sflux_out(:,:,bid)
dust%hflux_in(:,:,bid) = dust%hflux_out(:,:,bid)
POC%sflux_in(:,:,bid) = POC%sflux_out(:,:,bid)
POC%hflux_in(:,:,bid) = POC%hflux_out(:,:,bid)
P_iron%sflux_in(:,:,bid) = P_iron%sflux_out(:,:,bid)
P_iron%hflux_in(:,:,bid) = P_iron%hflux_out(:,:,bid)
!-----------------------------------------------------------------------
! compute decay factors
!-----------------------------------------------------------------------
if (partial_bottom_cells) then
DECAY_CaCO3 = exp(-DZT(:,:,k,bid) / P_CaCO3%diss)
DECAY_dust = exp(-DZT(:,:,k,bid) / dust%diss)
DECAY_Hard = exp(-DZT(:,:,k,bid) / 4.0e6_r8)
DECAY_HardDust = exp(-DZT(:,:,k,bid) / 1.2e7_r8)
else
DECAY_CaCO3 = exp(-dz(k) / P_CaCO3%diss)
DECAY_dust = exp(-dz(k) / dust%diss)
DECAY_Hard = exp(-dz(k) / 4.0e6_r8)
DECAY_HardDust = exp(-dz(k) / 1.2e7_r8)
endif
!----------------------------------------------------------------------
! Tref = 30.0 reference temperature (deg. C)
!
! Using q10 formulation with Q10 value of 1.1 soft POM (TfuncP) and
! a Q10 value of 2.5 soft bSi (TfuncS)
!-----------------------------------------------------------------------
!
! NOTE: Turning off temperature effect on POM lengthscale, three instances
! of TfuncP below have been removed, see comment lines.
!-------------------------------------------------------------------------
! TfuncP = 1.1_r8**(((TEMP + T0_Kelvin) - (Tref + T0_Kelvin)) / c10)
TfuncS = 2.5_r8**(((TEMP + T0_Kelvin) - (Tref + T0_Kelvin)) / c10)
poc_error = .false.
do j = 1,ny_block
do i = 1,nx_block
if (LAND_MASK(i,j,bid) .and. k <= KMT(i,j,bid)) then
if (partial_bottom_cells) then
dz_loc = DZT(i,j,k,bid)
else
dz_loc = dz(k)
endif
dzr_loc = c1 / dz_loc
!-----------------------------------------------------------------------
! increase POC diss where O2 concentrations are low
!-----------------------------------------------------------------------
! if (O2_loc(i,j) > c0) then
! poc_diss = min((POC%diss * 2.0_r8), (POC%diss / &
! (O2_loc(i,j) / (O2_loc(i,j) + 5.0_r8))))
! else
! poc_diss = POC%diss * 2.0_r8
! endif
poc_diss = POC%diss
!-----------------------------------------------------------------------
! decay_POC_E and decay_SiO2 set locally, modified by Tfunc
!-----------------------------------------------------------------------
! decay_POC_E = exp(-dz_loc / (poc_diss / TfuncP(i,j)))
decay_POC_E = exp(-dz_loc / poc_diss)
decay_SiO2 = exp(-dz_loc / (P_SiO2%diss / TfuncS(i,j)))
!-----------------------------------------------------------------------
! Set outgoing fluxes for non-iron pools.
! The outoing fluxes for ballast materials are from the
! solution of the coresponding continuous ODE across the model
! level. The ODE has a constant source term and linear decay.
! It is assumed that there is no sub-surface dust production.
!-----------------------------------------------------------------------
P_CaCO3%sflux_out(i,j,bid) = P_CaCO3%sflux_in(i,j,bid) * DECAY_CaCO3(i,j) + &
P_CaCO3%prod(i,j,bid) * ((c1 - P_CaCO3%gamma) * (c1 - DECAY_CaCO3(i,j)) &
* P_CaCO3%diss)
P_CaCO3%hflux_out(i,j,bid) = P_CaCO3%hflux_in(i,j,bid) * DECAY_Hard(i,j) + &
P_CaCO3%prod(i,j,bid) * (P_CaCO3%gamma * dz_loc)
P_SiO2%sflux_out(i,j,bid) = P_SiO2%sflux_in(i,j,bid) * decay_SiO2 + &
P_SiO2%prod(i,j,bid) * ((c1 - P_SiO2%gamma) * (c1 - decay_SiO2) &
* (P_SiO2%diss / TfuncS(i,j)))
P_SiO2%hflux_out(i,j,bid) = P_SiO2%hflux_in(i,j,bid) * DECAY_Hard(i,j) + &
P_SiO2%prod(i,j,bid) * (P_SiO2%gamma * dz_loc)
dust%sflux_out(i,j,bid) = dust%sflux_in(i,j,bid) * DECAY_dust(i,j)
dust%hflux_out(i,j,bid) = dust%hflux_in(i,j,bid) * DECAY_HardDust(i,j)
!-----------------------------------------------------------------------
! Compute how much POC_PROD is available for deficit reduction
! and excess POC flux after subtracting off fraction of non-dust
! ballast production from net POC_PROD.
!-----------------------------------------------------------------------
POC_PROD_avail = POC%prod(i,j,bid) - &
P_CaCO3%rho * P_CaCO3%prod(i,j,bid) - &
P_SiO2%rho * P_SiO2%prod(i,j,bid)
!-----------------------------------------------------------------------
! Check for POC production bounds violations
!-----------------------------------------------------------------------
if (POC_PROD_avail < c0) then
poc_error = .true.
endif
!-----------------------------------------------------------------------
! Compute 1st approximation to new QA_dust_def, the QA_dust
! deficit leaving the cell. Ignore POC_PROD_avail at this stage.
!-----------------------------------------------------------------------
if (QA_dust_def(i,j,bid) > 0) then
new_QA_dust_def = QA_dust_def(i,j,bid) * &
(dust%sflux_out(i,j,bid) + dust%hflux_out(i,j,bid)) / &
(dust%sflux_in(i,j,bid) + dust%hflux_in(i,j,bid))
else
new_QA_dust_def = c0
endif
!-----------------------------------------------------------------------
! Use POC_PROD_avail to reduce new_QA_dust_def.
!-----------------------------------------------------------------------
if (new_QA_dust_def > c0) then
new_QA_dust_def = new_QA_dust_def - POC_PROD_avail * dz_loc
if (new_QA_dust_def < c0) then
POC_PROD_avail = -new_QA_dust_def * dzr_loc
new_QA_dust_def = c0
else
POC_PROD_avail = c0
endif
endif
QA_dust_def(i,j,bid) = new_QA_dust_def
!-----------------------------------------------------------------------
! Compute outgoing POC fluxes. QA POC flux is computing using
! ballast fluxes and new_QA_dust_def. If no QA POC flux came in
! and no production occured, then no QA POC flux goes out. This
! shortcut is present to avoid roundoff cancellation errors from
! the dust%rho * dust_flux_out - QA_dust_def computation.
! Any POC_PROD_avail still remaining goes into excess POC flux.
!-----------------------------------------------------------------------
if (POC%hflux_in(i,j,bid) == c0 .and. POC%prod(i,j,bid) == c0) then
POC%hflux_out(i,j,bid) = c0
else
POC%hflux_out(i,j,bid) = P_CaCO3%rho * &
(P_CaCO3%sflux_out(i,j,bid) + P_CaCO3%hflux_out(i,j,bid)) + &
P_SiO2%rho * &
(P_SiO2%sflux_out(i,j,bid) + P_SiO2%hflux_out(i,j,bid)) + &
dust%rho * &
(dust%sflux_out(i,j,bid) + dust%hflux_out(i,j,bid)) - &
new_QA_dust_def
POC%hflux_out(i,j,bid) = max(POC%hflux_out(i,j,bid), c0)
endif
POC%sflux_out(i,j,bid) = POC%sflux_in(i,j,bid) * decay_POC_E + &
POC_PROD_avail *((c1 - decay_POC_E) * &
poc_diss)
! (poc_diss / TfuncP(i,j)))
!-----------------------------------------------------------------------
! Compute remineralization terms. It is assumed that there is no
! sub-surface dust production.
!-----------------------------------------------------------------------
P_CaCO3%remin(i,j,bid) = P_CaCO3%prod(i,j,bid) + &
((P_CaCO3%sflux_in(i,j,bid) - P_CaCO3%sflux_out(i,j,bid)) + &
(P_CaCO3%hflux_in(i,j,bid) - P_CaCO3%hflux_out(i,j,bid))) * dzr_loc
P_SiO2%remin(i,j,bid) = P_SiO2%prod(i,j,bid) + &
((P_SiO2%sflux_in(i,j,bid) - P_SiO2%sflux_out(i,j,bid)) + &
(P_SiO2%hflux_in(i,j,bid) - P_SiO2%hflux_out(i,j,bid))) * dzr_loc
POC%remin(i,j,bid) = POC%prod(i,j,bid) + &
((POC%sflux_in(i,j,bid) - POC%sflux_out(i,j,bid)) + &
(POC%hflux_in(i,j,bid) - POC%hflux_out(i,j,bid))) * dzr_loc
dust%remin(i,j,bid) = &
((dust%sflux_in(i,j,bid) - dust%sflux_out(i,j,bid)) + &
(dust%hflux_in(i,j,bid) - dust%hflux_out(i,j,bid))) * dzr_loc
!-----------------------------------------------------------------------
! Compute iron remineralization and flux out.
!-----------------------------------------------------------------------
if (POC%sflux_in(i,j,bid) + POC%hflux_in(i,j,bid) == c0) then
P_iron%remin(i,j,bid) = (POC%remin(i,j,bid) * parm_Red_Fe_C)
else
P_iron%remin(i,j,bid) = (POC%remin(i,j,bid) * &
(P_iron%sflux_in(i,j,bid) + P_iron%hflux_in(i,j,bid)) / &
(POC%sflux_in(i,j,bid) + POC%hflux_in(i,j,bid))) + &
(P_iron%sflux_in(i,j,bid) * 3.0e-6_r8)
endif
P_iron%sflux_out(i,j,bid) = P_iron%sflux_in(i,j,bid) + dz_loc * &
((c1 - P_iron%gamma) * P_iron%prod(i,j,bid) - P_iron%remin(i,j,bid))
if (P_iron%sflux_out(i,j,bid) < c0) then
P_iron%sflux_out(i,j,bid) = c0
P_iron%remin(i,j,bid) = P_iron%sflux_in(i,j,bid) * dzr_loc + &
(c1 - P_iron%gamma) * P_iron%prod(i,j,bid)
endif
!-----------------------------------------------------------------------
! Compute iron release from dust remin/dissolution
!
! dust remin gDust = 0.035 / 55.847 * 1.0e9 = 626712.0 nmolFe
! gFe molFe nmolFe
! Also add in Fe source from sediments if applicable to this cell.
!-----------------------------------------------------------------------
P_iron%remin(i,j,bid) = P_iron%remin(i,j,bid) &
+ dust%remin(i,j,bid) * dust_to_Fe &
+ (FESEDFLUX(i,j,k,bid) * dzr_loc)
!maltrud what is up here--jkm does not have this
P_iron%hflux_out(i,j,bid) = P_iron%hflux_in(i,j,bid)
else
P_CaCO3%sflux_out(i,j,bid) = c0
P_CaCO3%hflux_out(i,j,bid) = c0
P_CaCO3%remin(i,j,bid) = c0
P_SiO2%sflux_out(i,j,bid) = c0
P_SiO2%hflux_out(i,j,bid) = c0
P_SiO2%remin(i,j,bid) = c0
dust%sflux_out(i,j,bid) = c0
dust%hflux_out(i,j,bid) = c0
dust%remin(i,j,bid) = c0
POC%sflux_out(i,j,bid) = c0
POC%hflux_out(i,j,bid) = c0
POC%remin(i,j,bid) = c0
P_iron%sflux_out(i,j,bid) = c0
P_iron%hflux_out(i,j,bid) = c0
P_iron%remin(i,j,bid) = c0
endif
!-----------------------------------------------------------------------
! Remineralize everything in bottom cell.
!-----------------------------------------------------------------------
if (LAND_MASK(i,j,bid) .and. k == KMT(i,j,bid)) then
P_CaCO3%remin(i,j,bid) = P_CaCO3%remin(i,j,bid) + &
(P_CaCO3%sflux_out(i,j,bid) + P_CaCO3%hflux_out(i,j,bid)) * dzr_loc
P_CaCO3%sflux_out(i,j,bid) = c0
P_CaCO3%hflux_out(i,j,bid) = c0
P_SiO2%remin(i,j,bid) = P_SiO2%remin(i,j,bid) + &
(P_SiO2%sflux_out(i,j,bid) + P_SiO2%hflux_out(i,j,bid)) * dzr_loc
P_SiO2%sflux_out(i,j,bid) = c0
P_SiO2%hflux_out(i,j,bid) = c0
! dust%remin(i,j,bid) = dust%remin(i,j,bid) + &
! (dust%sflux_out(i,j,bid) + dust%hflux_out(i,j,bid)) * dzr_loc
! dust%sflux_out(i,j,bid) = c0
! dust%hflux_out(i,j,bid) = c0
POC%remin(i,j,bid) = POC%remin(i,j,bid) + &
(POC%sflux_out(i,j,bid) + POC%hflux_out(i,j,bid)) * dzr_loc
POC%sflux_out(i,j,bid) = c0
POC%hflux_out(i,j,bid) = c0
P_iron%remin(i,j,bid) = P_iron%remin(i,j,bid) + &
(P_iron%sflux_out(i,j,bid) + P_iron%hflux_out(i,j,bid)) * dzr_loc
P_iron%sflux_out(i,j,bid) = c0
P_iron%hflux_out(i,j,bid) = c0
endif
end do
end do
#ifdef CCSMCOUPLED
if (poc_error) then
call shr_sys_abort(subname /&
&/ ': mass ratio of ballast ' /&
&/ 'production exceeds POC production')
endif
#endif
!-----------------------------------------------------------------------
! Set tavg variables.
!-----------------------------------------------------------------------
if (tavg_requested(tavg_POC_FLUX_IN)) then
WORK = POC%sflux_in(:,:,bid) + POC%hflux_in(:,:,bid)
call accumulate_tavg_field(WORK, tavg_POC_FLUX_IN,bid,k)
endif
if (tavg_requested(tavg_POC_PROD)) then
call accumulate_tavg_field(POC%prod(:,:,bid), tavg_POC_PROD,bid,k)
endif
if (tavg_requested(tavg_POC_REMIN)) then
call accumulate_tavg_field(POC%remin(:,:,bid), tavg_POC_REMIN,bid,k)
endif
if (tavg_requested(tavg_CaCO3_FLUX_IN)) then
WORK = P_CaCO3%sflux_in(:,:,bid) + P_CaCO3%hflux_in(:,:,bid)
call accumulate_tavg_field(WORK, tavg_CaCO3_FLUX_IN,bid,k)
endif
if (tavg_requested(tavg_CaCO3_PROD)) then
call accumulate_tavg_field(P_CaCO3%prod(:,:,bid), tavg_CaCO3_PROD,bid,k)
endif
if (tavg_requested(tavg_CaCO3_REMIN)) then
call accumulate_tavg_field(P_CaCO3%remin(:,:,bid), tavg_CaCO3_REMIN,bid,k)
endif
if (tavg_requested(tavg_SiO2_FLUX_IN)) then
WORK = P_SiO2%sflux_in(:,:,bid) + P_SiO2%hflux_in(:,:,bid)
call accumulate_tavg_field(WORK, tavg_SiO2_FLUX_IN,bid,k)
endif
if (tavg_requested(tavg_SiO2_PROD)) then
call accumulate_tavg_field(P_SiO2%prod(:,:,bid), tavg_SiO2_PROD,bid,k)
endif
if (tavg_requested(tavg_SiO2_REMIN)) then
call accumulate_tavg_field(P_SiO2%remin(:,:,bid), tavg_SiO2_REMIN,bid,k)
endif
if (tavg_requested(tavg_dust_FLUX_IN)) then
WORK = dust%sflux_in(:,:,bid) + dust%hflux_in(:,:,bid)
call accumulate_tavg_field(WORK, tavg_dust_FLUX_IN,bid,k)
endif
if (tavg_requested(tavg_dust_REMIN)) then
call accumulate_tavg_field(dust%remin(:,:,bid), tavg_dust_REMIN,bid,k)
endif
if (tavg_requested(tavg_P_iron_FLUX_IN)) then
WORK = P_iron%sflux_in(:,:,bid) + P_iron%hflux_in(:,:,bid)
call accumulate_tavg_field(WORK, tavg_P_iron_FLUX_IN,bid,k)
endif
if (tavg_requested(tavg_P_iron_PROD)) then
call accumulate_tavg_field(P_iron%prod(:,:,bid), tavg_P_iron_PROD,bid,k)
endif
if (tavg_requested(tavg_P_iron_REMIN)) then
call accumulate_tavg_field(P_iron%remin(:,:,bid), tavg_P_iron_REMIN,bid,k)
endif
!-----------------------------------------------------------------------
!EOC
end subroutine compute_particulate_terms
!***********************************************************************
!BOP
! !IROUTINE: ecosys_init_sflux
! !INTERFACE:
subroutine ecosys_init_sflux 1,27
! !DESCRIPTION:
! Initialize surface flux computations for ecosys tracer module.
!
! !REVISION HISTORY:
! same as module
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
character(*), parameter :: subname = 'ecosys_mod:ecosys_init_sflux'
logical (log_kind) :: &
luse_INTERP_WORK ! does INTERP_WORK need to be allocated
integer (int_kind) :: &
n, & ! index for looping over tracers
iblock ! index for looping over blocks
real (r8), dimension (nx_block,ny_block) :: WORK
real (r8), dimension (nx_block,ny_block,12,max_blocks_clinic), target :: &
WORK_READ ! temporary space to read in fields
!-----------------------------------------------------------------------
luse_INTERP_WORK = .false.
!-----------------------------------------------------------------------
! read gas flux forcing (if required)
!-----------------------------------------------------------------------
if ((lflux_gas_o2 .or. lflux_gas_co2) .and. &
gas_flux_forcing_iopt == gas_flux_forcing_iopt_file) then
luse_INTERP_WORK = .true.
!-----------------------------------------------------------------------
! first, read ice file
!-----------------------------------------------------------------------
allocate(fice_file%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
if (trim(fice_file%input%filename) == 'unknown') &
fice_file%input%filename = gas_flux_forcing_file
call read_field(fice_file%input%file_fmt, &
fice_file%input%filename, &
fice_file%input%file_varname, &
WORK_READ)
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
fice_file%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
where (.not. LAND_MASK(:,:,iblock)) &
fice_file%DATA(:,:,iblock,1,n) = c0
fice_file%DATA(:,:,iblock,1,n) = &
fice_file%DATA(:,:,iblock,1,n) * fice_file%input%scale_factor
end do
end do
!$OMP END PARALLEL DO
call find_forcing_times(fice_file%data_time, &
fice_file%data_inc, fice_file%interp_type, &
fice_file%data_next, fice_file%data_time_min_loc, &
fice_file%data_update, fice_file%data_type)
!-----------------------------------------------------------------------
! next, read piston velocity file
!-----------------------------------------------------------------------
allocate(xkw_file%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
if (trim(xkw_file%input%filename) == 'unknown') &
xkw_file%input%filename = gas_flux_forcing_file
call read_field(xkw_file%input%file_fmt, &
xkw_file%input%filename, &
xkw_file%input%file_varname, &
WORK_READ)
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
xkw_file%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
where (.not. LAND_MASK(:,:,iblock)) &
xkw_file%DATA(:,:,iblock,1,n) = c0
xkw_file%DATA(:,:,iblock,1,n) = &
xkw_file%DATA(:,:,iblock,1,n) * xkw_file%input%scale_factor
end do
end do
!$OMP END PARALLEL DO
call find_forcing_times(xkw_file%data_time, &
xkw_file%data_inc, xkw_file%interp_type, &
xkw_file%data_next, xkw_file%data_time_min_loc, &
xkw_file%data_update, xkw_file%data_type)
!-----------------------------------------------------------------------
! last, read atmospheric pressure file
!-----------------------------------------------------------------------
allocate(ap_file%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
if (trim(ap_file%input%filename) == 'unknown') &
ap_file%input%filename = gas_flux_forcing_file
call read_field(ap_file%input%file_fmt, &
ap_file%input%filename, &
ap_file%input%file_varname, &
WORK_READ)
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
ap_file%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
where (.not. LAND_MASK(:,:,iblock)) &
ap_file%DATA(:,:,iblock,1,n) = c0
ap_file%DATA(:,:,iblock,1,n) = &
ap_file%DATA(:,:,iblock,1,n) * ap_file%input%scale_factor
end do
end do
!$OMP END PARALLEL DO
call find_forcing_times(ap_file%data_time, &
ap_file%data_inc, ap_file%interp_type, &
ap_file%data_next, ap_file%data_time_min_loc, &
ap_file%data_update, ap_file%data_type)
endif
!-----------------------------------------------------------------------
! load dust flux fields (if required)
!-----------------------------------------------------------------------
if (trim(dust_flux%input%filename) /= 'none' .and. &
trim(dust_flux%input%filename) /= 'unknown') then
luse_INTERP_WORK = .true.
allocate(dust_flux%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
if (trim(dust_flux%input%filename) == 'unknown') &
dust_flux%input%filename = gas_flux_forcing_file
call read_field(dust_flux%input%file_fmt, &
dust_flux%input%filename, &
dust_flux%input%file_varname, &
WORK_READ)
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
dust_flux%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
where (.not. LAND_MASK(:,:,iblock)) &
dust_flux%DATA(:,:,iblock,1,n) = c0
dust_flux%DATA(:,:,iblock,1,n) = &
dust_flux%DATA(:,:,iblock,1,n) * dust_flux%input%scale_factor
end do
end do
!$OMP END PARALLEL DO
call find_forcing_times(dust_flux%data_time, &
dust_flux%data_inc, dust_flux%interp_type, &
dust_flux%data_next, dust_flux%data_time_min_loc, &
dust_flux%data_update, dust_flux%data_type)
dust_flux%has_data = .true.
else
dust_flux%has_data = .false.
endif
!-----------------------------------------------------------------------
! load iron flux fields (if required)
!-----------------------------------------------------------------------
if (trim(iron_flux%input%filename) /= 'none' .and. &
trim(iron_flux%input%filename) /= 'unknown') then
luse_INTERP_WORK = .true.
allocate(iron_flux%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
if (trim(iron_flux%input%filename) == 'unknown') &
iron_flux%input%filename = gas_flux_forcing_file
call read_field(iron_flux%input%file_fmt, &
iron_flux%input%filename, &
iron_flux%input%file_varname, &
WORK_READ)
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
iron_flux%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
where (.not. LAND_MASK(:,:,iblock)) &
iron_flux%DATA(:,:,iblock,1,n) = c0
iron_flux%DATA(:,:,iblock,1,n) = &
iron_flux%DATA(:,:,iblock,1,n) * iron_flux%input%scale_factor
end do
end do
!$OMP END PARALLEL DO
call find_forcing_times(iron_flux%data_time, &
iron_flux%data_inc, iron_flux%interp_type, &
iron_flux%data_next, iron_flux%data_time_min_loc, &
iron_flux%data_update, iron_flux%data_type)
iron_flux%has_data = .true.
else
iron_flux%has_data = .false.
endif
!-----------------------------------------------------------------------
! load nox & noy flux fields (if required)
!-----------------------------------------------------------------------
if (trim(ndep_data_type) /= 'none' .and. &
trim(ndep_data_type) /= 'monthly-calendar' .and. &
trim(ndep_data_type) /= 'shr_stream') then
call document(subname, 'ndep_data_type', ndep_data_type)
call exit_POP(sigAbort, 'unknown ndep_data_type')
endif
if (trim(ndep_data_type) == 'monthly-calendar' .and. &
trim(nox_flux_monthly%input%filename) /= 'none' .and. &
trim(nox_flux_monthly%input%filename) /= 'unknown') then
luse_INTERP_WORK = .true.
allocate(nox_flux_monthly%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
if (trim(nox_flux_monthly%input%filename) == 'unknown') &
nox_flux_monthly%input%filename = gas_flux_forcing_file
call read_field(nox_flux_monthly%input%file_fmt, &
nox_flux_monthly%input%filename, &
nox_flux_monthly%input%file_varname, &
WORK_READ)
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
nox_flux_monthly%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
where (.not. LAND_MASK(:,:,iblock)) &
nox_flux_monthly%DATA(:,:,iblock,1,n) = c0
nox_flux_monthly%DATA(:,:,iblock,1,n) = &
nox_flux_monthly%DATA(:,:,iblock,1,n) * nox_flux_monthly%input%scale_factor
end do
end do
!$OMP END PARALLEL DO
call find_forcing_times(nox_flux_monthly%data_time, &
nox_flux_monthly%data_inc, nox_flux_monthly%interp_type, &
nox_flux_monthly%data_next, nox_flux_monthly%data_time_min_loc, &
nox_flux_monthly%data_update, nox_flux_monthly%data_type)
nox_flux_monthly%has_data = .true.
else
nox_flux_monthly%has_data = .false.
endif
if (trim(ndep_data_type) == 'monthly-calendar' .and. &
trim(nhy_flux_monthly%input%filename) /= 'none' .and. &
trim(nhy_flux_monthly%input%filename) /= 'unknown') then
luse_INTERP_WORK = .true.
allocate(nhy_flux_monthly%DATA(nx_block,ny_block,max_blocks_clinic,1,12))
if (trim(nhy_flux_monthly%input%filename) == 'unknown') &
nhy_flux_monthly%input%filename = gas_flux_forcing_file
call read_field(nhy_flux_monthly%input%file_fmt, &
nhy_flux_monthly%input%filename, &
nhy_flux_monthly%input%file_varname, &
WORK_READ)
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
nhy_flux_monthly%DATA(:,:,iblock,1,n) = WORK_READ(:,:,n,iblock)
where (.not. LAND_MASK(:,:,iblock)) &
nhy_flux_monthly%DATA(:,:,iblock,1,n) = c0
nhy_flux_monthly%DATA(:,:,iblock,1,n) = &
nhy_flux_monthly%DATA(:,:,iblock,1,n) * nhy_flux_monthly%input%scale_factor
end do
end do
!$OMP END PARALLEL DO
call find_forcing_times(nhy_flux_monthly%data_time, &
nhy_flux_monthly%data_inc, nhy_flux_monthly%interp_type, &
nhy_flux_monthly%data_next, nhy_flux_monthly%data_time_min_loc, &
nhy_flux_monthly%data_update, nhy_flux_monthly%data_type)
nhy_flux_monthly%has_data = .true.
else
nhy_flux_monthly%has_data = .false.
endif
!-----------------------------------------------------------------------
#ifndef CCSMCOUPLED
if (trim(ndep_data_type) == 'shr_stream') then
call document(subname, 'ndep_data_type', ndep_data_type)
call exit_POP(sigAbort, &
'shr_stream option only supported when CCSMCOUPLED is defined')
endif
#endif
!-----------------------------------------------------------------------
! allocate space for interpolate_forcing
!-----------------------------------------------------------------------
if (luse_INTERP_WORK) &
allocate(INTERP_WORK(nx_block,ny_block,max_blocks_clinic,1))
!-----------------------------------------------------------------------
! load iron PATCH flux fields (if required)
!-----------------------------------------------------------------------
if (liron_patch) then
!maltrud iron patch
! assume patch file has same normalization and format as deposition file
allocate(IRON_PATCH_FLUX(nx_block,ny_block,max_blocks_clinic))
if (trim(iron_flux%input%filename) == 'unknown') &
iron_flux%input%filename = gas_flux_forcing_file
call read_field(iron_flux%input%file_fmt, &
iron_flux%input%filename, &
iron_patch_flux_filename, &
IRON_PATCH_FLUX)
!$OMP PARALLEL DO PRIVATE(iblock, n)
do iblock=1,nblocks_clinic
do n=1,12
where (.not. LAND_MASK(:,:,iblock)) &
IRON_PATCH_FLUX(:,:,iblock) = c0
iron_flux%DATA(:,:,iblock,1,n) = &
IRON_PATCH_FLUX(:,:,iblock) * iron_flux%input%scale_factor
end do
end do
!$OMP END PARALLEL DO
endif
!-----------------------------------------------------------------------
! register and set
! fco2, the air-sea co2 gas flux
!-----------------------------------------------------------------------
call named_field_register('SFLUX_CO2', sflux_co2_nf_ind)
!$OMP PARALLEL DO PRIVATE(iblock,WORK)
do iblock=1,nblocks_clinic
WORK = c0
call named_field_set(sflux_co2_nf_ind, iblock, WORK)
end do
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
! verify running coupled if gas fluxes use coupler forcing
!-----------------------------------------------------------------------
if ((lflux_gas_o2 .or. lflux_gas_co2) .and. &
(gas_flux_forcing_iopt == gas_flux_forcing_iopt_drv .or. &
atm_co2_iopt == atm_co2_iopt_drv_prog .or. &
atm_co2_iopt == atm_co2_iopt_drv_diag) .and. &
.not. registry_match('lcoupled')) then
call exit_POP(sigAbort, 'ecosys_init: ecosys module requires the ' /&
&/ 'flux coupler when gas_flux_forcing_opt=drv')
endif
!-----------------------------------------------------------------------
! get named field index for atmospheric CO2, if required
!-----------------------------------------------------------------------
if (lflux_gas_co2 .and. atm_co2_iopt == atm_co2_iopt_drv_prog) then
call named_field_get_index('ATM_CO2_PROG', atm_co2_nf_ind, &
exit_on_err=.false.)
if (atm_co2_nf_ind == 0) then
call exit_POP(sigAbort, 'ecosys_init: ecosys module requires ' /&
&/ 'atmopsheric CO2 from the flux coupler ' /&
&/ 'and it is not present')
endif
endif
if (lflux_gas_co2 .and. atm_co2_iopt == atm_co2_iopt_drv_diag) then
call named_field_get_index('ATM_CO2_DIAG', atm_co2_nf_ind, &
exit_on_err=.false.)
if (atm_co2_nf_ind == 0) then
call exit_POP(sigAbort, 'ecosys_init: ecosys module requires ' /&
&/ 'atmopsheric CO2 from the flux coupler ' /&
&/ 'and it is not present')
endif
endif
!-----------------------------------------------------------------------
!EOC
end subroutine ecosys_init_sflux
!***********************************************************************
!BOP
! !IROUTINE: ecosys_init_interior_restore
! !INTERFACE:
subroutine ecosys_init_interior_restore 1,6
! !DESCRIPTION:
! Initialize interior restoring computations for ecosys tracer module.
!
! !REVISION HISTORY:
! same as module
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
integer (int_kind) :: &
n, & ! index for looping over tracers
k, & ! index for looping over levels
i,j, & ! index for looping over horiz. dims.
iblock ! index for looping over blocks
real (r8) :: &
subsurf_fesed ! sum of subsurface fesed values
!-----------------------------------------------------------------------
! initialize restoring timescale (if required)
!-----------------------------------------------------------------------
if (lrest_po4 .or. lrest_no3 .or. lrest_sio3) then
if (lnutr_variable_restore) then
allocate(NUTR_RESTORE_RTAU(nx_block,ny_block,max_blocks_clinic))
allocate(NUTR_RESTORE_MAX_LEVEL(nx_block,ny_block,max_blocks_clinic))
call read_field(nutr_variable_rest_file_fmt, &
nutr_variable_rest_file, &
'NUTR_RESTORE_MAX_LEVEL', &
NUTR_RESTORE_RTAU)
NUTR_RESTORE_MAX_LEVEL = nint(NUTR_RESTORE_RTAU)
call read_field(nutr_variable_rest_file_fmt, &
nutr_variable_rest_file, &
'NUTR_RESTORE_RTAU', &
NUTR_RESTORE_RTAU)
NUTR_RESTORE_RTAU = NUTR_RESTORE_RTAU/seconds_in_day ! convert days to secs
else
do k = 1,km
if (zt(k) < rest_z0) then
nutr_rest_time_inv(k) = rest_time_inv_surf
else if (zt(k) > rest_z1) then
nutr_rest_time_inv(k) = rest_time_inv_deep
else if (rest_z1 == rest_z0) then
nutr_rest_time_inv(k) = rest_time_inv_surf + p5 * &
(rest_time_inv_deep - rest_time_inv_surf)
else
nutr_rest_time_inv(k) = rest_time_inv_surf + &
(zt(k) - rest_z0) / (rest_z1 - rest_z0) * &
(rest_time_inv_deep - rest_time_inv_surf)
endif
end do
endif ! variable restoring
endif ! restoring
!-----------------------------------------------------------------------
! load restoring fields (if required)
!-----------------------------------------------------------------------
if (lrest_po4) then
allocate(PO4_CLIM(nx_block,ny_block,km,max_blocks_clinic))
call read_field(po4_rest%file_fmt, &
po4_rest%filename, &
po4_rest%file_varname, &
PO4_CLIM)
do iblock=1,nblocks_clinic
do k = 1, km
where (.not. LAND_MASK(:,:,iblock) .or. k > KMT(:,:,iblock)) &
PO4_CLIM(:,:,k,iblock) = c0
PO4_CLIM(:,:,k,iblock) = PO4_CLIM(:,:,k,iblock) * po4_rest%scale_factor
enddo
end do
endif
if (lrest_no3) then
allocate(NO3_CLIM(nx_block,ny_block,km,max_blocks_clinic))
call read_field(no3_rest%file_fmt, &
no3_rest%filename, &
no3_rest%file_varname, &
NO3_CLIM)
do iblock=1,nblocks_clinic
do k = 1, km
where (.not. LAND_MASK(:,:,iblock) .or. k > KMT(:,:,iblock)) &
NO3_CLIM(:,:,k,iblock) = c0
NO3_CLIM(:,:,k,iblock) = NO3_CLIM(:,:,k,iblock) * no3_rest%scale_factor
enddo
end do
endif
if (lrest_sio3) then
allocate(SiO3_CLIM(nx_block,ny_block,km,max_blocks_clinic))
call read_field(sio3_rest%file_fmt, &
sio3_rest%filename, &
sio3_rest%file_varname, &
SiO3_CLIM)
do iblock=1,nblocks_clinic
do k = 1, km
where (.not. LAND_MASK(:,:,iblock) .or. k > KMT(:,:,iblock)) &
SiO3_CLIM(:,:,k,iblock) = c0
SiO3_CLIM(:,:,k,iblock) = SiO3_CLIM(:,:,k,iblock) * sio3_rest%scale_factor
enddo
end do
endif
!-----------------------------------------------------------------------
! load fesedflux
! add subsurface positives to 1 level shallower, to accomodate overflow pop-ups
!-----------------------------------------------------------------------
allocate(FESEDFLUX(nx_block,ny_block,km,max_blocks_clinic))
call read_field(fesedflux_input%file_fmt, &
fesedflux_input%filename, &
fesedflux_input%file_varname, &
FESEDFLUX)
do iblock=1,nblocks_clinic
do j=1,ny_block
do i=1,nx_block
if (KMT(i,j,iblock) > 0 .and. KMT(i,j,iblock) < km) then
subsurf_fesed = c0
do k=KMT(i,j,iblock)+1,km
subsurf_fesed = subsurf_fesed + FESEDFLUX(i,j,k,iblock)
enddo
FESEDFLUX(i,j,KMT(i,j,iblock),iblock) = FESEDFLUX(i,j,KMT(i,j,iblock),iblock) + subsurf_fesed
endif
enddo
enddo
do k = 1, km
where (.not. LAND_MASK(:,:,iblock) .or. k > KMT(:,:,iblock)) &
FESEDFLUX(:,:,k,iblock) = c0
FESEDFLUX(:,:,k,iblock) = FESEDFLUX(:,:,k,iblock) * fesedflux_input%scale_factor
enddo
end do
!-----------------------------------------------------------------------
!EOC
end subroutine ecosys_init_interior_restore
!***********************************************************************
!BOP
! !IROUTINE: ecosys_set_sflux
! !INTERFACE:
subroutine ecosys_set_sflux(SHF_QSW_RAW, SHF_QSW, & 1,37
U10_SQR,IFRAC,PRESS,SST,SSS, &
SURF_VALS_OLD,SURF_VALS_CUR,STF_MODULE)
! !DESCRIPTION:
! Compute surface fluxes for ecosys tracer module.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,max_blocks_clinic), intent(in) :: &
SHF_QSW_RAW, &! penetrative solar heat flux, from coupler (degC*cm/s)
SHF_QSW, &! SHF_QSW used by physics, may have diurnal cylce imposed (degC*cm/s)
U10_SQR, &! 10m wind speed squared (cm/s)**2
IFRAC, &! sea ice fraction (non-dimensional)
PRESS, &! sea level atmospheric pressure (dyne/cm**2)
SST, &! sea surface temperature (C)
SSS ! sea surface salinity (psu)
real (r8), dimension(nx_block,ny_block,ecosys_tracer_cnt,max_blocks_clinic), &
intent(in) :: SURF_VALS_OLD, SURF_VALS_CUR ! module tracers
! !INPUT/OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,ecosys_tracer_cnt,max_blocks_clinic), &
intent(inout) :: STF_MODULE
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
character(*), parameter :: subname = 'ecosys_mod:ecosys_set_sflux'
logical (log_kind) :: first_call = .true.
type (block) :: &
this_block ! block info for the current block
integer (int_kind) :: &
i,j,iblock,n, & ! loop indices
mcdate,sec, & ! date vals for shr_strdata_advance
errorCode ! errorCode from HaloUpdate call
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
IFRAC_USED, & ! used ice fraction (non-dimensional)
XKW_USED, & ! portion of piston velocity (cm/s)
AP_USED, & ! used atm pressure (converted from dyne/cm**2 to atm)
IRON_FLUX_IN ! iron flux
real (r8), dimension(nx_block,ny_block,max_blocks_clinic) :: &
SHR_STREAM_WORK
real (r8), dimension(nx_block,ny_block) :: &
XKW_ICE, & ! common portion of piston vel., (1-fice)*xkw (cm/s)
SCHMIDT_USED, & ! used Schmidt number
PV, & ! piston velocity (cm/s)
O2SAT_1atm, & ! O2 saturation @ 1 atm (mmol/m^3)
O2SAT_USED, & ! used O2 saturation (mmol/m^3)
XCO2, & ! atmospheric co2 conc. (dry-air, 1 atm)
FLUX ! tracer flux (nmol/cm^2/s)
real (r8), dimension(nx_block) :: &
PHLO, & ! lower bound for ph in solver
PHHI, & ! upper bound for ph in solver
PH_NEW, & ! computed PH from solver
DIC_ROW, & ! row of DIC values for solver
ALK_ROW, & ! row of ALK values for solver
PO4_ROW, & ! row of PO4 values for solver
SiO3_ROW, & ! row of SiO3 values for solver
CO2STAR_ROW, & ! CO2STAR from solver
DCO2STAR_ROW, & ! DCO2STAR from solver
pCO2SURF_ROW, & ! pCO2SURF from solver
DpCO2_ROW ! DpCO2 from solver
character (char_len) :: &
tracer_data_label, & ! label for what is being updated
ndep_shr_stream_fldList
character (char_len), dimension(1) :: &
tracer_data_names ! short names for input data fields
integer (int_kind), dimension(1) :: &
tracer_bndy_loc, & ! location and field type for ghost
tracer_bndy_type ! cell updates
real (r8), dimension(nx_block,ny_block) :: &
WORK1, WORK2 ! temporaries for averages
real (r8) :: scalar_temp
!-----------------------------------------------------------------------
! local parameters
!-----------------------------------------------------------------------
real (r8), parameter :: &
xkw_coeff = 8.6e-9_r8 ! a = 0.31 cm/hr s^2/m^2 in (s/cm)
!-----------------------------------------------------------------------
call timer_start(ecosys_sflux_timer)
!-----------------------------------------------------------------------
if (check_time_flag(comp_surf_avg_flag)) &
call comp_surf_avg(SURF_VALS_OLD,SURF_VALS_CUR)
!-----------------------------------------------------------------------
! fluxes initially set to 0
! set Chl field for short-wave absorption
! store incoming shortwave in PAR_out field, converting to W/m^2
!-----------------------------------------------------------------------
scalar_temp = f_qsw_par / hflux_factor
!$OMP PARALLEL DO PRIVATE(iblock,WORK1)
do iblock = 1, nblocks_clinic
STF_MODULE(:,:,:,iblock) = c0
WORK1 = max(c0,p5*(SURF_VALS_OLD(:,:,spChl_ind,iblock) + &
SURF_VALS_CUR(:,:,spChl_ind,iblock))) + &
max(c0,p5*(SURF_VALS_OLD(:,:,diatChl_ind,iblock) + &
SURF_VALS_CUR(:,:,diatChl_ind,iblock))) + &
max(c0,p5*(SURF_VALS_OLD(:,:,diazChl_ind,iblock) + &
SURF_VALS_CUR(:,:,diazChl_ind,iblock)))
call named_field_set(totChl_surf_nf_ind, iblock, WORK1)
if (ecosys_qsw_distrb_const) then
PAR_out(:,:,iblock) = SHF_QSW_RAW(:,:,iblock)
else
PAR_out(:,:,iblock) = SHF_QSW(:,:,iblock)
endif
where (LAND_MASK(:,:,iblock))
PAR_out(:,:,iblock) = max(c0, scalar_temp * PAR_out(:,:,iblock))
elsewhere
PAR_out(:,:,iblock) = c0
end where
enddo
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
! Interpolate gas flux forcing data if necessary
!-----------------------------------------------------------------------
if ((lflux_gas_o2 .or. lflux_gas_co2) .and. &
gas_flux_forcing_iopt == gas_flux_forcing_iopt_file) then
if (thour00 >= fice_file%data_update) then
tracer_data_names = fice_file%input%file_varname
tracer_bndy_loc = field_loc_center
tracer_bndy_type = field_type_scalar
tracer_data_label = 'Ice Fraction'
call update_forcing_data(fice_file%data_time, &
fice_file%data_time_min_loc, fice_file%interp_type, &
fice_file%data_next, fice_file%data_update, &
fice_file%data_type, fice_file%data_inc, &
fice_file%DATA(:,:,:,:,1:12), fice_file%data_renorm, &
tracer_data_label, tracer_data_names, &
tracer_bndy_loc, tracer_bndy_type, &
fice_file%filename, fice_file%input%file_fmt)
endif
call interpolate_forcing(INTERP_WORK, &
fice_file%DATA(:,:,:,:,1:12), &
fice_file%data_time, fice_file%interp_type, &
fice_file%data_time_min_loc, fice_file%interp_freq, &
fice_file%interp_inc, fice_file%interp_next, &
fice_file%interp_last, 0)
IFRAC_USED = INTERP_WORK(:,:,:,1)
if (thour00 >= xkw_file%data_update) then
tracer_data_names = xkw_file%input%file_varname
tracer_bndy_loc = field_loc_center
tracer_bndy_type = field_type_scalar
tracer_data_label = 'Piston Velocity'
call update_forcing_data(xkw_file%data_time, &
xkw_file%data_time_min_loc, xkw_file%interp_type, &
xkw_file%data_next, xkw_file%data_update, &
xkw_file%data_type, xkw_file%data_inc, &
xkw_file%DATA(:,:,:,:,1:12), xkw_file%data_renorm, &
tracer_data_label, tracer_data_names, &
tracer_bndy_loc, tracer_bndy_type, &
xkw_file%filename, xkw_file%input%file_fmt)
endif
call interpolate_forcing(INTERP_WORK, &
xkw_file%DATA(:,:,:,:,1:12), &
xkw_file%data_time, xkw_file%interp_type, &
xkw_file%data_time_min_loc, xkw_file%interp_freq, &
xkw_file%interp_inc, xkw_file%interp_next, &
xkw_file%interp_last, 0)
XKW_USED = INTERP_WORK(:,:,:,1)
if (thour00 >= ap_file%data_update) then
tracer_data_names = ap_file%input%file_varname
tracer_bndy_loc = field_loc_center
tracer_bndy_type = field_type_scalar
tracer_data_label = 'Atmospheric Pressure'
call update_forcing_data(ap_file%data_time, &
ap_file%data_time_min_loc, ap_file%interp_type, &
ap_file%data_next, ap_file%data_update, &
ap_file%data_type, ap_file%data_inc, &
ap_file%DATA(:,:,:,:,1:12), ap_file%data_renorm, &
tracer_data_label, tracer_data_names, &
tracer_bndy_loc, tracer_bndy_type, &
ap_file%filename, ap_file%input%file_fmt)
endif
call interpolate_forcing(INTERP_WORK, &
ap_file%DATA(:,:,:,:,1:12), &
ap_file%data_time, ap_file%interp_type, &
ap_file%data_time_min_loc, ap_file%interp_freq, &
ap_file%interp_inc, ap_file%interp_next, &
ap_file%interp_last, 0)
AP_USED = INTERP_WORK(:,:,:,1)
endif
!-----------------------------------------------------------------------
! calculate gas flux quantities if necessary
!-----------------------------------------------------------------------
if (lflux_gas_o2 .or. lflux_gas_co2) then
!$OMP PARALLEL DO PRIVATE(iblock,j,XKW_ICE,SCHMIDT_USED,PV,O2SAT_USED, &
!$OMP O2SAT_1atm,FLUX,XCO2,PHLO,PHHI,DIC_ROW,ALK_ROW, &
!$OMP PO4_ROW,SiO3_ROW,PH_NEW,CO2STAR_ROW, &
!$OMP DCO2STAR_ROW,pCO2SURF_ROW,DpCO2_ROW)
do iblock = 1, nblocks_clinic
!-----------------------------------------------------------------------
! Apply OCMIP ice fraction mask when input is from a file.
!-----------------------------------------------------------------------
if (gas_flux_forcing_iopt == gas_flux_forcing_iopt_file) then
where (IFRAC_USED(:,:,iblock) < 0.2000_r8) &
IFRAC_USED(:,:,iblock) = 0.2000_r8
where (IFRAC_USED(:,:,iblock) > 0.9999_r8) &
IFRAC_USED(:,:,iblock) = 0.9999_r8
endif
if (gas_flux_forcing_iopt == gas_flux_forcing_iopt_drv) then
IFRAC_USED(:,:,iblock) = IFRAC(:,:,iblock)
where (IFRAC_USED(:,:,iblock) < c0) IFRAC_USED(:,:,iblock) = c0
where (IFRAC_USED(:,:,iblock) > c1) IFRAC_USED(:,:,iblock) = c1
XKW_USED(:,:,iblock) = xkw_coeff * U10_SQR(:,:,iblock)
AP_USED(:,:,iblock) = PRESS(:,:,iblock)
endif
!-----------------------------------------------------------------------
! assume PRESS is in cgs units (dyne/cm**2) since that is what is
! required for pressure forcing in barotropic
! want units to be atmospheres
! convertion from dyne/cm**2 to Pascals is P(mks) = P(cgs)/10.
! convertion from Pascals to atm is P(atm) = P(Pa)/101.325e+3_r8
!-----------------------------------------------------------------------
AP_USED(:,:,iblock) = PRESS(:,:,iblock) / 101.325e+4_r8
ECO_SFLUX_TAVG(:,:,1,iblock) = IFRAC_USED(:,:,iblock)
ECO_SFLUX_TAVG(:,:,2,iblock) = XKW_USED(:,:,iblock)
ECO_SFLUX_TAVG(:,:,3,iblock) = AP_USED(:,:,iblock)
!-----------------------------------------------------------------------
! Compute XKW_ICE. XKW is zero over land, so XKW_ICE is too.
!-----------------------------------------------------------------------
XKW_ICE = (c1 - IFRAC_USED(:,:,iblock)) * XKW_USED(:,:,iblock)
!-----------------------------------------------------------------------
! compute O2 flux
!-----------------------------------------------------------------------
if (lflux_gas_o2) then
SCHMIDT_USED = SCHMIDT_O2(SST(:,:,iblock), LAND_MASK(:,:,iblock))
O2SAT_1atm = O2SAT(SST(:,:,iblock),SSS(:,:,iblock), &
LAND_MASK(:,:,iblock))
where (LAND_MASK(:,:,iblock))
PV = XKW_ICE * SQRT(660.0_r8 / SCHMIDT_USED)
O2SAT_USED = AP_USED(:,:,iblock) * O2SAT_1atm
FLUX = PV * (O2SAT_USED - p5*(SURF_VALS_OLD(:,:,o2_ind,iblock) + &
SURF_VALS_CUR(:,:,o2_ind,iblock)))
STF_MODULE(:,:,o2_ind,iblock) = FLUX
elsewhere
O2SAT_USED = c0
STF_MODULE(:,:,o2_ind,iblock) = c0
end where
ECO_SFLUX_TAVG(:,:,4,iblock) = PV
ECO_SFLUX_TAVG(:,:,5,iblock) = SCHMIDT_USED
ECO_SFLUX_TAVG(:,:,6,iblock) = O2SAT_USED
ECO_SFLUX_TAVG(:,:,16,iblock) = STF_MODULE(:,:,o2_ind,iblock)
endif ! lflux_gas_o2
!-----------------------------------------------------------------------
! compute CO2 flux, computing disequilibrium one row at a time
!-----------------------------------------------------------------------
if (lflux_gas_co2) then
SCHMIDT_USED = SCHMIDT_CO2(SST(:,:,iblock), LAND_MASK(:,:,iblock))
where (LAND_MASK(:,:,iblock))
PV = XKW_ICE * SQRT(660.0_r8 / SCHMIDT_USED)
elsewhere
PV = c0
end where
!-----------------------------------------------------------------------
! Set XCO2
!-----------------------------------------------------------------------
select case (atm_co2_iopt)
case (atm_co2_iopt_const)
XCO2 = atm_co2_const
case (atm_co2_iopt_drv_prog, atm_co2_iopt_drv_diag)
call named_field_get(atm_co2_nf_ind, iblock, XCO2)
end select
do j = 1,ny_block
where (PH_PREV(:,j,iblock) /= c0)
PHLO = PH_PREV(:,j,iblock) - del_ph
PHHI = PH_PREV(:,j,iblock) + del_ph
elsewhere
PHLO = phlo_init
PHHI = phhi_init
end where
DIC_ROW = p5*(SURF_VALS_OLD(:,j,dic_ind,iblock) + &
SURF_VALS_CUR(:,j,dic_ind,iblock))
ALK_ROW = p5*(SURF_VALS_OLD(:,j,alk_ind,iblock) + &
SURF_VALS_CUR(:,j,alk_ind,iblock))
PO4_ROW = p5*(SURF_VALS_OLD(:,j,po4_ind,iblock) + &
SURF_VALS_CUR(:,j,po4_ind,iblock))
SiO3_ROW = p5*(SURF_VALS_OLD(:,j,sio3_ind,iblock) + &
SURF_VALS_CUR(:,j,sio3_ind,iblock))
call co2calc_row(iblock, LAND_MASK(:,j,iblock), locmip_k1_k2_bug_fix, &
SST(:,j,iblock), SSS(:,j,iblock), &
DIC_ROW, ALK_ROW, PO4_ROW, SiO3_ROW, &
PHLO, PHHI, PH_NEW, XCO2(:,j), &
AP_USED(:,j,iblock), CO2STAR_ROW, &
DCO2STAR_ROW, pCO2SURF_ROW, DpCO2_ROW)
PH_PREV(:,j,iblock) = PH_NEW
FLUX(:,j) = PV(:,j) * DCO2STAR_ROW
ECO_SFLUX_TAVG(:,j, 7,iblock) = CO2STAR_ROW
ECO_SFLUX_TAVG(:,j, 8,iblock) = DCO2STAR_ROW
ECO_SFLUX_TAVG(:,j, 9,iblock) = pCO2SURF_ROW
ECO_SFLUX_TAVG(:,j,10,iblock) = DpCO2_ROW
end do
!-----------------------------------------------------------------------
! set air-sea co2 gas flux named field, converting units from
! nmol/cm^2/s (positive down) to kg CO2/m^2/s (positive down)
!-----------------------------------------------------------------------
call named_field_set(sflux_co2_nf_ind, iblock, 44.0e-8_r8 * FLUX)
STF_MODULE(:,:,dic_ind,iblock) = STF_MODULE(:,:,dic_ind,iblock) + FLUX
ECO_SFLUX_TAVG(:,:,11,iblock) = PV
ECO_SFLUX_TAVG(:,:,12,iblock) = SCHMIDT_USED
ECO_SFLUX_TAVG(:,:,13,iblock) = FLUX
ECO_SFLUX_TAVG(:,:,14,iblock) = PH_PREV(:,:,iblock)
ECO_SFLUX_TAVG(:,:,15,iblock) = XCO2
endif ! lflux_gas_co2
enddo
!$OMP END PARALLEL DO
endif ! lflux_gas_o2 .or. lflux_gas_co2
!-----------------------------------------------------------------------
! calculate iron and dust fluxes if necessary
!-----------------------------------------------------------------------
if (iron_flux%has_data) then
if (thour00 >= iron_flux%data_update) then
tracer_data_names = iron_flux%input%file_varname
tracer_bndy_loc = field_loc_center
tracer_bndy_type = field_type_scalar
tracer_data_label = 'Iron Flux'
call update_forcing_data(iron_flux%data_time, &
iron_flux%data_time_min_loc, iron_flux%interp_type, &
iron_flux%data_next, iron_flux%data_update, &
iron_flux%data_type, iron_flux%data_inc, &
iron_flux%DATA(:,:,:,:,1:12), iron_flux%data_renorm, &
tracer_data_label, tracer_data_names, &
tracer_bndy_loc, tracer_bndy_type, &
iron_flux%filename, iron_flux%input%file_fmt)
endif
call interpolate_forcing(INTERP_WORK, &
iron_flux%DATA(:,:,:,:,1:12), &
iron_flux%data_time, iron_flux%interp_type, &
iron_flux%data_time_min_loc, iron_flux%interp_freq, &
iron_flux%interp_inc, iron_flux%interp_next, &
iron_flux%interp_last, 0)
if (liron_patch .and. imonth == iron_patch_month) then
IRON_FLUX_IN = INTERP_WORK(:,:,:,1) + IRON_PATCH_FLUX
else
IRON_FLUX_IN = INTERP_WORK(:,:,:,1)
endif
else
IRON_FLUX_IN = c0
endif
IRON_FLUX_IN = IRON_FLUX_IN * parm_Fe_bioavail
STF_MODULE(:,:,fe_ind,:) = IRON_FLUX_IN
if (dust_flux%has_data) then
if (thour00 >= dust_flux%data_update) then
tracer_data_names = dust_flux%input%file_varname
tracer_bndy_loc = field_loc_center
tracer_bndy_type = field_type_scalar
tracer_data_label = 'Dust Flux'
call update_forcing_data(dust_flux%data_time, &
dust_flux%data_time_min_loc, dust_flux%interp_type, &
dust_flux%data_next, dust_flux%data_update, &
dust_flux%data_type, dust_flux%data_inc, &
dust_flux%DATA(:,:,:,:,1:12), dust_flux%data_renorm, &
tracer_data_label, tracer_data_names, &
tracer_bndy_loc, tracer_bndy_type, &
dust_flux%filename, dust_flux%input%file_fmt)
endif
call interpolate_forcing(INTERP_WORK, &
dust_flux%DATA(:,:,:,:,1:12), &
dust_flux%data_time, dust_flux%interp_type, &
dust_flux%data_time_min_loc, dust_flux%interp_freq, &
dust_flux%interp_inc, dust_flux%interp_next, &
dust_flux%interp_last, 0)
dust_FLUX_IN = INTERP_WORK(:,:,:,1)
!-----------------------------------------------------------------------
! Reduce surface dust flux due to assumed instant surface dissolution
!-----------------------------------------------------------------------
dust_FLUX_IN = dust_FLUX_IN * (c1 - parm_Fe_bioavail)
else
dust_FLUX_IN = c0
endif
!-----------------------------------------------------------------------
! calculate nox and nhy fluxes if necessary
!-----------------------------------------------------------------------
if (nox_flux_monthly%has_data) then
if (thour00 >= nox_flux_monthly%data_update) then
tracer_data_names = nox_flux_monthly%input%file_varname
tracer_bndy_loc = field_loc_center
tracer_bndy_type = field_type_scalar
tracer_data_label = 'NOx Flux'
call update_forcing_data(nox_flux_monthly%data_time, &
nox_flux_monthly%data_time_min_loc, nox_flux_monthly%interp_type, &
nox_flux_monthly%data_next, nox_flux_monthly%data_update, &
nox_flux_monthly%data_type, nox_flux_monthly%data_inc, &
nox_flux_monthly%DATA(:,:,:,:,1:12), nox_flux_monthly%data_renorm, &
tracer_data_label, tracer_data_names, &
tracer_bndy_loc, tracer_bndy_type, &
nox_flux_monthly%filename, nox_flux_monthly%input%file_fmt)
endif
call interpolate_forcing(INTERP_WORK, &
nox_flux_monthly%DATA(:,:,:,:,1:12), &
nox_flux_monthly%data_time, nox_flux_monthly%interp_type, &
nox_flux_monthly%data_time_min_loc, nox_flux_monthly%interp_freq, &
nox_flux_monthly%interp_inc, nox_flux_monthly%interp_next, &
nox_flux_monthly%interp_last, 0)
STF_MODULE(:,:,no3_ind,:) = INTERP_WORK(:,:,:,1)
endif
if (nhy_flux_monthly%has_data) then
if (thour00 >= nhy_flux_monthly%data_update) then
tracer_data_names = nhy_flux_monthly%input%file_varname
tracer_bndy_loc = field_loc_center
tracer_bndy_type = field_type_scalar
tracer_data_label = 'NHy Flux'
call update_forcing_data(nhy_flux_monthly%data_time, &
nhy_flux_monthly%data_time_min_loc, nhy_flux_monthly%interp_type, &
nhy_flux_monthly%data_next, nhy_flux_monthly%data_update, &
nhy_flux_monthly%data_type, nhy_flux_monthly%data_inc, &
nhy_flux_monthly%DATA(:,:,:,:,1:12), nhy_flux_monthly%data_renorm, &
tracer_data_label, tracer_data_names, &
tracer_bndy_loc, tracer_bndy_type, &
nhy_flux_monthly%filename, nhy_flux_monthly%input%file_fmt)
endif
call interpolate_forcing(INTERP_WORK, &
nhy_flux_monthly%DATA(:,:,:,:,1:12), &
nhy_flux_monthly%data_time, nhy_flux_monthly%interp_type, &
nhy_flux_monthly%data_time_min_loc, nhy_flux_monthly%interp_freq, &
nhy_flux_monthly%interp_inc, nhy_flux_monthly%interp_next, &
nhy_flux_monthly%interp_last, 0)
STF_MODULE(:,:,nh4_ind,:) = INTERP_WORK(:,:,:,1)
endif
#ifdef CCSMCOUPLED
if (trim(ndep_data_type) == 'shr_stream') then
if (first_call) then
ndep_shr_stream_fldList = ' '
do n = 1, ndep_shr_stream_var_cnt
if (n == ndep_shr_stream_no_ind) &
ndep_shr_stream_fldList = trim(ndep_shr_stream_fldList) /&
&/ 'NOy_deposition'
if (n == ndep_shr_stream_nh_ind) &
ndep_shr_stream_fldList = trim(ndep_shr_stream_fldList) /&
&/ 'NHx_deposition'
if (n < ndep_shr_stream_var_cnt) &
ndep_shr_stream_fldList = trim(ndep_shr_stream_fldList) /&
&/ ':'
end do
call shr_strdata_create(ndep_sdat,name='ndep data', &
mpicom=POP_communicator, &
compid=POP_MCT_OCNID, &
gsmap=POP_MCT_gsMap_o, ggrid=POP_MCT_dom_o, &
nxg=nx_global, nyg=ny_global, &
yearFirst=ndep_shr_stream_year_first, &
yearLast=ndep_shr_stream_year_last, &
yearAlign=ndep_shr_stream_year_align, &
offset=0, &
domFilePath='', &
domFileName=ndep_shr_stream_file, &
domTvarName='time', &
domXvarName='TLONG', domYvarName='TLAT', &
domAreaName='TAREA', domMaskName='KMT', &
FilePath='', &
FileName=(/trim(ndep_shr_stream_file)/), &
fldListFile=ndep_shr_stream_fldList, &
fldListModel=ndep_shr_stream_fldList, &
fillalgo='none', mapalgo='none')
if (my_task == master_task) then
call shr_strdata_print(ndep_sdat)
endif
first_call = .false.
endif
mcdate = iyear*10000 + imonth*100 + iday
sec = isecond + 60 * (iminute + 60 * ihour)
call timer_start(ecosys_shr_strdata_advance_timer)
call shr_strdata_advance(ndep_sdat, mcdate, sec, &
POP_communicator, 'ndep')
call timer_stop(ecosys_shr_strdata_advance_timer)
call timer_print(ecosys_shr_strdata_advance_timer)
!
! process NO3 flux, store results in SHR_STREAM_WORK array
! instead of directly into STF_MODULE
! to avoid argument copies in HaloUpdate calls
!
n = 0
do iblock = 1, nblocks_clinic
this_block = get_block(blocks_clinic(iblock),iblock)
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
n = n + 1
SHR_STREAM_WORK(i,j,iblock) = &
ndep_sdat%avs(1)%rAttr(ndep_shr_stream_no_ind,n)
enddo
enddo
enddo
call POP_HaloUpdate(SHR_STREAM_WORK,POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode, &
fillValue = 0.0_POP_r8)
if (errorCode /= POP_Success) then
call exit_POP(sigAbort, subname /&
&/ ': error updating halo for Ndep fields')
endif
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock = 1, nblocks_clinic
where (LAND_MASK(:,:,iblock))
STF_MODULE(:,:,no3_ind,iblock) = &
ndep_shr_stream_scale_factor * SHR_STREAM_WORK(:,:,iblock)
elsewhere
STF_MODULE(:,:,no3_ind,iblock) = c0
endwhere
enddo
!$OMP END PARALLEL DO
!
! process NH4 flux, store results in SHR_STREAM_WORK array
! instead of directly into STF_MODULE
! to avoid argument copies in HaloUpdate calls
!
n = 0
do iblock = 1, nblocks_clinic
this_block = get_block(blocks_clinic(iblock),iblock)
do j=this_block%jb,this_block%je
do i=this_block%ib,this_block%ie
n = n + 1
SHR_STREAM_WORK(i,j,iblock) = &
ndep_sdat%avs(1)%rAttr(ndep_shr_stream_nh_ind,n)
enddo
enddo
enddo
call POP_HaloUpdate(SHR_STREAM_WORK,POP_haloClinic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode, &
fillValue = 0.0_POP_r8)
if (errorCode /= POP_Success) then
call exit_POP(sigAbort, subname /&
&/ ': error updating halo for Ndep fields')
endif
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock = 1, nblocks_clinic
where (LAND_MASK(:,:,iblock))
STF_MODULE(:,:,nh4_ind,iblock) = &
ndep_shr_stream_scale_factor * SHR_STREAM_WORK(:,:,iblock)
elsewhere
STF_MODULE(:,:,nh4_ind,iblock) = c0
endwhere
enddo
!$OMP END PARALLEL DO
endif
#endif
call timer_stop(ecosys_sflux_timer)
!-----------------------------------------------------------------------
!EOC
end subroutine ecosys_set_sflux
!*****************************************************************************
!BOP
! !IROUTINE: SCHMIDT_O2
! !INTERFACE:
function SCHMIDT_O2(SST, LAND_MASK) 1
! !DESCRIPTION:
! Compute Schmidt number of O2 in seawater as function of SST
! where LAND_MASK is true. Give zero where LAND_MASK is false.
!
! ref : Keeling et al, Global Biogeochem. Cycles, Vol. 12,
! No. 1, pp. 141-163, March 1998
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(in) :: SST
logical (log_kind), dimension(nx_block,ny_block), intent(in) :: &
LAND_MASK
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block) :: SCHMIDT_O2
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
real (r8), parameter :: &
a = 1638.0_r8, &
b = 81.83_r8, &
c = 1.483_r8, &
d = 0.008004_r8
!-----------------------------------------------------------------------
where (LAND_MASK)
SCHMIDT_O2 = a + SST * (-b + SST * (c + SST * (-d)))
elsewhere
SCHMIDT_O2 = c0
end where
!-----------------------------------------------------------------------
!EOC
end function SCHMIDT_O2
!*****************************************************************************
!BOP
! !IROUTINE: O2SAT
! !INTERFACE:
function O2SAT(SST, SSS, LAND_MASK) 2
! !DESCRIPTION:
!
! Computes oxygen saturation concentration at 1 atm total pressure
! in mmol/m^3 given the temperature (t, in deg C) and the salinity (s,
! in permil) where LAND_MASK is true. Give zero where LAND_MASK is false.
!
! FROM GARCIA AND GORDON (1992), LIMNOLOGY and OCEANOGRAPHY.
! THE FORMULA USED IS FROM PAGE 1310, EQUATION (8).
!
! *** NOTE: THE "A_3*TS^2" TERM (IN THE PAPER) IS INCORRECT. ***
! *** IT SHOULD NOT BE THERE. ***
!
! O2SAT IS DEFINED BETWEEN T(freezing) <= T <= 40(deg C) AND
! 0 permil <= S <= 42 permil
! CHECK VALUE: T = 10.0 deg C, S = 35.0 permil,
! O2SAT = 282.015 mmol/m^3
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(in) :: &
SST, & ! sea surface temperature (C)
SSS ! sea surface salinity (psu)
logical (log_kind), dimension(nx_block,ny_block), intent(in) :: &
LAND_MASK
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block) :: O2SAT
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
real (r8), dimension(nx_block,ny_block) :: TS
!-----------------------------------------------------------------------
! coefficients in expansion
!-----------------------------------------------------------------------
real (r8), parameter :: &
a_0 = 2.00907_r8, &
a_1 = 3.22014_r8, &
a_2 = 4.05010_r8, &
a_3 = 4.94457_r8, &
a_4 = -2.56847E-1_r8, &
a_5 = 3.88767_r8, &
b_0 = -6.24523E-3_r8, &
b_1 = -7.37614E-3_r8, &
b_2 = -1.03410E-2_r8, &
b_3 = -8.17083E-3_r8, &
c_0 = -4.88682E-7_r8
where (LAND_MASK)
TS = log( ((T0_Kelvin+25.0_r8) - SST) / (T0_Kelvin + SST) )
O2SAT = exp(a_0+TS*(a_1+TS*(a_2+TS*(a_3+TS*(a_4+TS*a_5)))) + &
SSS*( (b_0+TS*(b_1+TS*(b_2+TS*b_3))) + SSS*c_0 ))
elsewhere
O2SAT = c0
end where
!-----------------------------------------------------------------------
! Convert from ml/l to mmol/m^3
!-----------------------------------------------------------------------
O2SAT = O2SAT / 0.0223916_r8
!-----------------------------------------------------------------------
!EOC
end function O2SAT
!*****************************************************************************
!BOP
! !IROUTINE: SCHMIDT_CO2
! !INTERFACE:
function SCHMIDT_CO2(SST, LAND_MASK) 1
! !DESCRIPTION:
! Compute Schmidt number of CO2 in seawater as function of SST
! where LAND_MASK is true. Give zero where LAND_MASK is false.
!
! ref : Wanninkhof, J. Geophys. Res, Vol. 97, No. C5,
! pp. 7373-7382, May 15, 1992
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block), intent(in) :: SST
logical (log_kind), dimension(nx_block,ny_block), intent(in) :: &
LAND_MASK
! !OUTPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block) :: SCHMIDT_CO2
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
real (r8), parameter :: &
a = 2073.1_r8, &
b = 125.62_r8, &
c = 3.6276_r8, &
d = 0.043219_r8
where (LAND_MASK)
SCHMIDT_CO2 = a + SST * (-b + SST * (c + SST * (-d)))
elsewhere
SCHMIDT_CO2 = c0
end where
!-----------------------------------------------------------------------
!EOC
end function SCHMIDT_CO2
!*****************************************************************************
!BOP
! !IROUTINE: ecosys_tavg_forcing
! !INTERFACE:
subroutine ecosys_tavg_forcing(STF_MODULE) 1,48
! !DESCRIPTION:
! Accumulate non-standard forcing related tavg variables.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,ecosys_tracer_cnt,max_blocks_clinic), &
intent(in) :: STF_MODULE
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
integer (int_kind) :: &
iblock ! block loop index
!-----------------------------------------------------------------------
! multiply IRON, DUST fluxes by mpercm (.01) to convert from model
! units (cm/s)(mmol/m^3) to mmol/s/m^2
!-----------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock)
do iblock = 1,nblocks_clinic
if (tavg_requested(tavg_IRON_FLUX)) then
call accumulate_tavg_field(STF_MODULE(:,:,fe_ind,iblock)*mpercm, &
tavg_IRON_FLUX,iblock,1)
endif
if (tavg_requested(tavg_DUST_FLUX)) then
call accumulate_tavg_field(dust_FLUX_IN(:,:,iblock)*mpercm, &
tavg_DUST_FLUX,iblock,1)
endif
if (tavg_requested(tavg_NOx_FLUX)) then
call accumulate_tavg_field(STF_MODULE(:,:,no3_ind,iblock), &
tavg_NOx_FLUX,iblock,1)
endif
if (tavg_requested(tavg_NHy_FLUX)) then
call accumulate_tavg_field(STF_MODULE(:,:,nh4_ind,iblock), &
tavg_NHy_FLUX,iblock,1)
endif
!-----------------------------------------------------------------------
! now do components of flux calculations saved in hack array
!-----------------------------------------------------------------------
if (tavg_requested(tavg_ECOSYS_IFRAC)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,1,iblock), &
tavg_ECOSYS_IFRAC,iblock,1)
endif
if (tavg_requested(tavg_ECOSYS_IFRAC_2)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,1,iblock), &
tavg_ECOSYS_IFRAC_2,iblock,1)
endif
if (tavg_requested(tavg_ECOSYS_XKW)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,2,iblock), &
tavg_ECOSYS_XKW,iblock,1)
endif
if (tavg_requested(tavg_ECOSYS_XKW_2)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,2,iblock), &
tavg_ECOSYS_XKW_2,iblock,1)
endif
if (tavg_requested(tavg_ECOSYS_ATM_PRESS)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,3,iblock), &
tavg_ECOSYS_ATM_PRESS,iblock,1)
endif
if (tavg_requested(tavg_PV_O2)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,4,iblock), &
tavg_PV_O2,iblock,1)
endif
if (tavg_requested(tavg_SCHMIDT_O2)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,5,iblock), &
tavg_SCHMIDT_O2,iblock,1)
endif
if (tavg_requested(tavg_O2SAT)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,6,iblock), &
tavg_O2SAT,iblock,1)
endif
if (tavg_requested(tavg_CO2STAR)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,7,iblock), &
tavg_CO2STAR,iblock,1)
endif
if (tavg_requested(tavg_DCO2STAR)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,8,iblock), &
tavg_DCO2STAR,iblock,1)
endif
if (tavg_requested(tavg_pCO2SURF)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,9,iblock), &
tavg_pCO2SURF,iblock,1)
endif
if (tavg_requested(tavg_DpCO2)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,10,iblock), &
tavg_DpCO2,iblock,1)
endif
if (tavg_requested(tavg_DpCO2_2)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,10,iblock), &
tavg_DpCO2_2,iblock,1)
endif
if (tavg_requested(tavg_PV_CO2)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,11,iblock), &
tavg_PV_CO2,iblock,1)
endif
if (tavg_requested(tavg_SCHMIDT_CO2)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,12,iblock), &
tavg_SCHMIDT_CO2,iblock,1)
endif
if (tavg_requested(tavg_DIC_GAS_FLUX)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,13,iblock), &
tavg_DIC_GAS_FLUX,iblock,1)
endif
if (tavg_requested(tavg_DIC_GAS_FLUX_2)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,13,iblock), &
tavg_DIC_GAS_FLUX_2,iblock,1)
endif
if (tavg_requested(tavg_PH)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,14,iblock), &
tavg_PH,iblock,1)
endif
if (tavg_requested(tavg_ATM_CO2)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,15,iblock), &
tavg_ATM_CO2,iblock,1)
endif
if (tavg_requested(tavg_O2_GAS_FLUX_2)) then
call accumulate_tavg_field(ECO_SFLUX_TAVG(:,:,16,iblock), &
tavg_O2_GAS_FLUX_2,iblock,1)
endif
end do
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
!EOC
end subroutine ecosys_tavg_forcing
!*****************************************************************************
!BOP
! !IROUTINE: comp_surf_avg
! !INTERFACE:
subroutine comp_surf_avg(SURF_VALS_OLD,SURF_VALS_CUR) 3,1
! !DESCRIPTION:
! compute average surface tracer values
!
! avg = sum(SURF_VAL*TAREA) / sum(TAREA)
! with the sum taken over ocean points only
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (r8), dimension(nx_block,ny_block,ecosys_tracer_cnt,max_blocks_clinic), &
intent(in) :: SURF_VALS_OLD, SURF_VALS_CUR
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
integer (int_kind) :: &
n, & ! tracer index
iblock, & ! block index
dcount, & ! diag counter
ib,ie,jb,je
real (r8), dimension(max_blocks_clinic,ecosys_tracer_cnt) :: &
local_sums ! array for holding block sums of each diagnostic
real (r8) :: &
sum_tmp ! temp for local sum
real (r8), dimension(nx_block,ny_block) :: &
WORK1, &! local work space
TFACT ! factor for normalizing sums
type (block) :: &
this_block ! block information for current block
!-----------------------------------------------------------------------
local_sums = c0
!jw !$OMP PARALLEL DO PRIVATE(iblock,this_block,ib,ie,jb,je,TFACT,n,WORK1)
do iblock = 1,nblocks_clinic
this_block = get_block(blocks_clinic(iblock),iblock)
ib = this_block%ib
ie = this_block%ie
jb = this_block%jb
je = this_block%je
TFACT = TAREA(:,:,iblock)*RCALCT(:,:,iblock)
do n = 1, ecosys_tracer_cnt
if (vflux_flag(n)) then
WORK1 = p5*(SURF_VALS_OLD(:,:,n,iblock) + &
SURF_VALS_CUR(:,:,n,iblock))*TFACT
local_sums(iblock,n) = sum(WORK1(ib:ie,jb:je))
endif
end do
end do
!jw !$OMP END PARALLEL DO
do n = 1, ecosys_tracer_cnt
if (vflux_flag(n)) then
sum_tmp = sum(local_sums(:,n))
surf_avg(n) = global_sum(sum_tmp,distrb_clinic)/area_t
endif
end do
if(my_task == master_task) then
write(stdout,*)' Calculating surface tracer averages'
do n = 1, ecosys_tracer_cnt
if (vflux_flag(n)) then
write(stdout,*) n, surf_avg(n)
endif
end do
endif
!-----------------------------------------------------------------------
!EOC
end subroutine comp_surf_avg
!*****************************************************************************
!BOP
! !IROUTINE: ecosys_write_restart
! !INTERFACE:
subroutine ecosys_write_restart(restart_file, action) 1,6
! !DESCRIPTION:
! write auxiliary fields & scalars to restart files
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
character(*), intent(in) :: action
! !INPUT/OUTPUT PARAMETERS:
type (datafile), intent (inout) :: restart_file
!EOP
!BOC
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
character (char_len) :: &
short_name ! tracer name temporaries
type (io_dim) :: &
i_dim, j_dim ! dimension descriptors
integer (int_kind) :: n
type (io_field_desc), save :: PH_SURF
!-----------------------------------------------------------------------
if (trim(action) == 'add_attrib_file') then
short_name = char_blank
do n=1,ecosys_tracer_cnt
if (vflux_flag(n)) then
short_name = 'surf_avg_' /&
&/ ind_name_table(n)%name
call add_attrib_file(restart_file,trim(short_name),surf_avg(n))
endif
end do
endif
if (trim(action) == 'define') then
i_dim = construct_io_dim('i', nx_global)
j_dim = construct_io_dim('j', ny_global)
PH_SURF = construct_io_field('PH_SURF', i_dim, j_dim, &
long_name='surface pH at current time', &
units ='pH', grid_loc ='2110', &
field_loc = field_loc_center, &
field_type = field_type_scalar, &
d2d_array = PH_PREV)
call data_set (restart_file, 'define', PH_SURF)
endif
if (trim(action) == 'write') then
call data_set (restart_file, 'write', PH_SURF)
endif
!-----------------------------------------------------------------------
!EOC
end subroutine ecosys_write_restart
!*****************************************************************************
!BOP
! !IROUTINE: ecosys_tracer_ref_val
! !INTERFACE:
function ecosys_tracer_ref_val(ind) 1
! !DESCRIPTION:
! return reference value for tracers using virtual fluxes
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
integer (int_kind), intent(in) :: ind
! !OUTPUT PARAMETERS:
real (r8) :: ecosys_tracer_ref_val
!EOP
!BOC
!-----------------------------------------------------------------------
if (vflux_flag(ind)) then
ecosys_tracer_ref_val = surf_avg(ind)
else
ecosys_tracer_ref_val = c0
endif
!-----------------------------------------------------------------------
!EOC
end function ecosys_tracer_ref_val
!***********************************************************************
end module ecosys_mod
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||