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

Context Navigation

source: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 13191

Last change on this file since 13191 was 13191, checked in by jwhile, 4 years ago
Updates for 1d runnig
File size: 102.1 KB

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

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

Download in other formats: