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

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

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

Context Navigation

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

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

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

Download in other formats: