Context Navigation

asmbgc.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 11101

Download in other formats: