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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_v3/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 10077

Last change on this file since 10077 was 10077, checked in by dford, 6 years ago
Fix for compiling with key_medusa but not key_foam_medusa.
File size: 103.8 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 && defined key_foam_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 && defined key_foam_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 && defined key_foam_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 && defined key_foam_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 && defined key_foam_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 && defined key_foam_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 && defined key_foam_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 && defined key_foam_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 && defined key_foam_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	#if defined key_top && ( defined key_hadocc \|\| (defined key_medusa && defined key_foam_medusa) )
993	CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', &
994	& ' Mixed layer depth taken to be background maximum mld_max_bkg' )
995	zmld(:,:) = mld_max_bkg(:,:)
996	#else
997	CALL ctl_stop( ' Should not be doing biogeochemical assimilation without key_top' )
998	#endif
999	ELSE
1000	SELECT CASE( mld_choice_bgc )
1001	CASE ( 1 ) ! Turbocline/mixing depth [W points]
1002	zmld(:,:) = hmld(:,:)
1003	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
1004	zmld(:,:) = hmlp(:,:)
1005	CASE ( 3 ) ! Kara MLD [Interpolated]
1006	#if defined key_karaml
1007	IF ( ln_kara ) THEN
1008	zmld(:,:) = hmld_kara(:,:)
1009	ELSE
1010	CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
1011	& ' but ln_kara=.false.' )
1012	ENDIF
1013	#else
1014	CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
1015	& ' but is not defined' )
1016	#endif
1017	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
1018	!zmld(:,:) = hmld_tref(:,:)
1019	CALL ctl_stop( ' hmld_tref mixed layer requested for phyto2d assimilation,', &
1020	& ' but is not available in this version' )
1021	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
1022	zmld(:,:) = hmlpt(:,:)
1023	END SELECT
1024	ENDIF
1025
1026	zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
1027
1028	#if defined key_medusa && defined key_foam_medusa
1029	CALL asm_phyto2d_bal_medusa( (ln_slchltotinc .OR. ln_schltotinc), &
1030	& zinc_chltot, &
1031	& ln_slchldiainc, &
1032	& zinc_chldia, &
1033	& ln_slchlnoninc, &
1034	& zinc_chlnon, &
1035	& ln_slphytotinc, &
1036	& zinc_phytot, &
1037	& ln_slphydiainc, &
1038	& zinc_phydia, &
1039	& ln_slphynoninc, &
1040	& zinc_phynon, &
1041	& zincper, &
1042	& rn_maxchlinc, ln_phytobal, zmld, &
1043	& pgrow_avg_bkg, ploss_avg_bkg, &
1044	& phyt_avg_bkg, mld_max_bkg, &
1045	& tracer_bkg, phyto2d_balinc )
1046	#elif defined key_hadocc
1047	CALL asm_phyto2d_bal_hadocc( (ln_slchltotinc .OR. ln_schltotinc), &
1048	& zinc_chltot, &
1049	& ln_slphytotinc, &
1050	& zinc_phytot, &
1051	& zincper, &
1052	& rn_maxchlinc, ln_phytobal, zmld, &
1053	& pgrow_avg_bkg, ploss_avg_bkg, &
1054	& phyt_avg_bkg, mld_max_bkg, &
1055	& cchl_p_bkg(:,:,1), &
1056	& tracer_bkg, phyto2d_balinc )
1057	#else
1058	CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', &
1059	& 'but not defined a biogeochemical model' )
1060	#endif
1061
1062	ENDIF
1063
1064	IF ( ll_asmiau ) THEN
1065
1066	!--------------------------------------------------------------------
1067	! Incremental Analysis Updating
1068	!--------------------------------------------------------------------
1069
1070	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1071
1072	it = kt - nit000 + 1
1073	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1074	! note this is not a tendency so should not be divided by rdt
1075
1076	IF(lwp) THEN
1077	WRITE(numout,*)
1078	WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', &
1079	& kt,' with IAU weight = ', pwgtiau(it)
1080	WRITE(numout,*) '~~~~~~~~~~~~'
1081	ENDIF
1082
1083	! Update the biogeochemical variables
1084	! Add directly to trn and trb, rather than to tra, because tra gets
1085	! reset to zero at the start of trc_stp, called after this routine
1086	#if defined key_medusa && defined key_foam_medusa
1087	WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1088	& trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1089	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1090	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1091	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1092	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1093	END WHERE
1094	#elif defined key_hadocc
1095	WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1096	& trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1097	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1098	& phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1099	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1100	& phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1101	END WHERE
1102	#endif
1103
1104	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1105	! which is called at end of model run
1106	ENDIF
1107
1108	ELSEIF ( ll_asmdin ) THEN
1109
1110	!--------------------------------------------------------------------
1111	! Direct Initialization
1112	!--------------------------------------------------------------------
1113
1114	IF ( kt == nitdin_r ) THEN
1115
1116	neuler = 0 ! Force Euler forward step
1117
1118	! Initialize the now fields with the background + increment
1119	! Background currently is what the model is initialised with
1120	CALL ctl_warn( ' Doing direct initialisation with phyto2d assimilation', &
1121	& ' Background state is taken from model rather than background file' )
1122	#if defined key_medusa && defined key_foam_medusa
1123	WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1124	& trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1125	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1126	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1)
1127	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1128	END WHERE
1129	#elif defined key_hadocc
1130	WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1131	& trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1132	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1133	& phyto2d_balinc(:,:,:,jp_had0:jp_had1)
1134	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1135	END WHERE
1136	#endif
1137
1138	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1139	! which is called at end of model run
1140	ENDIF
1141	!
1142	ENDIF
1143	!
1144	END SUBROUTINE phyto2d_asm_inc
1145
1146	!!===========================================================================
1147	!!===========================================================================
1148	!!===========================================================================
1149
1150	SUBROUTINE phyto3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1151	!!------------------------------------------------------------------------
1152	!! * ROUTINE phyto3d_asm_inc *
1153	!!
1154	!! ** Purpose : Apply the profile chlorophyll assimilation increments.
1155	!!
1156	!! ** Method : Calculate increments to state variables.
1157	!! Direct initialization or Incremental Analysis Updating.
1158	!!
1159	!! ** Action :
1160	!!------------------------------------------------------------------------
1161	INTEGER, INTENT(IN) :: kt ! Current time step
1162	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1163	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1164	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1165	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1166	!
1167	INTEGER :: ji, jj, jk ! Loop counters
1168	INTEGER :: it ! Index
1169	REAL(wp) :: zincwgt ! IAU weight for timestep
1170	REAL(wp) :: zfrac_chn ! Fraction of jpchn
1171	REAL(wp) :: zfrac_chd ! Fraction of jpchd
1172	REAL(wp) :: zrat_phn_chn ! jpphn:jpchn ratio
1173	REAL(wp) :: zrat_phd_chd ! jpphd:jpchd ratio
1174	REAL(wp) :: zrat_pds_chd ! jppds:jpchd ratio
1175	REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc ! Chlorophyll increments
1176	REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl ! Chlorophyll background
1177	!!------------------------------------------------------------------------
1178
1179	IF ( kt <= nit000 ) THEN
1180
1181	IF ( ln_plchltotinc ) THEN
1182	! Convert log10(chlorophyll) increment back to a chlorophyll increment
1183	! In order to transform logchl incs to chl incs, need to account for model
1184	! background, cannot simply do 10^logchl_bkginc. Need to:
1185	! 1) Add logchl inc to log10(background) to get log10(analysis)
1186	! 2) Take 10^log10(analysis) to get analysis
1187	! 3) Subtract background from analysis to get chl incs
1188	! If rn_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
1189	#if defined key_medusa && defined key_foam_medusa
1190	bkg_chl(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd)
1191	#elif defined key_hadocc
1192	bkg_chl(:,:,:) = chl_bkg(:,:,:)
1193	#endif
1194	DO jk = 1, jpk
1195	DO jj = 1, jpj
1196	DO ji = 1, jpi
1197	IF ( bkg_chl(ji,jj,jk) > 0.0 ) THEN
1198	chl_inc(ji,jj,jk) = 10**( LOG10( bkg_chl(ji,jj,jk) ) + plchltot_bkginc(ji,jj,jk) ) - bkg_chl(ji,jj,jk)
1199	IF ( rn_maxchlinc > 0.0 ) THEN
1200	chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( chl_inc(ji,jj,jk), rn_maxchlinc ) )
1201	ENDIF
1202	ELSE
1203	chl_inc(ji,jj,jk) = 0.0
1204	ENDIF
1205	END DO
1206	END DO
1207	END DO
1208	ELSE IF ( ln_pchltotinc ) THEN
1209	DO jk = 1, jpk
1210	DO jj = 1, jpj
1211	DO ji = 1, jpi
1212	IF ( rn_maxchlinc > 0.0 ) THEN
1213	chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( pchltot_bkginc(ji,jj,jk), rn_maxchlinc ) )
1214	ELSE
1215	chl_inc(ji,jj,jk) = pchltot_bkginc(ji,jj,jk)
1216	ENDIF
1217	END DO
1218	END DO
1219	END DO
1220	ENDIF
1221
1222	#if defined key_medusa && defined key_foam_medusa
1223	! Loop over each grid point partioning the increments based on existing ratios
1224	DO jk = 1, jpk
1225	DO jj = 1, jpj
1226	DO ji = 1, jpi
1227	IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN
1228	zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd))
1229	zfrac_chd = 1.0 - zfrac_chn
1230	phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn
1231	phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd
1232	zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn)
1233	zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd)
1234	zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd)
1235	phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn
1236	phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd
1237	phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd
1238	ENDIF
1239	END DO
1240	END DO
1241	END DO
1242	#elif defined key_hadocc
1243	phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:)
1244	#else
1245	CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', &
1246	& 'but not defined a biogeochemical model' )
1247	#endif
1248
1249	ENDIF
1250
1251	IF ( ll_asmiau ) THEN
1252
1253	!--------------------------------------------------------------------
1254	! Incremental Analysis Updating
1255	!--------------------------------------------------------------------
1256
1257	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1258
1259	it = kt - nit000 + 1
1260	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1261	! note this is not a tendency so should not be divided by rdt
1262
1263	IF(lwp) THEN
1264	WRITE(numout,*)
1265	WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', &
1266	& kt,' with IAU weight = ', pwgtiau(it)
1267	WRITE(numout,*) '~~~~~~~~~~~~'
1268	ENDIF
1269
1270	! Update the biogeochemical variables
1271	! Add directly to trn and trb, rather than to tra, because tra gets
1272	! reset to zero at the start of trc_stp, called after this routine
1273	#if defined key_medusa && defined key_foam_medusa
1274	WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1275	& trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1276	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1277	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1278	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1279	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1280	END WHERE
1281	#elif defined key_hadocc
1282	WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1283	& trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1284	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1285	& phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1286	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1287	& phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1288	END WHERE
1289	#endif
1290
1291	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1292	! which is called at end of model run
1293	ENDIF
1294
1295	ELSEIF ( ll_asmdin ) THEN
1296
1297	!--------------------------------------------------------------------
1298	! Direct Initialization
1299	!--------------------------------------------------------------------
1300
1301	IF ( kt == nitdin_r ) THEN
1302
1303	neuler = 0 ! Force Euler forward step
1304
1305	! Initialize the now fields with the background + increment
1306	! Background currently is what the model is initialised with
1307	CALL ctl_warn( ' Doing direct initialisation with phyto3d assimilation', &
1308	& ' Background state is taken from model rather than background file' )
1309	#if defined key_medusa && defined key_foam_medusa
1310	WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1311	& trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1312	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1313	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1)
1314	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1315	END WHERE
1316	#elif defined key_hadocc
1317	WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1318	& trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1319	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1320	& phyto3d_balinc(:,:,:,jp_had0:jp_had1)
1321	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1322	END WHERE
1323	#endif
1324
1325	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1326	! which is called at end of model run
1327	ENDIF
1328	!
1329	ENDIF
1330	!
1331	END SUBROUTINE phyto3d_asm_inc
1332
1333	!!===========================================================================
1334	!!===========================================================================
1335	!!===========================================================================
1336
1337	SUBROUTINE pco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
1338	& ll_trainc, pt_bkginc, ps_bkginc )
1339	!!------------------------------------------------------------------------
1340	!! * ROUTINE pco2_asm_inc *
1341	!!
1342	!! ** Purpose : Apply the pco2/fco2 assimilation increments.
1343	!!
1344	!! ** Method : Calculate increments to state variables using carbon
1345	!! balancing.
1346	!! Direct initialization or Incremental Analysis Updating.
1347	!!
1348	!! ** Action :
1349	!!------------------------------------------------------------------------
1350	INTEGER, INTENT(IN) :: kt ! Current time step
1351	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1352	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1353	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1354	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1355	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
1356	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
1357	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
1358	!
1359	INTEGER :: ji, jj, jk ! Loop counters
1360	INTEGER :: it ! Index
1361	INTEGER :: jkmax ! Index of mixed layer depth
1362	REAL(wp) :: zincwgt ! IAU weight for current time step
1363	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth for balancing
1364	REAL(wp), DIMENSION(jpi,jpj) :: pco2_bkginc_temp ! For fCO2 conversion
1365	REAL(wp), DIMENSION(jpi,jpj) :: dic_bkg_temp ! DIC background
1366	REAL(wp), DIMENSION(jpi,jpj) :: alk_bkg_temp ! Alkalinity background
1367	REAL(wp), DIMENSION(jpi,jpj) :: tem_bkg_temp ! Temperature background
1368	REAL(wp), DIMENSION(jpi,jpj) :: sal_bkg_temp ! Salinity background
1369	REAL(wp), DIMENSION(1) :: patm ! Atmospheric pressure
1370	REAL(wp), DIMENSION(1) :: phyd ! Hydrostatic pressure
1371
1372	! Coefficients for fCO2 to pCO2 conversion
1373	REAL(wp) :: zcoef_fco2_1 = -1636.75
1374	REAL(wp) :: zcoef_fco2_2 = 12.0408
1375	REAL(wp) :: zcoef_fco2_3 = 0.0327957
1376	REAL(wp) :: zcoef_fco2_4 = 0.0000316528
1377	REAL(wp) :: zcoef_fco2_5 = 2.0
1378	REAL(wp) :: zcoef_fco2_6 = 57.7
1379	REAL(wp) :: zcoef_fco2_7 = 0.118
1380	REAL(wp) :: zcoef_fco2_8 = 82.0578
1381	!!------------------------------------------------------------------------
1382
1383	IF ( kt <= nit000 ) THEN
1384
1385	!----------------------------------------------------------------------
1386	! DIC and alkalinity balancing
1387	!----------------------------------------------------------------------
1388
1389	IF ( ln_sfco2inc ) THEN
1390	#if defined key_medusa && defined key_foam_medusa && defined key_roam
1391	! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
1392	patm(1) = 1.0
1393	phyd(1) = 0.0
1394	DO jj = 1, jpj
1395	DO ji = 1, jpi
1396	CALL f2pCO2( sfco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
1397	END DO
1398	END DO
1399	#else
1400	! If assimilating fCO2, then convert to pCO2 using temperature
1401	! See flux_gas.F90 within HadOCC for details of calculation
1402	pco2_bkginc_temp(:,:) = sfco2_bkginc(:,:) / &
1403	& EXP((zcoef_fco2_1 + &
1404	& zcoef_fco2_2 * (tsn(:,:,1,1)+rt0) - &
1405	& zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
1406	& zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0) + &
1407	& zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0))) / &
1408	& (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
1409	#endif
1410	ELSE
1411	pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
1412	ENDIF
1413
1414	! Account for physics assimilation if required
1415	IF ( ll_trainc ) THEN
1416	tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
1417	sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
1418	ELSE
1419	tem_bkg_temp(:,:) = tsn(:,:,1,1)
1420	sal_bkg_temp(:,:) = tsn(:,:,1,2)
1421	ENDIF
1422
1423	#if defined key_medusa && defined key_foam_medusa
1424	! Account for phytoplankton balancing if required
1425	IF ( ln_phytobal ) THEN
1426	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
1427	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk)
1428	ELSE
1429	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
1430	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
1431	ENDIF
1432
1433	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1434	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
1435	& pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
1436
1437	#elif defined key_hadocc
1438	! Account for phytoplankton balancing if required
1439	IF ( ln_phytobal ) THEN
1440	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
1441	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk)
1442	ELSE
1443	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
1444	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
1445	ENDIF
1446
1447	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1448	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
1449	& pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
1450
1451	#else
1452	CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
1453	& 'but not defined a biogeochemical model' )
1454	#endif
1455
1456	! Select mixed layer
1457	IF ( ll_asmdin ) THEN
1458	#if defined key_hadocc \|\| (defined key_medusa && defined key_foam_medusa)
1459	CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
1460	& ' Mixed layer depth taken to be background maximum mld_max_bkg' )
1461	zmld(:,:) = mld_max_bkg(:,:)
1462	#else
1463	CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
1464	#endif
1465	ELSE
1466	SELECT CASE( mld_choice_bgc )
1467	CASE ( 1 ) ! Turbocline/mixing depth [W points]
1468	zmld(:,:) = hmld(:,:)
1469	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
1470	zmld(:,:) = hmlp(:,:)
1471	CASE ( 3 ) ! Kara MLD [Interpolated]
1472	#if defined key_karaml
1473	IF ( ln_kara ) THEN
1474	zmld(:,:) = hmld_kara(:,:)
1475	ELSE
1476	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1477	& ' but ln_kara=.false.' )
1478	ENDIF
1479	#else
1480	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1481	& ' but is not defined' )
1482	#endif
1483	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
1484	!zmld(:,:) = hmld_tref(:,:)
1485	CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
1486	& ' but is not available in this version' )
1487	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
1488	zmld(:,:) = hmlpt(:,:)
1489	END SELECT
1490	ENDIF
1491
1492	! Propagate through mixed layer
1493	DO jj = 1, jpj
1494	DO ji = 1, jpi
1495	!
1496	jkmax = jpk-1
1497	DO jk = jpk-1, 1, -1
1498	IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. &
1499	& ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
1500	zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
1501	jkmax = jk
1502	ENDIF
1503	END DO
1504	!
1505	#if defined key_top
1506	DO jk = 2, jkmax
1507	pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
1508	END DO
1509	#endif
1510	!
1511	END DO
1512	END DO
1513
1514	ENDIF
1515
1516	IF ( ll_asmiau ) THEN
1517
1518	!--------------------------------------------------------------------
1519	! Incremental Analysis Updating
1520	!--------------------------------------------------------------------
1521
1522	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1523
1524	it = kt - nit000 + 1
1525	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1526	! note this is not a tendency so should not be divided by rdt
1527
1528	IF(lwp) THEN
1529	WRITE(numout,*)
1530	IF ( ln_spco2inc ) THEN
1531	WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
1532	& kt,' with IAU weight = ', pwgtiau(it)
1533	ELSE IF ( ln_sfco2inc ) THEN
1534	WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
1535	& kt,' with IAU weight = ', pwgtiau(it)
1536	ENDIF
1537	WRITE(numout,*) '~~~~~~~~~~~~'
1538	ENDIF
1539
1540	! Update the biogeochemical variables
1541	! Add directly to trn and trb, rather than to tra, because tra gets
1542	! reset to zero at the start of trc_stp, called after this routine
1543	#if defined key_medusa && defined key_foam_medusa
1544	WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1545	& trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1546	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1547	& pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1548	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1549	& pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1550	END WHERE
1551	#elif defined key_hadocc
1552	WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1553	& trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1554	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1555	& pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1556	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1557	& pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1558	END WHERE
1559	#endif
1560
1561	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1562	! which is called at end of model run
1563
1564	ENDIF
1565
1566	ELSEIF ( ll_asmdin ) THEN
1567
1568	!--------------------------------------------------------------------
1569	! Direct Initialization
1570	!--------------------------------------------------------------------
1571
1572	IF ( kt == nitdin_r ) THEN
1573
1574	neuler = 0 ! Force Euler forward step
1575
1576	! Initialize the now fields with the background + increment
1577	! Background currently is what the model is initialised with
1578	CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
1579	& ' Background state is taken from model rather than background file' )
1580	#if defined key_medusa && defined key_foam_medusa
1581	WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1582	& trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1583	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1584	& pco2_balinc(:,:,:,jp_msa0:jp_msa1)
1585	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1586	END WHERE
1587	#elif defined key_hadocc
1588	WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1589	& trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1590	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1591	& pco2_balinc(:,:,:,jp_had0:jp_had1)
1592	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1593	END WHERE
1594	#endif
1595
1596	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1597	! which is called at end of model run
1598	ENDIF
1599	!
1600	ENDIF
1601	!
1602	END SUBROUTINE pco2_asm_inc
1603
1604	!!===========================================================================
1605	!!===========================================================================
1606	!!===========================================================================
1607
1608	SUBROUTINE ph_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
1609	& ll_trainc, pt_bkginc, ps_bkginc )
1610	!!------------------------------------------------------------------------
1611	!! * ROUTINE ph_asm_inc *
1612	!!
1613	!! ** Purpose : Apply the pH assimilation increments.
1614	!!
1615	!! ** Method : Calculate increments to DIC and alkalinity from pH
1616	!! Use a similar approach to the pCO2 scheme
1617	!! Direct initialization or Incremental Analysis Updating.
1618	!!
1619	!! ** Action :
1620	!!------------------------------------------------------------------------
1621	INTEGER, INTENT(IN) :: kt ! Current time step
1622	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1623	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1624	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1625	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1626	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
1627	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
1628	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
1629
1630	REAL(wp) :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation
1631	REAL(wp) :: DIC_IN, ALK_IN ! DIC/alk in pH calculation
1632	REAL(wp) :: PH_OUT ! pH from pH calculation
1633	REAL(wp) :: ph_bkg ! pH from pH calculation
1634	REAL(wp) :: ph_dic ! pH from pH calculation
1635	REAL(wp) :: ph_alk ! pH from pH calculation
1636	REAL(wp) :: dph_ddic ! Change in pH wrt DIC
1637	REAL(wp) :: dph_dalk ! Change in pH wrt alk
1638	REAL(wp) :: weight ! Increment weighting
1639	REAL(wp) :: zincwgt ! IAU weight for current time step
1640	REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp ! DIC background
1641	REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp ! Alkalinity background
1642	REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp ! DIN background
1643	REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp ! Silicate background
1644	REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp ! Temperature background
1645	REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp ! Salinity background
1646	REAL(wp), DIMENSION(19) :: dummy_out ! Dummy variables
1647	INTEGER :: ji, jj, jk, jx ! Loop counters
1648	INTEGER :: it ! Index
1649	!!------------------------------------------------------------------------
1650
1651	#if ! (defined key_medusa && defined key_foam_medusa)
1652	CALL ctl_stop( ' pH balancing only implemented for MEDUSA' )
1653	#else
1654
1655	IF ( kt <= nit000 ) THEN
1656
1657	!----------------------------------------------------------------------
1658	! DIC and alkalinity balancing
1659	!----------------------------------------------------------------------
1660
1661	! Account for physics assimilation if required
1662	IF ( ll_trainc ) THEN
1663	tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:)
1664	sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:)
1665	ELSE
1666	tem_bkg_temp(:,:,:) = tsn(:,:,:,1)
1667	sal_bkg_temp(:,:,:) = tsn(:,:,:,2)
1668	ENDIF
1669
1670	! Account for phytoplankton balancing if required
1671	IF ( ln_phytobal ) THEN
1672	dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic)
1673	alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk)
1674	din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin)
1675	sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil)
1676	ELSE
1677	dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic)
1678	alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk)
1679	din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin)
1680	sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil)
1681	ENDIF
1682
1683	! Account for pCO2 balancing if required
1684	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1685	dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic)
1686	alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk)
1687	ENDIF
1688
1689	! Loop over grid points and calculate dpH/dDIC and dpH/dAlk
1690	! This requires three calls to the MOCSY carbonate package
1691	! One to find pH at background DIC and alkalinity
1692	! One to find pH when DIC is increased by 10
1693	! One to find pH when alkalinity is increased by 10
1694	! Then calculate the gradients
1695
1696	DO jk = 1, jpk
1697	DO jj = 1, jpj
1698	DO ji = 1, jpi
1699
1700	IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN
1701
1702	DO jx = 1, 3
1703	IF ( jx == 1 ) THEN
1704	DIC_IN = dic_bkg_temp(ji,jj,jk)
1705	ALK_IN = alk_bkg_temp(ji,jj,jk)
1706	ELSE IF ( jx == 2 ) THEN
1707	DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch
1708	ALK_IN = alk_bkg_temp(ji,jj,jk)
1709	ELSE IF ( jx == 3 ) THEN
1710	DIC_IN = dic_bkg_temp(ji,jj,jk)
1711	ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch
1712	ENDIF
1713
1714	CALL mocsy_interface(tsn(ji,jj,jk,jp_tem), & ! IN: temperature
1715	& tsn(ji,jj,jk,jp_sal), & ! salinity
1716	& ALK_IN, & ! alkalinity
1717	& DIC_IN, & ! DIC
1718	& sil_bkg_temp(ji,jj,jk), & ! silicate
1719	& din_bkg_temp(ji,jj,jk)/16.0, & ! phosphate (Redfield from DIN)
1720	& 1.0, & ! atmospheric pressure (dummy)
1721	& fsdept(ji,jj,jk), & ! depth
1722	& gphit(ji,jj), & ! latitude
1723	& 1.0, & ! gas transfer (dummy)
1724	& 1.0, & ! atmospheric xCO2 (dummy)
1725	& 1, & ! number of points
1726	& PH_OUT, & ! OUT: pH
1727	& dummy_out(1), dummy_out(2), & ! stuff we don't care about
1728	& dummy_out(3), dummy_out(4), &
1729	& dummy_out(5), dummy_out(6), &
1730	& dummy_out(7), dummy_out(8), &
1731	& dummy_out(9), dummy_out(10), &
1732	& dummy_out(11), dummy_out(12), &
1733	& dummy_out(13), dummy_out(14), &
1734	& dummy_out(15), dummy_out(16), &
1735	& dummy_out(17), dummy_out(18), &
1736	& dummy_out(19) )
1737
1738	IF ( jx == 1 ) THEN
1739	ph_bkg = PH_OUT
1740	ELSE IF ( jx == 2 ) THEN
1741	ph_dic = PH_OUT
1742	ELSE IF ( jx == 3 ) THEN
1743	ph_alk = PH_OUT
1744	ENDIF
1745	END DO
1746
1747	dph_ddic = (ph_dic - ph_bkg) / zsearch
1748	dph_dalk = (ph_alk - ph_bkg) / zsearch
1749	weight = pph_bkginc(ji,jj,jk) / ( (dph_ddicdph_ddic) + (dph_dalkdph_dalk) )
1750
1751	ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic
1752	ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk
1753
1754	ENDIF
1755
1756	END DO
1757	END DO
1758	END DO
1759
1760	ENDIF
1761
1762	IF ( ll_asmiau ) THEN
1763
1764	!--------------------------------------------------------------------
1765	! Incremental Analysis Updating
1766	!--------------------------------------------------------------------
1767
1768	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1769
1770	it = kt - nit000 + 1
1771	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1772	! note this is not a tendency so should not be divided by rdt
1773
1774	IF(lwp) THEN
1775	WRITE(numout,*)
1776	WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', &
1777	& kt,' with IAU weight = ', pwgtiau(it)
1778	WRITE(numout,*) '~~~~~~~~~~~~'
1779	ENDIF
1780
1781	! Update the biogeochemical variables
1782	! Add directly to trn and trb, rather than to tra, because tra gets
1783	! reset to zero at the start of trc_stp, called after this routine
1784	WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1785	& trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1786	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1787	& ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1788	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1789	& ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1790	END WHERE
1791
1792	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1793	! which is called at end of model run
1794
1795	ENDIF
1796
1797	ELSEIF ( ll_asmdin ) THEN
1798
1799	!--------------------------------------------------------------------
1800	! Direct Initialization
1801	!--------------------------------------------------------------------
1802
1803	IF ( kt == nitdin_r ) THEN
1804
1805	neuler = 0 ! Force Euler forward step
1806
1807	! Initialize the now fields with the background + increment
1808	! Background currently is what the model is initialised with
1809	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pH assimilation', &
1810	& ' Background state is taken from model rather than background file' )
1811	WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1812	& trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1813	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1814	& ph_balinc(:,:,:,jp_msa0:jp_msa1)
1815	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1816	END WHERE
1817
1818	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1819	! which is called at end of model run
1820	ENDIF
1821	!
1822	ENDIF
1823	#endif
1824	!
1825	END SUBROUTINE ph_asm_inc
1826
1827	!!===========================================================================
1828	!!===========================================================================
1829	!!===========================================================================
1830
1831	SUBROUTINE bgc3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1832	!!----------------------------------------------------------------------
1833	!! * ROUTINE dyn_asm_inc *
1834	!!
1835	!! ** Purpose : Apply generic 3D biogeochemistry assimilation increments.
1836	!!
1837	!! ** Method : Direct initialization or Incremental Analysis Updating.
1838	!!
1839	!! ** Action :
1840	!!----------------------------------------------------------------------
1841	INTEGER, INTENT(IN) :: kt ! Current time step
1842	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1843	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1844	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1845	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1846	!
1847	INTEGER :: jk ! Loop counter
1848	INTEGER :: it ! Index
1849	REAL(wp) :: zincwgt ! IAU weight for current time step
1850	REAL(wp) :: zincper ! IAU interval in seconds
1851	!!----------------------------------------------------------------------
1852
1853	IF ( kt <= nit000 ) THEN
1854
1855	!----------------------------------------------------------------------
1856	! Remove any other balancing increments
1857	!----------------------------------------------------------------------
1858
1859	IF ( ln_pno3inc ) THEN
1860	#if defined key_hadocc \|\| ( defined key_medusa && defined key_foam_medusa )
1861	#if defined key_hadocc
1862	it = jp_had_nut
1863	#elif defined key_medusa && defined key_foam_medusa
1864	it = jpdin
1865	#endif
1866	IF ( ln_phytobal ) THEN
1867	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1868	ENDIF
1869	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1870	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1871	ENDIF
1872	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1873	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1874	ENDIF
1875	IF ( ln_pphinc ) THEN
1876	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1877	ENDIF
1878	#else
1879	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1880	#endif
1881	ENDIF
1882
1883	IF ( ln_psi4inc ) THEN
1884	#if defined key_medusa && defined key_foam_medusa
1885	it = jpsil
1886	IF ( ln_phytobal ) THEN
1887	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1888	ENDIF
1889	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1890	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1891	ENDIF
1892	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1893	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1894	ENDIF
1895	IF ( ln_pphinc ) THEN
1896	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1897	ENDIF
1898	#else
1899	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1900	#endif
1901	ENDIF
1902
1903	IF ( ln_pdicinc ) THEN
1904	#if defined key_hadocc \|\| ( defined key_medusa && defined key_foam_medusa )
1905	#if defined key_hadocc
1906	it = jp_had_dic
1907	#elif defined key_medusa && defined key_foam_medusa
1908	it = jpdic
1909	#endif
1910	IF ( ln_phytobal ) THEN
1911	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1912	ENDIF
1913	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1914	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1915	ENDIF
1916	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1917	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1918	ENDIF
1919	IF ( ln_pphinc ) THEN
1920	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1921	ENDIF
1922	#else
1923	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1924	#endif
1925	ENDIF
1926
1927	IF ( ln_palkinc ) THEN
1928	#if defined key_hadocc \|\| ( defined key_medusa && defined key_foam_medusa )
1929	#if defined key_hadocc
1930	it = jp_had_alk
1931	#elif defined key_medusa && defined key_foam_medusa
1932	it = jpalk
1933	#endif
1934	IF ( ln_phytobal ) THEN
1935	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1936	ENDIF
1937	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1938	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1939	ENDIF
1940	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1941	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1942	ENDIF
1943	IF ( ln_pphinc ) THEN
1944	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1945	ENDIF
1946	#else
1947	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1948	#endif
1949	ENDIF
1950
1951	IF ( ln_po2inc ) THEN
1952	#if defined key_medusa && defined key_foam_medusa
1953	it = jpoxy
1954	IF ( ln_phytobal ) THEN
1955	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1956	ENDIF
1957	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1958	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1959	ENDIF
1960	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1961	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1962	ENDIF
1963	IF ( ln_pphinc ) THEN
1964	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1965	ENDIF
1966	#else
1967	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1968	#endif
1969	ENDIF
1970
1971	ENDIF
1972
1973	IF ( ll_asmiau ) THEN
1974
1975	!--------------------------------------------------------------------
1976	! Incremental Analysis Updating
1977	!--------------------------------------------------------------------
1978
1979	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1980
1981	it = kt - nit000 + 1
1982	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1983	! note this is not a tendency so should not be divided by rdt
1984
1985	IF(lwp) THEN
1986	WRITE(numout,*)
1987	WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', &
1988	& kt,' with IAU weight = ', pwgtiau(it)
1989	WRITE(numout,*) '~~~~~~~~~~~~'
1990	ENDIF
1991
1992	! Update the 3D BGC variables
1993	! Add directly to trn and trb, rather than to tra, because tra gets
1994	! reset to zero at the start of trc_stp, called after this routine
1995	! Don't apply increments if they'll take concentrations negative
1996
1997	IF ( ln_pno3inc ) THEN
1998	#if defined key_hadocc
1999	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2000	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2001	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
2002	trb(:,:,:,jp_had_nut) = trb(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
2003	END WHERE
2004	#elif defined key_medusa && defined key_foam_medusa
2005	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2006	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2007	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
2008	trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
2009	END WHERE
2010	#else
2011	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2012	#endif
2013	ENDIF
2014
2015	IF ( ln_psi4inc ) THEN
2016	#if defined key_medusa && defined key_foam_medusa
2017	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2018	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
2019	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
2020	trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
2021	END WHERE
2022	#else
2023	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2024	#endif
2025	ENDIF
2026
2027	IF ( ln_pdicinc ) THEN
2028	#if defined key_hadocc
2029	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2030	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2031	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
2032	trb(:,:,:,jp_had_dic) = trb(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
2033	END WHERE
2034	#elif defined key_medusa && defined key_foam_medusa
2035	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2036	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2037	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
2038	trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
2039	END WHERE
2040	#else
2041	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2042	#endif
2043	ENDIF
2044
2045	IF ( ln_palkinc ) THEN
2046	#if defined key_hadocc
2047	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2048	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2049	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
2050	trb(:,:,:,jp_had_alk) = trb(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
2051	END WHERE
2052	#elif defined key_medusa && defined key_foam_medusa
2053	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2054	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2055	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
2056	trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
2057	END WHERE
2058	#else
2059	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2060	#endif
2061	ENDIF
2062
2063	IF ( ln_po2inc ) THEN
2064	#if defined key_medusa && defined key_foam_medusa
2065	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2066	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
2067	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
2068	trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
2069	END WHERE
2070	#else
2071	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2072	#endif
2073	ENDIF
2074
2075	IF ( kt == nitiaufin_r ) THEN
2076	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
2077	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
2078	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
2079	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
2080	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
2081	ENDIF
2082
2083	ENDIF
2084
2085	ELSEIF ( ll_asmdin ) THEN
2086
2087	!--------------------------------------------------------------------
2088	! Direct Initialization
2089	!--------------------------------------------------------------------
2090
2091	IF ( kt == nitdin_r ) THEN
2092
2093	neuler = 0 ! Force Euler forward step
2094
2095	! Initialize the now fields with the background + increment
2096	! Background currently is what the model is initialised with
2097	#if defined key_hadocc
2098	CALL ctl_warn( ' Doing direct initialisation of HadOCC with 3D BGC assimilation', &
2099	& ' Background state is taken from model rather than background file' )
2100	#elif defined key_medusa && defined key_foam_medusa
2101	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with 3D BGC assimilation', &
2102	& ' Background state is taken from model rather than background file' )
2103	#endif
2104
2105	IF ( ln_pno3inc ) THEN
2106	#if defined key_hadocc
2107	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2108	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) > 0.0_wp )
2109	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:)
2110	trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut)
2111	END WHERE
2112	#elif defined key_medusa && defined key_foam_medusa
2113	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2114	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) > 0.0_wp )
2115	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:)
2116	trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
2117	END WHERE
2118	#else
2119	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2120	#endif
2121	ENDIF
2122
2123	IF ( ln_psi4inc ) THEN
2124	#if defined key_medusa && defined key_foam_medusa
2125	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2126	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) > 0.0_wp )
2127	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:)
2128	trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
2129	END WHERE
2130	#else
2131	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2132	#endif
2133	ENDIF
2134
2135	IF ( ln_pdicinc ) THEN
2136	#if defined key_hadocc
2137	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2138	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) > 0.0_wp )
2139	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:)
2140	trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic)
2141	END WHERE
2142	#elif defined key_medusa && defined key_foam_medusa
2143	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2144	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) > 0.0_wp )
2145	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:)
2146	trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
2147	END WHERE
2148	#else
2149	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2150	#endif
2151	ENDIF
2152
2153	IF ( ln_palkinc ) THEN
2154	#if defined key_hadocc
2155	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2156	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) > 0.0_wp )
2157	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:)
2158	trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk)
2159	END WHERE
2160	#elif defined key_medusa && defined key_foam_medusa
2161	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2162	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) > 0.0_wp )
2163	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:)
2164	trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
2165	END WHERE
2166	#else
2167	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2168	#endif
2169	ENDIF
2170
2171	IF ( ln_po2inc ) THEN
2172	#if defined key_medusa && defined key_foam_medusa
2173	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2174	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) > 0.0_wp )
2175	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:)
2176	trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
2177	END WHERE
2178	#else
2179	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2180	#endif
2181	ENDIF
2182
2183	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
2184	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
2185	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
2186	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
2187	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
2188	ENDIF
2189	!
2190	ENDIF
2191	!
2192	END SUBROUTINE bgc3d_asm_inc
2193
2194	!!===========================================================================
2195
2196	END MODULE asmbgc

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

Download in other formats: