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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_port2021/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 15332

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

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

Download in other formats: