New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

asmphyto2dbal_medusa.F90 in branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_medusa.F90 @ 9435

Last change on this file since 9435 was 9435, checked in by dford, 6 years ago
Allow assimilation of PFT chlorophyll.
File size: 28.9 KB

Rev	Line
[9432]	1	MODULE asmphyto2dbal_medusa
[8428]	2	!!======================================================================
[9432]	3	!! * MODULE asmphyto2dbal_medusa *
	4	!! Calculate increments to MEDUSA based on surface phyto2d increments
[8428]	5	!!
	6	!! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
	7	!! For licensing reasons this is kept in its own internal Met Office
	8	!! branch (dev/frdf/vn3.6_nitrogen_balancing) rather than in the Paris
	9	!! repository, and must be merged in when building.
	10	!!
	11	!!======================================================================
[9432]	12	!! History : 3.6 ! 2017-08 (D. Ford) Adapted from asmphyto2dbal_hadocc
[8428]	13	!!----------------------------------------------------------------------
[8436]	14	#if defined key_asminc && defined key_medusa && defined key_foam_medusa
[8428]	15	!!----------------------------------------------------------------------
	16	!! 'key_asminc' : assimilation increment interface
[8436]	17	!! 'key_medusa' : MEDUSA model
	18	!! 'key_foam_medusa' : MEDUSA extras for FOAM OBS and ASM
[8428]	19	!!----------------------------------------------------------------------
[9432]	20	!! asm_phyto2d_bal_medusa : routine to calculate increments to MEDUSA
[8428]	21	!!----------------------------------------------------------------------
	22	USE par_kind, ONLY: wp ! kind parameters
	23	USE par_oce, ONLY: jpi, jpj, jpk ! domain array sizes
	24	USE dom_oce, ONLY: gdepw_n ! domain information
[8485]	25	USE zdftmx, ONLY: ln_tmx_itf, & ! Indonesian Throughflow
	26	& mask_itf ! tidal mixing mask
[8428]	27	USE iom ! i/o
[8436]	28	USE sms_medusa ! MEDUSA parameters
	29	USE par_medusa ! MEDUSA parameters
[8428]	30	USE par_trc, ONLY: jptra ! Tracer parameters
	31	USE bioanalysis ! Nitrogen balancing
	32
	33	IMPLICIT NONE
	34	PRIVATE
	35
[9432]	36	PUBLIC asm_phyto2d_bal_medusa
[8428]	37
	38	! Default values for biological assimilation parameters
	39	! Should match Hemmings et al. (2008)
	40	REAL(wp), PARAMETER :: balnutext = 0.6 !: Default nutrient balancing factor
	41	REAL(wp), PARAMETER :: balnutmin = 0.1 !: Fraction of phytoplankton loss to nutrient
	42	REAL(wp), PARAMETER :: r = 1 !: Reliability of model specific growth rate
	43	REAL(wp), PARAMETER :: beta_g = 0.05 !: Low rate bias correction for growth rate estimator
	44	REAL(wp), PARAMETER :: beta_l = 0.05 !: Low rate bias correction for primary loss rate estimator
	45	REAL(wp), PARAMETER :: beta_m = 0.05 !: Low rate bias correction for secondary loss rate estimator
	46	REAL(wp), PARAMETER :: a_g = 0.2 !: Error s.d. for log10 of growth rate estimator
	47	REAL(wp), PARAMETER :: a_l = 0.4 !: Error s.d. for log10 of primary loss rate estimator
	48	REAL(wp), PARAMETER :: a_m = 0.7 !: Error s.d. for log10 of secondary loss rate estimator
	49	REAL(wp), PARAMETER :: zfracb0 = 0.7 !: Base zooplankton fraction of loss to Z & D
	50	REAL(wp), PARAMETER :: zfracb1 = 0 !: Phytoplankton sensitivity of zooplankton fraction
	51	REAL(wp), PARAMETER :: qrfmax = 1.1 !: Maximum nutrient limitation reduction factor
	52	REAL(wp), PARAMETER :: qafmax = 1.1 !: Maximum nutrient limitation amplification factor
	53	REAL(wp), PARAMETER :: zrfmax = 2 !: Maximum zooplankton reduction factor
	54	REAL(wp), PARAMETER :: zafmax = 2 !: Maximum zooplankton amplification factor
	55	REAL(wp), PARAMETER :: prfmax = 10 !: Maximum phytoplankton reduction factor (secondary)
	56	REAL(wp), PARAMETER :: incphymin = 0.0001 !: Minimum size of non-zero phytoplankton increment
	57	REAL(wp), PARAMETER :: integnstep = 20 !: Number of steps for p.d.f. integral evaluation
	58	REAL(wp), PARAMETER :: pthreshold = 0.01 !: Fractional threshold level for setting p.d.f.
	59	!
	60	LOGICAL, PARAMETER :: diag_active = .TRUE. !: Depth-independent diagnostics
	61	LOGICAL, PARAMETER :: diag_fulldepth_active = .TRUE. !: Full-depth diagnostics
	62	LOGICAL, PARAMETER :: gl_active = .TRUE. !: Growth/loss-based balancing
	63	LOGICAL, PARAMETER :: nbal_active = .TRUE. !: Nitrogen balancing
	64	LOGICAL, PARAMETER :: subsurf_active = .TRUE. !: Increments below MLD
	65	LOGICAL, PARAMETER :: deepneg_active = .FALSE. !: Negative primary increments below MLD
	66	LOGICAL, PARAMETER :: deeppos_active = .FALSE. !: Positive primary increments below MLD
	67	LOGICAL, PARAMETER :: nutprof_active = .TRUE. !: Secondary increments
	68
	69	CONTAINS
	70
[9432]	71	SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot, &
	72	& pinc_chltot, &
	73	& ld_chldia, &
	74	& pinc_chldia, &
	75	& ld_chlnon, &
	76	& pinc_chlnon, &
	77	& ld_phytot, &
	78	& pinc_phytot, &
	79	& ld_phydia, &
	80	& pinc_phydia, &
	81	& ld_phynon, &
	82	& pinc_phynon, &
	83	& pincper, &
	84	& p_maxchlinc, ld_phytobal, pmld, &
	85	& pgrow_avg_bkg, ploss_avg_bkg, &
	86	& phyt_avg_bkg, mld_max_bkg, &
	87	& tracer_bkg, phyto2d_balinc )
[8428]	88	!!---------------------------------------------------------------------------
[9432]	89	!! * ROUTINE asm_phyto2d_bal_medusa *
[8428]	90	!!
[9431]	91	!! ** Purpose : calculate increments to MEDUSA from 2d phytoplankton increments
[8428]	92	!!
[9431]	93	!! ** Method : average up MEDUSA to look like HadOCC
[8428]	94	!! call nitrogen balancing scheme
[8436]	95	!! separate back out to MEDUSA
[8428]	96	!!
[9431]	97	!! ** Action : populate phyto2d_balinc
[8428]	98	!!
	99	!! References : Hemmings et al., 2008, J. Mar. Res.
	100	!! Ford et al., 2012, Ocean Sci.
	101	!!---------------------------------------------------------------------------
	102	!!
[9431]	103	LOGICAL, INTENT(in ) :: ld_chltot ! Assim chltot y/n
	104	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: pinc_chltot ! chltot increments
	105	LOGICAL, INTENT(in ) :: ld_chldia ! Assim chldia y/n
	106	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: pinc_chldia ! chldia increments
	107	LOGICAL, INTENT(in ) :: ld_chlnon ! Assim chlnon y/n
	108	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: pinc_chlnon ! chlnon increments
	109	LOGICAL, INTENT(in ) :: ld_phytot ! Assim phytot y/n
	110	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: pinc_phytot ! phytot increments
	111	LOGICAL, INTENT(in ) :: ld_phydia ! Assim phydia y/n
	112	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: pinc_phydia ! phydia increments
	113	LOGICAL, INTENT(in ) :: ld_phynon ! Assim phynon y/n
	114	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: pinc_phynon ! phynon increments
	115	REAL(wp), INTENT(in ) :: pincper ! Assimilation period
	116	REAL(wp), INTENT(in ) :: p_maxchlinc ! Max chl increment
	117	LOGICAL, INTENT(in ) :: ld_phytobal ! Balancing y/n
	118	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: pmld ! Mixed layer depth
[8436]	119	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pgrow_avg_bkg ! Avg phyto growth
	120	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ploss_avg_bkg ! Avg phyto loss
	121	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: phyt_avg_bkg ! Avg phyto
	122	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: mld_max_bkg ! Max MLD
[8440]	123	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg ! State variables
[9431]	124	REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments
[8428]	125	!!
[8436]	126	INTEGER :: ji, jj, jk, jn ! Loop counters
[8428]	127	INTEGER :: jkmax ! Loop index
	128	INTEGER, DIMENSION(6) :: i_tracer ! Tracer indices
[8436]	129	REAL(wp) :: n2be_p ! N:biomass for total phy
	130	REAL(wp) :: n2be_z ! N:biomass for total zoo
	131	REAL(wp) :: n2be_d ! N:biomass for detritus
[9435]	132	REAL(wp) :: zfrac ! Fraction
[8436]	133	REAL(wp) :: zfrac_chn ! Fraction of jpchn
	134	REAL(wp) :: zfrac_chd ! Fraction of jpchd
	135	REAL(wp) :: zfrac_phn ! Fraction of jpphn
	136	REAL(wp) :: zfrac_phd ! Fraction of jpphd
	137	REAL(wp) :: zfrac_zmi ! Fraction of jpzmi
	138	REAL(wp) :: zfrac_zme ! Fraction of jpzme
	139	REAL(wp) :: zrat_pds_phd ! Ratio of jppds:jpphd
	140	REAL(wp) :: zrat_chd_phd ! Ratio of jpchd:jpphd
	141	REAL(wp) :: zrat_chn_phn ! Ratio of jpchn:jpphn
[9435]	142	REAL(wp) :: zrat_phn_chn ! Ratio of jpphn:jpchn
	143	REAL(wp) :: zrat_phd_chd ! Ratio of jpphd:jpchd
	144	REAL(wp) :: zrat_pds_chd ! Ratio of jppds:jpchd
[8436]	145	REAL(wp) :: zrat_dtc_det ! Ratio of jpdtc:jpdet
	146	REAL(wp), DIMENSION(jpi,jpj) :: cchl_p ! C:Chl for total phy
[8428]	147	REAL(wp), DIMENSION(16) :: modparm ! Model parameters
	148	REAL(wp), DIMENSION(20) :: assimparm ! Assimilation parameters
	149	REAL(wp), DIMENSION(jpi,jpj,jpk,6) :: bstate ! Background state
	150	REAL(wp), DIMENSION(jpi,jpj,jpk,6) :: outincs ! Balancing increments
	151	REAL(wp), DIMENSION(jpi,jpj,22) :: diag ! Depth-indep diagnostics
	152	REAL(wp), DIMENSION(jpi,jpj,jpk,22) :: diag_fulldepth ! Full-depth diagnostics
	153	!!---------------------------------------------------------------------------
	154
[9435]	155	! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
	156	IF ( p_maxchlinc > 0.0 ) THEN
	157	IF ( ld_chltot ) THEN
	158	DO jj = 1, jpj
	159	DO ji = 1, jpi
[9431]	160	pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
[9435]	161	END DO
[9431]	162	END DO
[9435]	163	ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN
	164	DO jj = 1, jpj
	165	DO ji = 1, jpi
	166	pinc_chltot(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnon(ji,jj)
	167	pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
	168	IF ( pinc_chltot(ji,jj) .NE. ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) ) ) THEN
	169	zfrac = pinc_chltot(ji,jj) / ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) )
	170	pinc_chldia(ji,jj) = pinc_chldia(ji,jj) * zfrac
	171	pinc_chlnon(ji,jj) = pinc_chlnon(ji,jj) * zfrac
	172	ENDIF
	173	END DO
	174	END DO
	175	ELSE IF ( ld_chldia ) THEN
	176	DO jj = 1, jpj
	177	DO ji = 1, jpi
	178	pinc_chldia(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chldia(ji,jj), p_maxchlinc ) )
	179	pinc_chltot(ji,jj) = pinc_chldia(ji,jj)
	180	END DO
	181	END DO
	182	ELSE IF ( ld_chlnon ) THEN
	183	DO jj = 1, jpj
	184	DO ji = 1, jpi
	185	pinc_chlnon(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chlnon(ji,jj), p_maxchlinc ) )
	186	pinc_chltot(ji,jj) = pinc_chlnon(ji,jj)
	187	END DO
	188	END DO
	189	ENDIF
[8648]	190	ENDIF
[9435]	191
	192	IF ( ld_phytot .OR. ld_phydia .OR. ld_phynon ) THEN
	193	CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' )
	194	ENDIF
[8428]	195
[9431]	196	IF ( ld_phytobal ) THEN ! Nitrogen balancing
[8428]	197
[8436]	198	! Set up model parameters to be passed into Hemmings balancing routine.
	199	! For now these are hardwired to the standard HadOCC parameter values
	200	! (except C:N ratios) as this is what the scheme was developed for.
	201	! Obviously, HadOCC and MEDUSA are rather different models, so this
	202	! isn't ideal, but there's not always direct analogues between the two
	203	! parameter sets, so it's the easiest way to get something running.
	204	! In the longer term, some serious MarMOT-based development is required.
	205	modparm(1) = 0.1 ! grow_sat
	206	modparm(2) = 2.0 ! psmax
	207	modparm(3) = 0.845 ! par
	208	modparm(4) = 0.02 ! alpha
	209	modparm(5) = 0.05 ! resp_rate
	210	modparm(6) = 0.05 ! pmort_rate
	211	modparm(7) = 0.01 ! phyto_min
	212	modparm(8) = 0.05 ! z_mort_1
	213	modparm(9) = 1.0 ! z_mort_2
	214	modparm(10) = ( xthetapn + xthetapd ) / 2.0 ! c2n_p
	215	modparm(11) = ( xthetazmi + xthetazme ) / 2.0 ! c2n_z
	216	modparm(12) = xthetad ! c2n_d
	217	modparm(13) = 0.01 ! graze_threshold
	218	modparm(14) = 2.0 ! holling_coef
	219	modparm(15) = 0.5 ! graze_sat
	220	modparm(16) = 2.0 ! graze_max
[8428]	221
	222	! Set up assimilation parameters to be passed into balancing routine
	223	! Not sure what assimparm(1) is meant to be, but it doesn't get used
	224	assimparm(2) = balnutext
	225	assimparm(3) = balnutmin
	226	assimparm(4) = r
	227	assimparm(5) = beta_g
	228	assimparm(6) = beta_l
	229	assimparm(7) = beta_m
	230	assimparm(8) = a_g
	231	assimparm(9) = a_l
	232	assimparm(10) = a_m
	233	assimparm(11) = zfracb0
	234	assimparm(12) = zfracb1
	235	assimparm(13) = qrfmax
	236	assimparm(14) = qafmax
	237	assimparm(15) = zrfmax
	238	assimparm(16) = zafmax
	239	assimparm(17) = prfmax
	240	assimparm(18) = incphymin
	241	assimparm(19) = integnstep
	242	assimparm(20) = pthreshold
	243
	244	! Set up external tracer indices array bstate
	245	i_tracer(1) = 1 ! nutrient
	246	i_tracer(2) = 2 ! phytoplankton
	247	i_tracer(3) = 3 ! zooplankton
	248	i_tracer(4) = 4 ! detritus
	249	i_tracer(5) = 5 ! DIC
	250	i_tracer(6) = 6 ! Alkalinity
	251
	252	! Set background state
[8440]	253	bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jpdin)
	254	bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jpphn) + tracer_bkg(:,:,:,jpphd)
	255	bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jpzmi) + tracer_bkg(:,:,:,jpzme)
	256	bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jpdet)
	257	bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jpdic)
	258	bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jpalk)
[8428]	259
[8436]	260	! Calculate carbon to chlorophyll ratio for combined phytoplankton
	261	! and nitrogen to biomass equivalent for PZD
	262	! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA
	263	cchl_p(:,:) = 0.0
	264	DO jj = 1, jpj
	265	DO ji = 1, jpi
[8440]	266	IF ( ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) ) .GT. 0.0 ) THEN
[8497]	267	cchl_p(ji,jj) = xmassc * ( ( tracer_bkg(ji,jj,1,jpphn) * xthetapn ) + &
	268	& ( tracer_bkg(ji,jj,1,jpphd) * xthetapd ) ) / &
[8440]	269	& ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) )
[8436]	270	ENDIF
	271	END DO
	272	END DO
	273	n2be_p = 14.01 + ( xmassc * ( ( xthetapn + xthetapd ) / 2.0 ) )
	274	n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) )
	275	n2be_d = 14.01 + ( xmassc * xthetad )
	276
[8428]	277	! Call nitrogen balancing routine
[8436]	278	CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm, &
	279	& n2be_p, n2be_z, n2be_d, assimparm, &
[9431]	280	& INT(pincper), 1, INT(SUM(tmask,3)), tmask(:,:,:), &
	281	& pmld(:,:), mld_max_bkg(:,:), pinc_chltot(:,:), cchl_p(:,:), &
[8436]	282	& nbal_active, phyt_avg_bkg(:,:), &
	283	& gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:), &
	284	& subsurf_active, deepneg_active, &
	285	& deeppos_active, nutprof_active, &
	286	& bstate, outincs, &
	287	& diag_active, diag, &
[8428]	288	& diag_fulldepth_active, diag_fulldepth )
[8436]	289
	290	! Loop over each grid point partioning the increments
[9431]	291	phyto2d_balinc(:,:,:,:) = 0.0
[8436]	292	DO jk = 1, jpk
	293	DO jj = 1, jpj
	294	DO ji = 1, jpi
	295
[9435]	296	! Phytoplankton
	297	IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. &
	298	& ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) .AND. &
	299	& ( pinc_chltot(ji,jj) /= 0.0 ) ) THEN
	300	IF ( ld_chltot ) THEN
	301	! Phytoplankton nitrogen split up based on existing ratios
	302	zfrac_phn = tracer_bkg(ji,jj,jk,jpphn) / &
	303	& (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd))
	304	zfrac_phd = tracer_bkg(ji,jj,jk,jpphd) / &
	305	& (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd))
	306	ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN
	307	! Phytoplankton nitrogen split up based on assimilation increments
	308	zfrac_phn = pinc_chlnon(ji,jj) / pinc_chltot(ji,jj)
	309	zfrac_phd = pinc_chldia(ji,jj) / pinc_chltot(ji,jj)
	310	ENDIF
	311
	312	! Phytoplankton silicate split up based on existing ratios
[8440]	313	zrat_pds_phd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpphd)
[9435]	314
[8436]	315	! Chlorophyll split up based on existing ratios to phytoplankton nitrogen
[9431]	316	! Not using pinc_chltot directly as it's only 2D
	317	! This method should give same results at surface as splitting pinc_chltot would
[8440]	318	zrat_chn_phn = tracer_bkg(ji,jj,jk,jpchn) / tracer_bkg(ji,jj,jk,jpphn)
	319	zrat_chd_phd = tracer_bkg(ji,jj,jk,jpchd) / tracer_bkg(ji,jj,jk,jpphd)
[9435]	320
	321	phyto2d_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn
	322	phyto2d_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd
	323	phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_pds_phd
[9431]	324	phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,jk,jpphn) * zrat_chn_phn
	325	phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_chd_phd
[8436]	326	ENDIF
	327
[9435]	328	! Zooplankton nitrogen split up based on existing ratios
[8440]	329	IF ( ( tracer_bkg(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpzme) > 0.0 ) ) THEN
[9435]	330	zfrac_zmi = tracer_bkg(ji,jj,jk,jpzmi) / &
	331	& (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme))
	332	zfrac_zme = tracer_bkg(ji,jj,jk,jpzme) / &
	333	& (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme))
[9431]	334	phyto2d_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi
	335	phyto2d_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme
[8436]	336	ENDIF
	337
	338	! Nitrogen nutrient straight from balancing scheme
[9431]	339	phyto2d_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1))
[8436]	340
	341	! Nitrogen detritus straight from balancing scheme
[9431]	342	phyto2d_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4))
[8436]	343
	344	! DIC straight from balancing scheme
[9431]	345	phyto2d_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5))
[8436]	346
	347	! Alkalinity straight from balancing scheme
[9431]	348	phyto2d_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6))
[8436]	349
	350	! Remove diatom silicate increment from nutrient silicate to conserve mass
[9431]	351	IF ( ( tracer_bkg(ji,jj,jk,jpsil) - phyto2d_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN
	352	phyto2d_balinc(ji,jj,jk,jpsil) = phyto2d_balinc(ji,jj,jk,jppds) * (-1.0)
[8436]	353	ENDIF
	354
[9435]	355	! Carbon detritus based on existing ratios
[8440]	356	IF ( ( tracer_bkg(ji,jj,jk,jpdet) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpdtc) > 0.0 ) ) THEN
	357	zrat_dtc_det = tracer_bkg(ji,jj,jk,jpdtc) / tracer_bkg(ji,jj,jk,jpdet)
[9431]	358	phyto2d_balinc(ji,jj,jk,jpdtc) = phyto2d_balinc(ji,jj,jk,jpdet) * zrat_dtc_det
[8436]	359	ENDIF
	360
	361	! Do nothing with iron or oxygen for the time being
[9431]	362	phyto2d_balinc(ji,jj,jk,jpfer) = 0.0
	363	phyto2d_balinc(ji,jj,jk,jpoxy) = 0.0
[8436]	364
	365	END DO
	366	END DO
	367	END DO
[8428]	368
	369	ELSE ! No nitrogen balancing
	370
[8436]	371	! Initialise individual chlorophyll increments to zero
[9431]	372	phyto2d_balinc(:,:,:,jpchn) = 0.0
	373	phyto2d_balinc(:,:,:,jpchd) = 0.0
[8428]	374
[8436]	375	! Split up total surface chlorophyll increments
	376	DO jj = 1, jpj
	377	DO ji = 1, jpi
[9435]	378	IF ( ( tracer_bkg(ji,jj,1,jpchn) > 0.0 ) .AND. &
	379	& ( tracer_bkg(ji,jj,1,jpchd) > 0.0 ) ) THEN
	380	IF ( ld_chltot ) THEN
	381	! Chlorophyll split up based on existing ratios
	382	zfrac_chn = tracer_bkg(ji,jj,1,jpchn) / &
	383	& ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) )
	384	zfrac_chd = tracer_bkg(ji,jj,1,jpchd) / &
	385	& ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) )
	386	phyto2d_balinc(ji,jj,1,jpchn) = pinc_chltot(ji,jj) * zfrac_chn
	387	phyto2d_balinc(ji,jj,1,jpchd) = pinc_chltot(ji,jj) * zfrac_chd
	388	ENDIF
	389	IF( ld_chldia ) THEN
	390	phyto2d_balinc(ji,jj,1,jpchd) = pinc_chldia(ji,jj)
	391	ENDIF
	392	IF( ld_chlnon ) THEN
	393	phyto2d_balinc(ji,jj,1,jpchn) = pinc_chlnon(ji,jj)
	394	ENDIF
	395
	396	! Maintain stoichiometric ratios of nitrogen and silicate
	397	IF ( ld_chltot .OR. ld_chlnon ) THEN
	398	zrat_phn_chn = tracer_bkg(ji,jj,1,jpphn) / tracer_bkg(ji,jj,1,jpchn)
	399	phyto2d_balinc(ji,jj,1,jpphn) = phyto2d_balinc(ji,jj,1,jpchn) * zrat_phn_chn
	400	ENDIF
	401	IF ( ld_chltot .OR. ld_chldia ) THEN
	402	zrat_phd_chd = tracer_bkg(ji,jj,1,jpphd) / tracer_bkg(ji,jj,1,jpchd)
	403	phyto2d_balinc(ji,jj,1,jpphd) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_phd_chd
	404	zrat_pds_chd = tracer_bkg(ji,jj,1,jppds) / tracer_bkg(ji,jj,1,jpchd)
	405	phyto2d_balinc(ji,jj,1,jppds) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_pds_chd
	406	ENDIF
[8436]	407	ENDIF
	408	END DO
	409	END DO
[8428]	410
	411	! Propagate through mixed layer
	412	DO jj = 1, jpj
	413	DO ji = 1, jpi
	414	!
	415	jkmax = jpk-1
	416	DO jk = jpk-1, 1, -1
[9431]	417	IF ( ( pmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. &
	418	& ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
	419	pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
[8428]	420	jkmax = jk
	421	ENDIF
	422	END DO
	423	!
	424	DO jk = 2, jkmax
[9431]	425	phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,1,jpchn)
	426	phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,1,jpchd)
[9435]	427	phyto2d_balinc(ji,jj,jk,jpphn) = phyto2d_balinc(ji,jj,1,jpphn)
	428	phyto2d_balinc(ji,jj,jk,jpphd) = phyto2d_balinc(ji,jj,1,jpphd)
	429	phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,1,jppds)
[8428]	430	END DO
	431	!
	432	END DO
	433	END DO
	434
	435	! Set other balancing increments to zero
[9431]	436	phyto2d_balinc(:,:,:,jpzmi) = 0.0
	437	phyto2d_balinc(:,:,:,jpzme) = 0.0
	438	phyto2d_balinc(:,:,:,jpdin) = 0.0
	439	phyto2d_balinc(:,:,:,jpsil) = 0.0
	440	phyto2d_balinc(:,:,:,jpfer) = 0.0
	441	phyto2d_balinc(:,:,:,jpdet) = 0.0
	442	phyto2d_balinc(:,:,:,jpdtc) = 0.0
	443	phyto2d_balinc(:,:,:,jpdic) = 0.0
	444	phyto2d_balinc(:,:,:,jpalk) = 0.0
	445	phyto2d_balinc(:,:,:,jpoxy) = 0.0
[8436]	446
[8428]	447	ENDIF
[8485]	448
	449	! If performing extra tidal mixing in the Indonesian Throughflow,
	450	! increments have been found to make the carbon cycle unstable
	451	! Therefore, mask these out
	452	IF ( ln_tmx_itf ) THEN
	453	DO jn = 1, jptra
	454	DO jk = 1, jpk
[9431]	455	phyto2d_balinc(:,:,jk,jn) = phyto2d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
[8485]	456	END DO
	457	END DO
	458	ENDIF
[8428]	459
[9432]	460	END SUBROUTINE asm_phyto2d_bal_medusa
[8428]	461
	462	#else
	463	!!----------------------------------------------------------------------
	464	!! Default option : Empty routine
	465	!!----------------------------------------------------------------------
	466	CONTAINS
[9432]	467	SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot, &
[9431]	468	& pinc_chltot, &
	469	& ld_chldia, &
	470	& pinc_chldia, &
	471	& ld_chlnon, &
	472	& pinc_chlnon, &
	473	& ld_phytot, &
	474	& pinc_phytot, &
	475	& ld_phydia, &
	476	& pinc_phydia, &
	477	& ld_phynon, &
	478	& pinc_phynon, &
	479	& pincper, &
	480	& p_maxchlinc, ld_phytobal, pmld, &
	481	& pgrow_avg_bkg, ploss_avg_bkg, &
	482	& phyt_avg_bkg, mld_max_bkg, &
	483	& tracer_bkg, phyto2d_balinc )
	484	LOGICAL :: ld_chltot
	485	REAL :: pinc_chltot(:,:)
	486	LOGICAL :: ld_chldia
	487	REAL :: pinc_chldia(:,:)
	488	LOGICAL :: ld_chlnon
	489	REAL :: pinc_chlnon(:,:)
	490	LOGICAL :: ld_phytot
	491	REAL :: pinc_phytot(:,:)
	492	LOGICAL :: ld_phydia
	493	REAL :: pinc_phydia(:,:)
	494	LOGICAL :: ld_phynon
	495	REAL :: pinc_phynon(:,:)
	496	REAL :: pincper
	497	REAL :: p_maxchlinc
	498	LOGICAL :: ld_phytobal
	499	REAL :: pmld(:,:)
[8440]	500	REAL :: pgrow_avg_bkg(:,:)
	501	REAL :: ploss_avg_bkg(:,:)
	502	REAL :: phyt_avg_bkg(:,:)
	503	REAL :: mld_max_bkg(:,:)
	504	REAL :: tracer_bkg(:,:,:,:)
[9431]	505	REAL :: phyto2d_balinc(:,:,:,:)
[9432]	506	WRITE(,) 'asm_phyto2d_bal_medusa: You should not have seen this print! error?'
	507	END SUBROUTINE asm_phyto2d_bal_medusa
[8428]	508	#endif
	509
	510	!!======================================================================
[9432]	511	END MODULE asmphyto2dbal_medusa

Note: See TracBrowser for help on using the repository browser.

Download in other formats: