Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

asmlogchlbal_ersem.F90 @ 9331

Last change on this file since 9331 was 9331, checked in by dford, 6 years ago
Add balancing code.
File size: 36.9 KB

Rev	Line
[6983]	1	MODULE asmlogchlbal_ersem
	2	!!======================================================================
	3	!! * MODULE asmlogchlbal_ersem *
	4	!! Calculate increments to ERSEM based on surface logchl increments
	5	!!======================================================================
	6	!! History : 3.6 ! 2016-09 (D. Ford) Original code
	7	!!----------------------------------------------------------------------
	8	#if defined key_asminc && defined key_fabm
	9	!!----------------------------------------------------------------------
	10	!! 'key_asminc' : assimilation increment interface
	11	!! 'key_fabm' : FABM-ERSEM coupling
	12	!!----------------------------------------------------------------------
	13	!! asm_logchl_bal_ersem : routine to calculate increments to ERSEM
	14	!!----------------------------------------------------------------------
	15	USE par_kind, ONLY: wp ! kind parameters
	16	USE par_oce, ONLY: jpi, jpj, jpk ! domain array sizes
	17	USE dom_oce, ONLY: gdepw_n ! domain information
	18	USE zdfmxl ! mixed layer depth
	19	USE iom ! i/o
	20	USE trc, ONLY: trn ! ERSEM state variables
	21	USE par_fabm ! FABM parameters
[7716]	22	USE par_trc, ONLY: jptra ! Tracer parameters
[9331]	23	USE bioanalysis
[6983]	24
	25	IMPLICIT NONE
	26	PRIVATE
	27
	28	PUBLIC asm_logchl_bal_ersem
	29
[9331]	30	! Default values for biological assimilation parameters
	31	! Should match Hemmings et al. (2008)
	32	REAL(wp), PARAMETER :: balnutext = 0.6 !: Default nutrient balancing factor
	33	REAL(wp), PARAMETER :: balnutmin = 0.1 !: Fraction of phytoplankton loss to nutrient
	34	REAL(wp), PARAMETER :: r = 1 !: Reliability of model specific growth rate
	35	REAL(wp), PARAMETER :: beta_g = 0.05 !: Low rate bias correction for growth rate estimator
	36	REAL(wp), PARAMETER :: beta_l = 0.05 !: Low rate bias correction for primary loss rate estimator
	37	REAL(wp), PARAMETER :: beta_m = 0.05 !: Low rate bias correction for secondary loss rate estimator
	38	REAL(wp), PARAMETER :: a_g = 0.2 !: Error s.d. for log10 of growth rate estimator
	39	REAL(wp), PARAMETER :: a_l = 0.4 !: Error s.d. for log10 of primary loss rate estimator
	40	REAL(wp), PARAMETER :: a_m = 0.7 !: Error s.d. for log10 of secondary loss rate estimator
	41	REAL(wp), PARAMETER :: zfracb0 = 0.7 !: Base zooplankton fraction of loss to Z & D
	42	REAL(wp), PARAMETER :: zfracb1 = 0 !: Phytoplankton sensitivity of zooplankton fraction
	43	REAL(wp), PARAMETER :: qrfmax = 1.1 !: Maximum nutrient limitation reduction factor
	44	REAL(wp), PARAMETER :: qafmax = 1.1 !: Maximum nutrient limitation amplification factor
	45	REAL(wp), PARAMETER :: zrfmax = 2 !: Maximum zooplankton reduction factor
	46	REAL(wp), PARAMETER :: zafmax = 2 !: Maximum zooplankton amplification factor
	47	REAL(wp), PARAMETER :: prfmax = 10 !: Maximum phytoplankton reduction factor (secondary)
	48	REAL(wp), PARAMETER :: incphymin = 0.0001 !: Minimum size of non-zero phytoplankton increment
	49	REAL(wp), PARAMETER :: integnstep = 20 !: Number of steps for p.d.f. integral evaluation
	50	REAL(wp), PARAMETER :: pthreshold = 0.01 !: Fractional threshold level for setting p.d.f.
	51	!
	52	LOGICAL, PARAMETER :: diag_active = .TRUE. !: Depth-independent diagnostics
	53	LOGICAL, PARAMETER :: diag_fulldepth_active = .TRUE. !: Full-depth diagnostics
	54	LOGICAL, PARAMETER :: gl_active = .TRUE. !: Growth/loss-based balancing
	55	LOGICAL, PARAMETER :: nbal_active = .TRUE. !: Nitrogen balancing
	56	LOGICAL, PARAMETER :: subsurf_active = .TRUE. !: Increments below MLD
	57	LOGICAL, PARAMETER :: deepneg_active = .FALSE. !: Negative primary increments below MLD
	58	LOGICAL, PARAMETER :: deeppos_active = .FALSE. !: Positive primary increments below MLD
	59	LOGICAL, PARAMETER :: nutprof_active = .TRUE. !: Secondary increments
	60
[6983]	61	CONTAINS
	62
[7803]	63	SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, &
[9331]	64	& k_maxchlinc, logchl_bkginc, logchl_balinc, &
	65	& pgrow_avg_bkg, ploss_avg_bkg, &
	66	& phyt_avg_bkg, mld_max_bkg, &
	67	& ld_logchlbal, aincper )
[6983]	68	!!---------------------------------------------------------------------------
	69	!! * ROUTINE asm_logchl_bal_ersem *
	70	!!
	71	!! ** Purpose : calculate increments to ERSEM from logchl increments
	72	!!
[7719]	73	!! ** Method : convert logchl increments to chl increments
	74	!! split between the ERSEM PFTs
	75	!! spread through the mixed layer
	76	!! [forthcoming: calculate increments to nutrients and zooplankton]
[6983]	77	!!
[7719]	78	!! ** Action : populate logchl_balinc
[6983]	79	!!
	80	!! References : forthcoming...
	81	!!---------------------------------------------------------------------------
	82	!!
[7803]	83	LOGICAL, INTENT(in ) :: ld_logchlpftinc
[7716]	84	INTEGER, INTENT(in ) :: npfts
	85	INTEGER, INTENT(in ) :: mld_choice_bgc
[7803]	86	REAL(wp), INTENT(in ) :: k_maxchlinc
[7716]	87	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,npfts) :: logchl_bkginc
	88	REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc
[9331]	89	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pgrow_avg_bkg ! Avg phyto growth
	90	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ploss_avg_bkg ! Avg phyto loss
	91	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: phyt_avg_bkg ! Avg phyto
	92	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: mld_max_bkg ! Max MLD
	93	LOGICAL, INTENT(in ) :: ld_logchlbal ! Balancing y/n
	94	REAL(wp), INTENT(in ) :: aincper
[6983]	95	!!
	96	INTEGER :: ji, jj, jk
	97	INTEGER :: jkmax
	98	REAL(wp) :: chl_tot, chl_inc
	99	REAL(wp), DIMENSION(jpi,jpj) :: zmld
[9331]	100	REAL(wp), DIMENSION(16) :: modparm ! Model parameters
	101	REAL(wp), DIMENSION(20) :: assimparm ! Assimilation parameters
	102	REAL(wp), DIMENSION(jpi,jpj,jpk,6) :: bstate ! Background state
	103	REAL(wp), DIMENSION(jpi,jpj,jpk,6) :: outincs ! Balancing increments
	104	REAL(wp), DIMENSION(jpi,jpj,22) :: diag ! Depth-indep diagnostics
	105	REAL(wp), DIMENSION(jpi,jpj,jpk,22) :: diag_fulldepth ! Full-depth diagnostics
	106	INTEGER, DIMENSION(6) :: i_tracer ! Tracer indices
	107	REAL(wp), DIMENSION(jpi,jpj) :: cchl_p ! C:Chl for total phy
	108	REAL(wp), DIMENSION(jpi,jpj) :: chlinctot
	109	REAL(wp) :: n2be_p ! N:biomass for total phy
	110	REAL(wp) :: n2be_z ! N:biomass for total zoo
	111	REAL(wp) :: n2be_d ! N:biomass for detritus
	112	REAL(wp) :: zfrac_n3n ! Fraction of jp_fabm_n3n
	113	REAL(wp) :: zfrac_n4n ! Fraction of jp_fabm_n4n
	114	REAL(wp) :: zfrac_r4n ! Fraction of jp_fabm_r4n
	115	REAL(wp) :: zfrac_r6n ! Fraction of jp_fabm_r6n
	116	REAL(wp) :: zfrac_r8n ! Fraction of jp_fabm_r8n
	117	REAL(wp) :: zfrac_z4n ! Fraction of jp_fabm_z4n
	118	REAL(wp) :: zfrac_z5n ! Fraction of jp_fabm_z5n
	119	REAL(wp) :: zfrac_z6n ! Fraction of jp_fabm_z6n
	120	REAL(wp) :: zfrac_p1n ! Fraction of jp_fabm_p1n
	121	REAL(wp) :: zfrac_p2n ! Fraction of jp_fabm_p2n
	122	REAL(wp) :: zfrac_p3n ! Fraction of jp_fabm_p3n
	123	REAL(wp) :: zfrac_p4n ! Fraction of jp_fabm_p4n
	124	REAL(wp) :: zrat_z4c_z4n
	125	REAL(wp) :: zrat_z5c_z5n
	126	REAL(wp) :: zrat_z5p_z5n
	127	REAL(wp) :: zrat_z6c_z6n
	128	REAL(wp) :: zrat_z6p_z6n
	129	REAL(wp) :: zrat_p1c_p1n
	130	REAL(wp) :: zrat_p1p_p1n
	131	REAL(wp) :: zrat_p1s_p1n
	132	REAL(wp) :: zrat_p2c_p2n
	133	REAL(wp) :: zrat_p2p_p2n
	134	REAL(wp) :: zrat_p3c_p3n
	135	REAL(wp) :: zrat_p3p_p3n
	136	REAL(wp) :: zrat_p4c_p4n
	137	REAL(wp) :: zrat_p4p_p4n
[6983]	138	!!---------------------------------------------------------------------------
	139
	140	! Split surface logchl incs into surface Chl1-4 incs
	141	!
[7716]	142	! In order to transform logchl incs to chl incs, need to account for the background,
	143	! cannot simply do 10^logchl_bkginc. Need to:
	144	! 1) Add logchl inc to log10(background) to get log10(analysis)
	145	! 2) Take 10^log10(analysis) to get analysis
	146	! 3) Subtract background from analysis to get chl incs
[7803]	147	! If k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
[7716]	148	!
[7719]	149	! Only apply increments if all of Chl1-4 background values are > 0
	150	! In theory, they always will be, and if any are not that's a sign
	151	! that something's going wrong which the assimilation might make worse
	152	!
[7803]	153	IF ( ld_logchlpftinc ) THEN
[6983]	154	!
	155	! Assimilating separate PFTs, so separately transform each from LogChl to Chl
	156	!
	157	IF ( npfts /= 4 ) THEN
	158	CALL ctl_stop( 'If assimilating PFTs into ERSEM, nn_asmpfts must be 4' )
	159	ENDIF
	160	DO jj = 1, jpj
	161	DO ji = 1, jpi
[7716]	162	IF ( ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
	163	& ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
	164	& ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
	165	& ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) ) THEN
[7719]	166	IF ( logchl_bkginc(ji,jj,1) /= 0.0 ) THEN
	167	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) ) + &
	168	& logchl_bkginc(ji,jj,1) ) - &
	169	& trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
	170	ENDIF
	171	IF ( logchl_bkginc(ji,jj,2) /= 0.0 ) THEN
	172	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) ) + &
	173	& logchl_bkginc(ji,jj,2) ) - &
	174	& trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
	175	ENDIF
	176	IF ( logchl_bkginc(ji,jj,3) /= 0.0 ) THEN
	177	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) ) + &
	178	& logchl_bkginc(ji,jj,3) ) - &
	179	& trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
	180	ENDIF
	181	IF ( logchl_bkginc(ji,jj,4) /= 0.0 ) THEN
	182	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) ) + &
	183	& logchl_bkginc(ji,jj,4) ) - &
	184	& trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
	185	ENDIF
[7803]	186	IF (k_maxchlinc > 0.0) THEN
	187	chl_inc = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
	188	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
	189	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
	190	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
	191	IF ( ABS(chl_inc) > k_maxchlinc ) THEN
	192	chl_tot = ABS(chl_inc) / k_maxchlinc
	193	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / chl_tot
	194	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / chl_tot
	195	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / chl_tot
	196	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / chl_tot
	197	ENDIF
	198	ENDIF
[6983]	199	ENDIF
	200	END DO
	201	END DO
	202	ELSE
	203	!
	204	! Assimilating total Chl, so transform total from LogChl to Chl
	205	! and split between PFTs according to the existing background ratios
	206	!
[7716]	207	IF ( npfts /= 1 ) THEN
	208	CALL ctl_stop( 'If assimilating total chlorophyll, nn_asmpfts must be 1' )
	209	ENDIF
[6983]	210	DO jj = 1, jpj
	211	DO ji = 1, jpi
[7719]	212	IF ( ( logchl_bkginc(ji,jj,1) /= 0.0 ) .AND. &
	213	& ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
	214	& ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
	215	& ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
	216	& ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) ) THEN
[6983]	217	chl_tot = trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
	218	& trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
[7716]	219	chl_inc = 10**( LOG10( chl_tot ) + logchl_bkginc(ji,jj,1) ) - chl_tot
[7803]	220	IF (k_maxchlinc > 0.0) THEN
	221	chl_inc = MAX( -1.0 * k_maxchlinc, MIN( chl_inc, k_maxchlinc ) )
	222	ENDIF
[7716]	223	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / chl_tot
	224	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / chl_tot
	225	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / chl_tot
	226	logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / chl_tot
[6983]	227	ENDIF
	228	END DO
	229	END DO
	230	ENDIF
	231
	232	! Propagate surface Chl1-4 incs through mixed layer
	233	! First, choose mixed layer definition
	234	!
	235	SELECT CASE( mld_choice_bgc )
	236	CASE ( 1 ) ! Turbocline/mixing depth [W points]
	237	zmld(:,:) = hmld(:,:)
	238	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
	239	zmld(:,:) = hmlp(:,:)
	240	CASE ( 3 ) ! Kara MLD [Interpolated]
	241	#if defined key_karaml
	242	IF ( ln_kara ) THEN
	243	zmld(:,:) = hmld_kara(:,:)
	244	ELSE
	245	CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
	246	& ' but ln_kara=.false.' )
	247	ENDIF
	248	#else
	249	CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
	250	& ' but is not defined' )
	251	#endif
	252	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
	253	zmld(:,:) = hmld_tref(:,:)
	254	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
	255	zmld(:,:) = hmlpt(:,:)
	256	END SELECT
[9331]	257
	258	IF ( ld_logchlbal ) THEN
	259
	260	! Set up model parameters to be passed into Hemmings balancing routine.
	261	! For now these are hardwired to the standard HadOCC parameter values
	262	! as this is what the scheme was developed for.
	263	! Obviously, HadOCC and ERSEM are rather different models, so this
	264	! isn't ideal, but there's not always direct analogues between the two
	265	! parameter sets, so it's the easiest way to get something running.
	266	! In the longer term, some serious MarMOT-based development is required.
	267	modparm(1) = 0.1 ! grow_sat
	268	modparm(2) = 2.0 ! psmax
	269	modparm(3) = 0.845 ! par
	270	modparm(4) = 0.02 ! alpha
	271	modparm(5) = 0.05 ! resp_rate
	272	modparm(6) = 0.05 ! pmort_rate
	273	modparm(7) = 0.01 ! phyto_min
	274	modparm(8) = 0.05 ! z_mort_1
	275	modparm(9) = 1.0 ! z_mort_2
	276	modparm(10) = 6.625 ! c2n_p
	277	modparm(11) = 5.625 ! c2n_z
	278	modparm(12) = 7.5 ! c2n_d
	279	modparm(13) = 0.01 ! graze_threshold
	280	modparm(14) = 2.0 ! holling_coef
	281	modparm(15) = 0.5 ! graze_sat
	282	modparm(16) = 2.0 ! graze_max
	283
	284	! Set up assimilation parameters to be passed into balancing routine
	285	! Not sure what assimparm(1) is meant to be, but it doesn't get used
	286	assimparm(2) = balnutext
	287	assimparm(3) = balnutmin
	288	assimparm(4) = r
	289	assimparm(5) = beta_g
	290	assimparm(6) = beta_l
	291	assimparm(7) = beta_m
	292	assimparm(8) = a_g
	293	assimparm(9) = a_l
	294	assimparm(10) = a_m
	295	assimparm(11) = zfracb0
	296	assimparm(12) = zfracb1
	297	assimparm(13) = qrfmax
	298	assimparm(14) = qafmax
	299	assimparm(15) = zrfmax
	300	assimparm(16) = zafmax
	301	assimparm(17) = prfmax
	302	assimparm(18) = incphymin
	303	assimparm(19) = integnstep
	304	assimparm(20) = pthreshold
	305
	306	! Set up external tracer indices array bstate
	307	i_tracer(1) = 1 ! nutrient
	308	i_tracer(2) = 2 ! phytoplankton
	309	i_tracer(3) = 3 ! zooplankton
	310	i_tracer(4) = 4 ! detritus
	311	i_tracer(5) = 5 ! DIC
	312	i_tracer(6) = 6 ! Alkalinity
	313
	314	! Set background state
	315	bstate(:,:,:,i_tracer(1)) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) + trn(:,:,:,jp_fabm_m1+jp_fabm_n4n)
	316	bstate(:,:,:,i_tracer(2)) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1n) + trn(:,:,:,jp_fabm_m1+jp_fabm_p2n) + &
	317	& trn(:,:,:,jp_fabm_m1+jp_fabm_p3n) + trn(:,:,:,jp_fabm_m1+jp_fabm_p4n)
	318	! Z4c needs converting by qnc, hardwire for now
	319	bstate(:,:,:,i_tracer(3)) = (trn(:,:,:,jp_fabm_m1+jp_fabm_z4c) * 0.0126 ) + &
	320	& trn(:,:,:,jp_fabm_m1+jp_fabm_z5n) + trn(:,:,:,jp_fabm_m1+jp_fabm_z6n)
	321	bstate(:,:,:,i_tracer(4)) = trn(:,:,:,jp_fabm_m1+jp_fabm_r4n) + trn(:,:,:,jp_fabm_m1+jp_fabm_r6n) + &
	322	& trn(:,:,:,jp_fabm_m1+jp_fabm_r8n)
	323	bstate(:,:,:,i_tracer(5)) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c)
	324	bstate(:,:,:,i_tracer(6)) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3a)
	325
	326	! Calculate carbon to chlorophyll ratio for combined phytoplankton
	327	! and nitrogen to biomass equivalent for PZD
	328	! Need a single number, so base on HadOCC
	329	cchl_p(:,:) = 0.0
	330	DO jj = 1, jpj
	331	DO ji = 1, jpi
	332	IF ( ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
	333	& trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) ) .GT. 0.0 ) THEN
	334	cchl_p(ji,jj) = ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_p1c) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) + &
	335	& trn(ji,jj,1,jp_fabm_m1+jp_fabm_p3c) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_p4c) ) / &
	336	& ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
	337	& trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) )
	338	ENDIF
	339	END DO
	340	END DO
	341	n2be_p = ( 14.01 + ( 12.01 * 6.625 ) ) ! / ( 14.01 + ( 12.01 * 6.625 ) )
	342	n2be_z = ( 14.01 + ( 12.01 * 5.625 ) ) ! / ( 14.01 + ( 12.01 * 6.625 ) )
	343	n2be_d = ( 14.01 + ( 12.01 * 7.5 ) ) ! / ( 14.01 + ( 12.01 * 6.625 ) )
	344
	345	chlinctot(:,:) = logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl1) + logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
	346	& logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl3) + logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl4)
	347
	348	! Call nitrogen balancing routine
	349	CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm, &
	350	& n2be_p, n2be_z, n2be_d, assimparm, &
	351	& INT(aincper), 1, INT(SUM(tmask,3)), tmask(:,:,:), &
	352	& zmld(:,:), mld_max_bkg(:,:), chlinctot(:,:), cchl_p(:,:), &
	353	& nbal_active, phyt_avg_bkg(:,:), &
	354	& gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:), &
	355	& subsurf_active, deepneg_active, &
	356	& deeppos_active, nutprof_active, &
	357	& bstate, outincs, &
	358	& diag_active, diag, &
	359	& diag_fulldepth_active, diag_fulldepth )
	360
	361	! Loop over each grid point partioning the increments
	362	DO jk = 1, jpk
	363	DO jj = 1, jpj
	364	DO ji = 1, jpi
	365
	366	! Nitrogen phytoplankton from balancing scheme
	367	! Split according to current ratios [ChlTot] or assimilation [PFTs]
	368	! Update carbon, phosphorus and silicon according to current ratios
	369	! Already have chlorophyll
	370	IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) > 0.0 ) .AND. &
	371	& ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) > 0.0 ) .AND. &
	372	& ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) > 0.0 ) .AND. &
	373	& ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) > 0.0 ) ) THEN
	374
	375	IF ( ld_logchlpftinc ) THEN
	376	IF ( (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
	377	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
	378	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
	379	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)) > 0.0 ) THEN
	380	zfrac_p1n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / &
	381	& (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
	382	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
	383	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
	384	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
	385	zfrac_p2n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / &
	386	& (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
	387	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
	388	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
	389	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
	390	zfrac_p3n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / &
	391	& (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
	392	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
	393	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
	394	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
	395	zfrac_p4n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / &
	396	& (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
	397	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
	398	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
	399	& logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
	400	ELSE
	401	zfrac_p1n = 0.25
	402	zfrac_p2n = 0.25
	403	zfrac_p3n = 0.25
	404	zfrac_p4n = 0.25
	405	ENDIF
	406	ELSE
	407	zfrac_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) / &
	408	& (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
	409	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
	410	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
	411	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
	412	zfrac_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) / &
	413	& (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
	414	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
	415	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
	416	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
	417	zfrac_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) / &
	418	& (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
	419	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
	420	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
	421	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
	422	zfrac_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) / &
	423	& (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
	424	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
	425	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
	426	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
	427	ENDIF
	428
	429	zrat_p1c_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
	430	zrat_p1p_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
	431	zrat_p1s_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
	432	zrat_p2c_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)
	433	zrat_p2p_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)
	434	zrat_p3c_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)
	435	zrat_p3p_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)
	436	zrat_p4c_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)
	437	zrat_p4p_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)
	438
	439	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p1n
	440	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p2n
	441	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p3n
	442	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p4n
	443	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1c_p1n
	444	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1p_p1n
	445	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1s_p1n
	446	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) * zrat_p2c_p2n
	447	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) * zrat_p2p_p2n
	448	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) * zrat_p3c_p3n
	449	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) * zrat_p3p_p3n
	450	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) * zrat_p4c_p4n
	451	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) * zrat_p4p_p4n
	452
	453	ENDIF
	454
	455	! Nitrogen nutrient from balancing scheme
	456	! Split between nitrate and ammonium according to current ratios
	457	IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) > 0.0 ) .AND. &
	458	& ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) > 0.0 ) ) THEN
	459	zfrac_n3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) / &
	460	& (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) + trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n))
	461	zfrac_n4n = 1.0 - zfrac_n3n
	462	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n3n
	463	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n4n
	464	ENDIF
	465
	466	! Nitrogen zooplankton from balancing scheme
	467	! Split according to current ratios
	468	! Update carbon and phosphorus according to current ratios
	469	IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) > 0.0 ) .AND. &
	470	& ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) > 0.0 ) .AND. &
	471	& ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) > 0.0 ) ) THEN
	472	zfrac_z4n = (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) / &
	473	& ((trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) + &
	474	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
	475	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n))
	476	zfrac_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) / &
	477	& ((trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) + &
	478	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
	479	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n))
	480	zfrac_z6n = 1.0 - zfrac_z4n - zfrac_z5n
	481	zrat_z4c_z4n = 1.0 / 0.0126
	482	zrat_z5c_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n)
	483	zrat_z5p_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n)
	484	zrat_z6c_z6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n)
	485	zrat_z6p_z6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n)
	486	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z5n
	487	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z6n
	488	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z4n * zrat_z4c_z4n
	489	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) * zrat_z5c_z5n
	490	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) * zrat_z6c_z6n
	491	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) * zrat_z5p_z5n
	492	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) * zrat_z6p_z6n
	493	ENDIF
	494
	495	! Nitrogen detritus from balancing scheme
	496	! Split according to current ratios
	497	IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) > 0.0 ) .AND. &
	498	& ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) > 0.0 ) .AND. &
	499	& ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) > 0.0 ) ) THEN
	500	zfrac_r4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) / &
	501	& (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
	502	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
	503	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
	504	zfrac_r6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) / &
	505	& (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
	506	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
	507	& trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
	508	zfrac_r8n = 1.0 - zfrac_r4n - zfrac_r6n
	509	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r4n
	510	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r6n
	511	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r8n
	512	ENDIF
	513
	514	! DIC straight from balancing scheme
	515	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3c) = outincs(ji,jj,jk,i_tracer(5))
	516
	517	! Alkalinity straight from balancing scheme
	518	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3a) = outincs(ji,jj,jk,i_tracer(6))
	519
	520	END DO
	521	END DO
	522	END DO
	523
	524	ENDIF
[6983]	525	!
[9331]	526	! Now set MLD to bottom of a level and propagate chlorophyll incs equally through mixed layer
	527	! If balancing, should really relate this back to phytoplankton, but stick with this for now
[6983]	528	!
	529	DO jj = 1, jpj
	530	DO ji = 1, jpi
	531	!
	532	jkmax = jpk-1
	533	DO jk = jpk-1, 1, -1
	534	IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. &
	535	& ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
	536	zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
	537	jkmax = jk
	538	ENDIF
	539	END DO
	540	!
	541	DO jk = 2, jkmax
[7716]	542	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
	543	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
	544	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
	545	logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
[6983]	546	END DO
	547	!
	548	END DO
	549	END DO
	550
	551	END SUBROUTINE asm_logchl_bal_ersem
	552
	553	#else
	554	!!----------------------------------------------------------------------
	555	!! Default option : Empty routine
	556	!!----------------------------------------------------------------------
	557	CONTAINS
[7803]	558	SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, &
	559	& k_maxchlinc, logchl_bkginc, logchl_balinc )
	560	LOGICAL :: ld_logchlpftinc
[6983]	561	INTEGER :: npfts
	562	INTEGER :: mld_choice_bgc
[7803]	563	REAL :: k_maxchlinc
[6983]	564	REAL :: logchl_bkginc(:,:,:)
[7716]	565	REAL :: logchl_balinc(:,:,:,:)
[6983]	566	WRITE(,) 'asm_logchl_bal_ersem: You should not have seen this print! error?'
	567	END SUBROUTINE asm_logchl_bal_ersem
	568	#endif
	569
	570	!!======================================================================
	571	END MODULE asmlogchlbal_ersem

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

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

Context Navigation

source: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90 @ 9331

Download in other formats: