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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 9862

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

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

Download in other formats: