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.

asmbgc.F90 in branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 10622

Last change on this file since 10622 was 10622, checked in by dford, 5 years ago
Implement biogeochemistry assimilation for FABM-ERSEM.
File size: 139.3 KB

Rev	Line
[9296]	1	MODULE asmbgc
[9297]	2	!!===========================================================================
[9296]	3	!! * MODULE asmbgc *
	4	!! Routines and declarations for biogeochemical assimilation
[9297]	5	!!===========================================================================
[9296]	6	!! History : 3.6 ! 2018-02 (D. Ford) Adapted from existing modules
[9297]	7	!!---------------------------------------------------------------------------
[9296]	8	!! 'key_asminc' : assimilation increment interface
	9	!! 'key_top' : passive tracers
	10	!! 'key_hadocc' : HadOCC model
	11	!! 'key_medusa' : MEDUSA model
	12	!! 'key_roam' : MEDUSA extras for carbonate cycle
[10622]	13	!! 'key_fabm' : ERSEM model coupled via FABM
[9296]	14	!! 'key_karaml' : Kara mixed layer depth
[9297]	15	!!---------------------------------------------------------------------------
	16	!! asm_bgc_check_options : check bgc assimilation options
[9296]	17	!! asm_bgc_init_incs : initialise bgc increments
	18	!! asm_bgc_init_bkg : initialise bgc background
[9297]	19	!! asm_bgc_bal_wri : write out bgc balancing increments
[9296]	20	!! asm_bgc_bkg_wri : write out bgc background
[9326]	21	!! phyto2d_asm_inc : apply the ocean colour increments
	22	!! phyto3d_asm_inc : apply the 3D phytoplankton increments
[9322]	23	!! pco2_asm_inc : apply the pCO2/fCO2 increments
	24	!! ph_asm_inc : apply the pH increments
	25	!! bgc3d_asm_inc : apply the generic 3D BGC increments
[9297]	26	!!---------------------------------------------------------------------------
	27	USE par_kind, ONLY: & ! kind parameters
	28	& wp
	29	USE par_oce, ONLY: & ! domain array sizes
	30	& jpi, &
	31	& jpj, &
	32	& jpk
	33	USE iom ! i/o
	34	USE oce, ONLY: & ! active tracer variables
	35	& tsn
	36	USE zdfmxl, ONLY : & ! mixed layer depth
[9296]	37	#if defined key_karaml
[9297]	38	& hmld_kara, &
	39	& ln_kara, &
[9296]	40	#endif
[10622]	41	& hmld, &
	42	& hmld_tref, &
[9297]	43	& hmlp, &
	44	& hmlpt
	45	USE asmpar, ONLY: & ! assimilation parameters
	46	& c_asmbkg, &
	47	& c_asmbal, &
	48	& nitbkg_r, &
[10622]	49	& nitavgbkg_r, &
[9297]	50	& nitdin_r, &
	51	& nitiaustr_r, &
	52	& nitiaufin_r
[10622]	53	USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA
	54	& asm_phyto2d_bal_medusa
	55	USE asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC
	56	& asm_phyto2d_bal_hadocc
	57	USE asmphyto2dbal_ersem, ONLY: & ! phyto2d balancing for ERSEM
	58	& asm_phyto2d_bal_ersem
	59	USE asmpco2bal, ONLY: & ! pCO2 balancing
	60	& asm_pco2_bal
[9296]	61	#if defined key_top
[9297]	62	USE trc, ONLY: & ! passive tracer variables
	63	& trn, &
[10622]	64	& trb, &
	65	& nittrc000
[9297]	66	USE par_trc, ONLY: & ! passive tracer parameters
[9296]	67	& jptra
	68	#endif
[10141]	69	#if defined key_medusa
[9431]	70	USE sms_medusa ! MEDUSA parameters
[9297]	71	USE par_medusa ! MEDUSA parameters
[9381]	72	USE mocsy_wrapper, ONLY: & ! MEDUSA carbonate system
	73	& mocsy_interface
	74	USE mocsy_mainmod, ONLY: & ! fCO2 to pCO2 conversion
[9297]	75	& f2pCO2
	76	USE sms_medusa, ONLY: & ! MEDUSA diagnostic variables
	77	& pgrow_avg, &
	78	& ploss_avg, &
	79	& phyt_avg, &
	80	& mld_max
[9296]	81	#elif defined key_hadocc
[9297]	82	USE par_hadocc ! HadOCC parameters
[9326]	83	USE had_bgc_const ! HadOCC parameters
[9297]	84	USE trc, ONLY: & ! HadOCC diagnostic variables
	85	& pgrow_avg, &
	86	& ploss_avg, &
	87	& phyt_avg, &
	88	& mld_max, &
	89	& HADOCC_CHL
	90	USE had_bgc_const, ONLY: & ! HadOCC C:Chl ratio
	91	& cchl_p
[10622]	92	#elif defined key_fabm
	93	USE par_fabm ! FABM-ERSEM parameters
	94	USE trc, ONLY: & ! FABM-ERSEM diagnostic variables
	95	& pgrow_avg, &
	96	& ploss_avg, &
	97	& phyt_avg, &
	98	& mld_max
[9296]	99	#endif
	100
	101	IMPLICIT NONE
	102	PRIVATE
	103
[9322]	104	PUBLIC asm_bgc_check_options ! called by asm_inc_init in asminc.F90
	105	PUBLIC asm_bgc_init_incs ! called by asm_inc_init in asminc.F90
	106	PRIVATE asm_bgc_read_incs_2d ! called by asm_bgc_init_incs
	107	PRIVATE asm_bgc_read_incs_3d ! called by asm_bgc_init_incs
	108	PUBLIC asm_bgc_init_bkg ! called by asm_inc_init in asminc.F90
	109	PUBLIC asm_bgc_bal_wri ! called by nemo_gcm in nemogcm.F90
[10622]	110	PUBLIC asm_bgc_bkg_alloc ! called by asm_bkg_wri in asmbkg.F90
	111	PUBLIC asm_bgc_bkg_tavg ! called by asm_bkg_wri in asmbkg.F90
[9322]	112	PUBLIC asm_bgc_bkg_wri ! called by asm_bkg_wri in asmbkg.F90
[9326]	113	PUBLIC phyto2d_asm_inc ! called by bgc_asm_inc in asminc.F90
	114	PUBLIC phyto3d_asm_inc ! called by bgc_asm_inc in asminc.F90
[9322]	115	PUBLIC pco2_asm_inc ! called by bgc_asm_inc in asminc.F90
	116	PUBLIC ph_asm_inc ! called by bgc_asm_inc in asminc.F90
	117	PUBLIC bgc3d_asm_inc ! called by bgc_asm_inc in asminc.F90
[9296]	118
[9322]	119	LOGICAL, PUBLIC :: ln_balwri = .FALSE. !: No output of balancing incs
	120	LOGICAL, PUBLIC :: ln_phytobal = .FALSE. !: No phytoplankton balancing
	121	LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE. !: No surface total log10(chlorophyll) increment
	122	LOGICAL, PUBLIC :: ln_slchldiainc = .FALSE. !: No surface diatom log10(chlorophyll) increment
	123	LOGICAL, PUBLIC :: ln_slchlnoninc = .FALSE. !: No surface non-diatom log10(chlorophyll) increment
[10622]	124	LOGICAL, PUBLIC :: ln_slchlnaninc = .FALSE. !: No surface nanophyto log10(chlorophyll) increment
	125	LOGICAL, PUBLIC :: ln_slchlpicinc = .FALSE. !: No surface picophyto log10(chlorophyll) increment
	126	LOGICAL, PUBLIC :: ln_slchldininc = .FALSE. !: No surface dinoflag log10(chlorophyll) increment
[9322]	127	LOGICAL, PUBLIC :: ln_schltotinc = .FALSE. !: No surface total chlorophyll increment
	128	LOGICAL, PUBLIC :: ln_slphytotinc = .FALSE. !: No surface total log10(phyto C) increment
	129	LOGICAL, PUBLIC :: ln_slphydiainc = .FALSE. !: No surface diatom log10(phyto C) increment
	130	LOGICAL, PUBLIC :: ln_slphynoninc = .FALSE. !: No surface non-diatom log10(phyto C) increment
	131	LOGICAL, PUBLIC :: ln_spco2inc = .FALSE. !: No surface pCO2 increment
	132	LOGICAL, PUBLIC :: ln_sfco2inc = .FALSE. !: No surface fCO2 increment
	133	LOGICAL, PUBLIC :: ln_plchltotinc = .FALSE. !: No profile total log10(chlorophyll) increment
	134	LOGICAL, PUBLIC :: ln_pchltotinc = .FALSE. !: No profile total chlorophyll increment
	135	LOGICAL, PUBLIC :: ln_pno3inc = .FALSE. !: No profile nitrate increment
	136	LOGICAL, PUBLIC :: ln_psi4inc = .FALSE. !: No profile silicate increment
	137	LOGICAL, PUBLIC :: ln_pdicinc = .FALSE. !: No profile dissolved inorganic carbon increment
	138	LOGICAL, PUBLIC :: ln_palkinc = .FALSE. !: No profile alkalinity increment
	139	LOGICAL, PUBLIC :: ln_pphinc = .FALSE. !: No profile pH increment
	140	LOGICAL, PUBLIC :: ln_po2inc = .FALSE. !: No profile oxygen increment
[9296]	141
[9297]	142	INTEGER, PUBLIC :: mld_choice_bgc = 1 !: choice of mld for bgc assimilation
	143	! 1) hmld - Turbocline/mixing depth
	144	! [W points]
	145	! 2) hmlp - Density criterion
	146	! (0.01 kg/m^3 change from 10m)
	147	! [W points]
	148	! 3) hmld_kara - Kara MLD
	149	! [Interpolated]
	150	! 4) hmld_tref - Temperature criterion
	151	! (0.2 K change from surface)
	152	! [T points]
	153	! 5) hmlpt - Density criterion
	154	! (0.01 kg/m^3 change from 10m)
	155	! [T points]
	156
	157	REAL(wp), PUBLIC :: rn_maxchlinc = -999.0 !: maximum absolute non-log chlorophyll inc
	158	! <= 0 no maximum applied (switch turned off)
	159	! > 0 absolute chl inc capped at this value
	160
	161	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchltot_bkginc ! slchltot inc
[9322]	162	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchldia_bkginc ! slchldia inc
	163	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchlnon_bkginc ! slchlnon inc
[10622]	164	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchlnan_bkginc ! slchlnan inc
	165	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchlpic_bkginc ! slchlpic inc
	166	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchldin_bkginc ! slchldin inc
[9322]	167	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: schltot_bkginc ! schltot inc
	168	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphytot_bkginc ! slphytot inc
	169	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphydia_bkginc ! slphydia inc
	170	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphynon_bkginc ! slphynon inc
	171	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sfco2_bkginc ! sfco2 inc
	172	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: spco2_bkginc ! spco2 inc
	173	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: plchltot_bkginc ! plchltot inc
	174	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pchltot_bkginc ! pchltot inc
	175	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pno3_bkginc ! pno3 inc
	176	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: psi4_bkginc ! psi4 inc
	177	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pdic_bkginc ! pdic inc
	178	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: palk_bkginc ! palk inc
	179	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pph_bkginc ! pph inc
	180	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: po2_bkginc ! po2 inc
[9326]	181	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto2d_balinc ! Balancing incs from ocean colour
	182	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto3d_balinc ! Balancing incs from p(l)chltot
[9322]	183	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc ! Balancing incs from spco2/sfco2
	184	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ph_balinc ! Balancing incs from pph
[9297]	185
	186	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pgrow_avg_bkg ! Background phyto growth
	187	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ploss_avg_bkg ! Background phyto loss
	188	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: phyt_avg_bkg ! Background phyto conc
	189	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mld_max_bkg ! Background max MLD
	190	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg ! Background tracer state
[9322]	191	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: chl_bkg ! Background chl
	192	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cchl_p_bkg ! Background c:chl
[10622]	193	#if defined key_fabm
	194	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: totalk_bkg ! Background total alkalinity
	195	#endif
[9296]	196
[10622]	197	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: pgrow_avg_tavg ! Time average of pgrow_avg
	198	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: ploss_avg_tavg ! Time average of ploss_avg
	199	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: phyt_avg_tavg ! Time average of phyt_avg
	200	REAL(wp), DIMENSION(:,:,:,:), SAVE, ALLOCATABLE :: trn_tavg ! Time average of trn
	201	#if defined key_hadocc
	202	REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE :: HADOCC_CHL_tavg ! Time average of HADOCC_CHL
	203	REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE :: cchl_p_tavg ! Time average of cchl_p
	204	#elif defined key_fabm
	205	REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE :: totalk_tavg ! Time average of O3_TA
	206	#endif
	207
[9381]	208	# include "domzgr_substitute.h90"
	209
[9296]	210	CONTAINS
	211
	212	SUBROUTINE asm_bgc_check_options
[9297]	213	!!------------------------------------------------------------------------
[9296]	214	!! * ROUTINE asm_bgc_check_options *
	215	!!
	216	!! ** Purpose : check validity of biogeochemical assimilation options
	217	!!
	218	!! ** Method : check settings of logicals
	219	!!
	220	!! ** Action : call ctl_stop/ctl_warn if required
	221	!!
	222	!! References : asm_inc_init
[9297]	223	!!------------------------------------------------------------------------
[9296]	224
[10622]	225	#if ! defined key_top \|\| ( ! defined key_hadocc && ! defined key_medusa && ! defined key_fabm )
[9322]	226	CALL ctl_stop( ' Attempting to assimilate biogeochemical observations', &
	227	& ' but no compatible biogeochemical model is available' )
	228	#endif
	229
	230	#if defined key_hadocc
	231	IF ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slphydiainc .OR. &
	232	& ln_slphynoninc .OR. ln_psi4inc .OR. ln_pphinc .OR. ln_po2inc ) THEN
	233	CALL ctl_stop( ' Cannot assimilate PFTs, Si4, pH or O2 into HadOCC' )
[9296]	234	ENDIF
[9322]	235	#endif
[9296]	236
[10141]	237	#if defined key_medusa
	238	IF ( .NOT. ln_foam_medusa ) THEN
	239	CALL ctl_stop( ' MEDUSA data assimilation options not turned on: set ln_foam_medusa' )
	240	ENDIF
	241	#endif
	242
[10622]	243	#if defined key_fabm
	244	IF ( ln_pphinc ) THEN
	245	CALL ctl_stop( ' Cannot currently assimilate pH into FABM-ERSEM' )
	246	ENDIF
	247	#endif
	248
[9322]	249	IF ( ( ln_phytobal ).AND. &
	250	& ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
[10622]	251	& ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. &
	252	& ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. &
	253	& ( .NOT. ln_schltotinc ).AND.( .NOT. ln_slphytotinc ).AND. &
	254	& ( .NOT. ln_slphydiainc ).AND.( .NOT. ln_slphynoninc ) ) THEN
[9322]	255	CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', &
	256	& ' if not assimilating ocean colour,', &
	257	& ' so ln_phytobal will be set to .false.')
	258	ln_phytobal = .FALSE.
[9296]	259	ENDIF
	260
[9381]	261	IF ( ( ln_balwri ).AND. &
	262	& ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
[10622]	263	& ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. &
	264	& ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. &
	265	& ( .NOT. ln_schltotinc ).AND. &
[9381]	266	& ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
	267	& ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. &
	268	& ( .NOT. ln_pchltotinc ).AND.( .NOT. ln_pphinc ).AND. &
	269	& ( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ) ) THEN
	270	CALL ctl_warn( ' Balancing increments not being calculated', &
	271	& ' so ln_balwri will be set to .false.')
[9322]	272	ln_balwri = .FALSE.
	273	ENDIF
	274
	275	IF ( ln_spco2inc .AND. ln_sfco2inc ) THEN
[9296]	276	CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
	277	ENDIF
	278
[9322]	279	IF ( ln_slchltotinc .AND. ln_schltotinc ) THEN
	280	CALL ctl_stop( ' Can only assimilate surface log10(chlorophyll) or chlorophyll, not both' )
	281	ENDIF
	282
	283	IF ( ln_plchltotinc .AND. ln_pchltotinc ) THEN
	284	CALL ctl_stop( ' Can only assimilate profile log10(chlorophyll) or chlorophyll, not both' )
	285	ENDIF
	286
	287	IF ( ( ln_slchltotinc .OR. ln_schltotinc ) .AND. &
[10622]	288	& ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slchlnaninc .OR. &
	289	& ln_slchlpicinc .OR. ln_slchldininc ) ) THEN
[9322]	290	CALL ctl_stop( ' Can only assimilate total or PFT surface chlorophyll, not both' )
	291	ENDIF
	292
[9435]	293	IF ( ln_phytobal .AND. &
	294	& ( ( ln_slchlnoninc .AND.( .NOT. ln_slchldiainc ) ).OR. &
[10622]	295	#if defined key_fabm
	296	& ( ( ln_slchldiainc .OR. ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc ) .AND. &
	297	& ( ( .NOT. ln_slchldiainc ) .OR. ( .NOT. ln_slchlnaninc ) .OR. &
	298	& ( .NOT. ln_slchlpicinc ) .OR. ( .NOT. ln_slchldininc ) ) ) .OR. &
	299	#endif
[9435]	300	& ( ln_slchldiainc .AND.( .NOT. ln_slchlnoninc ) ) ) ) THEN
	301	CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
	302	& ' unless assimilating all model PFTs,')
	303	ENDIF
	304
[9322]	305	IF ( ln_slphytotinc .AND. ( ln_slphydiainc .OR. ln_slphynoninc ) ) THEN
	306	CALL ctl_stop( ' Can only assimilate total or PFT surface phytoplankton carbon, not both' )
	307	ENDIF
	308
[9435]	309	IF ( ln_phytobal .AND. &
	310	& ( ( ln_slphynoninc .AND.( .NOT. ln_slphydiainc ) ).OR. &
	311	& ( ln_slphydiainc .AND.( .NOT. ln_slphynoninc ) ) ) ) THEN
	312	CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
	313	& ' unless assimilating all model PFTs,')
	314	ENDIF
	315
[9296]	316	END SUBROUTINE asm_bgc_check_options
	317
[9297]	318	!!===========================================================================
	319	!!===========================================================================
	320	!!===========================================================================
[9296]	321
	322	SUBROUTINE asm_bgc_init_incs( knum )
[9297]	323	!!------------------------------------------------------------------------
[9296]	324	!! * ROUTINE asm_bgc_init_incs *
	325	!!
	326	!! ** Purpose : initialise bgc increments
	327	!!
	328	!! ** Method : allocate arrays and read increments from file
	329	!!
	330	!! ** Action : allocate arrays and read increments from file
	331	!!
	332	!! References : asm_inc_init
[9297]	333	!!------------------------------------------------------------------------
[9296]	334	!!
	335	INTEGER, INTENT(in ) :: knum ! i/o unit of increments file
	336	!!
[9297]	337	!!------------------------------------------------------------------------
[9296]	338
[9322]	339	! Allocate and read increments
	340
[9296]	341	IF ( ln_slchltotinc ) THEN
	342	ALLOCATE( slchltot_bkginc(jpi,jpj) )
[9322]	343	CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc )
	344	ENDIF
	345
	346	IF ( ln_slchldiainc ) THEN
	347	ALLOCATE( slchldia_bkginc(jpi,jpj) )
	348	CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc )
	349	ENDIF
	350
	351	IF ( ln_slchlnoninc ) THEN
	352	ALLOCATE( slchlnon_bkginc(jpi,jpj) )
	353	CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc )
	354	ENDIF
	355
[10622]	356	IF ( ln_slchlnaninc ) THEN
	357	ALLOCATE( slchlnan_bkginc(jpi,jpj) )
	358	CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnan', slchlnan_bkginc )
	359	ENDIF
	360
	361	IF ( ln_slchlpicinc ) THEN
	362	ALLOCATE( slchlpic_bkginc(jpi,jpj) )
	363	CALL asm_bgc_read_incs_2d( knum, 'bckinslchlpic', slchlpic_bkginc )
	364	ENDIF
	365
	366	IF ( ln_slchldininc ) THEN
	367	ALLOCATE( slchldin_bkginc(jpi,jpj) )
	368	CALL asm_bgc_read_incs_2d( knum, 'bckinslchldin', slchldin_bkginc )
	369	ENDIF
	370
[9322]	371	IF ( ln_schltotinc ) THEN
	372	ALLOCATE( schltot_bkginc(jpi,jpj) )
	373	CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc )
	374	ENDIF
	375
	376	IF ( ln_slphytotinc ) THEN
	377	ALLOCATE( slphytot_bkginc(jpi,jpj) )
	378	CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc )
	379	ENDIF
	380
	381	IF ( ln_slphydiainc ) THEN
	382	ALLOCATE( slphydia_bkginc(jpi,jpj) )
	383	CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc )
	384	ENDIF
	385
	386	IF ( ln_slphynoninc ) THEN
	387	ALLOCATE( slphynon_bkginc(jpi,jpj) )
	388	CALL asm_bgc_read_incs_2d( knum, 'bckinslphynon', slphynon_bkginc )
	389	ENDIF
	390
	391	IF ( ln_sfco2inc ) THEN
	392	ALLOCATE( sfco2_bkginc(jpi,jpj) )
	393	CALL asm_bgc_read_incs_2d( knum, 'bckinsfco2', sfco2_bkginc )
	394	ENDIF
	395
	396	IF ( ln_spco2inc ) THEN
	397	ALLOCATE( spco2_bkginc(jpi,jpj) )
	398	CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc )
	399	ENDIF
	400
	401	IF ( ln_plchltotinc ) THEN
	402	ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) )
	403	CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc )
	404	ENDIF
	405
	406	IF ( ln_pchltotinc ) THEN
	407	ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) )
	408	CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc )
	409	ENDIF
	410
	411	IF ( ln_pno3inc ) THEN
	412	ALLOCATE( pno3_bkginc(jpi,jpj,jpk) )
	413	CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc )
	414	ENDIF
	415
	416	IF ( ln_psi4inc ) THEN
	417	ALLOCATE( psi4_bkginc(jpi,jpj,jpk) )
	418	CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc )
	419	ENDIF
	420
	421	IF ( ln_pdicinc ) THEN
	422	ALLOCATE( pdic_bkginc(jpi,jpj,jpk) )
	423	CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc )
	424	ENDIF
	425
	426	IF ( ln_palkinc ) THEN
	427	ALLOCATE( palk_bkginc(jpi,jpj,jpk) )
	428	CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc )
	429	ENDIF
	430
	431	IF ( ln_pphinc ) THEN
	432	ALLOCATE( pph_bkginc(jpi,jpj,jpk) )
	433	CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc )
	434	ENDIF
	435
	436	IF ( ln_po2inc ) THEN
	437	ALLOCATE( po2_bkginc(jpi,jpj,jpk) )
	438	CALL asm_bgc_read_incs_3d( knum, 'bckinpo2', po2_bkginc )
	439	ENDIF
	440
	441	! Allocate balancing increments
	442
	443	IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
[10622]	444	& ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
[9322]	445	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
	446	& ln_slphynoninc ) THEN
[9296]	447	#if defined key_top
[9326]	448	ALLOCATE( phyto2d_balinc(jpi,jpj,jpk,jptra) )
	449	phyto2d_balinc(:,:,:,:) = 0.0
[9322]	450	#else
	451	CALL ctl_stop( ' key_top must be set for balancing increments' )
[9296]	452	#endif
	453	ENDIF
[9326]	454
	455	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
	456	#if defined key_top
	457	ALLOCATE( phyto3d_balinc(jpi,jpj,jpk,jptra) )
	458	phyto3d_balinc(:,:,:,:) = 0.0
	459	#else
	460	CALL ctl_stop( ' key_top must be set for balancing increments' )
	461	#endif
	462	ENDIF
[9322]	463
[9296]	464	IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
	465	#if defined key_top
[9322]	466	ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) )
	467	pco2_balinc(:,:,:,:) = 0.0
	468	#else
	469	CALL ctl_stop( ' key_top must be set for balancing increments' )
[9296]	470	#endif
	471	ENDIF
	472
[9322]	473	IF ( ln_pphinc ) THEN
	474	#if defined key_top
	475	ALLOCATE( ph_balinc(jpi,jpj,jpk,jptra) )
	476	ph_balinc(:,:,:,:) = 0.0
	477	#else
	478	CALL ctl_stop( ' key_top must be set for balancing increments' )
	479	#endif
[9296]	480	ENDIF
	481
	482	END SUBROUTINE asm_bgc_init_incs
	483
[9297]	484	!!===========================================================================
	485	!!===========================================================================
	486	!!===========================================================================
[9296]	487
[9322]	488	SUBROUTINE asm_bgc_read_incs_2d( knum, cd_bgcname, p_incs )
	489	!!------------------------------------------------------------------------
	490	!! * ROUTINE asm_bgc_init_incs *
	491	!!
	492	!! ** Purpose : read 2d bgc increments
	493	!!
	494	!! ** Method : read increments from file
	495	!!
	496	!! ** Action : read increments from file
	497	!!
	498	!! References : asm_inc_init
	499	!!------------------------------------------------------------------------
	500	!!
	501	INTEGER, INTENT(in ) :: knum ! i/o unit
[10574]	502	CHARACTER(LEN=*), INTENT(in ) :: cd_bgcname ! variable
[9322]	503	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: p_incs ! increments
	504	!!
	505	!!------------------------------------------------------------------------
	506
	507	! Initialise
	508	p_incs(:,:) = 0.0
	509
	510	! read from file
	511	CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 )
	512
	513	! Apply the masks
	514	p_incs(:,:) = p_incs(:,:) * tmask(:,:,1)
	515
	516	! Set missing increments to 0.0 rather than 1e+20
	517	! to allow for differences in masks
	518	WHERE( ABS( p_incs(:,:) ) > 1.0e+10 ) p_incs(:,:) = 0.0
	519
	520	END SUBROUTINE asm_bgc_read_incs_2d
	521
	522	!!===========================================================================
	523	!!===========================================================================
	524	!!===========================================================================
	525
	526	SUBROUTINE asm_bgc_read_incs_3d( knum, cd_bgcname, p_incs )
	527	!!------------------------------------------------------------------------
	528	!! * ROUTINE asm_bgc_init_incs *
	529	!!
	530	!! ** Purpose : read 3d bgc increments
	531	!!
	532	!! ** Method : read increments from file
	533	!!
	534	!! ** Action : read increments from file
	535	!!
	536	!! References : asm_inc_init
	537	!!------------------------------------------------------------------------
	538	!!
	539	INTEGER, INTENT(in ) :: knum ! i/o unit
[10574]	540	CHARACTER(LEN=*), INTENT(in ) :: cd_bgcname ! variable
[9322]	541	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: p_incs ! increments
	542	!!
	543	!!------------------------------------------------------------------------
	544
	545	! Initialise
	546	p_incs(:,:,:) = 0.0
	547
	548	! read from file
	549	CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 )
	550
	551	! Apply the masks
	552	p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:)
	553
	554	! Set missing increments to 0.0 rather than 1e+20
	555	! to allow for differences in masks
	556	WHERE( ABS( p_incs(:,:,:) ) > 1.0e+10 ) p_incs(:,:,:) = 0.0
	557
	558	END SUBROUTINE asm_bgc_read_incs_3d
	559
	560	!!===========================================================================
	561	!!===========================================================================
	562	!!===========================================================================
	563
[9296]	564	SUBROUTINE asm_bgc_init_bkg
[9297]	565	!!------------------------------------------------------------------------
[9296]	566	!! * ROUTINE asm_bgc_init_bkg *
	567	!!
	568	!! ** Purpose : initialise bgc background
	569	!!
	570	!! ** Method : allocate arrays and read background from file
	571	!!
	572	!! ** Action : allocate arrays and read background from file
	573	!!
	574	!! References : asm_inc_init
[9297]	575	!!------------------------------------------------------------------------
[9296]	576	!!
	577	INTEGER :: inum ! i/o unit of background file
	578	INTEGER :: jt ! loop counter
	579	!!
[9297]	580	!!------------------------------------------------------------------------
[9296]	581
[10622]	582	#if defined key_top && ( defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm )
[9322]	583	IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
[10622]	584	& ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
[9322]	585	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
	586	& ln_slphynoninc .OR. ln_plchltotinc .OR. ln_pchltotinc ) THEN
[9296]	587
	588	ALLOCATE( pgrow_avg_bkg(jpi,jpj) )
	589	ALLOCATE( ploss_avg_bkg(jpi,jpj) )
	590	ALLOCATE( phyt_avg_bkg(jpi,jpj) )
	591	ALLOCATE( mld_max_bkg(jpi,jpj) )
	592	ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
	593	pgrow_avg_bkg(:,:) = 0.0
	594	ploss_avg_bkg(:,:) = 0.0
	595	phyt_avg_bkg(:,:) = 0.0
	596	mld_max_bkg(:,:) = 0.0
	597	tracer_bkg(:,:,:,:) = 0.0
	598
	599	#if defined key_hadocc
[9322]	600	ALLOCATE( chl_bkg(jpi,jpj,jpk) )
	601	ALLOCATE( cchl_p_bkg(jpi,jpj,jpk) )
	602	chl_bkg(:,:,:) = 0.0
	603	cchl_p_bkg(:,:,:) = 0.0
[9296]	604	#endif
	605
	606	!--------------------------------------------------------------------
[9322]	607	! Read background variables for phytoplankton assimilation
[9296]	608	! Some only required if performing balancing
	609	!--------------------------------------------------------------------
	610
	611	CALL iom_open( c_asmbkg, inum )
	612
	613	#if defined key_hadocc
	614	CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl', chl_bkg )
	615	CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg )
[9322]	616	chl_bkg(:,:,:) = chl_bkg(:,:,:) * tmask(:,:,:)
	617	cchl_p_bkg(:,:,:) = cchl_p_bkg(:,:,:) * tmask(:,:,:)
[10141]	618	#elif defined key_medusa
[9296]	619	CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) )
	620	CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) )
[9861]	621	CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) )
	622	CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
	623	CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
[10622]	624	#elif defined key_fabm
	625	CALL iom_get( inum, jpdom_autoglo, 'ersem_chl1', tracer_bkg(:,:,:,jp_fabm_chl1) )
	626	CALL iom_get( inum, jpdom_autoglo, 'ersem_chl2', tracer_bkg(:,:,:,jp_fabm_chl2) )
	627	CALL iom_get( inum, jpdom_autoglo, 'ersem_chl3', tracer_bkg(:,:,:,jp_fabm_chl3) )
	628	CALL iom_get( inum, jpdom_autoglo, 'ersem_chl4', tracer_bkg(:,:,:,jp_fabm_chl4) )
	629	CALL iom_get( inum, jpdom_autoglo, 'ersem_p1c', tracer_bkg(:,:,:,jp_fabm_p1c) )
	630	CALL iom_get( inum, jpdom_autoglo, 'ersem_p1n', tracer_bkg(:,:,:,jp_fabm_p1n) )
	631	CALL iom_get( inum, jpdom_autoglo, 'ersem_p1p', tracer_bkg(:,:,:,jp_fabm_p1p) )
	632	CALL iom_get( inum, jpdom_autoglo, 'ersem_p1s', tracer_bkg(:,:,:,jp_fabm_p1s) )
	633	CALL iom_get( inum, jpdom_autoglo, 'ersem_p2c', tracer_bkg(:,:,:,jp_fabm_p2c) )
	634	CALL iom_get( inum, jpdom_autoglo, 'ersem_p2n', tracer_bkg(:,:,:,jp_fabm_p2n) )
	635	CALL iom_get( inum, jpdom_autoglo, 'ersem_p2p', tracer_bkg(:,:,:,jp_fabm_p2p) )
	636	CALL iom_get( inum, jpdom_autoglo, 'ersem_p3c', tracer_bkg(:,:,:,jp_fabm_p3c) )
	637	CALL iom_get( inum, jpdom_autoglo, 'ersem_p3n', tracer_bkg(:,:,:,jp_fabm_p3n) )
	638	CALL iom_get( inum, jpdom_autoglo, 'ersem_p3p', tracer_bkg(:,:,:,jp_fabm_p3p) )
	639	CALL iom_get( inum, jpdom_autoglo, 'ersem_p4c', tracer_bkg(:,:,:,jp_fabm_p4c) )
	640	CALL iom_get( inum, jpdom_autoglo, 'ersem_p4n', tracer_bkg(:,:,:,jp_fabm_p4n) )
	641	CALL iom_get( inum, jpdom_autoglo, 'ersem_p4p', tracer_bkg(:,:,:,jp_fabm_p4p) )
[9296]	642	#endif
	643
[9322]	644	IF ( ln_phytobal ) THEN
[9296]	645
	646	CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg )
	647	CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg )
	648	CALL iom_get( inum, jpdom_autoglo, 'phyt_avg', phyt_avg_bkg )
	649	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
	650	pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1)
	651	ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1)
	652	phyt_avg_bkg(:,:) = phyt_avg_bkg(:,:) * tmask(:,:,1)
	653	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
	654
	655	#if defined key_hadocc
	656	CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) )
	657	CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) )
	658	CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) )
	659	CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) )
	660	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
	661	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
[10141]	662	#elif defined key_medusa
[9296]	663	CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) )
	664	CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) )
	665	CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) )
	666	CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) )
	667	CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) )
	668	CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) )
	669	CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) )
	670	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
	671	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
	672	CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
[10622]	673	#elif defined key_fabm
	674	CALL iom_get( inum, jpdom_autoglo, 'ersem_z4c', tracer_bkg(:,:,:,jp_fabm_z4c) )
	675	CALL iom_get( inum, jpdom_autoglo, 'ersem_z5c', tracer_bkg(:,:,:,jp_fabm_z5c) )
	676	CALL iom_get( inum, jpdom_autoglo, 'ersem_z5n', tracer_bkg(:,:,:,jp_fabm_z5n) )
	677	CALL iom_get( inum, jpdom_autoglo, 'ersem_z5p', tracer_bkg(:,:,:,jp_fabm_z5p) )
	678	CALL iom_get( inum, jpdom_autoglo, 'ersem_z6c', tracer_bkg(:,:,:,jp_fabm_z6c) )
	679	CALL iom_get( inum, jpdom_autoglo, 'ersem_z6n', tracer_bkg(:,:,:,jp_fabm_z6n) )
	680	CALL iom_get( inum, jpdom_autoglo, 'ersem_z6p', tracer_bkg(:,:,:,jp_fabm_z6p) )
	681	CALL iom_get( inum, jpdom_autoglo, 'ersem_n1p', tracer_bkg(:,:,:,jp_fabm_n1p) )
	682	CALL iom_get( inum, jpdom_autoglo, 'ersem_n3n', tracer_bkg(:,:,:,jp_fabm_n3n) )
	683	CALL iom_get( inum, jpdom_autoglo, 'ersem_n4n', tracer_bkg(:,:,:,jp_fabm_n4n) )
	684	CALL iom_get( inum, jpdom_autoglo, 'ersem_n5s', tracer_bkg(:,:,:,jp_fabm_n5s) )
	685	CALL iom_get( inum, jpdom_autoglo, 'ersem_o2o', tracer_bkg(:,:,:,jp_fabm_o2o) )
	686	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c', tracer_bkg(:,:,:,jp_fabm_o3c) )
	687	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) )
	688	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:) )
	689	totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
[9296]	690	#endif
[9862]	691	ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
[9296]	692	#if defined key_hadocc
	693	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
	694	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
[10141]	695	#elif defined key_medusa
[9296]	696	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
	697	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
[10622]	698	#elif defined key_fabm
	699	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c', tracer_bkg(:,:,:,jp_fabm_o3c) )
	700	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) )
	701	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:) )
	702	totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
[9296]	703	#endif
	704	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
	705	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
	706	ENDIF
	707
	708	CALL iom_close( inum )
	709
	710	DO jt = 1, jptra
	711	tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
	712	END DO
	713
[9322]	714	ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
[9296]	715
	716	ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
	717	ALLOCATE( mld_max_bkg(jpi,jpj) )
	718	tracer_bkg(:,:,:,:) = 0.0
	719	mld_max_bkg(:,:) = 0.0
	720
	721	CALL iom_open( c_asmbkg, inum )
	722
	723	#if defined key_hadocc
	724	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
	725	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
[10141]	726	#elif defined key_medusa
[9296]	727	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
	728	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
[10622]	729	#elif defined key_fabm
	730	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c', tracer_bkg(:,:,:,jp_fabm_o3c) )
	731	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) )
	732	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:) )
	733	totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
[9296]	734	#endif
	735	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
	736
	737	CALL iom_close( inum )
	738
	739	DO jt = 1, jptra
	740	tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
	741	END DO
	742	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
	743
	744	ENDIF
[9322]	745	#else
	746	CALL ctl_stop( ' asm_bgc_init_bkg: key_top and a compatible biogeochemical model required' )
[9296]	747	#endif
	748
	749	END SUBROUTINE asm_bgc_init_bkg
	750
[9297]	751	!!===========================================================================
	752	!!===========================================================================
	753	!!===========================================================================
[9296]	754
	755	SUBROUTINE asm_bgc_bal_wri( kt )
[9297]	756	!!------------------------------------------------------------------------
[9296]	757	!!
	758	!! * ROUTINE asm_bgc_bal_wri *
	759	!!
	760	!! ** Purpose : Write to file the balancing increments
	761	!! calculated for biogeochemistry
	762	!!
	763	!! ** Method : Write to file the balancing increments
	764	!! calculated for biogeochemistry
	765	!!
	766	!! ** Action :
	767	!!
	768	!! References :
	769	!!
	770	!! History :
	771	!! ! 2014-08 (D. Ford) Initial version, based on asm_bkg_wri
[9297]	772	!!------------------------------------------------------------------------
[9296]	773	!! * Arguments
	774	INTEGER, INTENT( IN ) :: kt ! Current time-step
	775	!! * Local declarations
	776	CHARACTER(LEN=50) :: cl_asmbal ! Filename (with extension)
	777	LOGICAL :: llok ! Check if file exists
	778	INTEGER :: inum ! File unit number
	779	REAL(wp) :: zdate ! Date
[9297]	780	!!------------------------------------------------------------------------
[9296]	781
	782	! Set things up
	783	zdate = REAL( ndastp )
	784	WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal )
	785
	786	! Check if file exists
	787	INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok )
	788	IF( .NOT. llok ) THEN
	789	IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// &
	790	& TRIM( c_asmbal ) // ' at timestep = ', kt
	791
	792	! Define the output file
	793	CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
	794
	795	! Write the information
	796	CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate )
	797
[9381]	798	IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
[10622]	799	& ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
[9381]	800	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
	801	& ln_slphynoninc ) THEN
[10141]	802	#if defined key_medusa
[9381]	803	CALL iom_rstput( kt, kt, inum, 'phy2d_chn', phyto2d_balinc(:,:,:,jpchn) )
	804	CALL iom_rstput( kt, kt, inum, 'phy2d_chd', phyto2d_balinc(:,:,:,jpchd) )
	805	CALL iom_rstput( kt, kt, inum, 'phy2d_phn', phyto2d_balinc(:,:,:,jpphn) )
	806	CALL iom_rstput( kt, kt, inum, 'phy2d_phd', phyto2d_balinc(:,:,:,jpphd) )
	807	CALL iom_rstput( kt, kt, inum, 'phy2d_pds', phyto2d_balinc(:,:,:,jppds) )
[9322]	808	IF ( ln_phytobal ) THEN
[9381]	809	CALL iom_rstput( kt, kt, inum, 'phy2d_zmi', phyto2d_balinc(:,:,:,jpzmi) )
	810	CALL iom_rstput( kt, kt, inum, 'phy2d_zme', phyto2d_balinc(:,:,:,jpzme) )
	811	CALL iom_rstput( kt, kt, inum, 'phy2d_din', phyto2d_balinc(:,:,:,jpdin) )
	812	CALL iom_rstput( kt, kt, inum, 'phy2d_sil', phyto2d_balinc(:,:,:,jpsil) )
	813	CALL iom_rstput( kt, kt, inum, 'phy2d_fer', phyto2d_balinc(:,:,:,jpfer) )
	814	CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jpdet) )
	815	CALL iom_rstput( kt, kt, inum, 'phy2d_dtc', phyto2d_balinc(:,:,:,jpdtc) )
	816	CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jpdic) )
	817	CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jpalk) )
	818	CALL iom_rstput( kt, kt, inum, 'phy2d_oxy', phyto2d_balinc(:,:,:,jpoxy) )
[9296]	819	ENDIF
	820	#elif defined key_hadocc
[9381]	821	CALL iom_rstput( kt, kt, inum, 'phy2d_phy', phyto2d_balinc(:,:,:,jp_had_phy) )
[9322]	822	IF ( ln_phytobal ) THEN
[9381]	823	CALL iom_rstput( kt, kt, inum, 'phy2d_nut', phyto2d_balinc(:,:,:,jp_had_nut) )
	824	CALL iom_rstput( kt, kt, inum, 'phy2d_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) )
	825	CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jp_had_pdn) )
	826	CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jp_had_dic) )
	827	CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) )
[9296]	828	ENDIF
[10622]	829	#elif defined key_fabm
	830	CALL iom_rstput( kt, kt, inum, 'phy2d_chl1', phyto2d_balinc(:,:,:,jp_fabm_chl1) )
	831	CALL iom_rstput( kt, kt, inum, 'phy2d_chl2', phyto2d_balinc(:,:,:,jp_fabm_chl2) )
	832	CALL iom_rstput( kt, kt, inum, 'phy2d_chl3', phyto2d_balinc(:,:,:,jp_fabm_chl3) )
	833	CALL iom_rstput( kt, kt, inum, 'phy2d_chl4', phyto2d_balinc(:,:,:,jp_fabm_chl4) )
	834	CALL iom_rstput( kt, kt, inum, 'phy2d_p1c', phyto2d_balinc(:,:,:,jp_fabm_p1c) )
	835	CALL iom_rstput( kt, kt, inum, 'phy2d_p1n', phyto2d_balinc(:,:,:,jp_fabm_p1n) )
	836	CALL iom_rstput( kt, kt, inum, 'phy2d_p1p', phyto2d_balinc(:,:,:,jp_fabm_p1p) )
	837	CALL iom_rstput( kt, kt, inum, 'phy2d_p1s', phyto2d_balinc(:,:,:,jp_fabm_p1s) )
	838	CALL iom_rstput( kt, kt, inum, 'phy2d_p2c', phyto2d_balinc(:,:,:,jp_fabm_p2c) )
	839	CALL iom_rstput( kt, kt, inum, 'phy2d_p2n', phyto2d_balinc(:,:,:,jp_fabm_p2n) )
	840	CALL iom_rstput( kt, kt, inum, 'phy2d_p2p', phyto2d_balinc(:,:,:,jp_fabm_p2p) )
	841	CALL iom_rstput( kt, kt, inum, 'phy2d_p3c', phyto2d_balinc(:,:,:,jp_fabm_p3c) )
	842	CALL iom_rstput( kt, kt, inum, 'phy2d_p3n', phyto2d_balinc(:,:,:,jp_fabm_p3n) )
	843	CALL iom_rstput( kt, kt, inum, 'phy2d_p3p', phyto2d_balinc(:,:,:,jp_fabm_p3p) )
	844	CALL iom_rstput( kt, kt, inum, 'phy2d_p4c', phyto2d_balinc(:,:,:,jp_fabm_p4c) )
	845	CALL iom_rstput( kt, kt, inum, 'phy2d_p4n', phyto2d_balinc(:,:,:,jp_fabm_p4n) )
	846	CALL iom_rstput( kt, kt, inum, 'phy2d_p4p', phyto2d_balinc(:,:,:,jp_fabm_p4p) )
	847	IF ( ln_phytobal ) THEN
	848	CALL iom_rstput( kt, kt, inum, 'phy2d_z4c', phyto2d_balinc(:,:,:,jp_fabm_z4c) )
	849	CALL iom_rstput( kt, kt, inum, 'phy2d_z5c', phyto2d_balinc(:,:,:,jp_fabm_z5c) )
	850	CALL iom_rstput( kt, kt, inum, 'phy2d_z5n', phyto2d_balinc(:,:,:,jp_fabm_z5n) )
	851	CALL iom_rstput( kt, kt, inum, 'phy2d_z5p', phyto2d_balinc(:,:,:,jp_fabm_z5p) )
	852	CALL iom_rstput( kt, kt, inum, 'phy2d_z6c', phyto2d_balinc(:,:,:,jp_fabm_z6c) )
	853	CALL iom_rstput( kt, kt, inum, 'phy2d_z6n', phyto2d_balinc(:,:,:,jp_fabm_z6n) )
	854	CALL iom_rstput( kt, kt, inum, 'phy2d_z6p', phyto2d_balinc(:,:,:,jp_fabm_z6p) )
	855	CALL iom_rstput( kt, kt, inum, 'phy2d_n1p', phyto2d_balinc(:,:,:,jp_fabm_n1p) )
	856	CALL iom_rstput( kt, kt, inum, 'phy2d_n3n', phyto2d_balinc(:,:,:,jp_fabm_n3n) )
	857	CALL iom_rstput( kt, kt, inum, 'phy2d_n4n', phyto2d_balinc(:,:,:,jp_fabm_n4n) )
	858	CALL iom_rstput( kt, kt, inum, 'phy2d_n5s', phyto2d_balinc(:,:,:,jp_fabm_n5s) )
	859	CALL iom_rstput( kt, kt, inum, 'phy2d_o2o', phyto2d_balinc(:,:,:,jp_fabm_o2o) )
	860	CALL iom_rstput( kt, kt, inum, 'phy2d_o3c', phyto2d_balinc(:,:,:,jp_fabm_o3c) )
	861	CALL iom_rstput( kt, kt, inum, 'phy2d_o3ba', phyto2d_balinc(:,:,:,jp_fabm_o3ba) )
	862	ENDIF
[9296]	863	#endif
	864	ENDIF
	865
[9381]	866	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
[10141]	867	#if defined key_medusa
[9381]	868	CALL iom_rstput( kt, kt, inum, 'phy3d_chn', phyto3d_balinc(:,:,:,jpchn) )
	869	CALL iom_rstput( kt, kt, inum, 'phy3d_chd', phyto3d_balinc(:,:,:,jpchd) )
	870	CALL iom_rstput( kt, kt, inum, 'phy3d_phn', phyto3d_balinc(:,:,:,jpphn) )
	871	CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) )
	872	CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) )
	873	#elif defined key_hadocc
	874	CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) )
[10622]	875	#elif defined key_fabm
	876	CALL iom_rstput( kt, kt, inum, 'phy3d_chl1', phyto3d_balinc(:,:,:,jp_fabm_chl1) )
	877	CALL iom_rstput( kt, kt, inum, 'phy3d_chl2', phyto3d_balinc(:,:,:,jp_fabm_chl2) )
	878	CALL iom_rstput( kt, kt, inum, 'phy3d_chl3', phyto3d_balinc(:,:,:,jp_fabm_chl3) )
	879	CALL iom_rstput( kt, kt, inum, 'phy3d_chl4', phyto3d_balinc(:,:,:,jp_fabm_chl4) )
	880	CALL iom_rstput( kt, kt, inum, 'phy3d_p1c', phyto3d_balinc(:,:,:,jp_fabm_p1c) )
	881	CALL iom_rstput( kt, kt, inum, 'phy3d_p1n', phyto3d_balinc(:,:,:,jp_fabm_p1n) )
	882	CALL iom_rstput( kt, kt, inum, 'phy3d_p1p', phyto3d_balinc(:,:,:,jp_fabm_p1p) )
	883	CALL iom_rstput( kt, kt, inum, 'phy3d_p1s', phyto3d_balinc(:,:,:,jp_fabm_p1s) )
	884	CALL iom_rstput( kt, kt, inum, 'phy3d_p2c', phyto3d_balinc(:,:,:,jp_fabm_p2c) )
	885	CALL iom_rstput( kt, kt, inum, 'phy3d_p2n', phyto3d_balinc(:,:,:,jp_fabm_p2n) )
	886	CALL iom_rstput( kt, kt, inum, 'phy3d_p2p', phyto3d_balinc(:,:,:,jp_fabm_p2p) )
	887	CALL iom_rstput( kt, kt, inum, 'phy3d_p3c', phyto3d_balinc(:,:,:,jp_fabm_p3c) )
	888	CALL iom_rstput( kt, kt, inum, 'phy3d_p3n', phyto3d_balinc(:,:,:,jp_fabm_p3n) )
	889	CALL iom_rstput( kt, kt, inum, 'phy3d_p3p', phyto3d_balinc(:,:,:,jp_fabm_p3p) )
	890	CALL iom_rstput( kt, kt, inum, 'phy3d_p4c', phyto3d_balinc(:,:,:,jp_fabm_p4c) )
	891	CALL iom_rstput( kt, kt, inum, 'phy3d_p4n', phyto3d_balinc(:,:,:,jp_fabm_p4n) )
	892	CALL iom_rstput( kt, kt, inum, 'phy3d_p4p', phyto3d_balinc(:,:,:,jp_fabm_p4p) )
[9381]	893	#endif
	894	ENDIF
	895
[9296]	896	IF ( ln_spco2inc ) THEN
[10141]	897	#if defined key_medusa
[9381]	898	CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jpdic) )
	899	CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jpalk) )
[9296]	900	#elif defined key_hadocc
[9381]	901	CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
	902	CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
[10622]	903	#elif defined key_fabm
	904	CALL iom_rstput( kt, kt, inum, 'pco2_o3c', pco2_balinc(:,:,:,jp_fabm_o3c) )
	905	CALL iom_rstput( kt, kt, inum, 'pco2_o3ba', pco2_balinc(:,:,:,jp_fabm_o3ba) )
[9296]	906	#endif
	907	ELSE IF ( ln_sfco2inc ) THEN
[10141]	908	#if defined key_medusa
[9381]	909	CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jpdic) )
	910	CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jpalk) )
[9296]	911	#elif defined key_hadocc
[9381]	912	CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
	913	CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
[10622]	914	#elif defined key_fabm
	915	CALL iom_rstput( kt, kt, inum, 'fco2_o3c', pco2_balinc(:,:,:,jp_fabm_o3c) )
	916	CALL iom_rstput( kt, kt, inum, 'fco2_o3ba', pco2_balinc(:,:,:,jp_fabm_o3ba) )
[9296]	917	#endif
	918	ENDIF
	919
[9381]	920	IF ( ln_pphinc ) THEN
[10141]	921	#if defined key_medusa
[9381]	922	CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jpdic) )
	923	CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jpalk) )
	924	#elif defined key_hadocc
	925	CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) )
	926	CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) )
[10622]	927	#elif defined key_fabm
	928	CALL iom_rstput( kt, kt, inum, 'ph_o3c', ph_balinc(:,:,:,jp_fabm_o3c) )
	929	CALL iom_rstput( kt, kt, inum, 'ph_o3ba', ph_balinc(:,:,:,jp_fabm_o3ba) )
[9381]	930	#endif
	931	ENDIF
	932
[9296]	933	CALL iom_close( inum )
	934	ELSE
	935	CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', &
	936	& ' Therefore not writing out balancing increments at this timestep', &
	937	& ' Something has probably gone wrong somewhere' )
	938	IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// &
	939	& TRIM( c_asmbal ) // ' at timestep = ', kt
	940	ENDIF
	941
	942	END SUBROUTINE asm_bgc_bal_wri
	943
[9297]	944	!!===========================================================================
	945	!!===========================================================================
	946	!!===========================================================================
[9296]	947
[10622]	948	SUBROUTINE asm_bgc_bkg_alloc
[9297]	949	!!------------------------------------------------------------------------
[10622]	950	!! * ROUTINE asm_bgc_bkg_alloc *
	951	!!
	952	!! ** Purpose : allocate time-average arrays for background
	953	!!
	954	!! ** Method : allocate time-average arrays for background
	955	!!
	956	!! ** Action : allocate time-average arrays for background
	957	!!
	958	!! References : asm_bkg_wri
	959	!!------------------------------------------------------------------------
	960	!!
	961	INTEGER :: ierror
	962	!!
	963	!!------------------------------------------------------------------------
	964
	965	ALLOCATE( pgrow_avg_tavg(jpi,jpj), STAT=ierror )
	966	IF( ierror > 0 ) THEN
	967	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate pgrow_avg_tavg' )
	968	ENDIF
	969	pgrow_avg_tavg(:,:) = 0.0
	970
	971	ALLOCATE( ploss_avg_tavg(jpi,jpj), STAT=ierror )
	972	IF( ierror > 0 ) THEN
	973	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate ploss_avg_tavg' )
	974	ENDIF
	975	ploss_avg_tavg(:,:) = 0.0
	976
	977	ALLOCATE( phyt_avg_tavg(jpi,jpj), STAT=ierror )
	978	IF( ierror > 0 ) THEN
	979	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate phyt_avg_tavg' )
	980	ENDIF
	981	phyt_avg_tavg(:,:) = 0.0
	982
	983	ALLOCATE( trn_tavg(jpi,jpj,jpk,jptra), STAT=ierror )
	984	IF( ierror > 0 ) THEN
	985	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate trn_tavg' )
	986	ENDIF
	987	trn_tavg(:,:,:,:) = 0.0
	988
	989	#if defined key_hadocc
	990	ALLOCATE( HADOCC_CHL_tavg(jpi,jpj,jpk), STAT=ierror )
	991	IF( ierror > 0 ) THEN
	992	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate HADOCC_CHL_tavg' )
	993	ENDIF
	994	HADOCC_CHL_tavg(:,:,:) = 0.0
	995
	996	ALLOCATE( cchl_p_tavg(jpi,jpj,jpk), STAT=ierror )
	997	IF( ierror > 0 ) THEN
	998	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate cchl_p_tavg' )
	999	ENDIF
	1000	cchl_p_tavg(:,:,:) = 0.0
	1001
	1002	#elif defined key_fabm
	1003	ALLOCATE( totalk_tavg(jpi,jpj,jpk), STAT=ierror )
	1004	IF( ierror > 0 ) THEN
	1005	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate totalk_tavg' )
	1006	ENDIF
	1007	totalk_tavg(:,:,:) = 0.0
	1008	#endif
	1009
	1010	END SUBROUTINE asm_bgc_bkg_alloc
	1011
	1012	!!===========================================================================
	1013	!!===========================================================================
	1014	!!===========================================================================
	1015
	1016	SUBROUTINE asm_bgc_bkg_tavg(kt, pnumtimes_tavg)
	1017	!!------------------------------------------------------------------------
[9296]	1018	!! * ROUTINE asm_bgc_bkg_wri *
	1019	!!
	1020	!! ** Purpose : write out bgc background
	1021	!!
	1022	!! ** Method : write out bgc background
	1023	!!
	1024	!! ** Action : write out bgc background
	1025	!!
	1026	!! References : asm_bkg_wri
[9297]	1027	!!------------------------------------------------------------------------
[9296]	1028	!!
[10622]	1029	INTEGER, INTENT(in ) :: kt ! Current time-step
[9296]	1030	!!
[10622]	1031	REAL(wp), INTENT(in ) :: pnumtimes_tavg ! No of times to average over
	1032	!!
[9297]	1033	!!------------------------------------------------------------------------
[10622]	1034
	1035	IF (kt == nittrc000) THEN
	1036	pgrow_avg_tavg(:,:) = 0.0
	1037	ploss_avg_tavg(:,:) = 0.0
	1038	phyt_avg_tavg(:,:) = 0.0
	1039	trn_tavg(:,:,:,:) = 0.0
	1040	#if defined key_hadocc
	1041	HADOCC_CHL_tavg(:,:,:) = 0.0
	1042	cchl_p_tavg(:,:,:) = 0.0
	1043	#elif defined key_fabm
	1044	totalk_tavg(:,:,:) = 0.0
	1045	#endif
	1046	ENDIF
	1047
	1048	pgrow_avg_tavg(:,:) = pgrow_avg_tavg(:,:) + pgrow_avg(:,:) / pnumtimes_tavg
	1049	ploss_avg_tavg(:,:) = ploss_avg_tavg(:,:) + ploss_avg(:,:) / pnumtimes_tavg
	1050	phyt_avg_tavg(:,:) = phyt_avg_tavg(:,:) + phyt_avg(:,:) / pnumtimes_tavg
	1051	trn_tavg(:,:,:,:) = trn_tavg(:,:,:,:) + trn(:,:,:,:) / pnumtimes_tavg
	1052	#if defined key_hadocc
	1053	HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL_tavg(:,:,:) + HADOCC_CHL(:,:,:) / pnumtimes_tavg
	1054	cchl_p_tavg(:,:,:) = cchl_p_tavg(:,:,:) + cchl_p(:,:,:) / pnumtimes_tavg
	1055	#elif defined key_fabm
	1056	totalk_tavg(:,:,:) = totalk_tavg(:,:,:) + &
	1057	& fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) / pnumtimes_tavg
	1058	totalk_tavg(:,:,:) = totalk_tavg(:,:,:) * tmask(:,:,:)
	1059	#endif
	1060
	1061	END SUBROUTINE asm_bgc_bkg_tavg
[9296]	1062
[10622]	1063	!!===========================================================================
	1064	!!===========================================================================
	1065	!!===========================================================================
	1066
	1067	SUBROUTINE asm_bgc_bkg_wri( kt, knum, ld_avgbkg )
	1068	!!------------------------------------------------------------------------
	1069	!! * ROUTINE asm_bgc_bkg_wri *
	1070	!!
	1071	!! ** Purpose : write out bgc background
	1072	!!
	1073	!! ** Method : write out bgc background
	1074	!!
	1075	!! ** Action : write out bgc background
	1076	!!
	1077	!! References : asm_bkg_wri
	1078	!!------------------------------------------------------------------------
	1079	!!
	1080	INTEGER, INTENT(in ) :: kt ! Current time-step
	1081	INTEGER, INTENT(in ) :: knum ! i/o unit of increments file
	1082	LOGICAL, INTENT(in ) :: ld_avgbkg ! Averaged background?
	1083	!!
	1084	INTEGER :: nitbgcbkg_r ! Period referenced to nit000
	1085	!!
	1086	!!------------------------------------------------------------------------
	1087
	1088	IF (ld_avgbkg) THEN
	1089	nitbgcbkg_r = nitavgbkg_r
	1090	ELSE
	1091	nitbgcbkg_r = nitbkg_r
	1092	pgrow_avg_tavg(:,:) = pgrow_avg(:,:)
	1093	ploss_avg_tavg(:,:) = ploss_avg(:,:)
	1094	phyt_avg_tavg(:,:) = phyt_avg(:,:)
	1095	trn_tavg(:,:,:,:) = trn(:,:,:,:)
[9296]	1096	#if defined key_hadocc
[10622]	1097	HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL(:,:,:)
	1098	cchl_p_tavg(:,:,:) = cchl_p(:,:,:)
	1099	#elif defined key_fabm
	1100	totalk_tavg(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta)
	1101	totalk_tavg(:,:,:) = totalk_tavg(:,:,:) * tmask(:,:,:)
	1102	#endif
	1103	ENDIF
	1104
	1105	#if defined key_hadocc
	1106	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg' , pgrow_avg_tavg )
	1107	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg' , ploss_avg_tavg )
	1108	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg' , phyt_avg_tavg )
	1109	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max' , mld_max )
	1110	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_nut' , trn_tavg(:,:,:,jp_had_nut) )
	1111	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_phy' , trn_tavg(:,:,:,jp_had_phy) )
	1112	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_zoo' , trn_tavg(:,:,:,jp_had_zoo) )
	1113	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_pdn' , trn_tavg(:,:,:,jp_had_pdn) )
	1114	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_dic' , trn_tavg(:,:,:,jp_had_dic) )
	1115	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_alk' , trn_tavg(:,:,:,jp_had_alk) )
	1116	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_chl' , HADOCC_CHL(:,:,:) )
	1117	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:) )
[10141]	1118	#elif defined key_medusa
[10622]	1119	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg' , pgrow_avg_tavg )
	1120	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg' , ploss_avg_tavg )
	1121	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg' , phyt_avg_tavg )
	1122	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max' , mld_max )
	1123	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chn' , trn_tavg(:,:,:,jpchn) )
	1124	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chd' , trn_tavg(:,:,:,jpchd) )
	1125	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phn' , trn_tavg(:,:,:,jpphn) )
	1126	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phd' , trn_tavg(:,:,:,jpphd) )
	1127	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_pds' , trn_tavg(:,:,:,jppds) )
	1128	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zmi' , trn_tavg(:,:,:,jpzmi) )
	1129	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zme' , trn_tavg(:,:,:,jpzme) )
	1130	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_din' , trn_tavg(:,:,:,jpdin) )
	1131	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_sil' , trn_tavg(:,:,:,jpsil) )
	1132	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_fer' , trn_tavg(:,:,:,jpfer) )
	1133	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_det' , trn_tavg(:,:,:,jpdet) )
	1134	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dtc' , trn_tavg(:,:,:,jpdtc) )
	1135	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dic' , trn_tavg(:,:,:,jpdic) )
	1136	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_alk' , trn_tavg(:,:,:,jpalk) )
	1137	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_oxy' , trn_tavg(:,:,:,jpoxy) )
	1138	#elif defined key_fabm
	1139	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg' , pgrow_avg_tavg )
	1140	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg' , ploss_avg_tavg )
	1141	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg' , phyt_avg_tavg )
	1142	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max' , mld_max )
	1143	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl1' , trn_tavg(:,:,:,jp_fabm_chl1) )
	1144	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl2' , trn_tavg(:,:,:,jp_fabm_chl2) )
	1145	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl3' , trn_tavg(:,:,:,jp_fabm_chl3) )
	1146	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl4' , trn_tavg(:,:,:,jp_fabm_chl4) )
	1147	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1c' , trn_tavg(:,:,:,jp_fabm_p1c) )
	1148	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1n' , trn_tavg(:,:,:,jp_fabm_p1n) )
	1149	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1p' , trn_tavg(:,:,:,jp_fabm_p1p) )
	1150	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1s' , trn_tavg(:,:,:,jp_fabm_p1s) )
	1151	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2c' , trn_tavg(:,:,:,jp_fabm_p2c) )
	1152	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2n' , trn_tavg(:,:,:,jp_fabm_p2n) )
	1153	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2p' , trn_tavg(:,:,:,jp_fabm_p2p) )
	1154	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3c' , trn_tavg(:,:,:,jp_fabm_p3c) )
	1155	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3n' , trn_tavg(:,:,:,jp_fabm_p3n) )
	1156	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3p' , trn_tavg(:,:,:,jp_fabm_p3p) )
	1157	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4c' , trn_tavg(:,:,:,jp_fabm_p4c) )
	1158	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4n' , trn_tavg(:,:,:,jp_fabm_p4n) )
	1159	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4p' , trn_tavg(:,:,:,jp_fabm_p4p) )
	1160	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z4c' , trn_tavg(:,:,:,jp_fabm_z4c) )
	1161	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5c' , trn_tavg(:,:,:,jp_fabm_z5c) )
	1162	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5n' , trn_tavg(:,:,:,jp_fabm_z5n) )
	1163	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5p' , trn_tavg(:,:,:,jp_fabm_z5p) )
	1164	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6c' , trn_tavg(:,:,:,jp_fabm_z6c) )
	1165	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6n' , trn_tavg(:,:,:,jp_fabm_z6n) )
	1166	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6p' , trn_tavg(:,:,:,jp_fabm_z6p) )
	1167	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n1p' , trn_tavg(:,:,:,jp_fabm_n1p) )
	1168	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n3n' , trn_tavg(:,:,:,jp_fabm_n3n) )
	1169	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n4n' , trn_tavg(:,:,:,jp_fabm_n4n) )
	1170	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n5s' , trn_tavg(:,:,:,jp_fabm_n5s) )
	1171	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o2o' , trn_tavg(:,:,:,jp_fabm_o2o) )
	1172	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3c' , trn_tavg(:,:,:,jp_fabm_o3c) )
	1173	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ba' , trn_tavg(:,:,:,jp_fabm_o3ba) )
	1174	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ta' , totalk_tavg )
[9296]	1175	#endif
	1176
	1177	END SUBROUTINE asm_bgc_bkg_wri
	1178
[9297]	1179	!!===========================================================================
	1180	!!===========================================================================
	1181	!!===========================================================================
[9296]	1182
[9431]	1183	SUBROUTINE asm_bgc_unlog_2d( pbkg, pinc_log, pinc_nonlog )
	1184	!!------------------------------------------------------------------------
	1185	!! * ROUTINE asm_bgc_init_incs *
	1186	!!
	1187	!! ** Purpose : convert log increments to non-log
	1188	!!
	1189	!! ** Method : need to account for model background,
	1190	!! cannot simply do 10^log_inc. Need to:
	1191	!! 1) Add log_inc to log10(background) to get log10(analysis)
	1192	!! 2) Take 10^log10(analysis) to get analysis
	1193	!! 3) Subtract background from analysis to get non-log incs
	1194	!!
	1195	!! ** Action : return non-log increments
	1196	!!
	1197	!! References :
	1198	!!------------------------------------------------------------------------
	1199	!!
	1200	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pbkg ! Background
	1201	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pinc_log ! Log incs
	1202	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: pinc_nonlog ! Non-log incs
	1203	!
	1204	INTEGER :: ji, jj ! Loop counters
	1205	!!
	1206	!!------------------------------------------------------------------------
	1207
	1208	DO jj = 1, jpj
	1209	DO ji = 1, jpi
	1210	IF ( pbkg(ji,jj) > 0.0 ) THEN
	1211	pinc_nonlog(ji,jj) = 10**( LOG10( pbkg(ji,jj) ) + &
	1212	& pinc_log(ji,jj) ) &
	1213	& - pbkg(ji,jj)
	1214	ELSE
	1215	pinc_nonlog(ji,jj) = 0.0
	1216	ENDIF
	1217	END DO
	1218	END DO
	1219
	1220	END SUBROUTINE asm_bgc_unlog_2d
	1221
	1222	!!===========================================================================
	1223	!!===========================================================================
	1224	!!===========================================================================
	1225
[9326]	1226	SUBROUTINE phyto2d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
[9297]	1227	!!------------------------------------------------------------------------
[9326]	1228	!! * ROUTINE phyto2d_asm_inc *
[9296]	1229	!!
	1230	!! ** Purpose : Apply the chlorophyll assimilation increments.
	1231	!!
	1232	!! ** Method : Calculate increments to state variables using nitrogen
	1233	!! balancing.
	1234	!! Direct initialization or Incremental Analysis Updating.
	1235	!!
	1236	!! ** Action :
[9297]	1237	!!------------------------------------------------------------------------
[9296]	1238	INTEGER, INTENT(IN) :: kt ! Current time step
	1239	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
	1240	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
	1241	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
	1242	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
	1243	!
[9431]	1244	INTEGER :: it ! Index
	1245	REAL(wp) :: zincwgt ! IAU weight for current time step
	1246	REAL(wp) :: zincper ! IAU interval in seconds
	1247	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth
	1248	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chltot ! Local chltot incs
	1249	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chltot ! Local chltot bkg
	1250	REAL(wp), DIMENSION(jpi,jpj) :: zinc_phytot ! Local phytot incs
	1251	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phytot ! Local phytot bkg
[10141]	1252	#if defined key_medusa
[9431]	1253	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs
	1254	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
	1255	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnon ! Local chlnon incs
	1256	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnon ! Local chlnon bkg
	1257	REAL(wp), DIMENSION(jpi,jpj) :: zinc_phydia ! Local phydia incs
	1258	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phydia ! Local phydia bkg
	1259	REAL(wp), DIMENSION(jpi,jpj) :: zinc_phynon ! Local phynon incs
	1260	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg
[10622]	1261	#elif defined key_fabm
	1262	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs
	1263	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
	1264	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldin ! Local chldin incs
	1265	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldin ! Local chldin bkg
	1266	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnan ! Local chlnan incs
	1267	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnan ! Local chlnan bkg
	1268	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlpic ! Local chlpic incs
	1269	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlpic ! Local chlpic bkg
[9431]	1270	#endif
[9297]	1271	!!------------------------------------------------------------------------
[9322]	1272
[9296]	1273	IF ( kt <= nit000 ) THEN
	1274
[9431]	1275	! Un-log any log increments for passing to balancing routines
[9435]	1276	! Remember that two sets of non-log increments should not be
	1277	! expected to be in the same ratio as their log equivalents
	1278
[9431]	1279	! Total chlorophyll
	1280	IF ( ln_slchltotinc ) THEN
[10141]	1281	#if defined key_medusa
[9431]	1282	zbkg_chltot(:,:) = tracer_bkg(:,:,1,jpchn) + tracer_bkg(:,:,1,jpchd)
	1283	#elif defined key_hadocc
	1284	zbkg_chltot(:,:) = chl_bkg(:,:,1)
[10622]	1285	#elif defined key_fabm
	1286	zbkg_chltot(:,:) = tracer_bkg(:,:,1,jp_fabm_chl1) + &
	1287	& tracer_bkg(:,:,1,jp_fabm_chl2) + &
	1288	& tracer_bkg(:,:,1,jp_fabm_chl3) + &
	1289	& tracer_bkg(:,:,1,jp_fabm_chl4)
[9431]	1290	#endif
	1291	CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot )
	1292	ELSE IF ( ln_schltotinc ) THEN
	1293	zinc_chltot(:,:) = schltot_bkginc(:,:)
	1294	ELSE
	1295	zinc_chltot(:,:) = 0.0
	1296	ENDIF
	1297
[10622]	1298	#if defined key_medusa \|\| defined key_fabm
[9431]	1299	! Diatom chlorophyll
	1300	IF ( ln_slchldiainc ) THEN
[10622]	1301	#if defined key_medusa
[9431]	1302	zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd)
[10622]	1303	#elif defined key_fabm
	1304	zbkg_chldia(:,:) = tracer_bkg(:,:,1,jp_fabm_chl1)
	1305	#endif
[9431]	1306	CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia )
	1307	ELSE
	1308	zinc_chldia(:,:) = 0.0
	1309	ENDIF
	1310	#endif
	1311
[10141]	1312	#if defined key_medusa
[9431]	1313	! Non-diatom chlorophyll
	1314	IF ( ln_slchlnoninc ) THEN
	1315	zbkg_chlnon(:,:) = tracer_bkg(:,:,1,jpchn)
	1316	CALL asm_bgc_unlog_2d( zbkg_chlnon, slchlnon_bkginc, zinc_chlnon )
	1317	ELSE
	1318	zinc_chlnon(:,:) = 0.0
	1319	ENDIF
	1320	#endif
	1321
[10622]	1322	#if defined key_fabm
	1323	! Nanophytoplankton chlorophyll
	1324	IF ( ln_slchlnaninc ) THEN
	1325	zbkg_chlnan(:,:) = tracer_bkg(:,:,1,jp_fabm_chl2)
	1326	CALL asm_bgc_unlog_2d( zbkg_chlnan, slchlnan_bkginc, zinc_chlnan )
	1327	ELSE
	1328	zinc_chlnan(:,:) = 0.0
	1329	ENDIF
	1330
	1331	! Picophytoplankton chlorophyll
	1332	IF ( ln_slchlpicinc ) THEN
	1333	zbkg_chlpic(:,:) = tracer_bkg(:,:,1,jp_fabm_chl3)
	1334	CALL asm_bgc_unlog_2d( zbkg_chlpic, slchlpic_bkginc, zinc_chlpic )
	1335	ELSE
	1336	zinc_chlpic(:,:) = 0.0
	1337	ENDIF
	1338
	1339	! Dinoflagellate chlorophyll
	1340	IF ( ln_slchldininc ) THEN
	1341	zbkg_chldin(:,:) = tracer_bkg(:,:,1,jp_fabm_chl4)
	1342	CALL asm_bgc_unlog_2d( zbkg_chldin, slchldin_bkginc, zinc_chldin )
	1343	ELSE
	1344	zinc_chldin(:,:) = 0.0
	1345	ENDIF
	1346	#endif
	1347
[9431]	1348	! Total phytoplankton carbon
	1349	IF ( ln_slphytotinc ) THEN
[10141]	1350	#if defined key_medusa
[9431]	1351	zbkg_phytot(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
	1352	#elif defined key_hadocc
	1353	zbkg_phytot(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
	1354	#endif
	1355	CALL asm_bgc_unlog_2d( zbkg_phytot, slphytot_bkginc, zinc_phytot )
	1356	ELSE
	1357	zinc_phytot(:,:) = 0.0
	1358	ENDIF
	1359
[10141]	1360	#if defined key_medusa
[9431]	1361	! Diatom phytoplankton carbon
	1362	IF ( ln_slphydiainc ) THEN
	1363	zbkg_phydia(:,:) = trn(:,:,1,jpphd) * xthetapd
	1364	CALL asm_bgc_unlog_2d( zbkg_phydia, slphydia_bkginc, zinc_phydia )
	1365	ELSE
	1366	zinc_phydia(:,:) = 0.0
	1367	ENDIF
	1368	#endif
	1369
[10141]	1370	#if defined key_medusa
[9431]	1371	! Non-diatom phytoplankton carbon
	1372	IF ( ln_slphynoninc ) THEN
	1373	zbkg_phynon(:,:) = trn(:,:,1,jpphn) * xthetapn
	1374	CALL asm_bgc_unlog_2d( zbkg_phynon, slphynon_bkginc, zinc_phynon )
	1375	ELSE
	1376	zinc_phynon(:,:) = 0.0
	1377	ENDIF
	1378	#endif
	1379
	1380	! Select mixed layer
	1381	IF ( ll_asmdin ) THEN
[10622]	1382	#if defined key_top && ( defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm )
[9431]	1383	CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', &
	1384	& ' Mixed layer depth taken to be background maximum mld_max_bkg' )
	1385	zmld(:,:) = mld_max_bkg(:,:)
[10057]	1386	#else
	1387	CALL ctl_stop( ' Should not be doing biogeochemical assimilation without key_top' )
	1388	#endif
[9431]	1389	ELSE
	1390	SELECT CASE( mld_choice_bgc )
	1391	CASE ( 1 ) ! Turbocline/mixing depth [W points]
	1392	zmld(:,:) = hmld(:,:)
	1393	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
	1394	zmld(:,:) = hmlp(:,:)
	1395	CASE ( 3 ) ! Kara MLD [Interpolated]
	1396	#if defined key_karaml
	1397	IF ( ln_kara ) THEN
	1398	zmld(:,:) = hmld_kara(:,:)
	1399	ELSE
	1400	CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
	1401	& ' but ln_kara=.false.' )
	1402	ENDIF
	1403	#else
	1404	CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
	1405	& ' but is not defined' )
	1406	#endif
	1407	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
[10622]	1408	zmld(:,:) = hmld_tref(:,:)
[9431]	1409	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
	1410	zmld(:,:) = hmlpt(:,:)
	1411	END SELECT
	1412	ENDIF
	1413
[9296]	1414	zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
	1415
[10141]	1416	#if defined key_medusa
[9432]	1417	CALL asm_phyto2d_bal_medusa( (ln_slchltotinc .OR. ln_schltotinc), &
	1418	& zinc_chltot, &
	1419	& ln_slchldiainc, &
	1420	& zinc_chldia, &
	1421	& ln_slchlnoninc, &
	1422	& zinc_chlnon, &
	1423	& ln_slphytotinc, &
	1424	& zinc_phytot, &
	1425	& ln_slphydiainc, &
	1426	& zinc_phydia, &
	1427	& ln_slphynoninc, &
	1428	& zinc_phynon, &
	1429	& zincper, &
	1430	& rn_maxchlinc, ln_phytobal, zmld, &
	1431	& pgrow_avg_bkg, ploss_avg_bkg, &
	1432	& phyt_avg_bkg, mld_max_bkg, &
	1433	& tracer_bkg, phyto2d_balinc )
[9296]	1434	#elif defined key_hadocc
[9432]	1435	CALL asm_phyto2d_bal_hadocc( (ln_slchltotinc .OR. ln_schltotinc), &
	1436	& zinc_chltot, &
	1437	& ln_slphytotinc, &
	1438	& zinc_phytot, &
	1439	& zincper, &
	1440	& rn_maxchlinc, ln_phytobal, zmld, &
	1441	& pgrow_avg_bkg, ploss_avg_bkg, &
	1442	& phyt_avg_bkg, mld_max_bkg, &
	1443	& cchl_p_bkg(:,:,1), &
	1444	& tracer_bkg, phyto2d_balinc )
[10622]	1445	#elif defined key_fabm
	1446	CALL asm_phyto2d_bal_ersem( (ln_slchltotinc .OR. ln_schltotinc), &
	1447	& zinc_chltot, &
	1448	& ln_slchldiainc, &
	1449	& zinc_chldia, &
	1450	& ln_slchlnaninc, &
	1451	& zinc_chlnan, &
	1452	& ln_slchlpicinc, &
	1453	& zinc_chlpic, &
	1454	& ln_slchldininc, &
	1455	& zinc_chldin, &
	1456	& zincper, &
	1457	& rn_maxchlinc, ln_phytobal, zmld, &
	1458	& pgrow_avg_bkg, ploss_avg_bkg, &
	1459	& phyt_avg_bkg, mld_max_bkg, &
	1460	& totalk_bkg, &
	1461	& tracer_bkg, phyto2d_balinc )
[9296]	1462	#else
[9431]	1463	CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', &
[9296]	1464	& 'but not defined a biogeochemical model' )
	1465	#endif
	1466
	1467	ENDIF
	1468
	1469	IF ( ll_asmiau ) THEN
	1470
	1471	!--------------------------------------------------------------------
	1472	! Incremental Analysis Updating
	1473	!--------------------------------------------------------------------
	1474
	1475	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
	1476
	1477	it = kt - nit000 + 1
	1478	zincwgt = pwgtiau(it) ! IAU weight for the current time step
	1479	! note this is not a tendency so should not be divided by rdt
	1480
	1481	IF(lwp) THEN
	1482	WRITE(numout,*)
[9431]	1483	WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', &
[9296]	1484	& kt,' with IAU weight = ', pwgtiau(it)
	1485	WRITE(numout,*) '~~~~~~~~~~~~'
	1486	ENDIF
	1487
	1488	! Update the biogeochemical variables
	1489	! Add directly to trn and trb, rather than to tra, because tra gets
	1490	! reset to zero at the start of trc_stp, called after this routine
[10141]	1491	#if defined key_medusa
[9382]	1492	WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1493	& trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
[9333]	1494	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1495	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1496	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
	1497	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1498	END WHERE
[9296]	1499	#elif defined key_hadocc
[9382]	1500	WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	1501	& trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
[9333]	1502	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	1503	& phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	1504	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
	1505	& phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	1506	END WHERE
[10622]	1507	#elif defined key_fabm
	1508	WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
	1509	& trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
	1510	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
	1511	& phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
	1512	trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
	1513	& phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
	1514	END WHERE
[9296]	1515	#endif
	1516
	1517	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1518	! which is called at end of model run
	1519	ENDIF
	1520
	1521	ELSEIF ( ll_asmdin ) THEN
	1522
	1523	!--------------------------------------------------------------------
	1524	! Direct Initialization
	1525	!--------------------------------------------------------------------
	1526
	1527	IF ( kt == nitdin_r ) THEN
	1528
	1529	neuler = 0 ! Force Euler forward step
	1530
	1531	! Initialize the now fields with the background + increment
	1532	! Background currently is what the model is initialised with
[9333]	1533	CALL ctl_warn( ' Doing direct initialisation with phyto2d assimilation', &
[9296]	1534	& ' Background state is taken from model rather than background file' )
[10141]	1535	#if defined key_medusa
[9382]	1536	WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1537	& trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
[9333]	1538	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1539	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1)
	1540	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
	1541	END WHERE
[9296]	1542	#elif defined key_hadocc
[9382]	1543	WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	1544	& trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
[9333]	1545	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	1546	& phyto2d_balinc(:,:,:,jp_had0:jp_had1)
	1547	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
	1548	END WHERE
[10622]	1549	#elif defined key_fabm
	1550	WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
	1551	& trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
	1552	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
	1553	& phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1)
	1554	trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
	1555	END WHERE
[9296]	1556	#endif
	1557
	1558	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1559	! which is called at end of model run
	1560	ENDIF
	1561	!
	1562	ENDIF
	1563	!
[9326]	1564	END SUBROUTINE phyto2d_asm_inc
[9296]	1565
[9297]	1566	!!===========================================================================
	1567	!!===========================================================================
	1568	!!===========================================================================
[9296]	1569
[9326]	1570	SUBROUTINE phyto3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
	1571	!!------------------------------------------------------------------------
	1572	!! * ROUTINE phyto3d_asm_inc *
	1573	!!
	1574	!! ** Purpose : Apply the profile chlorophyll assimilation increments.
	1575	!!
	1576	!! ** Method : Calculate increments to state variables.
	1577	!! Direct initialization or Incremental Analysis Updating.
	1578	!!
	1579	!! ** Action :
	1580	!!------------------------------------------------------------------------
	1581	INTEGER, INTENT(IN) :: kt ! Current time step
	1582	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
	1583	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
	1584	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
	1585	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
	1586	!
[10622]	1587	INTEGER :: ji, jj, jk ! Loop counters
	1588	INTEGER :: it ! Index
	1589	REAL(wp) :: zincwgt ! IAU weight for timestep
	1590	REAL(wp) :: zfrac_chn ! Fraction of jpchn
	1591	REAL(wp) :: zfrac_chd ! Fraction of jpchd
	1592	REAL(wp) :: zfrac_chl1 ! Fraction of jp_fabm_chl1
	1593	REAL(wp) :: zfrac_chl2 ! Fraction of jp_fabm_chl2
	1594	REAL(wp) :: zfrac_chl3 ! Fraction of jp_fabm_chl3
	1595	REAL(wp) :: zfrac_chl4 ! Fraction of jp_fabm_chl4
	1596	REAL(wp) :: zrat_phn_chn ! jpphn:jpchn ratio
	1597	REAL(wp) :: zrat_phd_chd ! jpphd:jpchd ratio
	1598	REAL(wp) :: zrat_pds_chd ! jppds:jpchd ratio
	1599	REAL(wp) :: zrat_p1c_chl1 ! jp_fabm_p1c:jp_fabm_chl1 ratio
	1600	REAL(wp) :: zrat_p1n_chl1 ! jp_fabm_p1n:jp_fabm_chl1 ratio
	1601	REAL(wp) :: zrat_p1p_chl1 ! jp_fabm_p1p:jp_fabm_chl1 ratio
	1602	REAL(wp) :: zrat_p1s_chl1 ! jp_fabm_p1s:jp_fabm_chl1 ratio
	1603	REAL(wp) :: zrat_p2c_chl2 ! jp_fabm_p2c:jp_fabm_chl2 ratio
	1604	REAL(wp) :: zrat_p2n_chl2 ! jp_fabm_p2n:jp_fabm_chl2 ratio
	1605	REAL(wp) :: zrat_p2p_chl2 ! jp_fabm_p2p:jp_fabm_chl2 ratio
	1606	REAL(wp) :: zrat_p3c_chl3 ! jp_fabm_p3c:jp_fabm_chl3 ratio
	1607	REAL(wp) :: zrat_p3n_chl3 ! jp_fabm_p3n:jp_fabm_chl3 ratio
	1608	REAL(wp) :: zrat_p3p_chl3 ! jp_fabm_p3p:jp_fabm_chl3 ratio
	1609	REAL(wp) :: zrat_p4c_chl4 ! jp_fabm_p4c:jp_fabm_chl4 ratio
	1610	REAL(wp) :: zrat_p4n_chl4 ! jp_fabm_p4n:jp_fabm_chl4 ratio
	1611	REAL(wp) :: zrat_p4p_chl4 ! jp_fabm_p4p:jp_fabm_chl4 ratio
	1612	REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc ! Chlorophyll increments
	1613	REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl ! Chlorophyll background
[9326]	1614	!!------------------------------------------------------------------------
	1615
	1616	IF ( kt <= nit000 ) THEN
	1617
	1618	IF ( ln_plchltotinc ) THEN
	1619	! Convert log10(chlorophyll) increment back to a chlorophyll increment
	1620	! In order to transform logchl incs to chl incs, need to account for model
	1621	! background, cannot simply do 10^logchl_bkginc. Need to:
	1622	! 1) Add logchl inc to log10(background) to get log10(analysis)
	1623	! 2) Take 10^log10(analysis) to get analysis
	1624	! 3) Subtract background from analysis to get chl incs
	1625	! If rn_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
[10141]	1626	#if defined key_medusa
[9326]	1627	bkg_chl(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd)
	1628	#elif defined key_hadocc
	1629	bkg_chl(:,:,:) = chl_bkg(:,:,:)
[10622]	1630	#elif defined key_fabm
	1631	bkg_chl(:,:,:) = tracer_bkg(:,:,:,jp_fabm_chl1) + &
	1632	& tracer_bkg(:,:,:,jp_fabm_chl2) + &
	1633	& tracer_bkg(:,:,:,jp_fabm_chl3) + &
	1634	& tracer_bkg(:,:,:,jp_fabm_chl4)
[9326]	1635	#endif
	1636	DO jk = 1, jpk
	1637	DO jj = 1, jpj
	1638	DO ji = 1, jpi
	1639	IF ( bkg_chl(ji,jj,jk) > 0.0 ) THEN
	1640	chl_inc(ji,jj,jk) = 10**( LOG10( bkg_chl(ji,jj,jk) ) + plchltot_bkginc(ji,jj,jk) ) - bkg_chl(ji,jj,jk)
	1641	IF ( rn_maxchlinc > 0.0 ) THEN
	1642	chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( chl_inc(ji,jj,jk), rn_maxchlinc ) )
	1643	ENDIF
	1644	ELSE
	1645	chl_inc(ji,jj,jk) = 0.0
	1646	ENDIF
	1647	END DO
	1648	END DO
	1649	END DO
	1650	ELSE IF ( ln_pchltotinc ) THEN
	1651	DO jk = 1, jpk
	1652	DO jj = 1, jpj
	1653	DO ji = 1, jpi
	1654	IF ( rn_maxchlinc > 0.0 ) THEN
	1655	chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( pchltot_bkginc(ji,jj,jk), rn_maxchlinc ) )
	1656	ELSE
	1657	chl_inc(ji,jj,jk) = pchltot_bkginc(ji,jj,jk)
	1658	ENDIF
	1659	END DO
	1660	END DO
	1661	END DO
	1662	ENDIF
	1663
[10141]	1664	#if defined key_medusa
[9326]	1665	! Loop over each grid point partioning the increments based on existing ratios
	1666	DO jk = 1, jpk
	1667	DO jj = 1, jpj
	1668	DO ji = 1, jpi
	1669	IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN
	1670	zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd))
	1671	zfrac_chd = 1.0 - zfrac_chn
	1672	phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn
	1673	phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd
[9381]	1674	zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn)
	1675	zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd)
	1676	zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd)
	1677	phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn
	1678	phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd
	1679	phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd
[9326]	1680	ENDIF
	1681	END DO
	1682	END DO
	1683	END DO
	1684	#elif defined key_hadocc
	1685	phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:)
[10622]	1686	#elif defined key_fabm
	1687	! Loop over each grid point partioning the increments based on existing ratios
	1688	DO jk = 1, jpk
	1689	DO jj = 1, jpj
	1690	DO ji = 1, jpi
	1691	IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_chl1) > 0.0 ) .AND. &
	1692	& ( tracer_bkg(ji,jj,jk,jp_fabm_chl2) > 0.0 ) .AND. &
	1693	& ( tracer_bkg(ji,jj,jk,jp_fabm_chl3) > 0.0 ) .AND. &
	1694	& ( tracer_bkg(ji,jj,jk,jp_fabm_chl4) > 0.0 ) ) THEN
	1695	zfrac_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_chl1) / bkg_chl(ji,jj,jk)
	1696	zfrac_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_chl2) / bkg_chl(ji,jj,jk)
	1697	zfrac_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_chl3) / bkg_chl(ji,jj,jk)
	1698	zfrac_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_chl4) / bkg_chl(ji,jj,jk)
	1699	phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) = chl_inc(ji,jj,jk) * zfrac_chl1
	1700	phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) = chl_inc(ji,jj,jk) * zfrac_chl2
	1701	phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) = chl_inc(ji,jj,jk) * zfrac_chl3
	1702	phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) = chl_inc(ji,jj,jk) * zfrac_chl4
	1703	zrat_p1c_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1c) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
	1704	zrat_p1n_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1n) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
	1705	zrat_p1p_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1p) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
	1706	zrat_p1s_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1s) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
	1707	zrat_p2c_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2c) / tracer_bkg(ji,jj,jk,jp_fabm_chl2)
	1708	zrat_p2n_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2n) / tracer_bkg(ji,jj,jk,jp_fabm_chl2)
	1709	zrat_p2p_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2p) / tracer_bkg(ji,jj,jk,jp_fabm_chl2)
	1710	zrat_p3c_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3c) / tracer_bkg(ji,jj,jk,jp_fabm_chl3)
	1711	zrat_p3n_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3n) / tracer_bkg(ji,jj,jk,jp_fabm_chl3)
	1712	zrat_p3p_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3p) / tracer_bkg(ji,jj,jk,jp_fabm_chl3)
	1713	zrat_p4c_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4c) / tracer_bkg(ji,jj,jk,jp_fabm_chl4)
	1714	zrat_p4n_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4n) / tracer_bkg(ji,jj,jk,jp_fabm_chl4)
	1715	zrat_p4p_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4p) / tracer_bkg(ji,jj,jk,jp_fabm_chl4)
	1716	phyto3d_balinc(ji,jj,jk,jp_fabm_p1c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1c_chl1
	1717	phyto3d_balinc(ji,jj,jk,jp_fabm_p1n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1n_chl1
	1718	phyto3d_balinc(ji,jj,jk,jp_fabm_p1p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1p_chl1
	1719	phyto3d_balinc(ji,jj,jk,jp_fabm_p1s) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1s_chl1
	1720	phyto3d_balinc(ji,jj,jk,jp_fabm_p2c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2c_chl2
	1721	phyto3d_balinc(ji,jj,jk,jp_fabm_p2n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2n_chl2
	1722	phyto3d_balinc(ji,jj,jk,jp_fabm_p2p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2p_chl2
	1723	phyto3d_balinc(ji,jj,jk,jp_fabm_p3c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3c_chl3
	1724	phyto3d_balinc(ji,jj,jk,jp_fabm_p3n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3n_chl3
	1725	phyto3d_balinc(ji,jj,jk,jp_fabm_p3p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3p_chl3
	1726	phyto3d_balinc(ji,jj,jk,jp_fabm_p4c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4c_chl4
	1727	phyto3d_balinc(ji,jj,jk,jp_fabm_p4n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4n_chl4
	1728	phyto3d_balinc(ji,jj,jk,jp_fabm_p4p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4p_chl4
	1729	ENDIF
	1730	END DO
	1731	END DO
	1732	END DO
[9326]	1733	#else
	1734	CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', &
	1735	& 'but not defined a biogeochemical model' )
	1736	#endif
	1737
	1738	ENDIF
	1739
	1740	IF ( ll_asmiau ) THEN
	1741
	1742	!--------------------------------------------------------------------
	1743	! Incremental Analysis Updating
	1744	!--------------------------------------------------------------------
	1745
	1746	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
	1747
	1748	it = kt - nit000 + 1
	1749	zincwgt = pwgtiau(it) ! IAU weight for the current time step
	1750	! note this is not a tendency so should not be divided by rdt
	1751
	1752	IF(lwp) THEN
	1753	WRITE(numout,*)
	1754	WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', &
	1755	& kt,' with IAU weight = ', pwgtiau(it)
	1756	WRITE(numout,*) '~~~~~~~~~~~~'
	1757	ENDIF
	1758
	1759	! Update the biogeochemical variables
	1760	! Add directly to trn and trb, rather than to tra, because tra gets
	1761	! reset to zero at the start of trc_stp, called after this routine
[10141]	1762	#if defined key_medusa
[9382]	1763	WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1764	& trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
[9326]	1765	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1766	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1767	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
	1768	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1769	END WHERE
	1770	#elif defined key_hadocc
[9382]	1771	WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	1772	& trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
[9326]	1773	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	1774	& phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	1775	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
	1776	& phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	1777	END WHERE
[10622]	1778	#elif defined key_fabm
	1779	WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
	1780	& trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
	1781	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
	1782	& phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
	1783	trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
	1784	& phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
	1785	END WHERE
[9326]	1786	#endif
	1787
	1788	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1789	! which is called at end of model run
	1790	ENDIF
	1791
	1792	ELSEIF ( ll_asmdin ) THEN
	1793
	1794	!--------------------------------------------------------------------
	1795	! Direct Initialization
	1796	!--------------------------------------------------------------------
	1797
	1798	IF ( kt == nitdin_r ) THEN
	1799
	1800	neuler = 0 ! Force Euler forward step
	1801
	1802	! Initialize the now fields with the background + increment
	1803	! Background currently is what the model is initialised with
	1804	CALL ctl_warn( ' Doing direct initialisation with phyto3d assimilation', &
	1805	& ' Background state is taken from model rather than background file' )
[10141]	1806	#if defined key_medusa
[9382]	1807	WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1808	& trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
[9326]	1809	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1810	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1)
	1811	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
	1812	END WHERE
	1813	#elif defined key_hadocc
[9382]	1814	WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	1815	& trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
[9326]	1816	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	1817	& phyto3d_balinc(:,:,:,jp_had0:jp_had1)
	1818	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
	1819	END WHERE
[10622]	1820	#elif defined key_fabm
	1821	WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
	1822	& trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
	1823	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
	1824	& phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1)
	1825	trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
	1826	END WHERE
[9326]	1827	#endif
	1828
	1829	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1830	! which is called at end of model run
	1831	ENDIF
	1832	!
	1833	ENDIF
	1834	!
	1835	END SUBROUTINE phyto3d_asm_inc
	1836
	1837	!!===========================================================================
	1838	!!===========================================================================
	1839	!!===========================================================================
	1840
[9322]	1841	SUBROUTINE pco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
	1842	& ll_trainc, pt_bkginc, ps_bkginc )
[9297]	1843	!!------------------------------------------------------------------------
[9296]	1844	!! * ROUTINE pco2_asm_inc *
	1845	!!
	1846	!! ** Purpose : Apply the pco2/fco2 assimilation increments.
	1847	!!
	1848	!! ** Method : Calculate increments to state variables using carbon
	1849	!! balancing.
	1850	!! Direct initialization or Incremental Analysis Updating.
	1851	!!
	1852	!! ** Action :
[9297]	1853	!!------------------------------------------------------------------------
[9296]	1854	INTEGER, INTENT(IN) :: kt ! Current time step
	1855	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
	1856	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
	1857	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
	1858	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
	1859	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
	1860	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
	1861	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
	1862	!
	1863	INTEGER :: ji, jj, jk ! Loop counters
	1864	INTEGER :: it ! Index
	1865	INTEGER :: jkmax ! Index of mixed layer depth
	1866	REAL(wp) :: zincwgt ! IAU weight for current time step
	1867	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth for balancing
	1868	REAL(wp), DIMENSION(jpi,jpj) :: pco2_bkginc_temp ! For fCO2 conversion
	1869	REAL(wp), DIMENSION(jpi,jpj) :: dic_bkg_temp ! DIC background
	1870	REAL(wp), DIMENSION(jpi,jpj) :: alk_bkg_temp ! Alkalinity background
	1871	REAL(wp), DIMENSION(jpi,jpj) :: tem_bkg_temp ! Temperature background
	1872	REAL(wp), DIMENSION(jpi,jpj) :: sal_bkg_temp ! Salinity background
	1873	REAL(wp), DIMENSION(1) :: patm ! Atmospheric pressure
	1874	REAL(wp), DIMENSION(1) :: phyd ! Hydrostatic pressure
	1875
	1876	! Coefficients for fCO2 to pCO2 conversion
	1877	REAL(wp) :: zcoef_fco2_1 = -1636.75
	1878	REAL(wp) :: zcoef_fco2_2 = 12.0408
	1879	REAL(wp) :: zcoef_fco2_3 = 0.0327957
	1880	REAL(wp) :: zcoef_fco2_4 = 0.0000316528
	1881	REAL(wp) :: zcoef_fco2_5 = 2.0
	1882	REAL(wp) :: zcoef_fco2_6 = 57.7
	1883	REAL(wp) :: zcoef_fco2_7 = 0.118
	1884	REAL(wp) :: zcoef_fco2_8 = 82.0578
[9297]	1885	!!------------------------------------------------------------------------
[9296]	1886
	1887	IF ( kt <= nit000 ) THEN
	1888
[9297]	1889	!----------------------------------------------------------------------
[9296]	1890	! DIC and alkalinity balancing
[9297]	1891	!----------------------------------------------------------------------
[9296]	1892
	1893	IF ( ln_sfco2inc ) THEN
[10141]	1894	#if defined key_medusa && defined key_roam
[9296]	1895	! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
	1896	patm(1) = 1.0
	1897	phyd(1) = 0.0
	1898	DO jj = 1, jpj
	1899	DO ji = 1, jpi
[9322]	1900	CALL f2pCO2( sfco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
[9296]	1901	END DO
	1902	END DO
	1903	#else
	1904	! If assimilating fCO2, then convert to pCO2 using temperature
	1905	! See flux_gas.F90 within HadOCC for details of calculation
[9326]	1906	pco2_bkginc_temp(:,:) = sfco2_bkginc(:,:) / &
[9296]	1907	& EXP((zcoef_fco2_1 + &
	1908	& zcoef_fco2_2 * (tsn(:,:,1,1)+rt0) - &
	1909	& zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
	1910	& zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0) + &
	1911	& zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0))) / &
	1912	& (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
	1913	#endif
	1914	ELSE
	1915	pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
	1916	ENDIF
	1917
	1918	! Account for physics assimilation if required
	1919	IF ( ll_trainc ) THEN
	1920	tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
	1921	sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
	1922	ELSE
	1923	tem_bkg_temp(:,:) = tsn(:,:,1,1)
	1924	sal_bkg_temp(:,:) = tsn(:,:,1,2)
	1925	ENDIF
	1926
[10141]	1927	#if defined key_medusa
[9381]	1928	! Account for phytoplankton balancing if required
	1929	IF ( ln_phytobal ) THEN
[9326]	1930	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
	1931	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk)
[9296]	1932	ELSE
	1933	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
	1934	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
	1935	ENDIF
	1936
	1937	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
	1938	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
[9322]	1939	& pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
[9296]	1940
	1941	#elif defined key_hadocc
[9381]	1942	! Account for phytoplankton balancing if required
	1943	IF ( ln_phytobal ) THEN
[9326]	1944	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
	1945	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk)
[9296]	1946	ELSE
	1947	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
	1948	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
	1949	ENDIF
	1950
	1951	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
	1952	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
[9322]	1953	& pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
[9296]	1954
[10622]	1955	#elif defined key_fabm
	1956	! Account for phytoplankton balancing if required
	1957	IF ( ln_phytobal ) THEN
	1958	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_o3c) + phyto2d_balinc(:,:,1,jp_fabm_o3c)
	1959	alk_bkg_temp(:,:) = totalk_bkg(:,:,1) + phyto2d_balinc(:,:,1,jp_fabm_o3ba)
	1960	ELSE
	1961	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_o3c)
	1962	alk_bkg_temp(:,:) = totalk_bkg(:,:,1)
	1963	ENDIF
	1964
	1965	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
	1966	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
	1967	& pco2_balinc(:,:,1,jp_fabm_o3c), pco2_balinc(:,:,1,jp_fabm_o3ba) )
	1968
[9296]	1969	#else
	1970	CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
	1971	& 'but not defined a biogeochemical model' )
	1972	#endif
	1973
	1974	! Select mixed layer
	1975	IF ( ll_asmdin ) THEN
[10622]	1976	#if defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm
[9296]	1977	CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
	1978	& ' Mixed layer depth taken to be background maximum mld_max_bkg' )
	1979	zmld(:,:) = mld_max_bkg(:,:)
	1980	#else
	1981	CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
	1982	#endif
	1983	ELSE
	1984	SELECT CASE( mld_choice_bgc )
	1985	CASE ( 1 ) ! Turbocline/mixing depth [W points]
	1986	zmld(:,:) = hmld(:,:)
	1987	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
	1988	zmld(:,:) = hmlp(:,:)
	1989	CASE ( 3 ) ! Kara MLD [Interpolated]
	1990	#if defined key_karaml
	1991	IF ( ln_kara ) THEN
	1992	zmld(:,:) = hmld_kara(:,:)
	1993	ELSE
	1994	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
	1995	& ' but ln_kara=.false.' )
	1996	ENDIF
	1997	#else
	1998	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
	1999	& ' but is not defined' )
	2000	#endif
	2001	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
[10622]	2002	zmld(:,:) = hmld_tref(:,:)
[9296]	2003	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
	2004	zmld(:,:) = hmlpt(:,:)
	2005	END SELECT
	2006	ENDIF
	2007
	2008	! Propagate through mixed layer
	2009	DO jj = 1, jpj
	2010	DO ji = 1, jpi
	2011	!
	2012	jkmax = jpk-1
	2013	DO jk = jpk-1, 1, -1
	2014	IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. &
	2015	& ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
	2016	zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
	2017	jkmax = jk
	2018	ENDIF
	2019	END DO
	2020	!
	2021	#if defined key_top
	2022	DO jk = 2, jkmax
[9322]	2023	pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
[9296]	2024	END DO
	2025	#endif
	2026	!
	2027	END DO
	2028	END DO
	2029
	2030	ENDIF
	2031
	2032	IF ( ll_asmiau ) THEN
	2033
	2034	!--------------------------------------------------------------------
	2035	! Incremental Analysis Updating
	2036	!--------------------------------------------------------------------
	2037
	2038	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
	2039
	2040	it = kt - nit000 + 1
	2041	zincwgt = pwgtiau(it) ! IAU weight for the current time step
	2042	! note this is not a tendency so should not be divided by rdt
	2043
	2044	IF(lwp) THEN
	2045	WRITE(numout,*)
	2046	IF ( ln_spco2inc ) THEN
	2047	WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
	2048	& kt,' with IAU weight = ', pwgtiau(it)
	2049	ELSE IF ( ln_sfco2inc ) THEN
	2050	WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
	2051	& kt,' with IAU weight = ', pwgtiau(it)
	2052	ENDIF
	2053	WRITE(numout,*) '~~~~~~~~~~~~'
	2054	ENDIF
	2055
	2056	! Update the biogeochemical variables
	2057	! Add directly to trn and trb, rather than to tra, because tra gets
	2058	! reset to zero at the start of trc_stp, called after this routine
[10141]	2059	#if defined key_medusa
[9382]	2060	WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	2061	& trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
	2062	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	2063	& pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	2064	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
	2065	& pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	2066	END WHERE
[9296]	2067	#elif defined key_hadocc
[9382]	2068	WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	2069	& trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
	2070	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	2071	& pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	2072	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
	2073	& pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	2074	END WHERE
[10622]	2075	#elif defined key_fabm
	2076	WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
	2077	& trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
	2078	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
	2079	& pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
	2080	trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
	2081	& pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
	2082	END WHERE
[9296]	2083	#endif
	2084
	2085	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	2086	! which is called at end of model run
	2087
	2088	ENDIF
	2089
	2090	ELSEIF ( ll_asmdin ) THEN
	2091
	2092	!--------------------------------------------------------------------
	2093	! Direct Initialization
	2094	!--------------------------------------------------------------------
	2095
	2096	IF ( kt == nitdin_r ) THEN
	2097
	2098	neuler = 0 ! Force Euler forward step
	2099
	2100	! Initialize the now fields with the background + increment
	2101	! Background currently is what the model is initialised with
[9382]	2102	CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
[9296]	2103	& ' Background state is taken from model rather than background file' )
[10141]	2104	#if defined key_medusa
[9382]	2105	WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	2106	& trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
	2107	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	2108	& pco2_balinc(:,:,:,jp_msa0:jp_msa1)
	2109	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
	2110	END WHERE
[9296]	2111	#elif defined key_hadocc
[9382]	2112	WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	2113	& trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
	2114	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	2115	& pco2_balinc(:,:,:,jp_had0:jp_had1)
	2116	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
	2117	END WHERE
[10622]	2118	#elif defined key_fabm
	2119	WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
	2120	& trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
	2121	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
	2122	& pco2_balinc(:,:,:,jp_fabm0:jp_fabm1)
	2123	trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
	2124	END WHERE
[9296]	2125	#endif
	2126
	2127	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	2128	! which is called at end of model run
	2129	ENDIF
	2130	!
	2131	ENDIF
	2132	!
[9322]	2133	END SUBROUTINE pco2_asm_inc
[9296]	2134
[9297]	2135	!!===========================================================================
[9322]	2136	!!===========================================================================
	2137	!!===========================================================================
[9296]	2138
[9322]	2139	SUBROUTINE ph_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
	2140	& ll_trainc, pt_bkginc, ps_bkginc )
	2141	!!------------------------------------------------------------------------
	2142	!! * ROUTINE ph_asm_inc *
	2143	!!
	2144	!! ** Purpose : Apply the pH assimilation increments.
	2145	!!
[9381]	2146	!! ** Method : Calculate increments to DIC and alkalinity from pH
	2147	!! Use a similar approach to the pCO2 scheme
[9322]	2148	!! Direct initialization or Incremental Analysis Updating.
	2149	!!
	2150	!! ** Action :
	2151	!!------------------------------------------------------------------------
[9381]	2152	INTEGER, INTENT(IN) :: kt ! Current time step
	2153	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
	2154	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
	2155	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
	2156	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
	2157	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
[9322]	2158	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
	2159	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
[9381]	2160
	2161	REAL(wp) :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation
	2162	REAL(wp) :: DIC_IN, ALK_IN ! DIC/alk in pH calculation
	2163	REAL(wp) :: PH_OUT ! pH from pH calculation
	2164	REAL(wp) :: ph_bkg ! pH from pH calculation
	2165	REAL(wp) :: ph_dic ! pH from pH calculation
	2166	REAL(wp) :: ph_alk ! pH from pH calculation
	2167	REAL(wp) :: dph_ddic ! Change in pH wrt DIC
	2168	REAL(wp) :: dph_dalk ! Change in pH wrt alk
	2169	REAL(wp) :: weight ! Increment weighting
	2170	REAL(wp) :: zincwgt ! IAU weight for current time step
	2171	REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp ! DIC background
	2172	REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp ! Alkalinity background
	2173	REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp ! DIN background
	2174	REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp ! Silicate background
	2175	REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp ! Temperature background
	2176	REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp ! Salinity background
	2177	REAL(wp), DIMENSION(19) :: dummy_out ! Dummy variables
	2178	INTEGER :: ji, jj, jk, jx ! Loop counters
	2179	INTEGER :: it ! Index
[9322]	2180	!!------------------------------------------------------------------------
[9381]	2181
[10141]	2182	#if ! defined key_medusa
[9381]	2183	CALL ctl_stop( ' pH balancing only implemented for MEDUSA' )
	2184	#else
	2185
	2186	IF ( kt <= nit000 ) THEN
	2187
	2188	!----------------------------------------------------------------------
	2189	! DIC and alkalinity balancing
	2190	!----------------------------------------------------------------------
	2191
	2192	! Account for physics assimilation if required
	2193	IF ( ll_trainc ) THEN
	2194	tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:)
	2195	sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:)
	2196	ELSE
	2197	tem_bkg_temp(:,:,:) = tsn(:,:,:,1)
	2198	sal_bkg_temp(:,:,:) = tsn(:,:,:,2)
	2199	ENDIF
	2200
	2201	! Account for phytoplankton balancing if required
	2202	IF ( ln_phytobal ) THEN
	2203	dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic)
	2204	alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk)
	2205	din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin)
	2206	sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil)
	2207	ELSE
	2208	dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic)
	2209	alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk)
	2210	din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin)
	2211	sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil)
	2212	ENDIF
	2213
	2214	! Account for pCO2 balancing if required
	2215	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	2216	dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic)
	2217	alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk)
	2218	ENDIF
	2219
	2220	! Loop over grid points and calculate dpH/dDIC and dpH/dAlk
	2221	! This requires three calls to the MOCSY carbonate package
	2222	! One to find pH at background DIC and alkalinity
	2223	! One to find pH when DIC is increased by 10
	2224	! One to find pH when alkalinity is increased by 10
	2225	! Then calculate the gradients
	2226
	2227	DO jk = 1, jpk
	2228	DO jj = 1, jpj
	2229	DO ji = 1, jpi
	2230
	2231	IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN
	2232
	2233	DO jx = 1, 3
	2234	IF ( jx == 1 ) THEN
	2235	DIC_IN = dic_bkg_temp(ji,jj,jk)
	2236	ALK_IN = alk_bkg_temp(ji,jj,jk)
	2237	ELSE IF ( jx == 2 ) THEN
	2238	DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch
	2239	ALK_IN = alk_bkg_temp(ji,jj,jk)
	2240	ELSE IF ( jx == 3 ) THEN
	2241	DIC_IN = dic_bkg_temp(ji,jj,jk)
	2242	ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch
	2243	ENDIF
	2244
	2245	CALL mocsy_interface(tsn(ji,jj,jk,jp_tem), & ! IN: temperature
	2246	& tsn(ji,jj,jk,jp_sal), & ! salinity
	2247	& ALK_IN, & ! alkalinity
	2248	& DIC_IN, & ! DIC
	2249	& sil_bkg_temp(ji,jj,jk), & ! silicate
	2250	& din_bkg_temp(ji,jj,jk)/16.0, & ! phosphate (Redfield from DIN)
	2251	& 1.0, & ! atmospheric pressure (dummy)
	2252	& fsdept(ji,jj,jk), & ! depth
	2253	& gphit(ji,jj), & ! latitude
	2254	& 1.0, & ! gas transfer (dummy)
	2255	& 1.0, & ! atmospheric xCO2 (dummy)
	2256	& 1, & ! number of points
	2257	& PH_OUT, & ! OUT: pH
	2258	& dummy_out(1), dummy_out(2), & ! stuff we don't care about
	2259	& dummy_out(3), dummy_out(4), &
	2260	& dummy_out(5), dummy_out(6), &
	2261	& dummy_out(7), dummy_out(8), &
	2262	& dummy_out(9), dummy_out(10), &
	2263	& dummy_out(11), dummy_out(12), &
	2264	& dummy_out(13), dummy_out(14), &
	2265	& dummy_out(15), dummy_out(16), &
	2266	& dummy_out(17), dummy_out(18), &
	2267	& dummy_out(19) )
	2268
	2269	IF ( jx == 1 ) THEN
	2270	ph_bkg = PH_OUT
	2271	ELSE IF ( jx == 2 ) THEN
	2272	ph_dic = PH_OUT
	2273	ELSE IF ( jx == 3 ) THEN
	2274	ph_alk = PH_OUT
	2275	ENDIF
	2276	END DO
	2277
	2278	dph_ddic = (ph_dic - ph_bkg) / zsearch
	2279	dph_dalk = (ph_alk - ph_bkg) / zsearch
	2280	weight = pph_bkginc(ji,jj,jk) / ( (dph_ddicdph_ddic) + (dph_dalkdph_dalk) )
	2281
	2282	ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic
	2283	ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk
	2284
	2285	ENDIF
	2286
	2287	END DO
	2288	END DO
	2289	END DO
	2290
	2291	ENDIF
[9322]	2292
[9381]	2293	IF ( ll_asmiau ) THEN
	2294
	2295	!--------------------------------------------------------------------
	2296	! Incremental Analysis Updating
	2297	!--------------------------------------------------------------------
	2298
	2299	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
	2300
	2301	it = kt - nit000 + 1
	2302	zincwgt = pwgtiau(it) ! IAU weight for the current time step
	2303	! note this is not a tendency so should not be divided by rdt
	2304
	2305	IF(lwp) THEN
	2306	WRITE(numout,*)
	2307	WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', &
	2308	& kt,' with IAU weight = ', pwgtiau(it)
	2309	WRITE(numout,*) '~~~~~~~~~~~~'
	2310	ENDIF
	2311
	2312	! Update the biogeochemical variables
	2313	! Add directly to trn and trb, rather than to tra, because tra gets
	2314	! reset to zero at the start of trc_stp, called after this routine
[9382]	2315	WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	2316	& trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
	2317	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	2318	& ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	2319	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
	2320	& ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	2321	END WHERE
[9381]	2322
	2323	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	2324	! which is called at end of model run
	2325
	2326	ENDIF
	2327
	2328	ELSEIF ( ll_asmdin ) THEN
	2329
	2330	!--------------------------------------------------------------------
	2331	! Direct Initialization
	2332	!--------------------------------------------------------------------
	2333
	2334	IF ( kt == nitdin_r ) THEN
	2335
	2336	neuler = 0 ! Force Euler forward step
	2337
	2338	! Initialize the now fields with the background + increment
	2339	! Background currently is what the model is initialised with
	2340	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pH assimilation', &
	2341	& ' Background state is taken from model rather than background file' )
[9382]	2342	WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	2343	& trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
	2344	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	2345	& ph_balinc(:,:,:,jp_msa0:jp_msa1)
	2346	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
	2347	END WHERE
[9381]	2348
	2349	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	2350	! which is called at end of model run
	2351	ENDIF
	2352	!
	2353	ENDIF
	2354	#endif
[9322]	2355	!
	2356	END SUBROUTINE ph_asm_inc
	2357
	2358	!!===========================================================================
	2359	!!===========================================================================
	2360	!!===========================================================================
	2361
	2362	SUBROUTINE bgc3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
	2363	!!----------------------------------------------------------------------
	2364	!! * ROUTINE dyn_asm_inc *
	2365	!!
	2366	!! ** Purpose : Apply generic 3D biogeochemistry assimilation increments.
	2367	!!
	2368	!! ** Method : Direct initialization or Incremental Analysis Updating.
	2369	!!
	2370	!! ** Action :
	2371	!!----------------------------------------------------------------------
	2372	INTEGER, INTENT(IN) :: kt ! Current time step
	2373	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
	2374	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
	2375	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
	2376	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
	2377	!
	2378	INTEGER :: it ! Index
	2379	REAL(wp) :: zincwgt ! IAU weight for current time step
	2380	!!----------------------------------------------------------------------
	2381
[9382]	2382	IF ( kt <= nit000 ) THEN
	2383
	2384	!----------------------------------------------------------------------
	2385	! Remove any other balancing increments
	2386	!----------------------------------------------------------------------
	2387
	2388	IF ( ln_pno3inc ) THEN
[10622]	2389	#if defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm
[9382]	2390	#if defined key_hadocc
	2391	it = jp_had_nut
[10141]	2392	#elif defined key_medusa
[9382]	2393	it = jpdin
[10622]	2394	#elif defined key_fabm
	2395	it = jp_fabm_n3n
[9382]	2396	#endif
	2397	IF ( ln_phytobal ) THEN
	2398	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
	2399	ENDIF
	2400	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
	2401	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
	2402	ENDIF
	2403	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	2404	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
	2405	ENDIF
	2406	IF ( ln_pphinc ) THEN
	2407	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - ph_balinc(:,:,:,it)
	2408	ENDIF
[10057]	2409	#else
	2410	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2411	#endif
[9382]	2412	ENDIF
	2413
	2414	IF ( ln_psi4inc ) THEN
[10622]	2415	#if defined key_medusa \|\| defined key_fabm
[10141]	2416	#if defined key_medusa
[9382]	2417	it = jpsil
[10622]	2418	#elif defined key_fabm
	2419	it = jp_fabm_n5s
	2420	#endif
[9382]	2421	IF ( ln_phytobal ) THEN
	2422	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
	2423	ENDIF
	2424	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
	2425	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
	2426	ENDIF
	2427	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	2428	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
	2429	ENDIF
	2430	IF ( ln_pphinc ) THEN
	2431	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - ph_balinc(:,:,:,it)
	2432	ENDIF
[10057]	2433	#else
	2434	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2435	#endif
[9382]	2436	ENDIF
	2437
	2438	IF ( ln_pdicinc ) THEN
[10622]	2439	#if defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm
[9382]	2440	#if defined key_hadocc
	2441	it = jp_had_dic
[10141]	2442	#elif defined key_medusa
[9382]	2443	it = jpdic
[10622]	2444	#elif defined key_fabm
	2445	it = jp_fabm_o3c
[9382]	2446	#endif
	2447	IF ( ln_phytobal ) THEN
	2448	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
	2449	ENDIF
	2450	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
	2451	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
	2452	ENDIF
	2453	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	2454	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
	2455	ENDIF
	2456	IF ( ln_pphinc ) THEN
	2457	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - ph_balinc(:,:,:,it)
	2458	ENDIF
[10057]	2459	#else
	2460	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2461	#endif
[9382]	2462	ENDIF
	2463
	2464	IF ( ln_palkinc ) THEN
[10622]	2465	#if defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm
[9382]	2466	#if defined key_hadocc
	2467	it = jp_had_alk
[10141]	2468	#elif defined key_medusa
[9382]	2469	it = jpalk
[10622]	2470	#elif defined key_fabm
	2471	it = jp_fabm_o3ba
[9382]	2472	#endif
	2473	IF ( ln_phytobal ) THEN
	2474	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
	2475	ENDIF
	2476	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
	2477	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
	2478	ENDIF
	2479	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	2480	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
	2481	ENDIF
	2482	IF ( ln_pphinc ) THEN
	2483	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - ph_balinc(:,:,:,it)
	2484	ENDIF
[10057]	2485	#else
	2486	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2487	#endif
[9382]	2488	ENDIF
	2489
	2490	IF ( ln_po2inc ) THEN
[10622]	2491	#if defined key_medusa \|\| defined key_fabm
[10141]	2492	#if defined key_medusa
[9382]	2493	it = jpoxy
[10622]	2494	#elif defined key_fabm
	2495	it = jp_fabm_o2o
	2496	#endif
[9382]	2497	IF ( ln_phytobal ) THEN
	2498	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
	2499	ENDIF
	2500	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
	2501	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
	2502	ENDIF
	2503	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	2504	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
	2505	ENDIF
	2506	IF ( ln_pphinc ) THEN
	2507	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - ph_balinc(:,:,:,it)
	2508	ENDIF
[10057]	2509	#else
	2510	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2511	#endif
[9382]	2512	ENDIF
	2513
	2514	ENDIF
	2515
[9322]	2516	IF ( ll_asmiau ) THEN
	2517
	2518	!--------------------------------------------------------------------
	2519	! Incremental Analysis Updating
	2520	!--------------------------------------------------------------------
	2521
	2522	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
	2523
	2524	it = kt - nit000 + 1
[9326]	2525	zincwgt = pwgtiau(it) ! IAU weight for the current time step
	2526	! note this is not a tendency so should not be divided by rdt
[9322]	2527
	2528	IF(lwp) THEN
	2529	WRITE(numout,*)
	2530	WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', &
	2531	& kt,' with IAU weight = ', pwgtiau(it)
	2532	WRITE(numout,*) '~~~~~~~~~~~~'
	2533	ENDIF
	2534
	2535	! Update the 3D BGC variables
	2536	! Add directly to trn and trb, rather than to tra, because tra gets
	2537	! reset to zero at the start of trc_stp, called after this routine
	2538	! Don't apply increments if they'll take concentrations negative
	2539
	2540	IF ( ln_pno3inc ) THEN
	2541	#if defined key_hadocc
	2542	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
	2543	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2544	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
	2545	trb(:,:,:,jp_had_nut) = trb(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
	2546	END WHERE
[10141]	2547	#elif defined key_medusa
[9322]	2548	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
	2549	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2550	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
	2551	trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
	2552	END WHERE
[10622]	2553	#elif defined key_fabm
	2554	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
	2555	& trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2556	trn(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt
	2557	trb(:,:,:,jp_fabm_n3n) = trb(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt
	2558	END WHERE
[9322]	2559	#else
	2560	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2561	#endif
	2562	ENDIF
	2563
	2564	IF ( ln_psi4inc ) THEN
[10141]	2565	#if defined key_medusa
[9322]	2566	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
	2567	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2568	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
	2569	trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
	2570	END WHERE
[10622]	2571	#elif defined key_fabm
	2572	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
	2573	& trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2574	trn(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt
	2575	trb(:,:,:,jp_fabm_n5s) = trb(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt
	2576	END WHERE
[9322]	2577	#else
	2578	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2579	#endif
	2580	ENDIF
	2581
	2582	IF ( ln_pdicinc ) THEN
	2583	#if defined key_hadocc
	2584	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
	2585	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2586	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
	2587	trb(:,:,:,jp_had_dic) = trb(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
	2588	END WHERE
[10141]	2589	#elif defined key_medusa
[9322]	2590	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
	2591	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2592	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
	2593	trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
	2594	END WHERE
[10622]	2595	#elif defined key_fabm
	2596	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
	2597	& trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2598	trn(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt
	2599	trb(:,:,:,jp_fabm_o3c) = trb(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt
	2600	END WHERE
[9322]	2601	#else
	2602	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2603	#endif
	2604	ENDIF
	2605
	2606	IF ( ln_palkinc ) THEN
	2607	#if defined key_hadocc
	2608	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
	2609	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2610	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
	2611	trb(:,:,:,jp_had_alk) = trb(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
	2612	END WHERE
[10141]	2613	#elif defined key_medusa
[9322]	2614	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
	2615	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2616	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
	2617	trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
	2618	END WHERE
[10622]	2619	#elif defined key_fabm
	2620	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
	2621	& trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2622	trn(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt
	2623	trb(:,:,:,jp_fabm_o3ba) = trb(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt
	2624	END WHERE
[9322]	2625	#else
	2626	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2627	#endif
	2628	ENDIF
	2629
	2630	IF ( ln_po2inc ) THEN
[10141]	2631	#if defined key_medusa
[9322]	2632	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
	2633	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2634	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
	2635	trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
	2636	END WHERE
[10622]	2637	#elif defined key_fabm
	2638	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
	2639	& trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2640	trn(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt
	2641	trb(:,:,:,jp_fabm_o2o) = trb(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt
	2642	END WHERE
[9322]	2643	#else
	2644	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2645	#endif
	2646	ENDIF
	2647
	2648	IF ( kt == nitiaufin_r ) THEN
	2649	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
	2650	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
	2651	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
	2652	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
	2653	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
	2654	ENDIF
	2655
	2656	ENDIF
	2657
	2658	ELSEIF ( ll_asmdin ) THEN
	2659
	2660	!--------------------------------------------------------------------
	2661	! Direct Initialization
	2662	!--------------------------------------------------------------------
	2663
	2664	IF ( kt == nitdin_r ) THEN
	2665
	2666	neuler = 0 ! Force Euler forward step
	2667
	2668	! Initialize the now fields with the background + increment
	2669	! Background currently is what the model is initialised with
[10622]	2670	CALL ctl_warn( ' Doing direct initialisation with 3D BGC assimilation', &
[9322]	2671	& ' Background state is taken from model rather than background file' )
	2672
	2673	IF ( ln_pno3inc ) THEN
	2674	#if defined key_hadocc
	2675	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
	2676	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) > 0.0_wp )
	2677	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:)
	2678	trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut)
	2679	END WHERE
[10141]	2680	#elif defined key_medusa
[9322]	2681	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
	2682	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) > 0.0_wp )
	2683	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:)
	2684	trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
	2685	END WHERE
[10622]	2686	#elif defined key_fabm
	2687	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
	2688	& trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) > 0.0_wp )
	2689	trn(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:)
	2690	trb(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n)
	2691	END WHERE
[9322]	2692	#else
	2693	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2694	#endif
	2695	ENDIF
	2696
	2697	IF ( ln_psi4inc ) THEN
[10141]	2698	#if defined key_medusa
[9322]	2699	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
	2700	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) > 0.0_wp )
	2701	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:)
	2702	trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
	2703	END WHERE
[10622]	2704	#elif defined key_fabm
	2705	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
	2706	& trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) > 0.0_wp )
	2707	trn(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:)
	2708	trb(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s)
	2709	END WHERE
[9322]	2710	#else
	2711	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2712	#endif
	2713	ENDIF
	2714
	2715	IF ( ln_pdicinc ) THEN
	2716	#if defined key_hadocc
	2717	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
	2718	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) > 0.0_wp )
	2719	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:)
	2720	trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic)
	2721	END WHERE
[10141]	2722	#elif defined key_medusa
[9322]	2723	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
	2724	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) > 0.0_wp )
	2725	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:)
	2726	trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
	2727	END WHERE
[10622]	2728	#elif defined key_fabm
	2729	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
	2730	& trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) > 0.0_wp )
	2731	trn(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:)
	2732	trb(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c)
	2733	END WHERE
[9322]	2734	#else
	2735	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2736	#endif
	2737	ENDIF
	2738
	2739	IF ( ln_palkinc ) THEN
	2740	#if defined key_hadocc
	2741	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
	2742	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) > 0.0_wp )
	2743	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:)
	2744	trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk)
	2745	END WHERE
[10141]	2746	#elif defined key_medusa
[9322]	2747	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
	2748	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) > 0.0_wp )
	2749	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:)
	2750	trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
	2751	END WHERE
[10622]	2752	#elif defined key_fabm
	2753	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
	2754	& trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) > 0.0_wp )
	2755	trn(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:)
	2756	trb(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba)
	2757	END WHERE
[9322]	2758	#else
	2759	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2760	#endif
	2761	ENDIF
	2762
	2763	IF ( ln_po2inc ) THEN
[10141]	2764	#if defined key_medusa
[9322]	2765	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
	2766	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) > 0.0_wp )
	2767	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:)
	2768	trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
	2769	END WHERE
[10622]	2770	#elif defined key_fabm
	2771	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
	2772	& trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) > 0.0_wp )
	2773	trn(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:)
	2774	trb(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o)
	2775	END WHERE
[9322]	2776	#else
	2777	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2778	#endif
	2779	ENDIF
	2780
	2781	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
	2782	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
	2783	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
	2784	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
	2785	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
	2786	ENDIF
	2787	!
	2788	ENDIF
	2789	!
	2790	END SUBROUTINE bgc3d_asm_inc
	2791
	2792	!!===========================================================================
	2793
[9296]	2794	END MODULE asmbgc

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

Download in other formats: