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

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

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

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_phytobal3d/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 13106

Last change on this file since 13106 was 13106, checked in by dford, 4 years ago
Allow 3D PFT chlorophyll increments to be applied.
File size: 155.0 KB

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

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

Download in other formats: