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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

asmbgc.F90 in branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 10574

Last change on this file since 10574 was 10574, checked in by dford, 5 years ago
Merge in functionality added to GO6 at r10149.
File size: 102.1 KB

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

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

Download in other formats: