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/dev_r5518_GO6_package_FOAMv14_phytobal3d/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_phytobal3d/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 13097

Last change on this file since 13097 was 13097, checked in by dford, 4 years ago
Allow nitrogen balancing for both 2D and 3D chlorophyll increments.
File size: 108.2 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
	13	!! 'key_karaml' : Kara mixed layer depth
[9297]	14	!!---------------------------------------------------------------------------
	15	!! asm_bgc_check_options : check bgc assimilation options
[9296]	16	!! asm_bgc_init_incs : initialise bgc increments
	17	!! asm_bgc_init_bkg : initialise bgc background
[9297]	18	!! asm_bgc_bal_wri : write out bgc balancing increments
[9296]	19	!! asm_bgc_bkg_wri : write out bgc background
[13097]	20	!! asm_bgc_unlog_2d : calculate non-log versions of 2D log increments
	21	!! asm_bgc_unlog_3d : calculate non-log versions of 3D log increments
	22	!! phyto2d_asm_inc : apply the 2D phytoplankton increments
[9326]	23	!! phyto3d_asm_inc : apply the 3D phytoplankton increments
[9322]	24	!! pco2_asm_inc : apply the pCO2/fCO2 increments
	25	!! ph_asm_inc : apply the pH increments
	26	!! bgc3d_asm_inc : apply the generic 3D BGC increments
[9297]	27	!!---------------------------------------------------------------------------
	28	USE par_kind, ONLY: & ! kind parameters
	29	& wp
	30	USE par_oce, ONLY: & ! domain array sizes
	31	& jpi, &
	32	& jpj, &
	33	& jpk
	34	USE iom ! i/o
	35	USE oce, ONLY: & ! active tracer variables
	36	& tsn
	37	USE zdfmxl, ONLY : & ! mixed layer depth
[9296]	38	#if defined key_karaml
[9297]	39	& hmld_kara, &
	40	& ln_kara, &
[9296]	41	#endif
[9297]	42	& hmld, &
	43	& hmlp, &
	44	& hmlpt
	45	USE asmpar, ONLY: & ! assimilation parameters
	46	& c_asmbkg, &
	47	& c_asmbal, &
	48	& nitbkg_r, &
	49	& nitdin_r, &
	50	& nitiaustr_r, &
	51	& nitiaufin_r
[9296]	52	#if defined key_top
[9297]	53	USE trc, ONLY: & ! passive tracer variables
	54	& trn, &
[9296]	55	& trb
[9297]	56	USE par_trc, ONLY: & ! passive tracer parameters
[9296]	57	& jptra
	58	#endif
[10141]	59	#if defined key_medusa
[13097]	60	USE asmphytobal_medusa, ONLY: & ! phytoplankton balancing for MEDUSA
	61	& asm_phyto_bal_medusa
[9297]	62	USE asmpco2bal, ONLY: & ! pCO2 balancing for MEDUSA
[9296]	63	& asm_pco2_bal
[9431]	64	USE sms_medusa ! MEDUSA parameters
[9297]	65	USE par_medusa ! MEDUSA parameters
[9381]	66	USE mocsy_wrapper, ONLY: & ! MEDUSA carbonate system
	67	& mocsy_interface
	68	USE mocsy_mainmod, ONLY: & ! fCO2 to pCO2 conversion
[9297]	69	& f2pCO2
	70	USE sms_medusa, ONLY: & ! MEDUSA diagnostic variables
	71	& pgrow_avg, &
	72	& ploss_avg, &
	73	& phyt_avg, &
[13097]	74	& pgrow_avg_3d, &
	75	& ploss_avg_3d, &
	76	& phyt_avg_3d, &
[9297]	77	& mld_max
[9296]	78	#elif defined key_hadocc
[13097]	79	USE asmphytobal_hadocc, ONLY: & ! phytoplankton balancing for HadOCC
	80	& asm_phyto_bal_hadocc
[9297]	81	USE asmpco2bal, ONLY: & ! pCO2 balancing for HadOCC
[9296]	82	& asm_pco2_bal
[9297]	83	USE par_hadocc ! HadOCC parameters
[9326]	84	USE had_bgc_const ! HadOCC parameters
[9297]	85	USE trc, ONLY: & ! HadOCC diagnostic variables
	86	& pgrow_avg, &
	87	& ploss_avg, &
	88	& phyt_avg, &
[13097]	89	& pgrow_avg_3d, &
	90	& ploss_avg_3d, &
	91	& phyt_avg_3d, &
[9297]	92	& mld_max, &
	93	& HADOCC_CHL
	94	USE had_bgc_const, ONLY: & ! HadOCC C:Chl ratio
	95	& cchl_p
[9296]	96	#endif
	97
	98	IMPLICIT NONE
	99	PRIVATE
	100
[9322]	101	PUBLIC asm_bgc_check_options ! called by asm_inc_init in asminc.F90
	102	PUBLIC asm_bgc_init_incs ! called by asm_inc_init in asminc.F90
	103	PRIVATE asm_bgc_read_incs_2d ! called by asm_bgc_init_incs
	104	PRIVATE asm_bgc_read_incs_3d ! called by asm_bgc_init_incs
	105	PUBLIC asm_bgc_init_bkg ! called by asm_inc_init in asminc.F90
	106	PUBLIC asm_bgc_bal_wri ! called by nemo_gcm in nemogcm.F90
	107	PUBLIC asm_bgc_bkg_wri ! called by asm_bkg_wri in asmbkg.F90
[13097]	108	PRIVATE asm_bgc_unlog_2d ! called by phyto2d_asm_inc
	109	PRIVATE asm_bgc_unlog_3d ! called by phyto3d_asm_inc
[9326]	110	PUBLIC phyto2d_asm_inc ! called by bgc_asm_inc in asminc.F90
	111	PUBLIC phyto3d_asm_inc ! called by bgc_asm_inc in asminc.F90
[9322]	112	PUBLIC pco2_asm_inc ! called by bgc_asm_inc in asminc.F90
	113	PUBLIC ph_asm_inc ! called by bgc_asm_inc in asminc.F90
	114	PUBLIC bgc3d_asm_inc ! called by bgc_asm_inc in asminc.F90
[9296]	115
[9322]	116	LOGICAL, PUBLIC :: ln_balwri = .FALSE. !: No output of balancing incs
	117	LOGICAL, PUBLIC :: ln_phytobal = .FALSE. !: No phytoplankton balancing
	118	LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE. !: No surface total log10(chlorophyll) increment
	119	LOGICAL, PUBLIC :: ln_slchldiainc = .FALSE. !: No surface diatom log10(chlorophyll) increment
	120	LOGICAL, PUBLIC :: ln_slchlnoninc = .FALSE. !: No surface non-diatom log10(chlorophyll) increment
	121	LOGICAL, PUBLIC :: ln_schltotinc = .FALSE. !: No surface total chlorophyll increment
	122	LOGICAL, PUBLIC :: ln_slphytotinc = .FALSE. !: No surface total log10(phyto C) increment
	123	LOGICAL, PUBLIC :: ln_slphydiainc = .FALSE. !: No surface diatom log10(phyto C) increment
	124	LOGICAL, PUBLIC :: ln_slphynoninc = .FALSE. !: No surface non-diatom log10(phyto C) increment
	125	LOGICAL, PUBLIC :: ln_spco2inc = .FALSE. !: No surface pCO2 increment
	126	LOGICAL, PUBLIC :: ln_sfco2inc = .FALSE. !: No surface fCO2 increment
	127	LOGICAL, PUBLIC :: ln_plchltotinc = .FALSE. !: No profile total log10(chlorophyll) increment
	128	LOGICAL, PUBLIC :: ln_pchltotinc = .FALSE. !: No profile total chlorophyll increment
	129	LOGICAL, PUBLIC :: ln_pno3inc = .FALSE. !: No profile nitrate increment
	130	LOGICAL, PUBLIC :: ln_psi4inc = .FALSE. !: No profile silicate increment
	131	LOGICAL, PUBLIC :: ln_pdicinc = .FALSE. !: No profile dissolved inorganic carbon increment
	132	LOGICAL, PUBLIC :: ln_palkinc = .FALSE. !: No profile alkalinity increment
	133	LOGICAL, PUBLIC :: ln_pphinc = .FALSE. !: No profile pH increment
	134	LOGICAL, PUBLIC :: ln_po2inc = .FALSE. !: No profile oxygen increment
[9296]	135
[9297]	136	INTEGER, PUBLIC :: mld_choice_bgc = 1 !: choice of mld for bgc assimilation
	137	! 1) hmld - Turbocline/mixing depth
	138	! [W points]
	139	! 2) hmlp - Density criterion
	140	! (0.01 kg/m^3 change from 10m)
	141	! [W points]
	142	! 3) hmld_kara - Kara MLD
	143	! [Interpolated]
	144	! 4) hmld_tref - Temperature criterion
	145	! (0.2 K change from surface)
	146	! [T points]
	147	! 5) hmlpt - Density criterion
	148	! (0.01 kg/m^3 change from 10m)
	149	! [T points]
	150
	151	REAL(wp), PUBLIC :: rn_maxchlinc = -999.0 !: maximum absolute non-log chlorophyll inc
	152	! <= 0 no maximum applied (switch turned off)
	153	! > 0 absolute chl inc capped at this value
	154
	155	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchltot_bkginc ! slchltot inc
[9322]	156	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchldia_bkginc ! slchldia inc
	157	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchlnon_bkginc ! slchlnon inc
	158	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: schltot_bkginc ! schltot inc
	159	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphytot_bkginc ! slphytot inc
	160	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphydia_bkginc ! slphydia inc
	161	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphynon_bkginc ! slphynon inc
	162	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sfco2_bkginc ! sfco2 inc
	163	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: spco2_bkginc ! spco2 inc
	164	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: plchltot_bkginc ! plchltot inc
	165	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pchltot_bkginc ! pchltot inc
	166	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pno3_bkginc ! pno3 inc
	167	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: psi4_bkginc ! psi4 inc
	168	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pdic_bkginc ! pdic inc
	169	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: palk_bkginc ! palk inc
	170	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pph_bkginc ! pph inc
	171	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: po2_bkginc ! po2 inc
[9326]	172	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto2d_balinc ! Balancing incs from ocean colour
	173	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto3d_balinc ! Balancing incs from p(l)chltot
[9322]	174	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc ! Balancing incs from spco2/sfco2
	175	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ph_balinc ! Balancing incs from pph
[9297]	176
	177	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pgrow_avg_bkg ! Background phyto growth
	178	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ploss_avg_bkg ! Background phyto loss
	179	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: phyt_avg_bkg ! Background phyto conc
[13097]	180	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pgrow_avg_3d_bkg ! Background phyto growth
	181	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ploss_avg_3d_bkg ! Background phyto loss
	182	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: phyt_avg_3d_bkg ! Background phyto conc
[9297]	183	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mld_max_bkg ! Background max MLD
	184	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg ! Background tracer state
[9322]	185	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: chl_bkg ! Background chl
	186	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cchl_p_bkg ! Background c:chl
[9296]	187
[9381]	188	# include "domzgr_substitute.h90"
	189
[9296]	190	CONTAINS
	191
	192	SUBROUTINE asm_bgc_check_options
[9297]	193	!!------------------------------------------------------------------------
[9296]	194	!! * ROUTINE asm_bgc_check_options *
	195	!!
	196	!! ** Purpose : check validity of biogeochemical assimilation options
	197	!!
	198	!! ** Method : check settings of logicals
	199	!!
	200	!! ** Action : call ctl_stop/ctl_warn if required
	201	!!
	202	!! References : asm_inc_init
[9297]	203	!!------------------------------------------------------------------------
[9296]	204
[10141]	205	#if ! defined key_top \|\| ( ! defined key_hadocc && ! defined key_medusa )
[9322]	206	CALL ctl_stop( ' Attempting to assimilate biogeochemical observations', &
	207	& ' but no compatible biogeochemical model is available' )
	208	#endif
	209
	210	#if defined key_hadocc
	211	IF ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slphydiainc .OR. &
	212	& ln_slphynoninc .OR. ln_psi4inc .OR. ln_pphinc .OR. ln_po2inc ) THEN
	213	CALL ctl_stop( ' Cannot assimilate PFTs, Si4, pH or O2 into HadOCC' )
[9296]	214	ENDIF
[9322]	215	#endif
[9296]	216
[10141]	217	#if defined key_medusa
	218	IF ( .NOT. ln_foam_medusa ) THEN
	219	CALL ctl_stop( ' MEDUSA data assimilation options not turned on: set ln_foam_medusa' )
	220	ENDIF
	221	#endif
	222
[9322]	223	IF ( ( ln_phytobal ).AND. &
	224	& ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
	225	& ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc ).AND. &
	226	& ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
[13097]	227	& ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. &
	228	& ( .NOT. ln_pchltotinc ) ) THEN
[9322]	229	CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', &
[13097]	230	& ' if not assimilating phytoplankton,', &
[9322]	231	& ' so ln_phytobal will be set to .false.')
	232	ln_phytobal = .FALSE.
[9296]	233	ENDIF
	234
[9381]	235	IF ( ( ln_balwri ).AND. &
	236	& ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
	237	& ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc ).AND. &
	238	& ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
	239	& ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. &
	240	& ( .NOT. ln_pchltotinc ).AND.( .NOT. ln_pphinc ).AND. &
	241	& ( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ) ) THEN
	242	CALL ctl_warn( ' Balancing increments not being calculated', &
	243	& ' so ln_balwri will be set to .false.')
[9322]	244	ln_balwri = .FALSE.
	245	ENDIF
	246
	247	IF ( ln_spco2inc .AND. ln_sfco2inc ) THEN
[9296]	248	CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
	249	ENDIF
	250
[9322]	251	IF ( ln_slchltotinc .AND. ln_schltotinc ) THEN
	252	CALL ctl_stop( ' Can only assimilate surface log10(chlorophyll) or chlorophyll, not both' )
	253	ENDIF
	254
	255	IF ( ln_plchltotinc .AND. ln_pchltotinc ) THEN
	256	CALL ctl_stop( ' Can only assimilate profile log10(chlorophyll) or chlorophyll, not both' )
	257	ENDIF
	258
	259	IF ( ( ln_slchltotinc .OR. ln_schltotinc ) .AND. &
	260	& ( ln_slchldiainc .OR. ln_slchlnoninc ) ) THEN
	261	CALL ctl_stop( ' Can only assimilate total or PFT surface chlorophyll, not both' )
	262	ENDIF
	263
[9435]	264	IF ( ln_phytobal .AND. &
	265	& ( ( ln_slchlnoninc .AND.( .NOT. ln_slchldiainc ) ).OR. &
	266	& ( ln_slchldiainc .AND.( .NOT. ln_slchlnoninc ) ) ) ) THEN
	267	CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
	268	& ' unless assimilating all model PFTs,')
	269	ENDIF
	270
[9322]	271	IF ( ln_slphytotinc .AND. ( ln_slphydiainc .OR. ln_slphynoninc ) ) THEN
	272	CALL ctl_stop( ' Can only assimilate total or PFT surface phytoplankton carbon, not both' )
	273	ENDIF
	274
[9435]	275	IF ( ln_phytobal .AND. &
	276	& ( ( ln_slphynoninc .AND.( .NOT. ln_slphydiainc ) ).OR. &
	277	& ( ln_slphydiainc .AND.( .NOT. ln_slphynoninc ) ) ) ) THEN
	278	CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
	279	& ' unless assimilating all model PFTs,')
	280	ENDIF
	281
[9296]	282	END SUBROUTINE asm_bgc_check_options
	283
[9297]	284	!!===========================================================================
	285	!!===========================================================================
	286	!!===========================================================================
[9296]	287
	288	SUBROUTINE asm_bgc_init_incs( knum )
[9297]	289	!!------------------------------------------------------------------------
[9296]	290	!! * ROUTINE asm_bgc_init_incs *
	291	!!
	292	!! ** Purpose : initialise bgc increments
	293	!!
	294	!! ** Method : allocate arrays and read increments from file
	295	!!
	296	!! ** Action : allocate arrays and read increments from file
	297	!!
	298	!! References : asm_inc_init
[9297]	299	!!------------------------------------------------------------------------
[9296]	300	!!
	301	INTEGER, INTENT(in ) :: knum ! i/o unit of increments file
	302	!!
[9297]	303	!!------------------------------------------------------------------------
[9296]	304
[9322]	305	! Allocate and read increments
	306
[9296]	307	IF ( ln_slchltotinc ) THEN
	308	ALLOCATE( slchltot_bkginc(jpi,jpj) )
[9322]	309	CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc )
	310	ENDIF
	311
	312	IF ( ln_slchldiainc ) THEN
	313	ALLOCATE( slchldia_bkginc(jpi,jpj) )
	314	CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc )
	315	ENDIF
	316
	317	IF ( ln_slchlnoninc ) THEN
	318	ALLOCATE( slchlnon_bkginc(jpi,jpj) )
	319	CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc )
	320	ENDIF
	321
	322	IF ( ln_schltotinc ) THEN
	323	ALLOCATE( schltot_bkginc(jpi,jpj) )
	324	CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc )
	325	ENDIF
	326
	327	IF ( ln_slphytotinc ) THEN
	328	ALLOCATE( slphytot_bkginc(jpi,jpj) )
	329	CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc )
	330	ENDIF
	331
	332	IF ( ln_slphydiainc ) THEN
	333	ALLOCATE( slphydia_bkginc(jpi,jpj) )
	334	CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc )
	335	ENDIF
	336
	337	IF ( ln_slphynoninc ) THEN
	338	ALLOCATE( slphynon_bkginc(jpi,jpj) )
	339	CALL asm_bgc_read_incs_2d( knum, 'bckinslphynon', slphynon_bkginc )
	340	ENDIF
	341
	342	IF ( ln_sfco2inc ) THEN
	343	ALLOCATE( sfco2_bkginc(jpi,jpj) )
	344	CALL asm_bgc_read_incs_2d( knum, 'bckinsfco2', sfco2_bkginc )
	345	ENDIF
	346
	347	IF ( ln_spco2inc ) THEN
	348	ALLOCATE( spco2_bkginc(jpi,jpj) )
	349	CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc )
	350	ENDIF
	351
	352	IF ( ln_plchltotinc ) THEN
	353	ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) )
	354	CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc )
	355	ENDIF
	356
	357	IF ( ln_pchltotinc ) THEN
	358	ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) )
	359	CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc )
	360	ENDIF
	361
	362	IF ( ln_pno3inc ) THEN
	363	ALLOCATE( pno3_bkginc(jpi,jpj,jpk) )
	364	CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc )
	365	ENDIF
	366
	367	IF ( ln_psi4inc ) THEN
	368	ALLOCATE( psi4_bkginc(jpi,jpj,jpk) )
	369	CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc )
	370	ENDIF
	371
	372	IF ( ln_pdicinc ) THEN
	373	ALLOCATE( pdic_bkginc(jpi,jpj,jpk) )
	374	CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc )
	375	ENDIF
	376
	377	IF ( ln_palkinc ) THEN
	378	ALLOCATE( palk_bkginc(jpi,jpj,jpk) )
	379	CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc )
	380	ENDIF
	381
	382	IF ( ln_pphinc ) THEN
	383	ALLOCATE( pph_bkginc(jpi,jpj,jpk) )
	384	CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc )
	385	ENDIF
	386
	387	IF ( ln_po2inc ) THEN
	388	ALLOCATE( po2_bkginc(jpi,jpj,jpk) )
	389	CALL asm_bgc_read_incs_3d( knum, 'bckinpo2', po2_bkginc )
	390	ENDIF
	391
	392	! Allocate balancing increments
	393
	394	IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
	395	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
	396	& ln_slphynoninc ) THEN
[9296]	397	#if defined key_top
[9326]	398	ALLOCATE( phyto2d_balinc(jpi,jpj,jpk,jptra) )
	399	phyto2d_balinc(:,:,:,:) = 0.0
[9322]	400	#else
	401	CALL ctl_stop( ' key_top must be set for balancing increments' )
[9296]	402	#endif
	403	ENDIF
[9326]	404
	405	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
	406	#if defined key_top
	407	ALLOCATE( phyto3d_balinc(jpi,jpj,jpk,jptra) )
	408	phyto3d_balinc(:,:,:,:) = 0.0
	409	#else
	410	CALL ctl_stop( ' key_top must be set for balancing increments' )
	411	#endif
	412	ENDIF
[9322]	413
[9296]	414	IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
	415	#if defined key_top
[9322]	416	ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) )
	417	pco2_balinc(:,:,:,:) = 0.0
	418	#else
	419	CALL ctl_stop( ' key_top must be set for balancing increments' )
[9296]	420	#endif
	421	ENDIF
	422
[9322]	423	IF ( ln_pphinc ) THEN
	424	#if defined key_top
	425	ALLOCATE( ph_balinc(jpi,jpj,jpk,jptra) )
	426	ph_balinc(:,:,:,:) = 0.0
	427	#else
	428	CALL ctl_stop( ' key_top must be set for balancing increments' )
	429	#endif
[9296]	430	ENDIF
	431
	432	END SUBROUTINE asm_bgc_init_incs
	433
[9297]	434	!!===========================================================================
	435	!!===========================================================================
	436	!!===========================================================================
[9296]	437
[9322]	438	SUBROUTINE asm_bgc_read_incs_2d( knum, cd_bgcname, p_incs )
	439	!!------------------------------------------------------------------------
	440	!! * ROUTINE asm_bgc_init_incs *
	441	!!
	442	!! ** Purpose : read 2d bgc increments
	443	!!
	444	!! ** Method : read increments from file
	445	!!
	446	!! ** Action : read increments from file
	447	!!
	448	!! References : asm_inc_init
	449	!!------------------------------------------------------------------------
	450	!!
	451	INTEGER, INTENT(in ) :: knum ! i/o unit
	452	CHARACTER(LEN=13), INTENT(in ) :: cd_bgcname ! variable
	453	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: p_incs ! increments
	454	!!
	455	!!------------------------------------------------------------------------
	456
	457	! Initialise
	458	p_incs(:,:) = 0.0
	459
	460	! read from file
	461	CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 )
	462
	463	! Apply the masks
	464	p_incs(:,:) = p_incs(:,:) * tmask(:,:,1)
	465
	466	! Set missing increments to 0.0 rather than 1e+20
	467	! to allow for differences in masks
	468	WHERE( ABS( p_incs(:,:) ) > 1.0e+10 ) p_incs(:,:) = 0.0
	469
	470	END SUBROUTINE asm_bgc_read_incs_2d
	471
	472	!!===========================================================================
	473	!!===========================================================================
	474	!!===========================================================================
	475
	476	SUBROUTINE asm_bgc_read_incs_3d( knum, cd_bgcname, p_incs )
	477	!!------------------------------------------------------------------------
	478	!! * ROUTINE asm_bgc_init_incs *
	479	!!
	480	!! ** Purpose : read 3d bgc increments
	481	!!
	482	!! ** Method : read increments from file
	483	!!
	484	!! ** Action : read increments from file
	485	!!
	486	!! References : asm_inc_init
	487	!!------------------------------------------------------------------------
	488	!!
	489	INTEGER, INTENT(in ) :: knum ! i/o unit
	490	CHARACTER(LEN=13), INTENT(in ) :: cd_bgcname ! variable
	491	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: p_incs ! increments
	492	!!
	493	!!------------------------------------------------------------------------
	494
	495	! Initialise
	496	p_incs(:,:,:) = 0.0
	497
	498	! read from file
	499	CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 )
	500
	501	! Apply the masks
	502	p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:)
	503
	504	! Set missing increments to 0.0 rather than 1e+20
	505	! to allow for differences in masks
	506	WHERE( ABS( p_incs(:,:,:) ) > 1.0e+10 ) p_incs(:,:,:) = 0.0
	507
	508	END SUBROUTINE asm_bgc_read_incs_3d
	509
	510	!!===========================================================================
	511	!!===========================================================================
	512	!!===========================================================================
	513
[9296]	514	SUBROUTINE asm_bgc_init_bkg
[9297]	515	!!------------------------------------------------------------------------
[9296]	516	!! * ROUTINE asm_bgc_init_bkg *
	517	!!
	518	!! ** Purpose : initialise bgc background
	519	!!
	520	!! ** Method : allocate arrays and read background from file
	521	!!
	522	!! ** Action : allocate arrays and read background from file
	523	!!
	524	!! References : asm_inc_init
[9297]	525	!!------------------------------------------------------------------------
[9296]	526	!!
	527	INTEGER :: inum ! i/o unit of background file
	528	INTEGER :: jt ! loop counter
	529	!!
[9297]	530	!!------------------------------------------------------------------------
[9296]	531
[10141]	532	#if defined key_top && ( defined key_hadocc \|\| defined key_medusa )
[9322]	533	IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
	534	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
	535	& ln_slphynoninc .OR. ln_plchltotinc .OR. ln_pchltotinc ) THEN
[9296]	536
	537	ALLOCATE( pgrow_avg_bkg(jpi,jpj) )
	538	ALLOCATE( ploss_avg_bkg(jpi,jpj) )
	539	ALLOCATE( phyt_avg_bkg(jpi,jpj) )
[13097]	540	ALLOCATE( pgrow_avg_3d_bkg(jpi,jpj,jpk) )
	541	ALLOCATE( ploss_avg_3d_bkg(jpi,jpj,jpk) )
	542	ALLOCATE( phyt_avg_3d_bkg(jpi,jpj,jpk) )
[9296]	543	ALLOCATE( mld_max_bkg(jpi,jpj) )
	544	ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
	545	pgrow_avg_bkg(:,:) = 0.0
	546	ploss_avg_bkg(:,:) = 0.0
	547	phyt_avg_bkg(:,:) = 0.0
[13097]	548	pgrow_avg_3d_bkg(:,:,:) = 0.0
	549	ploss_avg_3d_bkg(:,:,:) = 0.0
	550	phyt_avg_3d_bkg(:,:,:) = 0.0
[9296]	551	mld_max_bkg(:,:) = 0.0
	552	tracer_bkg(:,:,:,:) = 0.0
	553
	554	#if defined key_hadocc
[9322]	555	ALLOCATE( chl_bkg(jpi,jpj,jpk) )
	556	ALLOCATE( cchl_p_bkg(jpi,jpj,jpk) )
	557	chl_bkg(:,:,:) = 0.0
	558	cchl_p_bkg(:,:,:) = 0.0
[9296]	559	#endif
	560
	561	!--------------------------------------------------------------------
[9322]	562	! Read background variables for phytoplankton assimilation
[9296]	563	! Some only required if performing balancing
	564	!--------------------------------------------------------------------
	565
	566	CALL iom_open( c_asmbkg, inum )
	567
	568	#if defined key_hadocc
	569	CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl', chl_bkg )
	570	CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg )
[9322]	571	chl_bkg(:,:,:) = chl_bkg(:,:,:) * tmask(:,:,:)
	572	cchl_p_bkg(:,:,:) = cchl_p_bkg(:,:,:) * tmask(:,:,:)
[10141]	573	#elif defined key_medusa
[9296]	574	CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) )
	575	CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) )
[9861]	576	CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) )
	577	CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
	578	CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
[9296]	579	#endif
	580
[9322]	581	IF ( ln_phytobal ) THEN
[9296]	582
	583	CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg )
	584	CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg )
	585	CALL iom_get( inum, jpdom_autoglo, 'phyt_avg', phyt_avg_bkg )
[13097]	586	CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg_3d', pgrow_avg_3d_bkg )
	587	CALL iom_get( inum, jpdom_autoglo, 'ploss_avg_3d', ploss_avg_3d_bkg )
	588	CALL iom_get( inum, jpdom_autoglo, 'phyt_avg_3d', phyt_avg_3d_bkg )
[9296]	589	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
	590	pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1)
	591	ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1)
	592	phyt_avg_bkg(:,:) = phyt_avg_bkg(:,:) * tmask(:,:,1)
[13097]	593	pgrow_avg_3d_bkg(:,:,:) = pgrow_avg_3d_bkg(:,:,:) * tmask(:,:,:)
	594	ploss_avg_3d_bkg(:,:,:) = ploss_avg_3d_bkg(:,:,:) * tmask(:,:,:)
	595	phyt_avg_3d_bkg(:,:,:) = phyt_avg_3d_bkg(:,:,:) * tmask(:,:,:)
[9296]	596	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
	597
	598	#if defined key_hadocc
	599	CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) )
	600	CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) )
	601	CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) )
	602	CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) )
	603	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
	604	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
[10141]	605	#elif defined key_medusa
[9296]	606	CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) )
	607	CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) )
	608	CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) )
	609	CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) )
	610	CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) )
	611	CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) )
	612	CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) )
	613	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
	614	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
	615	CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
	616	#endif
[9862]	617	ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
[9296]	618	#if defined key_hadocc
	619	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
	620	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
[10141]	621	#elif defined key_medusa
[9296]	622	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
	623	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
	624	#endif
	625	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
	626	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
	627	ENDIF
	628
	629	CALL iom_close( inum )
	630
	631	DO jt = 1, jptra
	632	tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
	633	END DO
	634
[9322]	635	ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
[9296]	636
	637	ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
	638	ALLOCATE( mld_max_bkg(jpi,jpj) )
	639	tracer_bkg(:,:,:,:) = 0.0
	640	mld_max_bkg(:,:) = 0.0
	641
	642	CALL iom_open( c_asmbkg, inum )
	643
	644	#if defined key_hadocc
	645	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
	646	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
[10141]	647	#elif defined key_medusa
[9296]	648	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
	649	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
	650	#endif
	651	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
	652
	653	CALL iom_close( inum )
	654
	655	DO jt = 1, jptra
	656	tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
	657	END DO
	658	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
	659
	660	ENDIF
[9322]	661	#else
	662	CALL ctl_stop( ' asm_bgc_init_bkg: key_top and a compatible biogeochemical model required' )
[9296]	663	#endif
	664
	665	END SUBROUTINE asm_bgc_init_bkg
	666
[9297]	667	!!===========================================================================
	668	!!===========================================================================
	669	!!===========================================================================
[9296]	670
	671	SUBROUTINE asm_bgc_bal_wri( kt )
[9297]	672	!!------------------------------------------------------------------------
[9296]	673	!!
	674	!! * ROUTINE asm_bgc_bal_wri *
	675	!!
	676	!! ** Purpose : Write to file the balancing increments
	677	!! calculated for biogeochemistry
	678	!!
	679	!! ** Method : Write to file the balancing increments
	680	!! calculated for biogeochemistry
	681	!!
	682	!! ** Action :
	683	!!
	684	!! References :
	685	!!
	686	!! History :
	687	!! ! 2014-08 (D. Ford) Initial version, based on asm_bkg_wri
[9297]	688	!!------------------------------------------------------------------------
[9296]	689	!! * Arguments
	690	INTEGER, INTENT( IN ) :: kt ! Current time-step
	691	!! * Local declarations
	692	CHARACTER(LEN=50) :: cl_asmbal ! Filename (with extension)
	693	LOGICAL :: llok ! Check if file exists
	694	INTEGER :: inum ! File unit number
	695	REAL(wp) :: zdate ! Date
[9297]	696	!!------------------------------------------------------------------------
[9296]	697
	698	! Set things up
	699	zdate = REAL( ndastp )
	700	WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal )
	701
	702	! Check if file exists
	703	INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok )
	704	IF( .NOT. llok ) THEN
	705	IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// &
	706	& TRIM( c_asmbal ) // ' at timestep = ', kt
	707
	708	! Define the output file
	709	CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
	710
	711	! Write the information
	712	CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate )
	713
[9381]	714	IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
	715	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
	716	& ln_slphynoninc ) THEN
[10141]	717	#if defined key_medusa
[9381]	718	CALL iom_rstput( kt, kt, inum, 'phy2d_chn', phyto2d_balinc(:,:,:,jpchn) )
	719	CALL iom_rstput( kt, kt, inum, 'phy2d_chd', phyto2d_balinc(:,:,:,jpchd) )
	720	CALL iom_rstput( kt, kt, inum, 'phy2d_phn', phyto2d_balinc(:,:,:,jpphn) )
	721	CALL iom_rstput( kt, kt, inum, 'phy2d_phd', phyto2d_balinc(:,:,:,jpphd) )
	722	CALL iom_rstput( kt, kt, inum, 'phy2d_pds', phyto2d_balinc(:,:,:,jppds) )
[9322]	723	IF ( ln_phytobal ) THEN
[9381]	724	CALL iom_rstput( kt, kt, inum, 'phy2d_zmi', phyto2d_balinc(:,:,:,jpzmi) )
	725	CALL iom_rstput( kt, kt, inum, 'phy2d_zme', phyto2d_balinc(:,:,:,jpzme) )
	726	CALL iom_rstput( kt, kt, inum, 'phy2d_din', phyto2d_balinc(:,:,:,jpdin) )
	727	CALL iom_rstput( kt, kt, inum, 'phy2d_sil', phyto2d_balinc(:,:,:,jpsil) )
	728	CALL iom_rstput( kt, kt, inum, 'phy2d_fer', phyto2d_balinc(:,:,:,jpfer) )
	729	CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jpdet) )
	730	CALL iom_rstput( kt, kt, inum, 'phy2d_dtc', phyto2d_balinc(:,:,:,jpdtc) )
	731	CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jpdic) )
	732	CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jpalk) )
	733	CALL iom_rstput( kt, kt, inum, 'phy2d_oxy', phyto2d_balinc(:,:,:,jpoxy) )
[9296]	734	ENDIF
	735	#elif defined key_hadocc
[9381]	736	CALL iom_rstput( kt, kt, inum, 'phy2d_phy', phyto2d_balinc(:,:,:,jp_had_phy) )
[9322]	737	IF ( ln_phytobal ) THEN
[9381]	738	CALL iom_rstput( kt, kt, inum, 'phy2d_nut', phyto2d_balinc(:,:,:,jp_had_nut) )
	739	CALL iom_rstput( kt, kt, inum, 'phy2d_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) )
	740	CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jp_had_pdn) )
	741	CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jp_had_dic) )
	742	CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) )
[9296]	743	ENDIF
	744	#endif
	745	ENDIF
	746
[9381]	747	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
[10141]	748	#if defined key_medusa
[9381]	749	CALL iom_rstput( kt, kt, inum, 'phy3d_chn', phyto3d_balinc(:,:,:,jpchn) )
	750	CALL iom_rstput( kt, kt, inum, 'phy3d_chd', phyto3d_balinc(:,:,:,jpchd) )
	751	CALL iom_rstput( kt, kt, inum, 'phy3d_phn', phyto3d_balinc(:,:,:,jpphn) )
	752	CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) )
	753	CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) )
[13097]	754	IF ( ln_phytobal ) THEN
	755	CALL iom_rstput( kt, kt, inum, 'phy3d_zmi', phyto3d_balinc(:,:,:,jpzmi) )
	756	CALL iom_rstput( kt, kt, inum, 'phy3d_zme', phyto3d_balinc(:,:,:,jpzme) )
	757	CALL iom_rstput( kt, kt, inum, 'phy3d_din', phyto3d_balinc(:,:,:,jpdin) )
	758	CALL iom_rstput( kt, kt, inum, 'phy3d_sil', phyto3d_balinc(:,:,:,jpsil) )
	759	CALL iom_rstput( kt, kt, inum, 'phy3d_fer', phyto3d_balinc(:,:,:,jpfer) )
	760	CALL iom_rstput( kt, kt, inum, 'phy3d_det', phyto3d_balinc(:,:,:,jpdet) )
	761	CALL iom_rstput( kt, kt, inum, 'phy3d_dtc', phyto3d_balinc(:,:,:,jpdtc) )
	762	CALL iom_rstput( kt, kt, inum, 'phy3d_dic', phyto3d_balinc(:,:,:,jpdic) )
	763	CALL iom_rstput( kt, kt, inum, 'phy3d_alk', phyto3d_balinc(:,:,:,jpalk) )
	764	CALL iom_rstput( kt, kt, inum, 'phy3d_oxy', phyto3d_balinc(:,:,:,jpoxy) )
	765	ENDIF
[9381]	766	#elif defined key_hadocc
	767	CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) )
[13097]	768	IF ( ln_phytobal ) THEN
	769	CALL iom_rstput( kt, kt, inum, 'phy3d_nut', phyto3d_balinc(:,:,:,jp_had_nut) )
	770	CALL iom_rstput( kt, kt, inum, 'phy3d_zoo', phyto3d_balinc(:,:,:,jp_had_zoo) )
	771	CALL iom_rstput( kt, kt, inum, 'phy3d_det', phyto3d_balinc(:,:,:,jp_had_pdn) )
	772	CALL iom_rstput( kt, kt, inum, 'phy3d_dic', phyto3d_balinc(:,:,:,jp_had_dic) )
	773	CALL iom_rstput( kt, kt, inum, 'phy3d_alk', phyto3d_balinc(:,:,:,jp_had_alk) )
	774	ENDIF
[9381]	775	#endif
	776	ENDIF
	777
[9296]	778	IF ( ln_spco2inc ) THEN
[10141]	779	#if defined key_medusa
[9381]	780	CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jpdic) )
	781	CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jpalk) )
[9296]	782	#elif defined key_hadocc
[9381]	783	CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
	784	CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
[9296]	785	#endif
	786	ELSE IF ( ln_sfco2inc ) THEN
[10141]	787	#if defined key_medusa
[9381]	788	CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jpdic) )
	789	CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jpalk) )
[9296]	790	#elif defined key_hadocc
[9381]	791	CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
	792	CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
[9296]	793	#endif
	794	ENDIF
	795
[9381]	796	IF ( ln_pphinc ) THEN
[10141]	797	#if defined key_medusa
[9381]	798	CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jpdic) )
	799	CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jpalk) )
	800	#elif defined key_hadocc
	801	CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) )
	802	CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) )
	803	#endif
	804	ENDIF
	805
[9296]	806	CALL iom_close( inum )
	807	ELSE
	808	CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', &
	809	& ' Therefore not writing out balancing increments at this timestep', &
	810	& ' Something has probably gone wrong somewhere' )
	811	IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// &
	812	& TRIM( c_asmbal ) // ' at timestep = ', kt
	813	ENDIF
	814
	815	END SUBROUTINE asm_bgc_bal_wri
	816
[9297]	817	!!===========================================================================
	818	!!===========================================================================
	819	!!===========================================================================
[9296]	820
	821	SUBROUTINE asm_bgc_bkg_wri( kt, knum )
[9297]	822	!!------------------------------------------------------------------------
[9296]	823	!! * ROUTINE asm_bgc_bkg_wri *
	824	!!
	825	!! ** Purpose : write out bgc background
	826	!!
	827	!! ** Method : write out bgc background
	828	!!
	829	!! ** Action : write out bgc background
	830	!!
	831	!! References : asm_bkg_wri
[9297]	832	!!------------------------------------------------------------------------
[9296]	833	!!
	834	INTEGER, INTENT(in ) :: kt ! Current time-step
	835	INTEGER, INTENT(in ) :: knum ! i/o unit of increments file
	836	!!
[9297]	837	!!------------------------------------------------------------------------
[9296]	838
	839	CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg' , pgrow_avg )
	840	CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg' , ploss_avg )
	841	CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg' , phyt_avg )
[13097]	842	CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg_3d', pgrow_avg_3d )
	843	CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg_3d', ploss_avg_3d )
	844	CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg_3d' , phyt_avg_3d )
[9296]	845	CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max' , mld_max )
[13097]	846	#if defined key_hadocc
[9296]	847	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_nut' , trn(:,:,:,jp_had_nut) )
	848	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_phy' , trn(:,:,:,jp_had_phy) )
	849	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_zoo' , trn(:,:,:,jp_had_zoo) )
	850	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_pdn' , trn(:,:,:,jp_had_pdn) )
	851	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_dic' , trn(:,:,:,jp_had_dic) )
	852	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_alk' , trn(:,:,:,jp_had_alk) )
[9322]	853	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_chl' , HADOCC_CHL(:,:,:) )
	854	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:) )
[10141]	855	#elif defined key_medusa
[9296]	856	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn' , trn(:,:,:,jpchn) )
	857	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chd' , trn(:,:,:,jpchd) )
	858	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phn' , trn(:,:,:,jpphn) )
	859	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phd' , trn(:,:,:,jpphd) )
	860	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_pds' , trn(:,:,:,jppds) )
	861	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zmi' , trn(:,:,:,jpzmi) )
	862	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zme' , trn(:,:,:,jpzme) )
	863	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_din' , trn(:,:,:,jpdin) )
	864	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_sil' , trn(:,:,:,jpsil) )
	865	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_fer' , trn(:,:,:,jpfer) )
	866	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_det' , trn(:,:,:,jpdet) )
	867	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dtc' , trn(:,:,:,jpdtc) )
	868	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dic' , trn(:,:,:,jpdic) )
	869	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_alk' , trn(:,:,:,jpalk) )
	870	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_oxy' , trn(:,:,:,jpoxy) )
	871	#endif
	872
	873	END SUBROUTINE asm_bgc_bkg_wri
	874
[9297]	875	!!===========================================================================
	876	!!===========================================================================
	877	!!===========================================================================
[9296]	878
[9431]	879	SUBROUTINE asm_bgc_unlog_2d( pbkg, pinc_log, pinc_nonlog )
	880	!!------------------------------------------------------------------------
	881	!! * ROUTINE asm_bgc_init_incs *
	882	!!
	883	!! ** Purpose : convert log increments to non-log
	884	!!
	885	!! ** Method : need to account for model background,
	886	!! cannot simply do 10^log_inc. Need to:
	887	!! 1) Add log_inc to log10(background) to get log10(analysis)
	888	!! 2) Take 10^log10(analysis) to get analysis
	889	!! 3) Subtract background from analysis to get non-log incs
	890	!!
	891	!! ** Action : return non-log increments
	892	!!
	893	!! References :
	894	!!------------------------------------------------------------------------
	895	!!
	896	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pbkg ! Background
	897	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pinc_log ! Log incs
	898	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: pinc_nonlog ! Non-log incs
	899	!
	900	INTEGER :: ji, jj ! Loop counters
	901	!!
	902	!!------------------------------------------------------------------------
	903
	904	DO jj = 1, jpj
	905	DO ji = 1, jpi
	906	IF ( pbkg(ji,jj) > 0.0 ) THEN
	907	pinc_nonlog(ji,jj) = 10**( LOG10( pbkg(ji,jj) ) + &
	908	& pinc_log(ji,jj) ) &
	909	& - pbkg(ji,jj)
	910	ELSE
	911	pinc_nonlog(ji,jj) = 0.0
	912	ENDIF
	913	END DO
	914	END DO
	915
	916	END SUBROUTINE asm_bgc_unlog_2d
	917
	918	!!===========================================================================
	919	!!===========================================================================
	920	!!===========================================================================
	921
[13097]	922	SUBROUTINE asm_bgc_unlog_3d( pbkg, pinc_log, pinc_nonlog )
	923	!!------------------------------------------------------------------------
	924	!! * ROUTINE asm_bgc_init_incs *
	925	!!
	926	!! ** Purpose : convert log increments to non-log
	927	!!
	928	!! ** Method : need to account for model background,
	929	!! cannot simply do 10^log_inc. Need to:
	930	!! 1) Add log_inc to log10(background) to get log10(analysis)
	931	!! 2) Take 10^log10(analysis) to get analysis
	932	!! 3) Subtract background from analysis to get non-log incs
	933	!!
	934	!! ** Action : return non-log increments
	935	!!
	936	!! References :
	937	!!------------------------------------------------------------------------
	938	!!
	939	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pbkg ! Background
	940	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pinc_log ! Log incs
	941	REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpk) :: pinc_nonlog ! Non-log incs
	942	!
	943	INTEGER :: ji, jj, jk ! Loop counters
	944	!!
	945	!!------------------------------------------------------------------------
	946
	947	DO jk = 1, jpk
	948	DO jj = 1, jpj
	949	DO ji = 1, jpi
	950	IF ( pbkg(ji,jj,jk) > 0.0 ) THEN
	951	pinc_nonlog(ji,jj,jk) = 10**( LOG10( pbkg(ji,jj,jk) ) + &
	952	& pinc_log(ji,jj,jk) ) &
	953	& - pbkg(ji,jj,jk)
	954	ELSE
	955	pinc_nonlog(ji,jj,jk) = 0.0
	956	ENDIF
	957	END DO
	958	END DO
	959	END DO
	960
	961	END SUBROUTINE asm_bgc_unlog_3d
	962
	963	!!===========================================================================
	964	!!===========================================================================
	965	!!===========================================================================
	966
[9326]	967	SUBROUTINE phyto2d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
[9297]	968	!!------------------------------------------------------------------------
[9326]	969	!! * ROUTINE phyto2d_asm_inc *
[9296]	970	!!
	971	!! ** Purpose : Apply the chlorophyll assimilation increments.
	972	!!
	973	!! ** Method : Calculate increments to state variables using nitrogen
	974	!! balancing.
	975	!! Direct initialization or Incremental Analysis Updating.
	976	!!
	977	!! ** Action :
[9297]	978	!!------------------------------------------------------------------------
[9296]	979	INTEGER, INTENT(IN) :: kt ! Current time step
	980	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
	981	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
	982	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
	983	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
	984	!
[13097]	985	INTEGER :: jk ! Loop counter
	986	INTEGER :: it ! Index
	987	REAL(wp) :: zincwgt ! IAU weight for current time step
	988	REAL(wp) :: zincper ! IAU interval in seconds
	989	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth
	990	REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_chltot ! Local chltot incs
	991	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chltot ! Local chltot bkg
	992	REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_phytot ! Local phytot incs
	993	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phytot ! Local phytot bkg
[10141]	994	#if defined key_medusa
[13097]	995	REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_chldia ! Local chldia incs
	996	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
	997	REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_chlnon ! Local chlnon incs
	998	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnon ! Local chlnon bkg
	999	REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_phydia ! Local phydia incs
	1000	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phydia ! Local phydia bkg
	1001	REAL(wp), DIMENSION(jpi,jpj,1) :: zinc_phynon ! Local phynon incs
	1002	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg
[9431]	1003	#endif
[13097]	1004	REAL(wp), DIMENSION(jpi,jpj,1) :: zpgrow_avg_bkg ! Local pgrow_avg_bkg
	1005	REAL(wp), DIMENSION(jpi,jpj,1) :: zploss_avg_bkg ! Local ploss_avg_bkg
	1006	REAL(wp), DIMENSION(jpi,jpj,1) :: zphyt_avg_bkg ! Local phyt_avg_bkg
[9297]	1007	!!------------------------------------------------------------------------
[9322]	1008
[9296]	1009	IF ( kt <= nit000 ) THEN
	1010
[9431]	1011	! Un-log any log increments for passing to balancing routines
[9435]	1012	! Remember that two sets of non-log increments should not be
	1013	! expected to be in the same ratio as their log equivalents
	1014
[9431]	1015	! Total chlorophyll
	1016	IF ( ln_slchltotinc ) THEN
[10141]	1017	#if defined key_medusa
[9431]	1018	zbkg_chltot(:,:) = tracer_bkg(:,:,1,jpchn) + tracer_bkg(:,:,1,jpchd)
	1019	#elif defined key_hadocc
	1020	zbkg_chltot(:,:) = chl_bkg(:,:,1)
	1021	#endif
[13097]	1022	CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot(:,:,1) )
[9431]	1023	ELSE IF ( ln_schltotinc ) THEN
[13097]	1024	zinc_chltot(:,:,1) = schltot_bkginc(:,:)
[9431]	1025	ELSE
[13097]	1026	zinc_chltot(:,:,:) = 0.0
[9431]	1027	ENDIF
	1028
[10141]	1029	#if defined key_medusa
[9431]	1030	! Diatom chlorophyll
	1031	IF ( ln_slchldiainc ) THEN
	1032	zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd)
[13097]	1033	CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia(:,:,1) )
[9431]	1034	ELSE
[13097]	1035	zinc_chldia(:,:,:) = 0.0
[9431]	1036	ENDIF
	1037	#endif
	1038
[10141]	1039	#if defined key_medusa
[9431]	1040	! Non-diatom chlorophyll
	1041	IF ( ln_slchlnoninc ) THEN
	1042	zbkg_chlnon(:,:) = tracer_bkg(:,:,1,jpchn)
[13097]	1043	CALL asm_bgc_unlog_2d( zbkg_chlnon, slchlnon_bkginc, zinc_chlnon(:,:,1) )
[9431]	1044	ELSE
[13097]	1045	zinc_chlnon(:,:,:) = 0.0
[9431]	1046	ENDIF
	1047	#endif
	1048
	1049	! Total phytoplankton carbon
	1050	IF ( ln_slphytotinc ) THEN
[10141]	1051	#if defined key_medusa
[9431]	1052	zbkg_phytot(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
	1053	#elif defined key_hadocc
	1054	zbkg_phytot(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
	1055	#endif
[13097]	1056	CALL asm_bgc_unlog_2d( zbkg_phytot, slphytot_bkginc, zinc_phytot(:,:,1) )
[9431]	1057	ELSE
[13097]	1058	zinc_phytot(:,:,:) = 0.0
[9431]	1059	ENDIF
	1060
[10141]	1061	#if defined key_medusa
[9431]	1062	! Diatom phytoplankton carbon
	1063	IF ( ln_slphydiainc ) THEN
	1064	zbkg_phydia(:,:) = trn(:,:,1,jpphd) * xthetapd
[13097]	1065	CALL asm_bgc_unlog_2d( zbkg_phydia, slphydia_bkginc, zinc_phydia(:,:,1) )
[9431]	1066	ELSE
[13097]	1067	zinc_phydia(:,:,:) = 0.0
[9431]	1068	ENDIF
	1069	#endif
	1070
[10141]	1071	#if defined key_medusa
[9431]	1072	! Non-diatom phytoplankton carbon
	1073	IF ( ln_slphynoninc ) THEN
	1074	zbkg_phynon(:,:) = trn(:,:,1,jpphn) * xthetapn
[13097]	1075	CALL asm_bgc_unlog_2d( zbkg_phynon, slphynon_bkginc, zinc_phynon(:,:,1) )
[9431]	1076	ELSE
[13097]	1077	zinc_phynon(:,:,:) = 0.0
[9431]	1078	ENDIF
	1079	#endif
	1080
	1081	! Select mixed layer
	1082	IF ( ll_asmdin ) THEN
[10141]	1083	#if defined key_top && ( defined key_hadocc \|\| defined key_medusa )
[9431]	1084	CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', &
	1085	& ' Mixed layer depth taken to be background maximum mld_max_bkg' )
	1086	zmld(:,:) = mld_max_bkg(:,:)
[10057]	1087	#else
	1088	CALL ctl_stop( ' Should not be doing biogeochemical assimilation without key_top' )
	1089	#endif
[9431]	1090	ELSE
	1091	SELECT CASE( mld_choice_bgc )
	1092	CASE ( 1 ) ! Turbocline/mixing depth [W points]
	1093	zmld(:,:) = hmld(:,:)
	1094	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
	1095	zmld(:,:) = hmlp(:,:)
	1096	CASE ( 3 ) ! Kara MLD [Interpolated]
	1097	#if defined key_karaml
	1098	IF ( ln_kara ) THEN
	1099	zmld(:,:) = hmld_kara(:,:)
	1100	ELSE
	1101	CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
	1102	& ' but ln_kara=.false.' )
	1103	ENDIF
	1104	#else
	1105	CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
	1106	& ' but is not defined' )
	1107	#endif
	1108	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
	1109	!zmld(:,:) = hmld_tref(:,:)
[9432]	1110	CALL ctl_stop( ' hmld_tref mixed layer requested for phyto2d assimilation,', &
[9431]	1111	& ' but is not available in this version' )
	1112	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
	1113	zmld(:,:) = hmlpt(:,:)
	1114	END SELECT
	1115	ENDIF
	1116
[9296]	1117	zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
[13097]	1118
	1119	zpgrow_avg_bkg(:,:,1) = pgrow_avg_bkg(:,:)
	1120	zploss_avg_bkg(:,:,1) = ploss_avg_bkg(:,:)
	1121	zphyt_avg_bkg(:,:,1) = phyt_avg_bkg(:,:)
[9296]	1122
[10141]	1123	#if defined key_medusa
[13097]	1124	CALL asm_phyto_bal_medusa( 1, &
	1125	& (ln_slchltotinc .OR. ln_schltotinc), &
	1126	& zinc_chltot, &
	1127	& ln_slchldiainc, &
	1128	& zinc_chldia, &
	1129	& ln_slchlnoninc, &
	1130	& zinc_chlnon, &
	1131	& ln_slphytotinc, &
	1132	& zinc_phytot, &
	1133	& ln_slphydiainc, &
	1134	& zinc_phydia, &
	1135	& ln_slphynoninc, &
	1136	& zinc_phynon, &
	1137	& zincper, &
	1138	& rn_maxchlinc, ln_phytobal, zmld, &
	1139	& zpgrow_avg_bkg, zploss_avg_bkg, &
	1140	& zphyt_avg_bkg, mld_max_bkg, &
	1141	& tracer_bkg, phyto2d_balinc )
[9296]	1142	#elif defined key_hadocc
[13097]	1143	CALL asm_phyto_bal_hadocc( 1, &
	1144	& (ln_slchltotinc .OR. ln_schltotinc), &
	1145	& zinc_chltot, &
	1146	& ln_slphytotinc, &
	1147	& zinc_phytot, &
	1148	& zincper, &
	1149	& rn_maxchlinc, ln_phytobal, zmld, &
	1150	& zpgrow_avg_bkg, zploss_avg_bkg, &
	1151	& zphyt_avg_bkg, mld_max_bkg, &
	1152	& cchl_p_bkg(:,:,1), &
	1153	& tracer_bkg, phyto2d_balinc )
[9296]	1154	#else
[9431]	1155	CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', &
[9296]	1156	& 'but not defined a biogeochemical model' )
	1157	#endif
	1158
	1159	ENDIF
	1160
	1161	IF ( ll_asmiau ) THEN
	1162
	1163	!--------------------------------------------------------------------
	1164	! Incremental Analysis Updating
	1165	!--------------------------------------------------------------------
	1166
	1167	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
	1168
	1169	it = kt - nit000 + 1
	1170	zincwgt = pwgtiau(it) ! IAU weight for the current time step
	1171	! note this is not a tendency so should not be divided by rdt
	1172
	1173	IF(lwp) THEN
	1174	WRITE(numout,*)
[9431]	1175	WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', &
[9296]	1176	& kt,' with IAU weight = ', pwgtiau(it)
	1177	WRITE(numout,*) '~~~~~~~~~~~~'
	1178	ENDIF
	1179
	1180	! Update the biogeochemical variables
	1181	! Add directly to trn and trb, rather than to tra, because tra gets
	1182	! reset to zero at the start of trc_stp, called after this routine
[10141]	1183	#if defined key_medusa
[9382]	1184	WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1185	& trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
[9333]	1186	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1187	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1188	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
	1189	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1190	END WHERE
[9296]	1191	#elif defined key_hadocc
[9382]	1192	WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	1193	& trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
[9333]	1194	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	1195	& phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	1196	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
	1197	& phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	1198	END WHERE
[9296]	1199	#endif
	1200
	1201	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1202	! which is called at end of model run
	1203	ENDIF
	1204
	1205	ELSEIF ( ll_asmdin ) THEN
	1206
	1207	!--------------------------------------------------------------------
	1208	! Direct Initialization
	1209	!--------------------------------------------------------------------
	1210
	1211	IF ( kt == nitdin_r ) THEN
	1212
	1213	neuler = 0 ! Force Euler forward step
	1214
	1215	! Initialize the now fields with the background + increment
	1216	! Background currently is what the model is initialised with
[9333]	1217	CALL ctl_warn( ' Doing direct initialisation with phyto2d assimilation', &
[9296]	1218	& ' Background state is taken from model rather than background file' )
[10141]	1219	#if defined key_medusa
[9382]	1220	WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1221	& trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
[9333]	1222	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1223	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1)
	1224	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
	1225	END WHERE
[9296]	1226	#elif defined key_hadocc
[9382]	1227	WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	1228	& trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
[9333]	1229	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	1230	& phyto2d_balinc(:,:,:,jp_had0:jp_had1)
	1231	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
	1232	END WHERE
[9296]	1233	#endif
	1234
	1235	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1236	! which is called at end of model run
	1237	ENDIF
	1238	!
	1239	ENDIF
	1240	!
[9326]	1241	END SUBROUTINE phyto2d_asm_inc
[9296]	1242
[9297]	1243	!!===========================================================================
	1244	!!===========================================================================
	1245	!!===========================================================================
[9296]	1246
[9326]	1247	SUBROUTINE phyto3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
	1248	!!------------------------------------------------------------------------
	1249	!! * ROUTINE phyto3d_asm_inc *
	1250	!!
	1251	!! ** Purpose : Apply the profile chlorophyll assimilation increments.
	1252	!!
	1253	!! ** Method : Calculate increments to state variables.
	1254	!! Direct initialization or Incremental Analysis Updating.
	1255	!!
	1256	!! ** Action :
	1257	!!------------------------------------------------------------------------
	1258	INTEGER, INTENT(IN) :: kt ! Current time step
	1259	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
	1260	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
	1261	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
	1262	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
	1263	!
[9381]	1264	INTEGER :: ji, jj, jk ! Loop counters
	1265	INTEGER :: it ! Index
[13097]	1266	REAL(wp) :: zincper ! IAU interval in seconds
[9381]	1267	REAL(wp) :: zincwgt ! IAU weight for timestep
[13097]	1268	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zinc_chltot ! Chlorophyll increments
	1269	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbkg_chltot ! Chlorophyll background
	1270	REAL(wp), DIMENSION(jpi,jpj) :: zdummy_2d ! Dummy array for call
	1271	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdummy_3d ! Dummy array for call
[9326]	1272	!!------------------------------------------------------------------------
	1273
	1274	IF ( kt <= nit000 ) THEN
	1275
	1276	IF ( ln_plchltotinc ) THEN
[10141]	1277	#if defined key_medusa
[13097]	1278	zbkg_chltot(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd)
[9326]	1279	#elif defined key_hadocc
[13097]	1280	zbkg_chltot(:,:,:) = chl_bkg(:,:,:)
[9326]	1281	#endif
[13097]	1282	CALL asm_bgc_unlog_3d( zbkg_chltot, plchltot_bkginc, zinc_chltot )
[9326]	1283	ELSE IF ( ln_pchltotinc ) THEN
[13097]	1284	zinc_chltot(:,:,:) = pchltot_bkginc(:,:,:)
[9326]	1285	ENDIF
	1286
[13097]	1287	zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
	1288
[10141]	1289	#if defined key_medusa
[13097]	1290	CALL asm_phyto_bal_medusa( jpk, &
	1291	& (ln_plchltotinc .OR. ln_pchltotinc), &
	1292	& zinc_chltot, &
	1293	& .FALSE., &
	1294	& zdummy_3d, &
	1295	& .FALSE., &
	1296	& zdummy_3d, &
	1297	& .FALSE., &
	1298	& zdummy_3d, &
	1299	& .FALSE., &
	1300	& zdummy_3d, &
	1301	& .FALSE., &
	1302	& zdummy_3d, &
	1303	& zincper, &
	1304	& rn_maxchlinc, ln_phytobal, zdummy_2d, &
	1305	& pgrow_avg_3d_bkg, ploss_avg_3d_bkg, &
	1306	& phyt_avg_3d_bkg, mld_max_bkg, &
	1307	& tracer_bkg, phyto3d_balinc )
[9326]	1308	#elif defined key_hadocc
[13097]	1309	CALL asm_phyto_bal_hadocc( jpk, &
	1310	& (ln_plchltotinc .OR. ln_pchltotinc), &
	1311	& zinc_chltot, &
	1312	& .FALSE., &
	1313	& zdummy_3d, &
	1314	& zincper, &
	1315	& rn_maxchlinc, ln_phytobal, zdummy_2d, &
	1316	& pgrow_avg_3d_bkg, ploss_avg_3d_bkg, &
	1317	& phyt_avg_3d_bkg, mld_max_bkg, &
	1318	& cchl_p_bkg, &
	1319	& tracer_bkg, phyto3d_balinc )
[9326]	1320	#else
[13097]	1321	CALL ctl_stop( 'Attempting to assimilate phyto3d data, ', &
[9326]	1322	& 'but not defined a biogeochemical model' )
	1323	#endif
	1324
	1325	ENDIF
	1326
	1327	IF ( ll_asmiau ) THEN
	1328
	1329	!--------------------------------------------------------------------
	1330	! Incremental Analysis Updating
	1331	!--------------------------------------------------------------------
	1332
	1333	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
	1334
	1335	it = kt - nit000 + 1
	1336	zincwgt = pwgtiau(it) ! IAU weight for the current time step
	1337	! note this is not a tendency so should not be divided by rdt
	1338
	1339	IF(lwp) THEN
	1340	WRITE(numout,*)
	1341	WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', &
	1342	& kt,' with IAU weight = ', pwgtiau(it)
	1343	WRITE(numout,*) '~~~~~~~~~~~~'
	1344	ENDIF
	1345
	1346	! Update the biogeochemical variables
	1347	! Add directly to trn and trb, rather than to tra, because tra gets
	1348	! reset to zero at the start of trc_stp, called after this routine
[10141]	1349	#if defined key_medusa
[9382]	1350	WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1351	& trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
[9326]	1352	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1353	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1354	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
	1355	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1356	END WHERE
	1357	#elif defined key_hadocc
[9382]	1358	WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	1359	& trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
[9326]	1360	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	1361	& phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	1362	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
	1363	& phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	1364	END WHERE
	1365	#endif
	1366
	1367	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1368	! which is called at end of model run
	1369	ENDIF
	1370
	1371	ELSEIF ( ll_asmdin ) THEN
	1372
	1373	!--------------------------------------------------------------------
	1374	! Direct Initialization
	1375	!--------------------------------------------------------------------
	1376
	1377	IF ( kt == nitdin_r ) THEN
	1378
	1379	neuler = 0 ! Force Euler forward step
	1380
	1381	! Initialize the now fields with the background + increment
	1382	! Background currently is what the model is initialised with
	1383	CALL ctl_warn( ' Doing direct initialisation with phyto3d assimilation', &
	1384	& ' Background state is taken from model rather than background file' )
[10141]	1385	#if defined key_medusa
[9382]	1386	WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1387	& trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
[9326]	1388	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1389	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1)
	1390	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
	1391	END WHERE
	1392	#elif defined key_hadocc
[9382]	1393	WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	1394	& trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
[9326]	1395	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	1396	& phyto3d_balinc(:,:,:,jp_had0:jp_had1)
	1397	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
	1398	END WHERE
	1399	#endif
	1400
	1401	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1402	! which is called at end of model run
	1403	ENDIF
	1404	!
	1405	ENDIF
	1406	!
	1407	END SUBROUTINE phyto3d_asm_inc
	1408
	1409	!!===========================================================================
	1410	!!===========================================================================
	1411	!!===========================================================================
	1412
[9322]	1413	SUBROUTINE pco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
	1414	& ll_trainc, pt_bkginc, ps_bkginc )
[9297]	1415	!!------------------------------------------------------------------------
[9296]	1416	!! * ROUTINE pco2_asm_inc *
	1417	!!
	1418	!! ** Purpose : Apply the pco2/fco2 assimilation increments.
	1419	!!
	1420	!! ** Method : Calculate increments to state variables using carbon
	1421	!! balancing.
	1422	!! Direct initialization or Incremental Analysis Updating.
	1423	!!
	1424	!! ** Action :
[9297]	1425	!!------------------------------------------------------------------------
[9296]	1426	INTEGER, INTENT(IN) :: kt ! Current time step
	1427	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
	1428	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
	1429	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
	1430	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
	1431	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
	1432	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
	1433	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
	1434	!
	1435	INTEGER :: ji, jj, jk ! Loop counters
	1436	INTEGER :: it ! Index
	1437	INTEGER :: jkmax ! Index of mixed layer depth
	1438	REAL(wp) :: zincwgt ! IAU weight for current time step
	1439	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth for balancing
	1440	REAL(wp), DIMENSION(jpi,jpj) :: pco2_bkginc_temp ! For fCO2 conversion
	1441	REAL(wp), DIMENSION(jpi,jpj) :: dic_bkg_temp ! DIC background
	1442	REAL(wp), DIMENSION(jpi,jpj) :: alk_bkg_temp ! Alkalinity background
	1443	REAL(wp), DIMENSION(jpi,jpj) :: tem_bkg_temp ! Temperature background
	1444	REAL(wp), DIMENSION(jpi,jpj) :: sal_bkg_temp ! Salinity background
	1445	REAL(wp), DIMENSION(1) :: patm ! Atmospheric pressure
	1446	REAL(wp), DIMENSION(1) :: phyd ! Hydrostatic pressure
	1447
	1448	! Coefficients for fCO2 to pCO2 conversion
	1449	REAL(wp) :: zcoef_fco2_1 = -1636.75
	1450	REAL(wp) :: zcoef_fco2_2 = 12.0408
	1451	REAL(wp) :: zcoef_fco2_3 = 0.0327957
	1452	REAL(wp) :: zcoef_fco2_4 = 0.0000316528
	1453	REAL(wp) :: zcoef_fco2_5 = 2.0
	1454	REAL(wp) :: zcoef_fco2_6 = 57.7
	1455	REAL(wp) :: zcoef_fco2_7 = 0.118
	1456	REAL(wp) :: zcoef_fco2_8 = 82.0578
[9297]	1457	!!------------------------------------------------------------------------
[9296]	1458
	1459	IF ( kt <= nit000 ) THEN
	1460
[9297]	1461	!----------------------------------------------------------------------
[9296]	1462	! DIC and alkalinity balancing
[9297]	1463	!----------------------------------------------------------------------
[9296]	1464
	1465	IF ( ln_sfco2inc ) THEN
[10141]	1466	#if defined key_medusa && defined key_roam
[9296]	1467	! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
	1468	patm(1) = 1.0
	1469	phyd(1) = 0.0
	1470	DO jj = 1, jpj
	1471	DO ji = 1, jpi
[9322]	1472	CALL f2pCO2( sfco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
[9296]	1473	END DO
	1474	END DO
	1475	#else
	1476	! If assimilating fCO2, then convert to pCO2 using temperature
	1477	! See flux_gas.F90 within HadOCC for details of calculation
[9326]	1478	pco2_bkginc_temp(:,:) = sfco2_bkginc(:,:) / &
[9296]	1479	& EXP((zcoef_fco2_1 + &
	1480	& zcoef_fco2_2 * (tsn(:,:,1,1)+rt0) - &
	1481	& zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
	1482	& zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0) + &
	1483	& zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0))) / &
	1484	& (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
	1485	#endif
	1486	ELSE
	1487	pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
	1488	ENDIF
	1489
	1490	! Account for physics assimilation if required
	1491	IF ( ll_trainc ) THEN
	1492	tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
	1493	sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
	1494	ELSE
	1495	tem_bkg_temp(:,:) = tsn(:,:,1,1)
	1496	sal_bkg_temp(:,:) = tsn(:,:,1,2)
	1497	ENDIF
	1498
[10141]	1499	#if defined key_medusa
[9381]	1500	! Account for phytoplankton balancing if required
	1501	IF ( ln_phytobal ) THEN
[13097]	1502	IF ( ALLOCATED(phyto2d_balinc) ) THEN
	1503	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
	1504	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk)
	1505	ENDIF
	1506	IF ( ALLOCATED(phyto3d_balinc) ) THEN
	1507	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto3d_balinc(:,:,1,jpdic)
	1508	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto3d_balinc(:,:,1,jpalk)
	1509	ENDIF
[9296]	1510	ELSE
	1511	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
	1512	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
	1513	ENDIF
	1514
	1515	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
	1516	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
[9322]	1517	& pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
[9296]	1518
	1519	#elif defined key_hadocc
[9381]	1520	! Account for phytoplankton balancing if required
	1521	IF ( ln_phytobal ) THEN
[13097]	1522	IF ( ALLOCATED(phyto2d_balinc) ) THEN
	1523	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
	1524	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk)
	1525	ENDIF
	1526	IF ( ALLOCATED(phyto3d_balinc) ) THEN
	1527	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto3d_balinc(:,:,1,jp_had_dic)
	1528	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto3d_balinc(:,:,1,jp_had_alk)
	1529	ENDIF
[9296]	1530	ELSE
	1531	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
	1532	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
	1533	ENDIF
	1534
	1535	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
	1536	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
[9322]	1537	& pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
[9296]	1538
	1539	#else
	1540	CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
	1541	& 'but not defined a biogeochemical model' )
	1542	#endif
	1543
	1544	! Select mixed layer
	1545	IF ( ll_asmdin ) THEN
[10141]	1546	#if defined key_hadocc \|\| defined key_medusa
[9296]	1547	CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
	1548	& ' Mixed layer depth taken to be background maximum mld_max_bkg' )
	1549	zmld(:,:) = mld_max_bkg(:,:)
	1550	#else
	1551	CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
	1552	#endif
	1553	ELSE
	1554	SELECT CASE( mld_choice_bgc )
	1555	CASE ( 1 ) ! Turbocline/mixing depth [W points]
	1556	zmld(:,:) = hmld(:,:)
	1557	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
	1558	zmld(:,:) = hmlp(:,:)
	1559	CASE ( 3 ) ! Kara MLD [Interpolated]
	1560	#if defined key_karaml
	1561	IF ( ln_kara ) THEN
	1562	zmld(:,:) = hmld_kara(:,:)
	1563	ELSE
	1564	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
	1565	& ' but ln_kara=.false.' )
	1566	ENDIF
	1567	#else
	1568	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
	1569	& ' but is not defined' )
	1570	#endif
	1571	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
	1572	!zmld(:,:) = hmld_tref(:,:)
	1573	CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
	1574	& ' but is not available in this version' )
	1575	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
	1576	zmld(:,:) = hmlpt(:,:)
	1577	END SELECT
	1578	ENDIF
	1579
	1580	! Propagate through mixed layer
	1581	DO jj = 1, jpj
	1582	DO ji = 1, jpi
	1583	!
	1584	jkmax = jpk-1
	1585	DO jk = jpk-1, 1, -1
	1586	IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. &
	1587	& ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
	1588	zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
	1589	jkmax = jk
	1590	ENDIF
	1591	END DO
	1592	!
	1593	#if defined key_top
	1594	DO jk = 2, jkmax
[9322]	1595	pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
[9296]	1596	END DO
	1597	#endif
	1598	!
	1599	END DO
	1600	END DO
	1601
	1602	ENDIF
	1603
	1604	IF ( ll_asmiau ) THEN
	1605
	1606	!--------------------------------------------------------------------
	1607	! Incremental Analysis Updating
	1608	!--------------------------------------------------------------------
	1609
	1610	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
	1611
	1612	it = kt - nit000 + 1
	1613	zincwgt = pwgtiau(it) ! IAU weight for the current time step
	1614	! note this is not a tendency so should not be divided by rdt
	1615
	1616	IF(lwp) THEN
	1617	WRITE(numout,*)
	1618	IF ( ln_spco2inc ) THEN
	1619	WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
	1620	& kt,' with IAU weight = ', pwgtiau(it)
	1621	ELSE IF ( ln_sfco2inc ) THEN
	1622	WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
	1623	& kt,' with IAU weight = ', pwgtiau(it)
	1624	ENDIF
	1625	WRITE(numout,*) '~~~~~~~~~~~~'
	1626	ENDIF
	1627
	1628	! Update the biogeochemical variables
	1629	! Add directly to trn and trb, rather than to tra, because tra gets
	1630	! reset to zero at the start of trc_stp, called after this routine
[10141]	1631	#if defined key_medusa
[9382]	1632	WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1633	& trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
	1634	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1635	& pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1636	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
	1637	& pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1638	END WHERE
[9296]	1639	#elif defined key_hadocc
[9382]	1640	WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	1641	& trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
	1642	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	1643	& pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	1644	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
	1645	& pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
	1646	END WHERE
[9296]	1647	#endif
	1648
	1649	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1650	! which is called at end of model run
	1651
	1652	ENDIF
	1653
	1654	ELSEIF ( ll_asmdin ) THEN
	1655
	1656	!--------------------------------------------------------------------
	1657	! Direct Initialization
	1658	!--------------------------------------------------------------------
	1659
	1660	IF ( kt == nitdin_r ) THEN
	1661
	1662	neuler = 0 ! Force Euler forward step
	1663
	1664	! Initialize the now fields with the background + increment
	1665	! Background currently is what the model is initialised with
[9382]	1666	CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
[9296]	1667	& ' Background state is taken from model rather than background file' )
[10141]	1668	#if defined key_medusa
[9382]	1669	WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1670	& trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
	1671	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1672	& pco2_balinc(:,:,:,jp_msa0:jp_msa1)
	1673	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
	1674	END WHERE
[9296]	1675	#elif defined key_hadocc
[9382]	1676	WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
	1677	& trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
	1678	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
	1679	& pco2_balinc(:,:,:,jp_had0:jp_had1)
	1680	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
	1681	END WHERE
[9296]	1682	#endif
	1683
	1684	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1685	! which is called at end of model run
	1686	ENDIF
	1687	!
	1688	ENDIF
	1689	!
[9322]	1690	END SUBROUTINE pco2_asm_inc
[9296]	1691
[9297]	1692	!!===========================================================================
[9322]	1693	!!===========================================================================
	1694	!!===========================================================================
[9296]	1695
[9322]	1696	SUBROUTINE ph_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
	1697	& ll_trainc, pt_bkginc, ps_bkginc )
	1698	!!------------------------------------------------------------------------
	1699	!! * ROUTINE ph_asm_inc *
	1700	!!
	1701	!! ** Purpose : Apply the pH assimilation increments.
	1702	!!
[9381]	1703	!! ** Method : Calculate increments to DIC and alkalinity from pH
	1704	!! Use a similar approach to the pCO2 scheme
[9322]	1705	!! Direct initialization or Incremental Analysis Updating.
	1706	!!
	1707	!! ** Action :
	1708	!!------------------------------------------------------------------------
[9381]	1709	INTEGER, INTENT(IN) :: kt ! Current time step
	1710	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
	1711	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
	1712	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
	1713	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
	1714	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
[9322]	1715	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
	1716	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
[9381]	1717
	1718	REAL(wp) :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation
	1719	REAL(wp) :: DIC_IN, ALK_IN ! DIC/alk in pH calculation
	1720	REAL(wp) :: PH_OUT ! pH from pH calculation
	1721	REAL(wp) :: ph_bkg ! pH from pH calculation
	1722	REAL(wp) :: ph_dic ! pH from pH calculation
	1723	REAL(wp) :: ph_alk ! pH from pH calculation
	1724	REAL(wp) :: dph_ddic ! Change in pH wrt DIC
	1725	REAL(wp) :: dph_dalk ! Change in pH wrt alk
	1726	REAL(wp) :: weight ! Increment weighting
	1727	REAL(wp) :: zincwgt ! IAU weight for current time step
	1728	REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp ! DIC background
	1729	REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp ! Alkalinity background
	1730	REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp ! DIN background
	1731	REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp ! Silicate background
	1732	REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp ! Temperature background
	1733	REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp ! Salinity background
	1734	REAL(wp), DIMENSION(19) :: dummy_out ! Dummy variables
	1735	INTEGER :: ji, jj, jk, jx ! Loop counters
	1736	INTEGER :: it ! Index
[9322]	1737	!!------------------------------------------------------------------------
[9381]	1738
[10141]	1739	#if ! defined key_medusa
[9381]	1740	CALL ctl_stop( ' pH balancing only implemented for MEDUSA' )
	1741	#else
	1742
	1743	IF ( kt <= nit000 ) THEN
	1744
	1745	!----------------------------------------------------------------------
	1746	! DIC and alkalinity balancing
	1747	!----------------------------------------------------------------------
	1748
	1749	! Account for physics assimilation if required
	1750	IF ( ll_trainc ) THEN
	1751	tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:)
	1752	sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:)
	1753	ELSE
	1754	tem_bkg_temp(:,:,:) = tsn(:,:,:,1)
	1755	sal_bkg_temp(:,:,:) = tsn(:,:,:,2)
	1756	ENDIF
	1757
	1758	! Account for phytoplankton balancing if required
	1759	IF ( ln_phytobal ) THEN
[13097]	1760	IF ( ALLOCATED(phyto2d_balinc) ) THEN
	1761	dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic)
	1762	alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk)
	1763	din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin)
	1764	sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil)
	1765	ENDIF
	1766	IF ( ALLOCATED(phyto3d_balinc) ) THEN
	1767	dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto3d_balinc(:,:,:,jpdic)
	1768	alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto3d_balinc(:,:,:,jpalk)
	1769	din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto3d_balinc(:,:,:,jpdin)
	1770	sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto3d_balinc(:,:,:,jpsil)
	1771	ENDIF
[9381]	1772	ELSE
	1773	dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic)
	1774	alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk)
	1775	din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin)
	1776	sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil)
	1777	ENDIF
	1778
	1779	! Account for pCO2 balancing if required
	1780	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	1781	dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic)
	1782	alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk)
	1783	ENDIF
	1784
	1785	! Loop over grid points and calculate dpH/dDIC and dpH/dAlk
	1786	! This requires three calls to the MOCSY carbonate package
	1787	! One to find pH at background DIC and alkalinity
	1788	! One to find pH when DIC is increased by 10
	1789	! One to find pH when alkalinity is increased by 10
	1790	! Then calculate the gradients
	1791
	1792	DO jk = 1, jpk
	1793	DO jj = 1, jpj
	1794	DO ji = 1, jpi
	1795
	1796	IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN
	1797
	1798	DO jx = 1, 3
	1799	IF ( jx == 1 ) THEN
	1800	DIC_IN = dic_bkg_temp(ji,jj,jk)
	1801	ALK_IN = alk_bkg_temp(ji,jj,jk)
	1802	ELSE IF ( jx == 2 ) THEN
	1803	DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch
	1804	ALK_IN = alk_bkg_temp(ji,jj,jk)
	1805	ELSE IF ( jx == 3 ) THEN
	1806	DIC_IN = dic_bkg_temp(ji,jj,jk)
	1807	ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch
	1808	ENDIF
	1809
	1810	CALL mocsy_interface(tsn(ji,jj,jk,jp_tem), & ! IN: temperature
	1811	& tsn(ji,jj,jk,jp_sal), & ! salinity
	1812	& ALK_IN, & ! alkalinity
	1813	& DIC_IN, & ! DIC
	1814	& sil_bkg_temp(ji,jj,jk), & ! silicate
	1815	& din_bkg_temp(ji,jj,jk)/16.0, & ! phosphate (Redfield from DIN)
	1816	& 1.0, & ! atmospheric pressure (dummy)
	1817	& fsdept(ji,jj,jk), & ! depth
	1818	& gphit(ji,jj), & ! latitude
	1819	& 1.0, & ! gas transfer (dummy)
	1820	& 1.0, & ! atmospheric xCO2 (dummy)
	1821	& 1, & ! number of points
	1822	& PH_OUT, & ! OUT: pH
	1823	& dummy_out(1), dummy_out(2), & ! stuff we don't care about
	1824	& dummy_out(3), dummy_out(4), &
	1825	& dummy_out(5), dummy_out(6), &
	1826	& dummy_out(7), dummy_out(8), &
	1827	& dummy_out(9), dummy_out(10), &
	1828	& dummy_out(11), dummy_out(12), &
	1829	& dummy_out(13), dummy_out(14), &
	1830	& dummy_out(15), dummy_out(16), &
	1831	& dummy_out(17), dummy_out(18), &
	1832	& dummy_out(19) )
	1833
	1834	IF ( jx == 1 ) THEN
	1835	ph_bkg = PH_OUT
	1836	ELSE IF ( jx == 2 ) THEN
	1837	ph_dic = PH_OUT
	1838	ELSE IF ( jx == 3 ) THEN
	1839	ph_alk = PH_OUT
	1840	ENDIF
	1841	END DO
	1842
	1843	dph_ddic = (ph_dic - ph_bkg) / zsearch
	1844	dph_dalk = (ph_alk - ph_bkg) / zsearch
	1845	weight = pph_bkginc(ji,jj,jk) / ( (dph_ddicdph_ddic) + (dph_dalkdph_dalk) )
	1846
	1847	ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic
	1848	ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk
	1849
	1850	ENDIF
	1851
	1852	END DO
	1853	END DO
	1854	END DO
	1855
	1856	ENDIF
[9322]	1857
[9381]	1858	IF ( ll_asmiau ) THEN
	1859
	1860	!--------------------------------------------------------------------
	1861	! Incremental Analysis Updating
	1862	!--------------------------------------------------------------------
	1863
	1864	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
	1865
	1866	it = kt - nit000 + 1
	1867	zincwgt = pwgtiau(it) ! IAU weight for the current time step
	1868	! note this is not a tendency so should not be divided by rdt
	1869
	1870	IF(lwp) THEN
	1871	WRITE(numout,*)
	1872	WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', &
	1873	& kt,' with IAU weight = ', pwgtiau(it)
	1874	WRITE(numout,*) '~~~~~~~~~~~~'
	1875	ENDIF
	1876
	1877	! Update the biogeochemical variables
	1878	! Add directly to trn and trb, rather than to tra, because tra gets
	1879	! reset to zero at the start of trc_stp, called after this routine
[9382]	1880	WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1881	& trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
	1882	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1883	& ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1884	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
	1885	& ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
	1886	END WHERE
[9381]	1887
	1888	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1889	! which is called at end of model run
	1890
	1891	ENDIF
	1892
	1893	ELSEIF ( ll_asmdin ) THEN
	1894
	1895	!--------------------------------------------------------------------
	1896	! Direct Initialization
	1897	!--------------------------------------------------------------------
	1898
	1899	IF ( kt == nitdin_r ) THEN
	1900
	1901	neuler = 0 ! Force Euler forward step
	1902
	1903	! Initialize the now fields with the background + increment
	1904	! Background currently is what the model is initialised with
	1905	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pH assimilation', &
	1906	& ' Background state is taken from model rather than background file' )
[9382]	1907	WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
	1908	& trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
	1909	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
	1910	& ph_balinc(:,:,:,jp_msa0:jp_msa1)
	1911	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
	1912	END WHERE
[9381]	1913
	1914	! Do not deallocate arrays - needed by asm_bgc_bal_wri
	1915	! which is called at end of model run
	1916	ENDIF
	1917	!
	1918	ENDIF
	1919	#endif
[9322]	1920	!
	1921	END SUBROUTINE ph_asm_inc
	1922
	1923	!!===========================================================================
	1924	!!===========================================================================
	1925	!!===========================================================================
	1926
	1927	SUBROUTINE bgc3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
	1928	!!----------------------------------------------------------------------
	1929	!! * ROUTINE dyn_asm_inc *
	1930	!!
	1931	!! ** Purpose : Apply generic 3D biogeochemistry assimilation increments.
	1932	!!
	1933	!! ** Method : Direct initialization or Incremental Analysis Updating.
	1934	!!
	1935	!! ** Action :
	1936	!!----------------------------------------------------------------------
	1937	INTEGER, INTENT(IN) :: kt ! Current time step
	1938	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
	1939	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
	1940	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
	1941	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
	1942	!
	1943	INTEGER :: jk ! Loop counter
	1944	INTEGER :: it ! Index
	1945	REAL(wp) :: zincwgt ! IAU weight for current time step
	1946	REAL(wp) :: zincper ! IAU interval in seconds
	1947	!!----------------------------------------------------------------------
	1948
[9382]	1949	IF ( kt <= nit000 ) THEN
	1950
	1951	!----------------------------------------------------------------------
	1952	! Remove any other balancing increments
	1953	!----------------------------------------------------------------------
	1954
	1955	IF ( ln_pno3inc ) THEN
[10141]	1956	#if defined key_hadocc \|\| defined key_medusa
[9382]	1957	#if defined key_hadocc
	1958	it = jp_had_nut
[10141]	1959	#elif defined key_medusa
[9382]	1960	it = jpdin
	1961	#endif
[13097]	1962	IF ( ALLOCATED(phyto2d_balinc) ) THEN
[9382]	1963	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
	1964	ENDIF
[13097]	1965	IF ( ALLOCATED(phyto3d_balinc) ) THEN
[9382]	1966	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
	1967	ENDIF
	1968	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	1969	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
	1970	ENDIF
	1971	IF ( ln_pphinc ) THEN
	1972	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - ph_balinc(:,:,:,it)
	1973	ENDIF
[10057]	1974	#else
	1975	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	1976	#endif
[9382]	1977	ENDIF
	1978
	1979	IF ( ln_psi4inc ) THEN
[10141]	1980	#if defined key_medusa
[9382]	1981	it = jpsil
[13097]	1982	IF ( ALLOCATED(phyto2d_balinc) ) THEN
[9382]	1983	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
	1984	ENDIF
[13097]	1985	IF ( ALLOCATED(phyto3d_balinc) ) THEN
[9382]	1986	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
	1987	ENDIF
	1988	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	1989	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
	1990	ENDIF
	1991	IF ( ln_pphinc ) THEN
	1992	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - ph_balinc(:,:,:,it)
	1993	ENDIF
[10057]	1994	#else
	1995	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	1996	#endif
[9382]	1997	ENDIF
	1998
	1999	IF ( ln_pdicinc ) THEN
[10141]	2000	#if defined key_hadocc \|\| defined key_medusa
[9382]	2001	#if defined key_hadocc
	2002	it = jp_had_dic
[10141]	2003	#elif defined key_medusa
[9382]	2004	it = jpdic
	2005	#endif
[13097]	2006	IF ( ALLOCATED(phyto2d_balinc) ) THEN
[9382]	2007	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
	2008	ENDIF
[13097]	2009	IF ( ALLOCATED(phyto3d_balinc) ) THEN
[9382]	2010	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
	2011	ENDIF
	2012	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	2013	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
	2014	ENDIF
	2015	IF ( ln_pphinc ) THEN
	2016	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - ph_balinc(:,:,:,it)
	2017	ENDIF
[10057]	2018	#else
	2019	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2020	#endif
[9382]	2021	ENDIF
	2022
	2023	IF ( ln_palkinc ) THEN
[10141]	2024	#if defined key_hadocc \|\| defined key_medusa
[9382]	2025	#if defined key_hadocc
	2026	it = jp_had_alk
[10141]	2027	#elif defined key_medusa
[9382]	2028	it = jpalk
	2029	#endif
[13097]	2030	IF ( ALLOCATED(phyto2d_balinc) ) THEN
[9382]	2031	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
	2032	ENDIF
[13097]	2033	IF ( ALLOCATED(phyto3d_balinc) ) THEN
[9382]	2034	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
	2035	ENDIF
	2036	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	2037	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
	2038	ENDIF
	2039	IF ( ln_pphinc ) THEN
	2040	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - ph_balinc(:,:,:,it)
	2041	ENDIF
[10057]	2042	#else
	2043	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2044	#endif
[9382]	2045	ENDIF
	2046
	2047	IF ( ln_po2inc ) THEN
[10141]	2048	#if defined key_medusa
[9382]	2049	it = jpoxy
[13097]	2050	IF ( ALLOCATED(phyto2d_balinc) ) THEN
[9382]	2051	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
	2052	ENDIF
[13097]	2053	IF ( ALLOCATED(phyto3d_balinc) ) THEN
[9382]	2054	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
	2055	ENDIF
	2056	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
	2057	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
	2058	ENDIF
	2059	IF ( ln_pphinc ) THEN
	2060	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - ph_balinc(:,:,:,it)
	2061	ENDIF
[10057]	2062	#else
	2063	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2064	#endif
[9382]	2065	ENDIF
	2066
	2067	ENDIF
	2068
[9322]	2069	IF ( ll_asmiau ) THEN
	2070
	2071	!--------------------------------------------------------------------
	2072	! Incremental Analysis Updating
	2073	!--------------------------------------------------------------------
	2074
	2075	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
	2076
	2077	it = kt - nit000 + 1
[9326]	2078	zincwgt = pwgtiau(it) ! IAU weight for the current time step
	2079	! note this is not a tendency so should not be divided by rdt
[9322]	2080
	2081	IF(lwp) THEN
	2082	WRITE(numout,*)
	2083	WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', &
	2084	& kt,' with IAU weight = ', pwgtiau(it)
	2085	WRITE(numout,*) '~~~~~~~~~~~~'
	2086	ENDIF
	2087
	2088	! Update the 3D BGC variables
	2089	! Add directly to trn and trb, rather than to tra, because tra gets
	2090	! reset to zero at the start of trc_stp, called after this routine
	2091	! Don't apply increments if they'll take concentrations negative
	2092
	2093	IF ( ln_pno3inc ) THEN
	2094	#if defined key_hadocc
	2095	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
	2096	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2097	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
	2098	trb(:,:,:,jp_had_nut) = trb(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
	2099	END WHERE
[10141]	2100	#elif defined key_medusa
[9322]	2101	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
	2102	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2103	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
	2104	trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
	2105	END WHERE
	2106	#else
	2107	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2108	#endif
	2109	ENDIF
	2110
	2111	IF ( ln_psi4inc ) THEN
[10141]	2112	#if defined key_medusa
[9322]	2113	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
	2114	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2115	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
	2116	trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
	2117	END WHERE
	2118	#else
	2119	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2120	#endif
	2121	ENDIF
	2122
	2123	IF ( ln_pdicinc ) THEN
	2124	#if defined key_hadocc
	2125	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
	2126	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2127	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
	2128	trb(:,:,:,jp_had_dic) = trb(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
	2129	END WHERE
[10141]	2130	#elif defined key_medusa
[9322]	2131	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
	2132	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2133	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
	2134	trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
	2135	END WHERE
	2136	#else
	2137	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2138	#endif
	2139	ENDIF
	2140
	2141	IF ( ln_palkinc ) THEN
	2142	#if defined key_hadocc
	2143	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
	2144	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2145	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
	2146	trb(:,:,:,jp_had_alk) = trb(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
	2147	END WHERE
[10141]	2148	#elif defined key_medusa
[9322]	2149	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
	2150	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2151	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
	2152	trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
	2153	END WHERE
	2154	#else
	2155	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2156	#endif
	2157	ENDIF
	2158
	2159	IF ( ln_po2inc ) THEN
[10141]	2160	#if defined key_medusa
[9322]	2161	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
	2162	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
	2163	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
	2164	trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
	2165	END WHERE
	2166	#else
	2167	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2168	#endif
	2169	ENDIF
	2170
	2171	IF ( kt == nitiaufin_r ) THEN
	2172	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
	2173	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
	2174	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
	2175	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
	2176	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
	2177	ENDIF
	2178
	2179	ENDIF
	2180
	2181	ELSEIF ( ll_asmdin ) THEN
	2182
	2183	!--------------------------------------------------------------------
	2184	! Direct Initialization
	2185	!--------------------------------------------------------------------
	2186
	2187	IF ( kt == nitdin_r ) THEN
	2188
	2189	neuler = 0 ! Force Euler forward step
	2190
	2191	! Initialize the now fields with the background + increment
	2192	! Background currently is what the model is initialised with
	2193	#if defined key_hadocc
	2194	CALL ctl_warn( ' Doing direct initialisation of HadOCC with 3D BGC assimilation', &
	2195	& ' Background state is taken from model rather than background file' )
[10141]	2196	#elif defined key_medusa
[9322]	2197	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with 3D BGC assimilation', &
	2198	& ' Background state is taken from model rather than background file' )
	2199	#endif
	2200
	2201	IF ( ln_pno3inc ) THEN
	2202	#if defined key_hadocc
	2203	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
	2204	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) > 0.0_wp )
	2205	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:)
	2206	trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut)
	2207	END WHERE
[10141]	2208	#elif defined key_medusa
[9322]	2209	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
	2210	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) > 0.0_wp )
	2211	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:)
	2212	trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
	2213	END WHERE
	2214	#else
	2215	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2216	#endif
	2217	ENDIF
	2218
	2219	IF ( ln_psi4inc ) THEN
[10141]	2220	#if defined key_medusa
[9322]	2221	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
	2222	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) > 0.0_wp )
	2223	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:)
	2224	trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
	2225	END WHERE
	2226	#else
	2227	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2228	#endif
	2229	ENDIF
	2230
	2231	IF ( ln_pdicinc ) THEN
	2232	#if defined key_hadocc
	2233	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
	2234	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) > 0.0_wp )
	2235	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:)
	2236	trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic)
	2237	END WHERE
[10141]	2238	#elif defined key_medusa
[9322]	2239	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
	2240	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) > 0.0_wp )
	2241	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:)
	2242	trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
	2243	END WHERE
	2244	#else
	2245	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2246	#endif
	2247	ENDIF
	2248
	2249	IF ( ln_palkinc ) THEN
	2250	#if defined key_hadocc
	2251	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
	2252	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) > 0.0_wp )
	2253	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:)
	2254	trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk)
	2255	END WHERE
[10141]	2256	#elif defined key_medusa
[9322]	2257	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
	2258	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) > 0.0_wp )
	2259	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:)
	2260	trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
	2261	END WHERE
	2262	#else
	2263	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2264	#endif
	2265	ENDIF
	2266
	2267	IF ( ln_po2inc ) THEN
[10141]	2268	#if defined key_medusa
[9322]	2269	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
	2270	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) > 0.0_wp )
	2271	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:)
	2272	trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
	2273	END WHERE
	2274	#else
	2275	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
	2276	#endif
	2277	ENDIF
	2278
	2279	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
	2280	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
	2281	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
	2282	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
	2283	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
	2284	ENDIF
	2285	!
	2286	ENDIF
	2287	!
	2288	END SUBROUTINE bgc3d_asm_inc
	2289
	2290	!!===========================================================================
	2291
[9296]	2292	END MODULE asmbgc

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

Download in other formats: