Context Navigation

asmbgc.F90 @ 15337

Last change on this file since 15337 was 15337, checked in by andmirek, 3 years ago
Ticket GO29: fix tests 7 and 8
File size: 102.6 KB

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

Download in other formats: