Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

asmbgc.F90 @ 10774

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 10774

Download in other formats: