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 @ 10661

Last change on this file since 10661 was 10661, checked in by dford, 5 years ago
Add option to apply phosphate increments.
File size: 141.9 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_fabm' : ERSEM model coupled via FABM
14	!! 'key_karaml' : Kara mixed layer depth
15	!!---------------------------------------------------------------------------
16	!! asm_bgc_check_options : check bgc assimilation options
17	!! asm_bgc_init_incs : initialise bgc increments
18	!! asm_bgc_init_bkg : initialise bgc background
19	!! asm_bgc_bal_wri : write out bgc balancing increments
20	!! asm_bgc_bkg_wri : write out bgc background
21	!! phyto2d_asm_inc : apply the ocean colour increments
22	!! phyto3d_asm_inc : apply the 3D phytoplankton increments
23	!! pco2_asm_inc : apply the pCO2/fCO2 increments
24	!! ph_asm_inc : apply the pH increments
25	!! bgc3d_asm_inc : apply the generic 3D BGC increments
26	!!---------------------------------------------------------------------------
27	USE par_kind, ONLY: & ! kind parameters
28	& wp
29	USE par_oce, ONLY: & ! domain array sizes
30	& jpi, &
31	& jpj, &
32	& jpk
33	USE iom ! i/o
34	USE oce, ONLY: & ! active tracer variables
35	& tsn
36	USE zdfmxl, ONLY : & ! mixed layer depth
37	#if defined key_karaml
38	& hmld_kara, &
39	& ln_kara, &
40	#endif
41	& hmld, &
42	& hmld_tref, &
43	& hmlp, &
44	& hmlpt
45	USE asmpar, ONLY: & ! assimilation parameters
46	& c_asmbkg, &
47	& c_asmbal, &
48	& nitbkg_r, &
49	& nitavgbkg_r, &
50	& nitdin_r, &
51	& nitiaustr_r, &
52	& nitiaufin_r
53	USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA
54	& asm_phyto2d_bal_medusa
55	USE asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC
56	& asm_phyto2d_bal_hadocc
57	USE asmphyto2dbal_ersem, ONLY: & ! phyto2d balancing for ERSEM
58	& asm_phyto2d_bal_ersem
59	USE asmpco2bal, ONLY: & ! pCO2 balancing
60	& asm_pco2_bal
61	#if defined key_top
62	USE trc, ONLY: & ! passive tracer variables
63	& trn, &
64	& trb, &
65	& nittrc000
66	USE par_trc, ONLY: & ! passive tracer parameters
67	& jptra
68	#endif
69	#if defined key_medusa
70	USE sms_medusa ! MEDUSA parameters
71	USE par_medusa ! MEDUSA parameters
72	USE mocsy_wrapper, ONLY: & ! MEDUSA carbonate system
73	& mocsy_interface
74	USE mocsy_mainmod, ONLY: & ! fCO2 to pCO2 conversion
75	& f2pCO2
76	USE sms_medusa, ONLY: & ! MEDUSA diagnostic variables
77	& pgrow_avg, &
78	& ploss_avg, &
79	& phyt_avg, &
80	& mld_max
81	#elif defined key_hadocc
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	#elif defined key_fabm
93	USE par_fabm ! FABM-ERSEM parameters
94	USE trc, ONLY: & ! FABM-ERSEM diagnostic variables
95	& pgrow_avg, &
96	& ploss_avg, &
97	& phyt_avg, &
98	& mld_max
99	#endif
100
101	IMPLICIT NONE
102	PRIVATE
103
104	PUBLIC asm_bgc_check_options ! called by asm_inc_init in asminc.F90
105	PUBLIC asm_bgc_init_incs ! called by asm_inc_init in asminc.F90
106	PRIVATE asm_bgc_read_incs_2d ! called by asm_bgc_init_incs
107	PRIVATE asm_bgc_read_incs_3d ! called by asm_bgc_init_incs
108	PUBLIC asm_bgc_init_bkg ! called by asm_inc_init in asminc.F90
109	PUBLIC asm_bgc_bal_wri ! called by nemo_gcm in nemogcm.F90
110	PUBLIC asm_bgc_bkg_alloc ! called by asm_bkg_wri in asmbkg.F90
111	PUBLIC asm_bgc_bkg_tavg ! called by asm_bkg_wri in asmbkg.F90
112	PUBLIC asm_bgc_bkg_wri ! called by asm_bkg_wri in asmbkg.F90
113	PUBLIC phyto2d_asm_inc ! called by bgc_asm_inc in asminc.F90
114	PUBLIC phyto3d_asm_inc ! called by bgc_asm_inc in asminc.F90
115	PUBLIC pco2_asm_inc ! called by bgc_asm_inc in asminc.F90
116	PUBLIC ph_asm_inc ! called by bgc_asm_inc in asminc.F90
117	PUBLIC bgc3d_asm_inc ! called by bgc_asm_inc in asminc.F90
118
119	LOGICAL, PUBLIC :: ln_balwri = .FALSE. !: No output of balancing incs
120	LOGICAL, PUBLIC :: ln_phytobal = .FALSE. !: No phytoplankton balancing
121	LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE. !: No surface total log10(chlorophyll) increment
122	LOGICAL, PUBLIC :: ln_slchldiainc = .FALSE. !: No surface diatom log10(chlorophyll) increment
123	LOGICAL, PUBLIC :: ln_slchlnoninc = .FALSE. !: No surface non-diatom log10(chlorophyll) increment
124	LOGICAL, PUBLIC :: ln_slchlnaninc = .FALSE. !: No surface nanophyto log10(chlorophyll) increment
125	LOGICAL, PUBLIC :: ln_slchlpicinc = .FALSE. !: No surface picophyto log10(chlorophyll) increment
126	LOGICAL, PUBLIC :: ln_slchldininc = .FALSE. !: No surface dinoflag log10(chlorophyll) increment
127	LOGICAL, PUBLIC :: ln_schltotinc = .FALSE. !: No surface total chlorophyll increment
128	LOGICAL, PUBLIC :: ln_slphytotinc = .FALSE. !: No surface total log10(phyto C) increment
129	LOGICAL, PUBLIC :: ln_slphydiainc = .FALSE. !: No surface diatom log10(phyto C) increment
130	LOGICAL, PUBLIC :: ln_slphynoninc = .FALSE. !: No surface non-diatom log10(phyto C) increment
131	LOGICAL, PUBLIC :: ln_spco2inc = .FALSE. !: No surface pCO2 increment
132	LOGICAL, PUBLIC :: ln_sfco2inc = .FALSE. !: No surface fCO2 increment
133	LOGICAL, PUBLIC :: ln_plchltotinc = .FALSE. !: No profile total log10(chlorophyll) increment
134	LOGICAL, PUBLIC :: ln_pchltotinc = .FALSE. !: No profile total chlorophyll increment
135	LOGICAL, PUBLIC :: ln_pno3inc = .FALSE. !: No profile nitrate increment
136	LOGICAL, PUBLIC :: ln_psi4inc = .FALSE. !: No profile silicate increment
137	LOGICAL, PUBLIC :: ln_ppo4inc = .FALSE. !: No profile phosphate increment
138	LOGICAL, PUBLIC :: ln_pdicinc = .FALSE. !: No profile dissolved inorganic carbon increment
139	LOGICAL, PUBLIC :: ln_palkinc = .FALSE. !: No profile alkalinity increment
140	LOGICAL, PUBLIC :: ln_pphinc = .FALSE. !: No profile pH increment
141	LOGICAL, PUBLIC :: ln_po2inc = .FALSE. !: No profile oxygen increment
142
143	INTEGER, PUBLIC :: mld_choice_bgc = 1 !: choice of mld for bgc assimilation
144	! 1) hmld - Turbocline/mixing depth
145	! [W points]
146	! 2) hmlp - Density criterion
147	! (0.01 kg/m^3 change from 10m)
148	! [W points]
149	! 3) hmld_kara - Kara MLD
150	! [Interpolated]
151	! 4) hmld_tref - Temperature criterion
152	! (0.2 K change from surface)
153	! [T points]
154	! 5) hmlpt - Density criterion
155	! (0.01 kg/m^3 change from 10m)
156	! [T points]
157
158	REAL(wp), PUBLIC :: rn_maxchlinc = -999.0 !: maximum absolute non-log chlorophyll inc
159	! <= 0 no maximum applied (switch turned off)
160	! > 0 absolute chl inc capped at this value
161
162	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchltot_bkginc ! slchltot inc
163	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchldia_bkginc ! slchldia inc
164	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchlnon_bkginc ! slchlnon inc
165	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchlnan_bkginc ! slchlnan inc
166	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchlpic_bkginc ! slchlpic inc
167	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchldin_bkginc ! slchldin inc
168	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: schltot_bkginc ! schltot inc
169	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphytot_bkginc ! slphytot inc
170	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphydia_bkginc ! slphydia inc
171	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphynon_bkginc ! slphynon inc
172	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sfco2_bkginc ! sfco2 inc
173	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: spco2_bkginc ! spco2 inc
174	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: plchltot_bkginc ! plchltot inc
175	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pchltot_bkginc ! pchltot inc
176	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pno3_bkginc ! pno3 inc
177	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: psi4_bkginc ! psi4 inc
178	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ppo4_bkginc ! ppo4 inc
179	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pdic_bkginc ! pdic inc
180	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: palk_bkginc ! palk inc
181	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pph_bkginc ! pph inc
182	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: po2_bkginc ! po2 inc
183	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto2d_balinc ! Balancing incs from ocean colour
184	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto3d_balinc ! Balancing incs from p(l)chltot
185	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc ! Balancing incs from spco2/sfco2
186	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ph_balinc ! Balancing incs from pph
187
188	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pgrow_avg_bkg ! Background phyto growth
189	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ploss_avg_bkg ! Background phyto loss
190	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: phyt_avg_bkg ! Background phyto conc
191	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mld_max_bkg ! Background max MLD
192	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg ! Background tracer state
193	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: chl_bkg ! Background chl
194	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cchl_p_bkg ! Background c:chl
195	#if defined key_fabm
196	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: totalk_bkg ! Background total alkalinity
197	#endif
198
199	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: pgrow_avg_tavg ! Time average of pgrow_avg
200	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: ploss_avg_tavg ! Time average of ploss_avg
201	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: phyt_avg_tavg ! Time average of phyt_avg
202	REAL(wp), DIMENSION(:,:,:,:), SAVE, ALLOCATABLE :: trn_tavg ! Time average of trn
203	#if defined key_hadocc
204	REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE :: HADOCC_CHL_tavg ! Time average of HADOCC_CHL
205	REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE :: cchl_p_tavg ! Time average of cchl_p
206	#elif defined key_fabm
207	REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE :: totalk_tavg ! Time average of O3_TA
208	#endif
209
210	# include "domzgr_substitute.h90"
211
212	CONTAINS
213
214	SUBROUTINE asm_bgc_check_options
215	!!------------------------------------------------------------------------
216	!! * ROUTINE asm_bgc_check_options *
217	!!
218	!! ** Purpose : check validity of biogeochemical assimilation options
219	!!
220	!! ** Method : check settings of logicals
221	!!
222	!! ** Action : call ctl_stop/ctl_warn if required
223	!!
224	!! References : asm_inc_init
225	!!------------------------------------------------------------------------
226
227	#if ! defined key_top \|\| ( ! defined key_hadocc && ! defined key_medusa && ! defined key_fabm )
228	CALL ctl_stop( ' Attempting to assimilate biogeochemical observations', &
229	& ' but no compatible biogeochemical model is available' )
230	#endif
231
232	#if defined key_hadocc
233	IF ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slchlnaninc .OR. &
234	& ln_slchldininc .OR. ln_slchlpicinc .OR. ln_slphydiainc .OR. &
235	& ln_slphynoninc .OR. ln_psi4inc .OR. ln_pphinc .OR. &
236	& ln_po2inc .OR. ln_ppo4inc ) THEN
237	CALL ctl_stop( ' Cannot assimilate PFTs, Si4, PO4, pH or O2 into HadOCC' )
238	ENDIF
239	#endif
240
241	#if defined key_medusa
242	IF ( .NOT. ln_foam_medusa ) THEN
243	CALL ctl_stop( ' MEDUSA data assimilation options not turned on: set ln_foam_medusa' )
244	ENDIF
245
246	IF ( ln_ppo4inc ) THEN
247	CALL ctl_stop( ' Cannot assimilate PO4 into MEDUSA' )
248	ENDIF
249	#endif
250
251	#if defined key_fabm
252	IF ( ln_pphinc ) THEN
253	CALL ctl_stop( ' Cannot currently assimilate pH into FABM-ERSEM' )
254	ENDIF
255	#endif
256
257	IF ( ( ln_phytobal ).AND. &
258	& ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
259	& ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. &
260	& ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. &
261	& ( .NOT. ln_schltotinc ).AND.( .NOT. ln_slphytotinc ).AND. &
262	& ( .NOT. ln_slphydiainc ).AND.( .NOT. ln_slphynoninc ) ) THEN
263	CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', &
264	& ' if not assimilating ocean colour,', &
265	& ' so ln_phytobal will be set to .false.')
266	ln_phytobal = .FALSE.
267	ENDIF
268
269	IF ( ( ln_balwri ).AND. &
270	& ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
271	& ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. &
272	& ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. &
273	& ( .NOT. ln_schltotinc ).AND. &
274	& ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
275	& ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. &
276	& ( .NOT. ln_pchltotinc ).AND.( .NOT. ln_pphinc ).AND. &
277	& ( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ) ) THEN
278	CALL ctl_warn( ' Balancing increments not being calculated', &
279	& ' so ln_balwri will be set to .false.')
280	ln_balwri = .FALSE.
281	ENDIF
282
283	IF ( ln_spco2inc .AND. ln_sfco2inc ) THEN
284	CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
285	ENDIF
286
287	IF ( ln_slchltotinc .AND. ln_schltotinc ) THEN
288	CALL ctl_stop( ' Can only assimilate surface log10(chlorophyll) or chlorophyll, not both' )
289	ENDIF
290
291	IF ( ln_plchltotinc .AND. ln_pchltotinc ) THEN
292	CALL ctl_stop( ' Can only assimilate profile log10(chlorophyll) or chlorophyll, not both' )
293	ENDIF
294
295	IF ( ( ln_slchltotinc .OR. ln_schltotinc ) .AND. &
296	& ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slchlnaninc .OR. &
297	& ln_slchlpicinc .OR. ln_slchldininc ) ) THEN
298	CALL ctl_stop( ' Can only assimilate total or PFT surface chlorophyll, not both' )
299	ENDIF
300
301	IF ( ln_phytobal .AND. &
302	& ( ( ln_slchlnoninc .AND.( .NOT. ln_slchldiainc ) ).OR. &
303	#if defined key_fabm
304	& ( ( ( ln_slchldiainc .OR. ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc ) .AND. &
305	& ( ( .NOT. ln_slchldiainc ) .OR. ( .NOT. ln_slchlnaninc ) .OR. &
306	& ( .NOT. ln_slchlpicinc ) .OR. ( .NOT. ln_slchldininc ) ) ) .AND. &
307	& ( ln_slchldiainc .AND.( .NOT. ln_slchlnoninc ) ) ) ) ) THEN
308	#elif defined key_medusa
309	& ( ln_slchldiainc .AND.( .NOT. ln_slchlnoninc ) ) ) ) THEN
310	#endif
311	CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
312	& ' unless assimilating all model PFTs,')
313	ENDIF
314
315	IF ( ln_slphytotinc .AND. ( ln_slphydiainc .OR. ln_slphynoninc ) ) THEN
316	CALL ctl_stop( ' Can only assimilate total or PFT surface phytoplankton carbon, not both' )
317	ENDIF
318
319	IF ( ln_phytobal .AND. &
320	& ( ( ln_slphynoninc .AND.( .NOT. ln_slphydiainc ) ).OR. &
321	& ( ln_slphydiainc .AND.( .NOT. ln_slphynoninc ) ) ) ) THEN
322	CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
323	& ' unless assimilating all model PFTs,')
324	ENDIF
325
326	END SUBROUTINE asm_bgc_check_options
327
328	!!===========================================================================
329	!!===========================================================================
330	!!===========================================================================
331
332	SUBROUTINE asm_bgc_init_incs( knum )
333	!!------------------------------------------------------------------------
334	!! * ROUTINE asm_bgc_init_incs *
335	!!
336	!! ** Purpose : initialise bgc increments
337	!!
338	!! ** Method : allocate arrays and read increments from file
339	!!
340	!! ** Action : allocate arrays and read increments from file
341	!!
342	!! References : asm_inc_init
343	!!------------------------------------------------------------------------
344	!!
345	INTEGER, INTENT(in ) :: knum ! i/o unit of increments file
346	!!
347	!!------------------------------------------------------------------------
348
349	! Allocate and read increments
350
351	IF ( ln_slchltotinc ) THEN
352	ALLOCATE( slchltot_bkginc(jpi,jpj) )
353	CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc )
354	ENDIF
355
356	IF ( ln_slchldiainc ) THEN
357	ALLOCATE( slchldia_bkginc(jpi,jpj) )
358	CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc )
359	ENDIF
360
361	IF ( ln_slchlnoninc ) THEN
362	ALLOCATE( slchlnon_bkginc(jpi,jpj) )
363	CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc )
364	ENDIF
365
366	IF ( ln_slchlnaninc ) THEN
367	ALLOCATE( slchlnan_bkginc(jpi,jpj) )
368	CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnan', slchlnan_bkginc )
369	ENDIF
370
371	IF ( ln_slchlpicinc ) THEN
372	ALLOCATE( slchlpic_bkginc(jpi,jpj) )
373	CALL asm_bgc_read_incs_2d( knum, 'bckinslchlpic', slchlpic_bkginc )
374	ENDIF
375
376	IF ( ln_slchldininc ) THEN
377	ALLOCATE( slchldin_bkginc(jpi,jpj) )
378	CALL asm_bgc_read_incs_2d( knum, 'bckinslchldin', slchldin_bkginc )
379	ENDIF
380
381	IF ( ln_schltotinc ) THEN
382	ALLOCATE( schltot_bkginc(jpi,jpj) )
383	CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc )
384	ENDIF
385
386	IF ( ln_slphytotinc ) THEN
387	ALLOCATE( slphytot_bkginc(jpi,jpj) )
388	CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc )
389	ENDIF
390
391	IF ( ln_slphydiainc ) THEN
392	ALLOCATE( slphydia_bkginc(jpi,jpj) )
393	CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc )
394	ENDIF
395
396	IF ( ln_slphynoninc ) THEN
397	ALLOCATE( slphynon_bkginc(jpi,jpj) )
398	CALL asm_bgc_read_incs_2d( knum, 'bckinslphynon', slphynon_bkginc )
399	ENDIF
400
401	IF ( ln_sfco2inc ) THEN
402	ALLOCATE( sfco2_bkginc(jpi,jpj) )
403	CALL asm_bgc_read_incs_2d( knum, 'bckinsfco2', sfco2_bkginc )
404	ENDIF
405
406	IF ( ln_spco2inc ) THEN
407	ALLOCATE( spco2_bkginc(jpi,jpj) )
408	CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc )
409	ENDIF
410
411	IF ( ln_plchltotinc ) THEN
412	ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) )
413	CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc )
414	ENDIF
415
416	IF ( ln_pchltotinc ) THEN
417	ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) )
418	CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc )
419	ENDIF
420
421	IF ( ln_pno3inc ) THEN
422	ALLOCATE( pno3_bkginc(jpi,jpj,jpk) )
423	CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc )
424	ENDIF
425
426	IF ( ln_psi4inc ) THEN
427	ALLOCATE( psi4_bkginc(jpi,jpj,jpk) )
428	CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc )
429	ENDIF
430
431	IF ( ln_ppo4inc ) THEN
432	ALLOCATE( ppo4_bkginc(jpi,jpj,jpk) )
433	CALL asm_bgc_read_incs_3d( knum, 'bckinppo4', ppo4_bkginc )
434	ENDIF
435
436	IF ( ln_pdicinc ) THEN
437	ALLOCATE( pdic_bkginc(jpi,jpj,jpk) )
438	CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc )
439	ENDIF
440
441	IF ( ln_palkinc ) THEN
442	ALLOCATE( palk_bkginc(jpi,jpj,jpk) )
443	CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc )
444	ENDIF
445
446	IF ( ln_pphinc ) THEN
447	ALLOCATE( pph_bkginc(jpi,jpj,jpk) )
448	CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc )
449	ENDIF
450
451	IF ( ln_po2inc ) THEN
452	ALLOCATE( po2_bkginc(jpi,jpj,jpk) )
453	CALL asm_bgc_read_incs_3d( knum, 'bckinpo2', po2_bkginc )
454	ENDIF
455
456	! Allocate balancing increments
457
458	IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
459	& ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
460	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
461	& ln_slphynoninc ) THEN
462	#if defined key_top
463	ALLOCATE( phyto2d_balinc(jpi,jpj,jpk,jptra) )
464	phyto2d_balinc(:,:,:,:) = 0.0
465	#else
466	CALL ctl_stop( ' key_top must be set for balancing increments' )
467	#endif
468	ENDIF
469
470	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
471	#if defined key_top
472	ALLOCATE( phyto3d_balinc(jpi,jpj,jpk,jptra) )
473	phyto3d_balinc(:,:,:,:) = 0.0
474	#else
475	CALL ctl_stop( ' key_top must be set for balancing increments' )
476	#endif
477	ENDIF
478
479	IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
480	#if defined key_top
481	ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) )
482	pco2_balinc(:,:,:,:) = 0.0
483	#else
484	CALL ctl_stop( ' key_top must be set for balancing increments' )
485	#endif
486	ENDIF
487
488	IF ( ln_pphinc ) THEN
489	#if defined key_top
490	ALLOCATE( ph_balinc(jpi,jpj,jpk,jptra) )
491	ph_balinc(:,:,:,:) = 0.0
492	#else
493	CALL ctl_stop( ' key_top must be set for balancing increments' )
494	#endif
495	ENDIF
496
497	END SUBROUTINE asm_bgc_init_incs
498
499	!!===========================================================================
500	!!===========================================================================
501	!!===========================================================================
502
503	SUBROUTINE asm_bgc_read_incs_2d( knum, cd_bgcname, p_incs )
504	!!------------------------------------------------------------------------
505	!! * ROUTINE asm_bgc_init_incs *
506	!!
507	!! ** Purpose : read 2d bgc increments
508	!!
509	!! ** Method : read increments from file
510	!!
511	!! ** Action : read increments from file
512	!!
513	!! References : asm_inc_init
514	!!------------------------------------------------------------------------
515	!!
516	INTEGER, INTENT(in ) :: knum ! i/o unit
517	CHARACTER(LEN=*), INTENT(in ) :: cd_bgcname ! variable
518	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: p_incs ! increments
519	!!
520	!!------------------------------------------------------------------------
521
522	! Initialise
523	p_incs(:,:) = 0.0
524
525	! read from file
526	CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 )
527
528	! Apply the masks
529	p_incs(:,:) = p_incs(:,:) * tmask(:,:,1)
530
531	! Set missing increments to 0.0 rather than 1e+20
532	! to allow for differences in masks
533	WHERE( ABS( p_incs(:,:) ) > 1.0e+10 ) p_incs(:,:) = 0.0
534
535	END SUBROUTINE asm_bgc_read_incs_2d
536
537	!!===========================================================================
538	!!===========================================================================
539	!!===========================================================================
540
541	SUBROUTINE asm_bgc_read_incs_3d( knum, cd_bgcname, p_incs )
542	!!------------------------------------------------------------------------
543	!! * ROUTINE asm_bgc_init_incs *
544	!!
545	!! ** Purpose : read 3d bgc increments
546	!!
547	!! ** Method : read increments from file
548	!!
549	!! ** Action : read increments from file
550	!!
551	!! References : asm_inc_init
552	!!------------------------------------------------------------------------
553	!!
554	INTEGER, INTENT(in ) :: knum ! i/o unit
555	CHARACTER(LEN=*), INTENT(in ) :: cd_bgcname ! variable
556	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: p_incs ! increments
557	!!
558	!!------------------------------------------------------------------------
559
560	! Initialise
561	p_incs(:,:,:) = 0.0
562
563	! read from file
564	CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 )
565
566	! Apply the masks
567	p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:)
568
569	! Set missing increments to 0.0 rather than 1e+20
570	! to allow for differences in masks
571	WHERE( ABS( p_incs(:,:,:) ) > 1.0e+10 ) p_incs(:,:,:) = 0.0
572
573	END SUBROUTINE asm_bgc_read_incs_3d
574
575	!!===========================================================================
576	!!===========================================================================
577	!!===========================================================================
578
579	SUBROUTINE asm_bgc_init_bkg
580	!!------------------------------------------------------------------------
581	!! * ROUTINE asm_bgc_init_bkg *
582	!!
583	!! ** Purpose : initialise bgc background
584	!!
585	!! ** Method : allocate arrays and read background from file
586	!!
587	!! ** Action : allocate arrays and read background from file
588	!!
589	!! References : asm_inc_init
590	!!------------------------------------------------------------------------
591	!!
592	INTEGER :: inum ! i/o unit of background file
593	INTEGER :: jt ! loop counter
594	!!
595	!!------------------------------------------------------------------------
596
597	#if defined key_top && ( defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm )
598	IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
599	& ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
600	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
601	& ln_slphynoninc .OR. ln_plchltotinc .OR. ln_pchltotinc ) THEN
602
603	ALLOCATE( pgrow_avg_bkg(jpi,jpj) )
604	ALLOCATE( ploss_avg_bkg(jpi,jpj) )
605	ALLOCATE( phyt_avg_bkg(jpi,jpj) )
606	ALLOCATE( mld_max_bkg(jpi,jpj) )
607	ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
608	pgrow_avg_bkg(:,:) = 0.0
609	ploss_avg_bkg(:,:) = 0.0
610	phyt_avg_bkg(:,:) = 0.0
611	mld_max_bkg(:,:) = 0.0
612	tracer_bkg(:,:,:,:) = 0.0
613
614	#if defined key_hadocc
615	ALLOCATE( chl_bkg(jpi,jpj,jpk) )
616	ALLOCATE( cchl_p_bkg(jpi,jpj,jpk) )
617	chl_bkg(:,:,:) = 0.0
618	cchl_p_bkg(:,:,:) = 0.0
619	#elif defined key_fabm
620	ALLOCATE( totalk_bkg(jpi,jpj,jpk) )
621	totalk_bkg(:,:,:) = 0.0
622	#endif
623
624	!--------------------------------------------------------------------
625	! Read background variables for phytoplankton assimilation
626	! Some only required if performing balancing
627	!--------------------------------------------------------------------
628
629	CALL iom_open( c_asmbkg, inum )
630
631	#if defined key_hadocc
632	CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl', chl_bkg )
633	CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg )
634	chl_bkg(:,:,:) = chl_bkg(:,:,:) * tmask(:,:,:)
635	cchl_p_bkg(:,:,:) = cchl_p_bkg(:,:,:) * tmask(:,:,:)
636	#elif defined key_medusa
637	CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) )
638	CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) )
639	CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) )
640	CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
641	CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
642	#elif defined key_fabm
643	CALL iom_get( inum, jpdom_autoglo, 'ersem_chl1', tracer_bkg(:,:,:,jp_fabm_chl1) )
644	CALL iom_get( inum, jpdom_autoglo, 'ersem_chl2', tracer_bkg(:,:,:,jp_fabm_chl2) )
645	CALL iom_get( inum, jpdom_autoglo, 'ersem_chl3', tracer_bkg(:,:,:,jp_fabm_chl3) )
646	CALL iom_get( inum, jpdom_autoglo, 'ersem_chl4', tracer_bkg(:,:,:,jp_fabm_chl4) )
647	CALL iom_get( inum, jpdom_autoglo, 'ersem_p1c', tracer_bkg(:,:,:,jp_fabm_p1c) )
648	CALL iom_get( inum, jpdom_autoglo, 'ersem_p1n', tracer_bkg(:,:,:,jp_fabm_p1n) )
649	CALL iom_get( inum, jpdom_autoglo, 'ersem_p1p', tracer_bkg(:,:,:,jp_fabm_p1p) )
650	CALL iom_get( inum, jpdom_autoglo, 'ersem_p1s', tracer_bkg(:,:,:,jp_fabm_p1s) )
651	CALL iom_get( inum, jpdom_autoglo, 'ersem_p2c', tracer_bkg(:,:,:,jp_fabm_p2c) )
652	CALL iom_get( inum, jpdom_autoglo, 'ersem_p2n', tracer_bkg(:,:,:,jp_fabm_p2n) )
653	CALL iom_get( inum, jpdom_autoglo, 'ersem_p2p', tracer_bkg(:,:,:,jp_fabm_p2p) )
654	CALL iom_get( inum, jpdom_autoglo, 'ersem_p3c', tracer_bkg(:,:,:,jp_fabm_p3c) )
655	CALL iom_get( inum, jpdom_autoglo, 'ersem_p3n', tracer_bkg(:,:,:,jp_fabm_p3n) )
656	CALL iom_get( inum, jpdom_autoglo, 'ersem_p3p', tracer_bkg(:,:,:,jp_fabm_p3p) )
657	CALL iom_get( inum, jpdom_autoglo, 'ersem_p4c', tracer_bkg(:,:,:,jp_fabm_p4c) )
658	CALL iom_get( inum, jpdom_autoglo, 'ersem_p4n', tracer_bkg(:,:,:,jp_fabm_p4n) )
659	CALL iom_get( inum, jpdom_autoglo, 'ersem_p4p', tracer_bkg(:,:,:,jp_fabm_p4p) )
660	#endif
661
662	IF ( ln_phytobal ) THEN
663
664	CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg )
665	CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg )
666	CALL iom_get( inum, jpdom_autoglo, 'phyt_avg', phyt_avg_bkg )
667	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
668	pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1)
669	ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1)
670	phyt_avg_bkg(:,:) = phyt_avg_bkg(:,:) * tmask(:,:,1)
671	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
672
673	#if defined key_hadocc
674	CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) )
675	CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) )
676	CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) )
677	CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) )
678	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
679	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
680	#elif defined key_medusa
681	CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) )
682	CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) )
683	CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) )
684	CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) )
685	CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) )
686	CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) )
687	CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) )
688	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
689	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
690	CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
691	#elif defined key_fabm
692	CALL iom_get( inum, jpdom_autoglo, 'ersem_z4c', tracer_bkg(:,:,:,jp_fabm_z4c) )
693	CALL iom_get( inum, jpdom_autoglo, 'ersem_z5c', tracer_bkg(:,:,:,jp_fabm_z5c) )
694	CALL iom_get( inum, jpdom_autoglo, 'ersem_z5n', tracer_bkg(:,:,:,jp_fabm_z5n) )
695	CALL iom_get( inum, jpdom_autoglo, 'ersem_z5p', tracer_bkg(:,:,:,jp_fabm_z5p) )
696	CALL iom_get( inum, jpdom_autoglo, 'ersem_z6c', tracer_bkg(:,:,:,jp_fabm_z6c) )
697	CALL iom_get( inum, jpdom_autoglo, 'ersem_z6n', tracer_bkg(:,:,:,jp_fabm_z6n) )
698	CALL iom_get( inum, jpdom_autoglo, 'ersem_z6p', tracer_bkg(:,:,:,jp_fabm_z6p) )
699	CALL iom_get( inum, jpdom_autoglo, 'ersem_n1p', tracer_bkg(:,:,:,jp_fabm_n1p) )
700	CALL iom_get( inum, jpdom_autoglo, 'ersem_n3n', tracer_bkg(:,:,:,jp_fabm_n3n) )
701	CALL iom_get( inum, jpdom_autoglo, 'ersem_n4n', tracer_bkg(:,:,:,jp_fabm_n4n) )
702	CALL iom_get( inum, jpdom_autoglo, 'ersem_n5s', tracer_bkg(:,:,:,jp_fabm_n5s) )
703	CALL iom_get( inum, jpdom_autoglo, 'ersem_o2o', tracer_bkg(:,:,:,jp_fabm_o2o) )
704	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c', tracer_bkg(:,:,:,jp_fabm_o3c) )
705	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) )
706	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:) )
707	totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
708	#endif
709	ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
710	#if defined key_hadocc
711	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
712	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
713	#elif defined key_medusa
714	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
715	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
716	#elif defined key_fabm
717	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c', tracer_bkg(:,:,:,jp_fabm_o3c) )
718	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) )
719	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:) )
720	totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
721	#endif
722	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
723	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
724	ENDIF
725
726	CALL iom_close( inum )
727
728	DO jt = 1, jptra
729	tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
730	END DO
731
732	ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
733
734	ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
735	ALLOCATE( mld_max_bkg(jpi,jpj) )
736	tracer_bkg(:,:,:,:) = 0.0
737	mld_max_bkg(:,:) = 0.0
738
739	CALL iom_open( c_asmbkg, inum )
740
741	#if defined key_hadocc
742	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
743	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
744	#elif defined key_medusa
745	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
746	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
747	#elif defined key_fabm
748	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c', tracer_bkg(:,:,:,jp_fabm_o3c) )
749	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) )
750	CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:) )
751	totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
752	#endif
753	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
754
755	CALL iom_close( inum )
756
757	DO jt = 1, jptra
758	tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
759	END DO
760	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
761
762	ENDIF
763	#else
764	CALL ctl_stop( ' asm_bgc_init_bkg: key_top and a compatible biogeochemical model required' )
765	#endif
766
767	END SUBROUTINE asm_bgc_init_bkg
768
769	!!===========================================================================
770	!!===========================================================================
771	!!===========================================================================
772
773	SUBROUTINE asm_bgc_bal_wri( kt )
774	!!------------------------------------------------------------------------
775	!!
776	!! * ROUTINE asm_bgc_bal_wri *
777	!!
778	!! ** Purpose : Write to file the balancing increments
779	!! calculated for biogeochemistry
780	!!
781	!! ** Method : Write to file the balancing increments
782	!! calculated for biogeochemistry
783	!!
784	!! ** Action :
785	!!
786	!! References :
787	!!
788	!! History :
789	!! ! 2014-08 (D. Ford) Initial version, based on asm_bkg_wri
790	!!------------------------------------------------------------------------
791	!! * Arguments
792	INTEGER, INTENT( IN ) :: kt ! Current time-step
793	!! * Local declarations
794	CHARACTER(LEN=50) :: cl_asmbal ! Filename (with extension)
795	LOGICAL :: llok ! Check if file exists
796	INTEGER :: inum ! File unit number
797	REAL(wp) :: zdate ! Date
798	!!------------------------------------------------------------------------
799
800	! Set things up
801	zdate = REAL( ndastp )
802	WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal )
803
804	! Check if file exists
805	INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok )
806	IF( .NOT. llok ) THEN
807	IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// &
808	& TRIM( c_asmbal ) // ' at timestep = ', kt
809
810	! Define the output file
811	CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
812
813	! Write the information
814	CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate )
815
816	IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
817	& ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
818	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
819	& ln_slphynoninc ) THEN
820	#if defined key_medusa
821	CALL iom_rstput( kt, kt, inum, 'phy2d_chn', phyto2d_balinc(:,:,:,jpchn) )
822	CALL iom_rstput( kt, kt, inum, 'phy2d_chd', phyto2d_balinc(:,:,:,jpchd) )
823	CALL iom_rstput( kt, kt, inum, 'phy2d_phn', phyto2d_balinc(:,:,:,jpphn) )
824	CALL iom_rstput( kt, kt, inum, 'phy2d_phd', phyto2d_balinc(:,:,:,jpphd) )
825	CALL iom_rstput( kt, kt, inum, 'phy2d_pds', phyto2d_balinc(:,:,:,jppds) )
826	IF ( ln_phytobal ) THEN
827	CALL iom_rstput( kt, kt, inum, 'phy2d_zmi', phyto2d_balinc(:,:,:,jpzmi) )
828	CALL iom_rstput( kt, kt, inum, 'phy2d_zme', phyto2d_balinc(:,:,:,jpzme) )
829	CALL iom_rstput( kt, kt, inum, 'phy2d_din', phyto2d_balinc(:,:,:,jpdin) )
830	CALL iom_rstput( kt, kt, inum, 'phy2d_sil', phyto2d_balinc(:,:,:,jpsil) )
831	CALL iom_rstput( kt, kt, inum, 'phy2d_fer', phyto2d_balinc(:,:,:,jpfer) )
832	CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jpdet) )
833	CALL iom_rstput( kt, kt, inum, 'phy2d_dtc', phyto2d_balinc(:,:,:,jpdtc) )
834	CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jpdic) )
835	CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jpalk) )
836	CALL iom_rstput( kt, kt, inum, 'phy2d_oxy', phyto2d_balinc(:,:,:,jpoxy) )
837	ENDIF
838	#elif defined key_hadocc
839	CALL iom_rstput( kt, kt, inum, 'phy2d_phy', phyto2d_balinc(:,:,:,jp_had_phy) )
840	IF ( ln_phytobal ) THEN
841	CALL iom_rstput( kt, kt, inum, 'phy2d_nut', phyto2d_balinc(:,:,:,jp_had_nut) )
842	CALL iom_rstput( kt, kt, inum, 'phy2d_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) )
843	CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jp_had_pdn) )
844	CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jp_had_dic) )
845	CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) )
846	ENDIF
847	#elif defined key_fabm
848	CALL iom_rstput( kt, kt, inum, 'phy2d_chl1', phyto2d_balinc(:,:,:,jp_fabm_chl1) )
849	CALL iom_rstput( kt, kt, inum, 'phy2d_chl2', phyto2d_balinc(:,:,:,jp_fabm_chl2) )
850	CALL iom_rstput( kt, kt, inum, 'phy2d_chl3', phyto2d_balinc(:,:,:,jp_fabm_chl3) )
851	CALL iom_rstput( kt, kt, inum, 'phy2d_chl4', phyto2d_balinc(:,:,:,jp_fabm_chl4) )
852	CALL iom_rstput( kt, kt, inum, 'phy2d_p1c', phyto2d_balinc(:,:,:,jp_fabm_p1c) )
853	CALL iom_rstput( kt, kt, inum, 'phy2d_p1n', phyto2d_balinc(:,:,:,jp_fabm_p1n) )
854	CALL iom_rstput( kt, kt, inum, 'phy2d_p1p', phyto2d_balinc(:,:,:,jp_fabm_p1p) )
855	CALL iom_rstput( kt, kt, inum, 'phy2d_p1s', phyto2d_balinc(:,:,:,jp_fabm_p1s) )
856	CALL iom_rstput( kt, kt, inum, 'phy2d_p2c', phyto2d_balinc(:,:,:,jp_fabm_p2c) )
857	CALL iom_rstput( kt, kt, inum, 'phy2d_p2n', phyto2d_balinc(:,:,:,jp_fabm_p2n) )
858	CALL iom_rstput( kt, kt, inum, 'phy2d_p2p', phyto2d_balinc(:,:,:,jp_fabm_p2p) )
859	CALL iom_rstput( kt, kt, inum, 'phy2d_p3c', phyto2d_balinc(:,:,:,jp_fabm_p3c) )
860	CALL iom_rstput( kt, kt, inum, 'phy2d_p3n', phyto2d_balinc(:,:,:,jp_fabm_p3n) )
861	CALL iom_rstput( kt, kt, inum, 'phy2d_p3p', phyto2d_balinc(:,:,:,jp_fabm_p3p) )
862	CALL iom_rstput( kt, kt, inum, 'phy2d_p4c', phyto2d_balinc(:,:,:,jp_fabm_p4c) )
863	CALL iom_rstput( kt, kt, inum, 'phy2d_p4n', phyto2d_balinc(:,:,:,jp_fabm_p4n) )
864	CALL iom_rstput( kt, kt, inum, 'phy2d_p4p', phyto2d_balinc(:,:,:,jp_fabm_p4p) )
865	IF ( ln_phytobal ) THEN
866	CALL iom_rstput( kt, kt, inum, 'phy2d_z4c', phyto2d_balinc(:,:,:,jp_fabm_z4c) )
867	CALL iom_rstput( kt, kt, inum, 'phy2d_z5c', phyto2d_balinc(:,:,:,jp_fabm_z5c) )
868	CALL iom_rstput( kt, kt, inum, 'phy2d_z5n', phyto2d_balinc(:,:,:,jp_fabm_z5n) )
869	CALL iom_rstput( kt, kt, inum, 'phy2d_z5p', phyto2d_balinc(:,:,:,jp_fabm_z5p) )
870	CALL iom_rstput( kt, kt, inum, 'phy2d_z6c', phyto2d_balinc(:,:,:,jp_fabm_z6c) )
871	CALL iom_rstput( kt, kt, inum, 'phy2d_z6n', phyto2d_balinc(:,:,:,jp_fabm_z6n) )
872	CALL iom_rstput( kt, kt, inum, 'phy2d_z6p', phyto2d_balinc(:,:,:,jp_fabm_z6p) )
873	CALL iom_rstput( kt, kt, inum, 'phy2d_n1p', phyto2d_balinc(:,:,:,jp_fabm_n1p) )
874	CALL iom_rstput( kt, kt, inum, 'phy2d_n3n', phyto2d_balinc(:,:,:,jp_fabm_n3n) )
875	CALL iom_rstput( kt, kt, inum, 'phy2d_n4n', phyto2d_balinc(:,:,:,jp_fabm_n4n) )
876	CALL iom_rstput( kt, kt, inum, 'phy2d_n5s', phyto2d_balinc(:,:,:,jp_fabm_n5s) )
877	CALL iom_rstput( kt, kt, inum, 'phy2d_o2o', phyto2d_balinc(:,:,:,jp_fabm_o2o) )
878	CALL iom_rstput( kt, kt, inum, 'phy2d_o3c', phyto2d_balinc(:,:,:,jp_fabm_o3c) )
879	CALL iom_rstput( kt, kt, inum, 'phy2d_o3ba', phyto2d_balinc(:,:,:,jp_fabm_o3ba) )
880	ENDIF
881	#endif
882	ENDIF
883
884	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
885	#if defined key_medusa
886	CALL iom_rstput( kt, kt, inum, 'phy3d_chn', phyto3d_balinc(:,:,:,jpchn) )
887	CALL iom_rstput( kt, kt, inum, 'phy3d_chd', phyto3d_balinc(:,:,:,jpchd) )
888	CALL iom_rstput( kt, kt, inum, 'phy3d_phn', phyto3d_balinc(:,:,:,jpphn) )
889	CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) )
890	CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) )
891	#elif defined key_hadocc
892	CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) )
893	#elif defined key_fabm
894	CALL iom_rstput( kt, kt, inum, 'phy3d_chl1', phyto3d_balinc(:,:,:,jp_fabm_chl1) )
895	CALL iom_rstput( kt, kt, inum, 'phy3d_chl2', phyto3d_balinc(:,:,:,jp_fabm_chl2) )
896	CALL iom_rstput( kt, kt, inum, 'phy3d_chl3', phyto3d_balinc(:,:,:,jp_fabm_chl3) )
897	CALL iom_rstput( kt, kt, inum, 'phy3d_chl4', phyto3d_balinc(:,:,:,jp_fabm_chl4) )
898	CALL iom_rstput( kt, kt, inum, 'phy3d_p1c', phyto3d_balinc(:,:,:,jp_fabm_p1c) )
899	CALL iom_rstput( kt, kt, inum, 'phy3d_p1n', phyto3d_balinc(:,:,:,jp_fabm_p1n) )
900	CALL iom_rstput( kt, kt, inum, 'phy3d_p1p', phyto3d_balinc(:,:,:,jp_fabm_p1p) )
901	CALL iom_rstput( kt, kt, inum, 'phy3d_p1s', phyto3d_balinc(:,:,:,jp_fabm_p1s) )
902	CALL iom_rstput( kt, kt, inum, 'phy3d_p2c', phyto3d_balinc(:,:,:,jp_fabm_p2c) )
903	CALL iom_rstput( kt, kt, inum, 'phy3d_p2n', phyto3d_balinc(:,:,:,jp_fabm_p2n) )
904	CALL iom_rstput( kt, kt, inum, 'phy3d_p2p', phyto3d_balinc(:,:,:,jp_fabm_p2p) )
905	CALL iom_rstput( kt, kt, inum, 'phy3d_p3c', phyto3d_balinc(:,:,:,jp_fabm_p3c) )
906	CALL iom_rstput( kt, kt, inum, 'phy3d_p3n', phyto3d_balinc(:,:,:,jp_fabm_p3n) )
907	CALL iom_rstput( kt, kt, inum, 'phy3d_p3p', phyto3d_balinc(:,:,:,jp_fabm_p3p) )
908	CALL iom_rstput( kt, kt, inum, 'phy3d_p4c', phyto3d_balinc(:,:,:,jp_fabm_p4c) )
909	CALL iom_rstput( kt, kt, inum, 'phy3d_p4n', phyto3d_balinc(:,:,:,jp_fabm_p4n) )
910	CALL iom_rstput( kt, kt, inum, 'phy3d_p4p', phyto3d_balinc(:,:,:,jp_fabm_p4p) )
911	#endif
912	ENDIF
913
914	IF ( ln_spco2inc ) THEN
915	#if defined key_medusa
916	CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jpdic) )
917	CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jpalk) )
918	#elif defined key_hadocc
919	CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
920	CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
921	#elif defined key_fabm
922	CALL iom_rstput( kt, kt, inum, 'pco2_o3c', pco2_balinc(:,:,:,jp_fabm_o3c) )
923	CALL iom_rstput( kt, kt, inum, 'pco2_o3ba', pco2_balinc(:,:,:,jp_fabm_o3ba) )
924	#endif
925	ELSE IF ( ln_sfco2inc ) THEN
926	#if defined key_medusa
927	CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jpdic) )
928	CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jpalk) )
929	#elif defined key_hadocc
930	CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
931	CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
932	#elif defined key_fabm
933	CALL iom_rstput( kt, kt, inum, 'fco2_o3c', pco2_balinc(:,:,:,jp_fabm_o3c) )
934	CALL iom_rstput( kt, kt, inum, 'fco2_o3ba', pco2_balinc(:,:,:,jp_fabm_o3ba) )
935	#endif
936	ENDIF
937
938	IF ( ln_pphinc ) THEN
939	#if defined key_medusa
940	CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jpdic) )
941	CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jpalk) )
942	#elif defined key_hadocc
943	CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) )
944	CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) )
945	#elif defined key_fabm
946	CALL iom_rstput( kt, kt, inum, 'ph_o3c', ph_balinc(:,:,:,jp_fabm_o3c) )
947	CALL iom_rstput( kt, kt, inum, 'ph_o3ba', ph_balinc(:,:,:,jp_fabm_o3ba) )
948	#endif
949	ENDIF
950
951	CALL iom_close( inum )
952	ELSE
953	CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', &
954	& ' Therefore not writing out balancing increments at this timestep', &
955	& ' Something has probably gone wrong somewhere' )
956	IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// &
957	& TRIM( c_asmbal ) // ' at timestep = ', kt
958	ENDIF
959
960	END SUBROUTINE asm_bgc_bal_wri
961
962	!!===========================================================================
963	!!===========================================================================
964	!!===========================================================================
965
966	SUBROUTINE asm_bgc_bkg_alloc
967	!!------------------------------------------------------------------------
968	!! * ROUTINE asm_bgc_bkg_alloc *
969	!!
970	!! ** Purpose : allocate time-average arrays for background
971	!!
972	!! ** Method : allocate time-average arrays for background
973	!!
974	!! ** Action : allocate time-average arrays for background
975	!!
976	!! References : asm_bkg_wri
977	!!------------------------------------------------------------------------
978	!!
979	INTEGER :: ierror
980	!!
981	!!------------------------------------------------------------------------
982
983	ALLOCATE( pgrow_avg_tavg(jpi,jpj), STAT=ierror )
984	IF( ierror > 0 ) THEN
985	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate pgrow_avg_tavg' )
986	ENDIF
987	pgrow_avg_tavg(:,:) = 0.0
988
989	ALLOCATE( ploss_avg_tavg(jpi,jpj), STAT=ierror )
990	IF( ierror > 0 ) THEN
991	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate ploss_avg_tavg' )
992	ENDIF
993	ploss_avg_tavg(:,:) = 0.0
994
995	ALLOCATE( phyt_avg_tavg(jpi,jpj), STAT=ierror )
996	IF( ierror > 0 ) THEN
997	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate phyt_avg_tavg' )
998	ENDIF
999	phyt_avg_tavg(:,:) = 0.0
1000
1001	ALLOCATE( trn_tavg(jpi,jpj,jpk,jptra), STAT=ierror )
1002	IF( ierror > 0 ) THEN
1003	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate trn_tavg' )
1004	ENDIF
1005	trn_tavg(:,:,:,:) = 0.0
1006
1007	#if defined key_hadocc
1008	ALLOCATE( HADOCC_CHL_tavg(jpi,jpj,jpk), STAT=ierror )
1009	IF( ierror > 0 ) THEN
1010	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate HADOCC_CHL_tavg' )
1011	ENDIF
1012	HADOCC_CHL_tavg(:,:,:) = 0.0
1013
1014	ALLOCATE( cchl_p_tavg(jpi,jpj,jpk), STAT=ierror )
1015	IF( ierror > 0 ) THEN
1016	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate cchl_p_tavg' )
1017	ENDIF
1018	cchl_p_tavg(:,:,:) = 0.0
1019
1020	#elif defined key_fabm
1021	ALLOCATE( totalk_tavg(jpi,jpj,jpk), STAT=ierror )
1022	IF( ierror > 0 ) THEN
1023	CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate totalk_tavg' )
1024	ENDIF
1025	totalk_tavg(:,:,:) = 0.0
1026	#endif
1027
1028	END SUBROUTINE asm_bgc_bkg_alloc
1029
1030	!!===========================================================================
1031	!!===========================================================================
1032	!!===========================================================================
1033
1034	SUBROUTINE asm_bgc_bkg_tavg(kt, pnumtimes_tavg)
1035	!!------------------------------------------------------------------------
1036	!! * ROUTINE asm_bgc_bkg_wri *
1037	!!
1038	!! ** Purpose : write out bgc background
1039	!!
1040	!! ** Method : write out bgc background
1041	!!
1042	!! ** Action : write out bgc background
1043	!!
1044	!! References : asm_bkg_wri
1045	!!------------------------------------------------------------------------
1046	!!
1047	INTEGER, INTENT(in ) :: kt ! Current time-step
1048	!!
1049	REAL(wp), INTENT(in ) :: pnumtimes_tavg ! No of times to average over
1050	!!
1051	!!------------------------------------------------------------------------
1052
1053	IF (kt == nittrc000) THEN
1054	pgrow_avg_tavg(:,:) = 0.0
1055	ploss_avg_tavg(:,:) = 0.0
1056	phyt_avg_tavg(:,:) = 0.0
1057	trn_tavg(:,:,:,:) = 0.0
1058	#if defined key_hadocc
1059	HADOCC_CHL_tavg(:,:,:) = 0.0
1060	cchl_p_tavg(:,:,:) = 0.0
1061	#elif defined key_fabm
1062	totalk_tavg(:,:,:) = 0.0
1063	#endif
1064	ENDIF
1065
1066	pgrow_avg_tavg(:,:) = pgrow_avg_tavg(:,:) + pgrow_avg(:,:) / pnumtimes_tavg
1067	ploss_avg_tavg(:,:) = ploss_avg_tavg(:,:) + ploss_avg(:,:) / pnumtimes_tavg
1068	phyt_avg_tavg(:,:) = phyt_avg_tavg(:,:) + phyt_avg(:,:) / pnumtimes_tavg
1069	trn_tavg(:,:,:,:) = trn_tavg(:,:,:,:) + trn(:,:,:,:) / pnumtimes_tavg
1070	#if defined key_hadocc
1071	HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL_tavg(:,:,:) + HADOCC_CHL(:,:,:) / pnumtimes_tavg
1072	cchl_p_tavg(:,:,:) = cchl_p_tavg(:,:,:) + cchl_p(:,:,:) / pnumtimes_tavg
1073	#elif defined key_fabm
1074	totalk_tavg(:,:,:) = totalk_tavg(:,:,:) + &
1075	& fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) / pnumtimes_tavg
1076	totalk_tavg(:,:,:) = totalk_tavg(:,:,:) * tmask(:,:,:)
1077	#endif
1078
1079	END SUBROUTINE asm_bgc_bkg_tavg
1080
1081	!!===========================================================================
1082	!!===========================================================================
1083	!!===========================================================================
1084
1085	SUBROUTINE asm_bgc_bkg_wri( kt, knum, ld_avgbkg )
1086	!!------------------------------------------------------------------------
1087	!! * ROUTINE asm_bgc_bkg_wri *
1088	!!
1089	!! ** Purpose : write out bgc background
1090	!!
1091	!! ** Method : write out bgc background
1092	!!
1093	!! ** Action : write out bgc background
1094	!!
1095	!! References : asm_bkg_wri
1096	!!------------------------------------------------------------------------
1097	!!
1098	INTEGER, INTENT(in ) :: kt ! Current time-step
1099	INTEGER, INTENT(in ) :: knum ! i/o unit of increments file
1100	LOGICAL, INTENT(in ) :: ld_avgbkg ! Averaged background?
1101	!!
1102	INTEGER :: nitbgcbkg_r ! Period referenced to nit000
1103	!!
1104	!!------------------------------------------------------------------------
1105
1106	IF (ld_avgbkg) THEN
1107	nitbgcbkg_r = nitavgbkg_r
1108	ELSE
1109	nitbgcbkg_r = nitbkg_r
1110	pgrow_avg_tavg(:,:) = pgrow_avg(:,:)
1111	ploss_avg_tavg(:,:) = ploss_avg(:,:)
1112	phyt_avg_tavg(:,:) = phyt_avg(:,:)
1113	trn_tavg(:,:,:,:) = trn(:,:,:,:)
1114	#if defined key_hadocc
1115	HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL(:,:,:)
1116	cchl_p_tavg(:,:,:) = cchl_p(:,:,:)
1117	#elif defined key_fabm
1118	totalk_tavg(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta)
1119	totalk_tavg(:,:,:) = totalk_tavg(:,:,:) * tmask(:,:,:)
1120	#endif
1121	ENDIF
1122
1123	#if defined key_hadocc
1124	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg' , pgrow_avg_tavg )
1125	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg' , ploss_avg_tavg )
1126	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg' , phyt_avg_tavg )
1127	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max' , mld_max )
1128	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_nut' , trn_tavg(:,:,:,jp_had_nut) )
1129	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_phy' , trn_tavg(:,:,:,jp_had_phy) )
1130	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_zoo' , trn_tavg(:,:,:,jp_had_zoo) )
1131	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_pdn' , trn_tavg(:,:,:,jp_had_pdn) )
1132	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_dic' , trn_tavg(:,:,:,jp_had_dic) )
1133	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_alk' , trn_tavg(:,:,:,jp_had_alk) )
1134	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_chl' , HADOCC_CHL(:,:,:) )
1135	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:) )
1136	#elif defined key_medusa
1137	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg' , pgrow_avg_tavg )
1138	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg' , ploss_avg_tavg )
1139	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg' , phyt_avg_tavg )
1140	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max' , mld_max )
1141	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chn' , trn_tavg(:,:,:,jpchn) )
1142	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chd' , trn_tavg(:,:,:,jpchd) )
1143	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phn' , trn_tavg(:,:,:,jpphn) )
1144	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phd' , trn_tavg(:,:,:,jpphd) )
1145	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_pds' , trn_tavg(:,:,:,jppds) )
1146	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zmi' , trn_tavg(:,:,:,jpzmi) )
1147	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zme' , trn_tavg(:,:,:,jpzme) )
1148	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_din' , trn_tavg(:,:,:,jpdin) )
1149	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_sil' , trn_tavg(:,:,:,jpsil) )
1150	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_fer' , trn_tavg(:,:,:,jpfer) )
1151	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_det' , trn_tavg(:,:,:,jpdet) )
1152	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dtc' , trn_tavg(:,:,:,jpdtc) )
1153	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dic' , trn_tavg(:,:,:,jpdic) )
1154	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_alk' , trn_tavg(:,:,:,jpalk) )
1155	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_oxy' , trn_tavg(:,:,:,jpoxy) )
1156	#elif defined key_fabm
1157	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg' , pgrow_avg_tavg )
1158	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg' , ploss_avg_tavg )
1159	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg' , phyt_avg_tavg )
1160	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max' , mld_max )
1161	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl1' , trn_tavg(:,:,:,jp_fabm_chl1) )
1162	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl2' , trn_tavg(:,:,:,jp_fabm_chl2) )
1163	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl3' , trn_tavg(:,:,:,jp_fabm_chl3) )
1164	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl4' , trn_tavg(:,:,:,jp_fabm_chl4) )
1165	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1c' , trn_tavg(:,:,:,jp_fabm_p1c) )
1166	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1n' , trn_tavg(:,:,:,jp_fabm_p1n) )
1167	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1p' , trn_tavg(:,:,:,jp_fabm_p1p) )
1168	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1s' , trn_tavg(:,:,:,jp_fabm_p1s) )
1169	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2c' , trn_tavg(:,:,:,jp_fabm_p2c) )
1170	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2n' , trn_tavg(:,:,:,jp_fabm_p2n) )
1171	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2p' , trn_tavg(:,:,:,jp_fabm_p2p) )
1172	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3c' , trn_tavg(:,:,:,jp_fabm_p3c) )
1173	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3n' , trn_tavg(:,:,:,jp_fabm_p3n) )
1174	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3p' , trn_tavg(:,:,:,jp_fabm_p3p) )
1175	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4c' , trn_tavg(:,:,:,jp_fabm_p4c) )
1176	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4n' , trn_tavg(:,:,:,jp_fabm_p4n) )
1177	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4p' , trn_tavg(:,:,:,jp_fabm_p4p) )
1178	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z4c' , trn_tavg(:,:,:,jp_fabm_z4c) )
1179	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5c' , trn_tavg(:,:,:,jp_fabm_z5c) )
1180	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5n' , trn_tavg(:,:,:,jp_fabm_z5n) )
1181	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5p' , trn_tavg(:,:,:,jp_fabm_z5p) )
1182	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6c' , trn_tavg(:,:,:,jp_fabm_z6c) )
1183	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6n' , trn_tavg(:,:,:,jp_fabm_z6n) )
1184	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6p' , trn_tavg(:,:,:,jp_fabm_z6p) )
1185	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n1p' , trn_tavg(:,:,:,jp_fabm_n1p) )
1186	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n3n' , trn_tavg(:,:,:,jp_fabm_n3n) )
1187	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n4n' , trn_tavg(:,:,:,jp_fabm_n4n) )
1188	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n5s' , trn_tavg(:,:,:,jp_fabm_n5s) )
1189	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o2o' , trn_tavg(:,:,:,jp_fabm_o2o) )
1190	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3c' , trn_tavg(:,:,:,jp_fabm_o3c) )
1191	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ba' , trn_tavg(:,:,:,jp_fabm_o3ba) )
1192	CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ta' , totalk_tavg )
1193	#endif
1194
1195	END SUBROUTINE asm_bgc_bkg_wri
1196
1197	!!===========================================================================
1198	!!===========================================================================
1199	!!===========================================================================
1200
1201	SUBROUTINE asm_bgc_unlog_2d( pbkg, pinc_log, pinc_nonlog )
1202	!!------------------------------------------------------------------------
1203	!! * ROUTINE asm_bgc_init_incs *
1204	!!
1205	!! ** Purpose : convert log increments to non-log
1206	!!
1207	!! ** Method : need to account for model background,
1208	!! cannot simply do 10^log_inc. Need to:
1209	!! 1) Add log_inc to log10(background) to get log10(analysis)
1210	!! 2) Take 10^log10(analysis) to get analysis
1211	!! 3) Subtract background from analysis to get non-log incs
1212	!!
1213	!! ** Action : return non-log increments
1214	!!
1215	!! References :
1216	!!------------------------------------------------------------------------
1217	!!
1218	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pbkg ! Background
1219	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: pinc_log ! Log incs
1220	REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: pinc_nonlog ! Non-log incs
1221	!
1222	INTEGER :: ji, jj ! Loop counters
1223	!!
1224	!!------------------------------------------------------------------------
1225
1226	DO jj = 1, jpj
1227	DO ji = 1, jpi
1228	IF ( pbkg(ji,jj) > 0.0 ) THEN
1229	pinc_nonlog(ji,jj) = 10**( LOG10( pbkg(ji,jj) ) + &
1230	& pinc_log(ji,jj) ) &
1231	& - pbkg(ji,jj)
1232	ELSE
1233	pinc_nonlog(ji,jj) = 0.0
1234	ENDIF
1235	END DO
1236	END DO
1237
1238	END SUBROUTINE asm_bgc_unlog_2d
1239
1240	!!===========================================================================
1241	!!===========================================================================
1242	!!===========================================================================
1243
1244	SUBROUTINE phyto2d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1245	!!------------------------------------------------------------------------
1246	!! * ROUTINE phyto2d_asm_inc *
1247	!!
1248	!! ** Purpose : Apply the chlorophyll assimilation increments.
1249	!!
1250	!! ** Method : Calculate increments to state variables using nitrogen
1251	!! balancing.
1252	!! Direct initialization or Incremental Analysis Updating.
1253	!!
1254	!! ** Action :
1255	!!------------------------------------------------------------------------
1256	INTEGER, INTENT(IN) :: kt ! Current time step
1257	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1258	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1259	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1260	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1261	!
1262	INTEGER :: it ! Index
1263	REAL(wp) :: zincwgt ! IAU weight for current time step
1264	REAL(wp) :: zincper ! IAU interval in seconds
1265	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth
1266	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chltot ! Local chltot incs
1267	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chltot ! Local chltot bkg
1268	REAL(wp), DIMENSION(jpi,jpj) :: zinc_phytot ! Local phytot incs
1269	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phytot ! Local phytot bkg
1270	#if defined key_medusa
1271	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs
1272	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
1273	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnon ! Local chlnon incs
1274	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnon ! Local chlnon bkg
1275	REAL(wp), DIMENSION(jpi,jpj) :: zinc_phydia ! Local phydia incs
1276	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phydia ! Local phydia bkg
1277	REAL(wp), DIMENSION(jpi,jpj) :: zinc_phynon ! Local phynon incs
1278	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg
1279	#elif defined key_fabm
1280	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs
1281	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
1282	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldin ! Local chldin incs
1283	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldin ! Local chldin bkg
1284	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnan ! Local chlnan incs
1285	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnan ! Local chlnan bkg
1286	REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlpic ! Local chlpic incs
1287	REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlpic ! Local chlpic bkg
1288	#endif
1289	!!------------------------------------------------------------------------
1290
1291	IF ( kt <= nit000 ) THEN
1292
1293	! Un-log any log increments for passing to balancing routines
1294	! Remember that two sets of non-log increments should not be
1295	! expected to be in the same ratio as their log equivalents
1296
1297	! Total chlorophyll
1298	IF ( ln_slchltotinc ) THEN
1299	#if defined key_medusa
1300	zbkg_chltot(:,:) = tracer_bkg(:,:,1,jpchn) + tracer_bkg(:,:,1,jpchd)
1301	#elif defined key_hadocc
1302	zbkg_chltot(:,:) = chl_bkg(:,:,1)
1303	#elif defined key_fabm
1304	zbkg_chltot(:,:) = tracer_bkg(:,:,1,jp_fabm_chl1) + &
1305	& tracer_bkg(:,:,1,jp_fabm_chl2) + &
1306	& tracer_bkg(:,:,1,jp_fabm_chl3) + &
1307	& tracer_bkg(:,:,1,jp_fabm_chl4)
1308	#endif
1309	CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot )
1310	ELSE IF ( ln_schltotinc ) THEN
1311	zinc_chltot(:,:) = schltot_bkginc(:,:)
1312	ELSE
1313	zinc_chltot(:,:) = 0.0
1314	ENDIF
1315
1316	#if defined key_medusa \|\| defined key_fabm
1317	! Diatom chlorophyll
1318	IF ( ln_slchldiainc ) THEN
1319	#if defined key_medusa
1320	zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd)
1321	#elif defined key_fabm
1322	zbkg_chldia(:,:) = tracer_bkg(:,:,1,jp_fabm_chl1)
1323	#endif
1324	CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia )
1325	ELSE
1326	zinc_chldia(:,:) = 0.0
1327	ENDIF
1328	#endif
1329
1330	#if defined key_medusa
1331	! Non-diatom chlorophyll
1332	IF ( ln_slchlnoninc ) THEN
1333	zbkg_chlnon(:,:) = tracer_bkg(:,:,1,jpchn)
1334	CALL asm_bgc_unlog_2d( zbkg_chlnon, slchlnon_bkginc, zinc_chlnon )
1335	ELSE
1336	zinc_chlnon(:,:) = 0.0
1337	ENDIF
1338	#endif
1339
1340	#if defined key_fabm
1341	! Nanophytoplankton chlorophyll
1342	IF ( ln_slchlnaninc ) THEN
1343	zbkg_chlnan(:,:) = tracer_bkg(:,:,1,jp_fabm_chl2)
1344	CALL asm_bgc_unlog_2d( zbkg_chlnan, slchlnan_bkginc, zinc_chlnan )
1345	ELSE
1346	zinc_chlnan(:,:) = 0.0
1347	ENDIF
1348
1349	! Picophytoplankton chlorophyll
1350	IF ( ln_slchlpicinc ) THEN
1351	zbkg_chlpic(:,:) = tracer_bkg(:,:,1,jp_fabm_chl3)
1352	CALL asm_bgc_unlog_2d( zbkg_chlpic, slchlpic_bkginc, zinc_chlpic )
1353	ELSE
1354	zinc_chlpic(:,:) = 0.0
1355	ENDIF
1356
1357	! Dinoflagellate chlorophyll
1358	IF ( ln_slchldininc ) THEN
1359	zbkg_chldin(:,:) = tracer_bkg(:,:,1,jp_fabm_chl4)
1360	CALL asm_bgc_unlog_2d( zbkg_chldin, slchldin_bkginc, zinc_chldin )
1361	ELSE
1362	zinc_chldin(:,:) = 0.0
1363	ENDIF
1364	#endif
1365
1366	! Total phytoplankton carbon
1367	IF ( ln_slphytotinc ) THEN
1368	#if defined key_medusa
1369	zbkg_phytot(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
1370	#elif defined key_hadocc
1371	zbkg_phytot(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
1372	#endif
1373	CALL asm_bgc_unlog_2d( zbkg_phytot, slphytot_bkginc, zinc_phytot )
1374	ELSE
1375	zinc_phytot(:,:) = 0.0
1376	ENDIF
1377
1378	#if defined key_medusa
1379	! Diatom phytoplankton carbon
1380	IF ( ln_slphydiainc ) THEN
1381	zbkg_phydia(:,:) = trn(:,:,1,jpphd) * xthetapd
1382	CALL asm_bgc_unlog_2d( zbkg_phydia, slphydia_bkginc, zinc_phydia )
1383	ELSE
1384	zinc_phydia(:,:) = 0.0
1385	ENDIF
1386	#endif
1387
1388	#if defined key_medusa
1389	! Non-diatom phytoplankton carbon
1390	IF ( ln_slphynoninc ) THEN
1391	zbkg_phynon(:,:) = trn(:,:,1,jpphn) * xthetapn
1392	CALL asm_bgc_unlog_2d( zbkg_phynon, slphynon_bkginc, zinc_phynon )
1393	ELSE
1394	zinc_phynon(:,:) = 0.0
1395	ENDIF
1396	#endif
1397
1398	! Select mixed layer
1399	IF ( ll_asmdin ) THEN
1400	#if defined key_top && ( defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm )
1401	CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', &
1402	& ' Mixed layer depth taken to be background maximum mld_max_bkg' )
1403	zmld(:,:) = mld_max_bkg(:,:)
1404	#else
1405	CALL ctl_stop( ' Should not be doing biogeochemical assimilation without key_top' )
1406	#endif
1407	ELSE
1408	SELECT CASE( mld_choice_bgc )
1409	CASE ( 1 ) ! Turbocline/mixing depth [W points]
1410	zmld(:,:) = hmld(:,:)
1411	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
1412	zmld(:,:) = hmlp(:,:)
1413	CASE ( 3 ) ! Kara MLD [Interpolated]
1414	#if defined key_karaml
1415	IF ( ln_kara ) THEN
1416	zmld(:,:) = hmld_kara(:,:)
1417	ELSE
1418	CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
1419	& ' but ln_kara=.false.' )
1420	ENDIF
1421	#else
1422	CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
1423	& ' but is not defined' )
1424	#endif
1425	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
1426	zmld(:,:) = hmld_tref(:,:)
1427	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
1428	zmld(:,:) = hmlpt(:,:)
1429	END SELECT
1430	ENDIF
1431
1432	zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
1433
1434	#if defined key_medusa
1435	CALL asm_phyto2d_bal_medusa( (ln_slchltotinc .OR. ln_schltotinc), &
1436	& zinc_chltot, &
1437	& ln_slchldiainc, &
1438	& zinc_chldia, &
1439	& ln_slchlnoninc, &
1440	& zinc_chlnon, &
1441	& ln_slphytotinc, &
1442	& zinc_phytot, &
1443	& ln_slphydiainc, &
1444	& zinc_phydia, &
1445	& ln_slphynoninc, &
1446	& zinc_phynon, &
1447	& zincper, &
1448	& rn_maxchlinc, ln_phytobal, zmld, &
1449	& pgrow_avg_bkg, ploss_avg_bkg, &
1450	& phyt_avg_bkg, mld_max_bkg, &
1451	& tracer_bkg, phyto2d_balinc )
1452	#elif defined key_hadocc
1453	CALL asm_phyto2d_bal_hadocc( (ln_slchltotinc .OR. ln_schltotinc), &
1454	& zinc_chltot, &
1455	& ln_slphytotinc, &
1456	& zinc_phytot, &
1457	& zincper, &
1458	& rn_maxchlinc, ln_phytobal, zmld, &
1459	& pgrow_avg_bkg, ploss_avg_bkg, &
1460	& phyt_avg_bkg, mld_max_bkg, &
1461	& cchl_p_bkg(:,:,1), &
1462	& tracer_bkg, phyto2d_balinc )
1463	#elif defined key_fabm
1464	CALL asm_phyto2d_bal_ersem( (ln_slchltotinc .OR. ln_schltotinc), &
1465	& zinc_chltot, &
1466	& ln_slchldiainc, &
1467	& zinc_chldia, &
1468	& ln_slchlnaninc, &
1469	& zinc_chlnan, &
1470	& ln_slchlpicinc, &
1471	& zinc_chlpic, &
1472	& ln_slchldininc, &
1473	& zinc_chldin, &
1474	& zincper, &
1475	& rn_maxchlinc, ln_phytobal, zmld, &
1476	& pgrow_avg_bkg, ploss_avg_bkg, &
1477	& phyt_avg_bkg, mld_max_bkg, &
1478	& totalk_bkg, &
1479	& tracer_bkg, phyto2d_balinc )
1480	#else
1481	CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', &
1482	& 'but not defined a biogeochemical model' )
1483	#endif
1484
1485	ENDIF
1486
1487	IF ( ll_asmiau ) THEN
1488
1489	!--------------------------------------------------------------------
1490	! Incremental Analysis Updating
1491	!--------------------------------------------------------------------
1492
1493	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1494
1495	it = kt - nit000 + 1
1496	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1497	! note this is not a tendency so should not be divided by rdt
1498
1499	IF(lwp) THEN
1500	WRITE(numout,*)
1501	WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', &
1502	& kt,' with IAU weight = ', pwgtiau(it)
1503	WRITE(numout,*) '~~~~~~~~~~~~'
1504	ENDIF
1505
1506	! Update the biogeochemical variables
1507	! Add directly to trn and trb, rather than to tra, because tra gets
1508	! reset to zero at the start of trc_stp, called after this routine
1509	#if defined key_medusa
1510	WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1511	& trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1512	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1513	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1514	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1515	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1516	END WHERE
1517	#elif defined key_hadocc
1518	WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1519	& trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1520	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1521	& phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1522	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1523	& phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1524	END WHERE
1525	#elif defined key_fabm
1526	WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
1527	& trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
1528	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
1529	& phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
1530	trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
1531	& phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
1532	END WHERE
1533	#endif
1534
1535	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1536	! which is called at end of model run
1537	ENDIF
1538
1539	ELSEIF ( ll_asmdin ) THEN
1540
1541	!--------------------------------------------------------------------
1542	! Direct Initialization
1543	!--------------------------------------------------------------------
1544
1545	IF ( kt == nitdin_r ) THEN
1546
1547	neuler = 0 ! Force Euler forward step
1548
1549	! Initialize the now fields with the background + increment
1550	! Background currently is what the model is initialised with
1551	CALL ctl_warn( ' Doing direct initialisation with phyto2d assimilation', &
1552	& ' Background state is taken from model rather than background file' )
1553	#if defined key_medusa
1554	WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1555	& trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1556	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1557	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1)
1558	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1559	END WHERE
1560	#elif defined key_hadocc
1561	WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1562	& trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1563	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1564	& phyto2d_balinc(:,:,:,jp_had0:jp_had1)
1565	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1566	END WHERE
1567	#elif defined key_fabm
1568	WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
1569	& trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
1570	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
1571	& phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1)
1572	trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
1573	END WHERE
1574	#endif
1575
1576	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1577	! which is called at end of model run
1578	ENDIF
1579	!
1580	ENDIF
1581	!
1582	END SUBROUTINE phyto2d_asm_inc
1583
1584	!!===========================================================================
1585	!!===========================================================================
1586	!!===========================================================================
1587
1588	SUBROUTINE phyto3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1589	!!------------------------------------------------------------------------
1590	!! * ROUTINE phyto3d_asm_inc *
1591	!!
1592	!! ** Purpose : Apply the profile chlorophyll assimilation increments.
1593	!!
1594	!! ** Method : Calculate increments to state variables.
1595	!! Direct initialization or Incremental Analysis Updating.
1596	!!
1597	!! ** Action :
1598	!!------------------------------------------------------------------------
1599	INTEGER, INTENT(IN) :: kt ! Current time step
1600	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1601	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1602	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1603	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1604	!
1605	INTEGER :: ji, jj, jk ! Loop counters
1606	INTEGER :: it ! Index
1607	REAL(wp) :: zincwgt ! IAU weight for timestep
1608	REAL(wp) :: zfrac_chn ! Fraction of jpchn
1609	REAL(wp) :: zfrac_chd ! Fraction of jpchd
1610	REAL(wp) :: zfrac_chl1 ! Fraction of jp_fabm_chl1
1611	REAL(wp) :: zfrac_chl2 ! Fraction of jp_fabm_chl2
1612	REAL(wp) :: zfrac_chl3 ! Fraction of jp_fabm_chl3
1613	REAL(wp) :: zfrac_chl4 ! Fraction of jp_fabm_chl4
1614	REAL(wp) :: zrat_phn_chn ! jpphn:jpchn ratio
1615	REAL(wp) :: zrat_phd_chd ! jpphd:jpchd ratio
1616	REAL(wp) :: zrat_pds_chd ! jppds:jpchd ratio
1617	REAL(wp) :: zrat_p1c_chl1 ! jp_fabm_p1c:jp_fabm_chl1 ratio
1618	REAL(wp) :: zrat_p1n_chl1 ! jp_fabm_p1n:jp_fabm_chl1 ratio
1619	REAL(wp) :: zrat_p1p_chl1 ! jp_fabm_p1p:jp_fabm_chl1 ratio
1620	REAL(wp) :: zrat_p1s_chl1 ! jp_fabm_p1s:jp_fabm_chl1 ratio
1621	REAL(wp) :: zrat_p2c_chl2 ! jp_fabm_p2c:jp_fabm_chl2 ratio
1622	REAL(wp) :: zrat_p2n_chl2 ! jp_fabm_p2n:jp_fabm_chl2 ratio
1623	REAL(wp) :: zrat_p2p_chl2 ! jp_fabm_p2p:jp_fabm_chl2 ratio
1624	REAL(wp) :: zrat_p3c_chl3 ! jp_fabm_p3c:jp_fabm_chl3 ratio
1625	REAL(wp) :: zrat_p3n_chl3 ! jp_fabm_p3n:jp_fabm_chl3 ratio
1626	REAL(wp) :: zrat_p3p_chl3 ! jp_fabm_p3p:jp_fabm_chl3 ratio
1627	REAL(wp) :: zrat_p4c_chl4 ! jp_fabm_p4c:jp_fabm_chl4 ratio
1628	REAL(wp) :: zrat_p4n_chl4 ! jp_fabm_p4n:jp_fabm_chl4 ratio
1629	REAL(wp) :: zrat_p4p_chl4 ! jp_fabm_p4p:jp_fabm_chl4 ratio
1630	REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc ! Chlorophyll increments
1631	REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl ! Chlorophyll background
1632	!!------------------------------------------------------------------------
1633
1634	IF ( kt <= nit000 ) THEN
1635
1636	IF ( ln_plchltotinc ) THEN
1637	! Convert log10(chlorophyll) increment back to a chlorophyll increment
1638	! In order to transform logchl incs to chl incs, need to account for model
1639	! background, cannot simply do 10^logchl_bkginc. Need to:
1640	! 1) Add logchl inc to log10(background) to get log10(analysis)
1641	! 2) Take 10^log10(analysis) to get analysis
1642	! 3) Subtract background from analysis to get chl incs
1643	! If rn_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
1644	#if defined key_medusa
1645	bkg_chl(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd)
1646	#elif defined key_hadocc
1647	bkg_chl(:,:,:) = chl_bkg(:,:,:)
1648	#elif defined key_fabm
1649	bkg_chl(:,:,:) = tracer_bkg(:,:,:,jp_fabm_chl1) + &
1650	& tracer_bkg(:,:,:,jp_fabm_chl2) + &
1651	& tracer_bkg(:,:,:,jp_fabm_chl3) + &
1652	& tracer_bkg(:,:,:,jp_fabm_chl4)
1653	#endif
1654	DO jk = 1, jpk
1655	DO jj = 1, jpj
1656	DO ji = 1, jpi
1657	IF ( bkg_chl(ji,jj,jk) > 0.0 ) THEN
1658	chl_inc(ji,jj,jk) = 10**( LOG10( bkg_chl(ji,jj,jk) ) + plchltot_bkginc(ji,jj,jk) ) - bkg_chl(ji,jj,jk)
1659	IF ( rn_maxchlinc > 0.0 ) THEN
1660	chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( chl_inc(ji,jj,jk), rn_maxchlinc ) )
1661	ENDIF
1662	ELSE
1663	chl_inc(ji,jj,jk) = 0.0
1664	ENDIF
1665	END DO
1666	END DO
1667	END DO
1668	ELSE IF ( ln_pchltotinc ) THEN
1669	DO jk = 1, jpk
1670	DO jj = 1, jpj
1671	DO ji = 1, jpi
1672	IF ( rn_maxchlinc > 0.0 ) THEN
1673	chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( pchltot_bkginc(ji,jj,jk), rn_maxchlinc ) )
1674	ELSE
1675	chl_inc(ji,jj,jk) = pchltot_bkginc(ji,jj,jk)
1676	ENDIF
1677	END DO
1678	END DO
1679	END DO
1680	ENDIF
1681
1682	#if defined key_medusa
1683	! Loop over each grid point partioning the increments based on existing ratios
1684	DO jk = 1, jpk
1685	DO jj = 1, jpj
1686	DO ji = 1, jpi
1687	IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN
1688	zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd))
1689	zfrac_chd = 1.0 - zfrac_chn
1690	phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn
1691	phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd
1692	zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn)
1693	zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd)
1694	zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd)
1695	phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn
1696	phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd
1697	phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd
1698	ENDIF
1699	END DO
1700	END DO
1701	END DO
1702	#elif defined key_hadocc
1703	phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:)
1704	#elif defined key_fabm
1705	! Loop over each grid point partioning the increments based on existing ratios
1706	DO jk = 1, jpk
1707	DO jj = 1, jpj
1708	DO ji = 1, jpi
1709	IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_chl1) > 0.0 ) .AND. &
1710	& ( tracer_bkg(ji,jj,jk,jp_fabm_chl2) > 0.0 ) .AND. &
1711	& ( tracer_bkg(ji,jj,jk,jp_fabm_chl3) > 0.0 ) .AND. &
1712	& ( tracer_bkg(ji,jj,jk,jp_fabm_chl4) > 0.0 ) ) THEN
1713	zfrac_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_chl1) / bkg_chl(ji,jj,jk)
1714	zfrac_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_chl2) / bkg_chl(ji,jj,jk)
1715	zfrac_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_chl3) / bkg_chl(ji,jj,jk)
1716	zfrac_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_chl4) / bkg_chl(ji,jj,jk)
1717	phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) = chl_inc(ji,jj,jk) * zfrac_chl1
1718	phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) = chl_inc(ji,jj,jk) * zfrac_chl2
1719	phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) = chl_inc(ji,jj,jk) * zfrac_chl3
1720	phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) = chl_inc(ji,jj,jk) * zfrac_chl4
1721	zrat_p1c_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1c) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
1722	zrat_p1n_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1n) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
1723	zrat_p1p_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1p) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
1724	zrat_p1s_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1s) / tracer_bkg(ji,jj,jk,jp_fabm_chl1)
1725	zrat_p2c_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2c) / tracer_bkg(ji,jj,jk,jp_fabm_chl2)
1726	zrat_p2n_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2n) / tracer_bkg(ji,jj,jk,jp_fabm_chl2)
1727	zrat_p2p_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2p) / tracer_bkg(ji,jj,jk,jp_fabm_chl2)
1728	zrat_p3c_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3c) / tracer_bkg(ji,jj,jk,jp_fabm_chl3)
1729	zrat_p3n_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3n) / tracer_bkg(ji,jj,jk,jp_fabm_chl3)
1730	zrat_p3p_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3p) / tracer_bkg(ji,jj,jk,jp_fabm_chl3)
1731	zrat_p4c_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4c) / tracer_bkg(ji,jj,jk,jp_fabm_chl4)
1732	zrat_p4n_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4n) / tracer_bkg(ji,jj,jk,jp_fabm_chl4)
1733	zrat_p4p_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4p) / tracer_bkg(ji,jj,jk,jp_fabm_chl4)
1734	phyto3d_balinc(ji,jj,jk,jp_fabm_p1c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1c_chl1
1735	phyto3d_balinc(ji,jj,jk,jp_fabm_p1n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1n_chl1
1736	phyto3d_balinc(ji,jj,jk,jp_fabm_p1p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1p_chl1
1737	phyto3d_balinc(ji,jj,jk,jp_fabm_p1s) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1s_chl1
1738	phyto3d_balinc(ji,jj,jk,jp_fabm_p2c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2c_chl2
1739	phyto3d_balinc(ji,jj,jk,jp_fabm_p2n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2n_chl2
1740	phyto3d_balinc(ji,jj,jk,jp_fabm_p2p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2p_chl2
1741	phyto3d_balinc(ji,jj,jk,jp_fabm_p3c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3c_chl3
1742	phyto3d_balinc(ji,jj,jk,jp_fabm_p3n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3n_chl3
1743	phyto3d_balinc(ji,jj,jk,jp_fabm_p3p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3p_chl3
1744	phyto3d_balinc(ji,jj,jk,jp_fabm_p4c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4c_chl4
1745	phyto3d_balinc(ji,jj,jk,jp_fabm_p4n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4n_chl4
1746	phyto3d_balinc(ji,jj,jk,jp_fabm_p4p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4p_chl4
1747	ENDIF
1748	END DO
1749	END DO
1750	END DO
1751	#else
1752	CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', &
1753	& 'but not defined a biogeochemical model' )
1754	#endif
1755
1756	ENDIF
1757
1758	IF ( ll_asmiau ) THEN
1759
1760	!--------------------------------------------------------------------
1761	! Incremental Analysis Updating
1762	!--------------------------------------------------------------------
1763
1764	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1765
1766	it = kt - nit000 + 1
1767	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1768	! note this is not a tendency so should not be divided by rdt
1769
1770	IF(lwp) THEN
1771	WRITE(numout,*)
1772	WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', &
1773	& kt,' with IAU weight = ', pwgtiau(it)
1774	WRITE(numout,*) '~~~~~~~~~~~~'
1775	ENDIF
1776
1777	! Update the biogeochemical variables
1778	! Add directly to trn and trb, rather than to tra, because tra gets
1779	! reset to zero at the start of trc_stp, called after this routine
1780	#if defined key_medusa
1781	WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1782	& trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1783	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1784	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1785	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1786	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1787	END WHERE
1788	#elif defined key_hadocc
1789	WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1790	& trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1791	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1792	& phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1793	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1794	& phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1795	END WHERE
1796	#elif defined key_fabm
1797	WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
1798	& trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
1799	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
1800	& phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
1801	trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
1802	& phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
1803	END WHERE
1804	#endif
1805
1806	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1807	! which is called at end of model run
1808	ENDIF
1809
1810	ELSEIF ( ll_asmdin ) THEN
1811
1812	!--------------------------------------------------------------------
1813	! Direct Initialization
1814	!--------------------------------------------------------------------
1815
1816	IF ( kt == nitdin_r ) THEN
1817
1818	neuler = 0 ! Force Euler forward step
1819
1820	! Initialize the now fields with the background + increment
1821	! Background currently is what the model is initialised with
1822	CALL ctl_warn( ' Doing direct initialisation with phyto3d assimilation', &
1823	& ' Background state is taken from model rather than background file' )
1824	#if defined key_medusa
1825	WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1826	& trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1827	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1828	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1)
1829	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1830	END WHERE
1831	#elif defined key_hadocc
1832	WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1833	& trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1834	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1835	& phyto3d_balinc(:,:,:,jp_had0:jp_had1)
1836	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1837	END WHERE
1838	#elif defined key_fabm
1839	WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
1840	& trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
1841	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
1842	& phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1)
1843	trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
1844	END WHERE
1845	#endif
1846
1847	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1848	! which is called at end of model run
1849	ENDIF
1850	!
1851	ENDIF
1852	!
1853	END SUBROUTINE phyto3d_asm_inc
1854
1855	!!===========================================================================
1856	!!===========================================================================
1857	!!===========================================================================
1858
1859	SUBROUTINE pco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
1860	& ll_trainc, pt_bkginc, ps_bkginc )
1861	!!------------------------------------------------------------------------
1862	!! * ROUTINE pco2_asm_inc *
1863	!!
1864	!! ** Purpose : Apply the pco2/fco2 assimilation increments.
1865	!!
1866	!! ** Method : Calculate increments to state variables using carbon
1867	!! balancing.
1868	!! Direct initialization or Incremental Analysis Updating.
1869	!!
1870	!! ** Action :
1871	!!------------------------------------------------------------------------
1872	INTEGER, INTENT(IN) :: kt ! Current time step
1873	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1874	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1875	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1876	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1877	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
1878	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
1879	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
1880	!
1881	INTEGER :: ji, jj, jk ! Loop counters
1882	INTEGER :: it ! Index
1883	INTEGER :: jkmax ! Index of mixed layer depth
1884	REAL(wp) :: zincwgt ! IAU weight for current time step
1885	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth for balancing
1886	REAL(wp), DIMENSION(jpi,jpj) :: pco2_bkginc_temp ! For fCO2 conversion
1887	REAL(wp), DIMENSION(jpi,jpj) :: dic_bkg_temp ! DIC background
1888	REAL(wp), DIMENSION(jpi,jpj) :: alk_bkg_temp ! Alkalinity background
1889	REAL(wp), DIMENSION(jpi,jpj) :: tem_bkg_temp ! Temperature background
1890	REAL(wp), DIMENSION(jpi,jpj) :: sal_bkg_temp ! Salinity background
1891	REAL(wp), DIMENSION(1) :: patm ! Atmospheric pressure
1892	REAL(wp), DIMENSION(1) :: phyd ! Hydrostatic pressure
1893
1894	! Coefficients for fCO2 to pCO2 conversion
1895	REAL(wp) :: zcoef_fco2_1 = -1636.75
1896	REAL(wp) :: zcoef_fco2_2 = 12.0408
1897	REAL(wp) :: zcoef_fco2_3 = 0.0327957
1898	REAL(wp) :: zcoef_fco2_4 = 0.0000316528
1899	REAL(wp) :: zcoef_fco2_5 = 2.0
1900	REAL(wp) :: zcoef_fco2_6 = 57.7
1901	REAL(wp) :: zcoef_fco2_7 = 0.118
1902	REAL(wp) :: zcoef_fco2_8 = 82.0578
1903	!!------------------------------------------------------------------------
1904
1905	IF ( kt <= nit000 ) THEN
1906
1907	!----------------------------------------------------------------------
1908	! DIC and alkalinity balancing
1909	!----------------------------------------------------------------------
1910
1911	IF ( ln_sfco2inc ) THEN
1912	#if defined key_medusa && defined key_roam
1913	! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
1914	patm(1) = 1.0
1915	phyd(1) = 0.0
1916	DO jj = 1, jpj
1917	DO ji = 1, jpi
1918	CALL f2pCO2( sfco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
1919	END DO
1920	END DO
1921	#else
1922	! If assimilating fCO2, then convert to pCO2 using temperature
1923	! See flux_gas.F90 within HadOCC for details of calculation
1924	pco2_bkginc_temp(:,:) = sfco2_bkginc(:,:) / &
1925	& EXP((zcoef_fco2_1 + &
1926	& zcoef_fco2_2 * (tsn(:,:,1,1)+rt0) - &
1927	& zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
1928	& zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0) + &
1929	& zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0))) / &
1930	& (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
1931	#endif
1932	ELSE
1933	pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
1934	ENDIF
1935
1936	! Account for physics assimilation if required
1937	IF ( ll_trainc ) THEN
1938	tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
1939	sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
1940	ELSE
1941	tem_bkg_temp(:,:) = tsn(:,:,1,1)
1942	sal_bkg_temp(:,:) = tsn(:,:,1,2)
1943	ENDIF
1944
1945	#if defined key_medusa
1946	! Account for phytoplankton balancing if required
1947	IF ( ln_phytobal ) THEN
1948	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
1949	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk)
1950	ELSE
1951	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
1952	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
1953	ENDIF
1954
1955	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1956	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
1957	& pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
1958
1959	#elif defined key_hadocc
1960	! Account for phytoplankton balancing if required
1961	IF ( ln_phytobal ) THEN
1962	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
1963	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk)
1964	ELSE
1965	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
1966	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
1967	ENDIF
1968
1969	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1970	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
1971	& pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
1972
1973	#elif defined key_fabm
1974	! Account for phytoplankton balancing if required
1975	IF ( ln_phytobal ) THEN
1976	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_o3c) + phyto2d_balinc(:,:,1,jp_fabm_o3c)
1977	alk_bkg_temp(:,:) = totalk_bkg(:,:,1) + phyto2d_balinc(:,:,1,jp_fabm_o3ba)
1978	ELSE
1979	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_o3c)
1980	alk_bkg_temp(:,:) = totalk_bkg(:,:,1)
1981	ENDIF
1982
1983	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1984	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
1985	& pco2_balinc(:,:,1,jp_fabm_o3c), pco2_balinc(:,:,1,jp_fabm_o3ba) )
1986
1987	#else
1988	CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
1989	& 'but not defined a biogeochemical model' )
1990	#endif
1991
1992	! Select mixed layer
1993	IF ( ll_asmdin ) THEN
1994	#if defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm
1995	CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
1996	& ' Mixed layer depth taken to be background maximum mld_max_bkg' )
1997	zmld(:,:) = mld_max_bkg(:,:)
1998	#else
1999	CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
2000	#endif
2001	ELSE
2002	SELECT CASE( mld_choice_bgc )
2003	CASE ( 1 ) ! Turbocline/mixing depth [W points]
2004	zmld(:,:) = hmld(:,:)
2005	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
2006	zmld(:,:) = hmlp(:,:)
2007	CASE ( 3 ) ! Kara MLD [Interpolated]
2008	#if defined key_karaml
2009	IF ( ln_kara ) THEN
2010	zmld(:,:) = hmld_kara(:,:)
2011	ELSE
2012	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
2013	& ' but ln_kara=.false.' )
2014	ENDIF
2015	#else
2016	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
2017	& ' but is not defined' )
2018	#endif
2019	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
2020	zmld(:,:) = hmld_tref(:,:)
2021	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
2022	zmld(:,:) = hmlpt(:,:)
2023	END SELECT
2024	ENDIF
2025
2026	! Propagate through mixed layer
2027	DO jj = 1, jpj
2028	DO ji = 1, jpi
2029	!
2030	jkmax = jpk-1
2031	DO jk = jpk-1, 1, -1
2032	IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. &
2033	& ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
2034	zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
2035	jkmax = jk
2036	ENDIF
2037	END DO
2038	!
2039	#if defined key_top
2040	DO jk = 2, jkmax
2041	pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
2042	END DO
2043	#endif
2044	!
2045	END DO
2046	END DO
2047
2048	ENDIF
2049
2050	IF ( ll_asmiau ) THEN
2051
2052	!--------------------------------------------------------------------
2053	! Incremental Analysis Updating
2054	!--------------------------------------------------------------------
2055
2056	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
2057
2058	it = kt - nit000 + 1
2059	zincwgt = pwgtiau(it) ! IAU weight for the current time step
2060	! note this is not a tendency so should not be divided by rdt
2061
2062	IF(lwp) THEN
2063	WRITE(numout,*)
2064	IF ( ln_spco2inc ) THEN
2065	WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
2066	& kt,' with IAU weight = ', pwgtiau(it)
2067	ELSE IF ( ln_sfco2inc ) THEN
2068	WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
2069	& kt,' with IAU weight = ', pwgtiau(it)
2070	ENDIF
2071	WRITE(numout,*) '~~~~~~~~~~~~'
2072	ENDIF
2073
2074	! Update the biogeochemical variables
2075	! Add directly to trn and trb, rather than to tra, because tra gets
2076	! reset to zero at the start of trc_stp, called after this routine
2077	#if defined key_medusa
2078	WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
2079	& trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
2080	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
2081	& pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
2082	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
2083	& pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
2084	END WHERE
2085	#elif defined key_hadocc
2086	WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
2087	& trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
2088	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
2089	& pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
2090	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
2091	& pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
2092	END WHERE
2093	#elif defined key_fabm
2094	WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
2095	& trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
2096	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
2097	& pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
2098	trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
2099	& pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
2100	END WHERE
2101	#endif
2102
2103	! Do not deallocate arrays - needed by asm_bgc_bal_wri
2104	! which is called at end of model run
2105
2106	ENDIF
2107
2108	ELSEIF ( ll_asmdin ) THEN
2109
2110	!--------------------------------------------------------------------
2111	! Direct Initialization
2112	!--------------------------------------------------------------------
2113
2114	IF ( kt == nitdin_r ) THEN
2115
2116	neuler = 0 ! Force Euler forward step
2117
2118	! Initialize the now fields with the background + increment
2119	! Background currently is what the model is initialised with
2120	CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
2121	& ' Background state is taken from model rather than background file' )
2122	#if defined key_medusa
2123	WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
2124	& trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
2125	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
2126	& pco2_balinc(:,:,:,jp_msa0:jp_msa1)
2127	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
2128	END WHERE
2129	#elif defined key_hadocc
2130	WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
2131	& trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
2132	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
2133	& pco2_balinc(:,:,:,jp_had0:jp_had1)
2134	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
2135	END WHERE
2136	#elif defined key_fabm
2137	WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
2138	& trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
2139	trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
2140	& pco2_balinc(:,:,:,jp_fabm0:jp_fabm1)
2141	trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
2142	END WHERE
2143	#endif
2144
2145	! Do not deallocate arrays - needed by asm_bgc_bal_wri
2146	! which is called at end of model run
2147	ENDIF
2148	!
2149	ENDIF
2150	!
2151	END SUBROUTINE pco2_asm_inc
2152
2153	!!===========================================================================
2154	!!===========================================================================
2155	!!===========================================================================
2156
2157	SUBROUTINE ph_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
2158	& ll_trainc, pt_bkginc, ps_bkginc )
2159	!!------------------------------------------------------------------------
2160	!! * ROUTINE ph_asm_inc *
2161	!!
2162	!! ** Purpose : Apply the pH assimilation increments.
2163	!!
2164	!! ** Method : Calculate increments to DIC and alkalinity from pH
2165	!! Use a similar approach to the pCO2 scheme
2166	!! Direct initialization or Incremental Analysis Updating.
2167	!!
2168	!! ** Action :
2169	!!------------------------------------------------------------------------
2170	INTEGER, INTENT(IN) :: kt ! Current time step
2171	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
2172	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
2173	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
2174	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
2175	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
2176	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
2177	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
2178
2179	REAL(wp) :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation
2180	REAL(wp) :: DIC_IN, ALK_IN ! DIC/alk in pH calculation
2181	REAL(wp) :: PH_OUT ! pH from pH calculation
2182	REAL(wp) :: ph_bkg ! pH from pH calculation
2183	REAL(wp) :: ph_dic ! pH from pH calculation
2184	REAL(wp) :: ph_alk ! pH from pH calculation
2185	REAL(wp) :: dph_ddic ! Change in pH wrt DIC
2186	REAL(wp) :: dph_dalk ! Change in pH wrt alk
2187	REAL(wp) :: weight ! Increment weighting
2188	REAL(wp) :: zincwgt ! IAU weight for current time step
2189	REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp ! DIC background
2190	REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp ! Alkalinity background
2191	REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp ! DIN background
2192	REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp ! Silicate background
2193	REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp ! Temperature background
2194	REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp ! Salinity background
2195	REAL(wp), DIMENSION(19) :: dummy_out ! Dummy variables
2196	INTEGER :: ji, jj, jk, jx ! Loop counters
2197	INTEGER :: it ! Index
2198	!!------------------------------------------------------------------------
2199
2200	#if ! defined key_medusa
2201	CALL ctl_stop( ' pH balancing only implemented for MEDUSA' )
2202	#else
2203
2204	IF ( kt <= nit000 ) THEN
2205
2206	!----------------------------------------------------------------------
2207	! DIC and alkalinity balancing
2208	!----------------------------------------------------------------------
2209
2210	! Account for physics assimilation if required
2211	IF ( ll_trainc ) THEN
2212	tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:)
2213	sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:)
2214	ELSE
2215	tem_bkg_temp(:,:,:) = tsn(:,:,:,1)
2216	sal_bkg_temp(:,:,:) = tsn(:,:,:,2)
2217	ENDIF
2218
2219	! Account for phytoplankton balancing if required
2220	IF ( ln_phytobal ) THEN
2221	dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic)
2222	alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk)
2223	din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin)
2224	sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil)
2225	ELSE
2226	dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic)
2227	alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk)
2228	din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin)
2229	sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil)
2230	ENDIF
2231
2232	! Account for pCO2 balancing if required
2233	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2234	dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic)
2235	alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk)
2236	ENDIF
2237
2238	! Loop over grid points and calculate dpH/dDIC and dpH/dAlk
2239	! This requires three calls to the MOCSY carbonate package
2240	! One to find pH at background DIC and alkalinity
2241	! One to find pH when DIC is increased by 10
2242	! One to find pH when alkalinity is increased by 10
2243	! Then calculate the gradients
2244
2245	DO jk = 1, jpk
2246	DO jj = 1, jpj
2247	DO ji = 1, jpi
2248
2249	IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN
2250
2251	DO jx = 1, 3
2252	IF ( jx == 1 ) THEN
2253	DIC_IN = dic_bkg_temp(ji,jj,jk)
2254	ALK_IN = alk_bkg_temp(ji,jj,jk)
2255	ELSE IF ( jx == 2 ) THEN
2256	DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch
2257	ALK_IN = alk_bkg_temp(ji,jj,jk)
2258	ELSE IF ( jx == 3 ) THEN
2259	DIC_IN = dic_bkg_temp(ji,jj,jk)
2260	ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch
2261	ENDIF
2262
2263	CALL mocsy_interface(tsn(ji,jj,jk,jp_tem), & ! IN: temperature
2264	& tsn(ji,jj,jk,jp_sal), & ! salinity
2265	& ALK_IN, & ! alkalinity
2266	& DIC_IN, & ! DIC
2267	& sil_bkg_temp(ji,jj,jk), & ! silicate
2268	& din_bkg_temp(ji,jj,jk)/16.0, & ! phosphate (Redfield from DIN)
2269	& 1.0, & ! atmospheric pressure (dummy)
2270	& fsdept(ji,jj,jk), & ! depth
2271	& gphit(ji,jj), & ! latitude
2272	& 1.0, & ! gas transfer (dummy)
2273	& 1.0, & ! atmospheric xCO2 (dummy)
2274	& 1, & ! number of points
2275	& PH_OUT, & ! OUT: pH
2276	& dummy_out(1), dummy_out(2), & ! stuff we don't care about
2277	& dummy_out(3), dummy_out(4), &
2278	& dummy_out(5), dummy_out(6), &
2279	& dummy_out(7), dummy_out(8), &
2280	& dummy_out(9), dummy_out(10), &
2281	& dummy_out(11), dummy_out(12), &
2282	& dummy_out(13), dummy_out(14), &
2283	& dummy_out(15), dummy_out(16), &
2284	& dummy_out(17), dummy_out(18), &
2285	& dummy_out(19) )
2286
2287	IF ( jx == 1 ) THEN
2288	ph_bkg = PH_OUT
2289	ELSE IF ( jx == 2 ) THEN
2290	ph_dic = PH_OUT
2291	ELSE IF ( jx == 3 ) THEN
2292	ph_alk = PH_OUT
2293	ENDIF
2294	END DO
2295
2296	dph_ddic = (ph_dic - ph_bkg) / zsearch
2297	dph_dalk = (ph_alk - ph_bkg) / zsearch
2298	weight = pph_bkginc(ji,jj,jk) / ( (dph_ddicdph_ddic) + (dph_dalkdph_dalk) )
2299
2300	ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic
2301	ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk
2302
2303	ENDIF
2304
2305	END DO
2306	END DO
2307	END DO
2308
2309	ENDIF
2310
2311	IF ( ll_asmiau ) THEN
2312
2313	!--------------------------------------------------------------------
2314	! Incremental Analysis Updating
2315	!--------------------------------------------------------------------
2316
2317	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
2318
2319	it = kt - nit000 + 1
2320	zincwgt = pwgtiau(it) ! IAU weight for the current time step
2321	! note this is not a tendency so should not be divided by rdt
2322
2323	IF(lwp) THEN
2324	WRITE(numout,*)
2325	WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', &
2326	& kt,' with IAU weight = ', pwgtiau(it)
2327	WRITE(numout,*) '~~~~~~~~~~~~'
2328	ENDIF
2329
2330	! Update the biogeochemical variables
2331	! Add directly to trn and trb, rather than to tra, because tra gets
2332	! reset to zero at the start of trc_stp, called after this routine
2333	WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
2334	& trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
2335	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
2336	& ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
2337	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
2338	& ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
2339	END WHERE
2340
2341	! Do not deallocate arrays - needed by asm_bgc_bal_wri
2342	! which is called at end of model run
2343
2344	ENDIF
2345
2346	ELSEIF ( ll_asmdin ) THEN
2347
2348	!--------------------------------------------------------------------
2349	! Direct Initialization
2350	!--------------------------------------------------------------------
2351
2352	IF ( kt == nitdin_r ) THEN
2353
2354	neuler = 0 ! Force Euler forward step
2355
2356	! Initialize the now fields with the background + increment
2357	! Background currently is what the model is initialised with
2358	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pH assimilation', &
2359	& ' Background state is taken from model rather than background file' )
2360	WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
2361	& trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
2362	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
2363	& ph_balinc(:,:,:,jp_msa0:jp_msa1)
2364	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
2365	END WHERE
2366
2367	! Do not deallocate arrays - needed by asm_bgc_bal_wri
2368	! which is called at end of model run
2369	ENDIF
2370	!
2371	ENDIF
2372	#endif
2373	!
2374	END SUBROUTINE ph_asm_inc
2375
2376	!!===========================================================================
2377	!!===========================================================================
2378	!!===========================================================================
2379
2380	SUBROUTINE bgc3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
2381	!!----------------------------------------------------------------------
2382	!! * ROUTINE dyn_asm_inc *
2383	!!
2384	!! ** Purpose : Apply generic 3D biogeochemistry assimilation increments.
2385	!!
2386	!! ** Method : Direct initialization or Incremental Analysis Updating.
2387	!!
2388	!! ** Action :
2389	!!----------------------------------------------------------------------
2390	INTEGER, INTENT(IN) :: kt ! Current time step
2391	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
2392	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
2393	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
2394	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
2395	!
2396	INTEGER :: it ! Index
2397	REAL(wp) :: zincwgt ! IAU weight for current time step
2398	!!----------------------------------------------------------------------
2399
2400	IF ( kt <= nit000 ) THEN
2401
2402	!----------------------------------------------------------------------
2403	! Remove any other balancing increments
2404	!----------------------------------------------------------------------
2405
2406	IF ( ln_pno3inc ) THEN
2407	#if defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm
2408	#if defined key_hadocc
2409	it = jp_had_nut
2410	#elif defined key_medusa
2411	it = jpdin
2412	#elif defined key_fabm
2413	it = jp_fabm_n3n
2414	#endif
2415	IF ( ln_phytobal ) THEN
2416	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2417	ENDIF
2418	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2419	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2420	ENDIF
2421	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2422	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2423	ENDIF
2424	IF ( ln_pphinc ) THEN
2425	pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2426	ENDIF
2427	#else
2428	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2429	#endif
2430	ENDIF
2431
2432	IF ( ln_psi4inc ) THEN
2433	#if defined key_medusa \|\| defined key_fabm
2434	#if defined key_medusa
2435	it = jpsil
2436	#elif defined key_fabm
2437	it = jp_fabm_n5s
2438	#endif
2439	IF ( ln_phytobal ) THEN
2440	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2441	ENDIF
2442	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2443	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2444	ENDIF
2445	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2446	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2447	ENDIF
2448	IF ( ln_pphinc ) THEN
2449	psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2450	ENDIF
2451	#else
2452	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2453	#endif
2454	ENDIF
2455
2456	IF ( ln_ppo4inc ) THEN
2457	#if defined key_fabm
2458	it = jp_fabm_n1p
2459
2460	IF ( ln_phytobal ) THEN
2461	ppo4_bkginc(:,:,:) = ppo4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2462	ENDIF
2463	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2464	ppo4_bkginc(:,:,:) = ppo4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2465	ENDIF
2466	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2467	ppo4_bkginc(:,:,:) = ppo4_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2468	ENDIF
2469	IF ( ln_pphinc ) THEN
2470	ppo4_bkginc(:,:,:) = ppo4_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2471	ENDIF
2472	#else
2473	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2474	#endif
2475	ENDIF
2476
2477	IF ( ln_pdicinc ) THEN
2478	#if defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm
2479	#if defined key_hadocc
2480	it = jp_had_dic
2481	#elif defined key_medusa
2482	it = jpdic
2483	#elif defined key_fabm
2484	it = jp_fabm_o3c
2485	#endif
2486	IF ( ln_phytobal ) THEN
2487	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2488	ENDIF
2489	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2490	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2491	ENDIF
2492	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2493	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2494	ENDIF
2495	IF ( ln_pphinc ) THEN
2496	pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2497	ENDIF
2498	#else
2499	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2500	#endif
2501	ENDIF
2502
2503	IF ( ln_palkinc ) THEN
2504	#if defined key_hadocc \|\| defined key_medusa \|\| defined key_fabm
2505	#if defined key_hadocc
2506	it = jp_had_alk
2507	#elif defined key_medusa
2508	it = jpalk
2509	#elif defined key_fabm
2510	it = jp_fabm_o3ba
2511	#endif
2512	IF ( ln_phytobal ) THEN
2513	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2514	ENDIF
2515	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2516	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2517	ENDIF
2518	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2519	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2520	ENDIF
2521	IF ( ln_pphinc ) THEN
2522	palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2523	ENDIF
2524	#else
2525	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2526	#endif
2527	ENDIF
2528
2529	IF ( ln_po2inc ) THEN
2530	#if defined key_medusa \|\| defined key_fabm
2531	#if defined key_medusa
2532	it = jpoxy
2533	#elif defined key_fabm
2534	it = jp_fabm_o2o
2535	#endif
2536	IF ( ln_phytobal ) THEN
2537	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2538	ENDIF
2539	IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2540	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2541	ENDIF
2542	IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2543	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2544	ENDIF
2545	IF ( ln_pphinc ) THEN
2546	po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2547	ENDIF
2548	#else
2549	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2550	#endif
2551	ENDIF
2552
2553	ENDIF
2554
2555	IF ( ll_asmiau ) THEN
2556
2557	!--------------------------------------------------------------------
2558	! Incremental Analysis Updating
2559	!--------------------------------------------------------------------
2560
2561	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
2562
2563	it = kt - nit000 + 1
2564	zincwgt = pwgtiau(it) ! IAU weight for the current time step
2565	! note this is not a tendency so should not be divided by rdt
2566
2567	IF(lwp) THEN
2568	WRITE(numout,*)
2569	WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', &
2570	& kt,' with IAU weight = ', pwgtiau(it)
2571	WRITE(numout,*) '~~~~~~~~~~~~'
2572	ENDIF
2573
2574	! Update the 3D BGC variables
2575	! Add directly to trn and trb, rather than to tra, because tra gets
2576	! reset to zero at the start of trc_stp, called after this routine
2577	! Don't apply increments if they'll take concentrations negative
2578
2579	IF ( ln_pno3inc ) THEN
2580	#if defined key_hadocc
2581	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2582	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2583	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
2584	trb(:,:,:,jp_had_nut) = trb(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
2585	END WHERE
2586	#elif defined key_medusa
2587	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2588	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2589	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
2590	trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
2591	END WHERE
2592	#elif defined key_fabm
2593	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2594	& trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2595	trn(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt
2596	trb(:,:,:,jp_fabm_n3n) = trb(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt
2597	END WHERE
2598	#else
2599	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2600	#endif
2601	ENDIF
2602
2603	IF ( ln_psi4inc ) THEN
2604	#if defined key_medusa
2605	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2606	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
2607	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
2608	trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
2609	END WHERE
2610	#elif defined key_fabm
2611	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2612	& trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
2613	trn(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt
2614	trb(:,:,:,jp_fabm_n5s) = trb(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt
2615	END WHERE
2616	#else
2617	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2618	#endif
2619	ENDIF
2620
2621	IF ( ln_ppo4inc ) THEN
2622	#if defined key_fabm
2623	WHERE( ppo4_bkginc(:,:,:) > 0.0_wp .OR. &
2624	& trn(:,:,:,jp_fabm_n1p) + ppo4_bkginc(:,:,:) * zincwgt > 0.0_wp )
2625	trn(:,:,:,jp_fabm_n1p) = trn(:,:,:,jp_fabm_n1p) + ppo4_bkginc(:,:,:) * zincwgt
2626	trb(:,:,:,jp_fabm_n1p) = trb(:,:,:,jp_fabm_n1p) + ppo4_bkginc(:,:,:) * zincwgt
2627	END WHERE
2628	#else
2629	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2630	#endif
2631	ENDIF
2632
2633	IF ( ln_pdicinc ) THEN
2634	#if defined key_hadocc
2635	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2636	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2637	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
2638	trb(:,:,:,jp_had_dic) = trb(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
2639	END WHERE
2640	#elif defined key_medusa
2641	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2642	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2643	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
2644	trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
2645	END WHERE
2646	#elif defined key_fabm
2647	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2648	& trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2649	trn(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt
2650	trb(:,:,:,jp_fabm_o3c) = trb(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt
2651	END WHERE
2652	#else
2653	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2654	#endif
2655	ENDIF
2656
2657	IF ( ln_palkinc ) THEN
2658	#if defined key_hadocc
2659	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2660	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2661	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
2662	trb(:,:,:,jp_had_alk) = trb(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
2663	END WHERE
2664	#elif defined key_medusa
2665	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2666	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2667	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
2668	trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
2669	END WHERE
2670	#elif defined key_fabm
2671	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2672	& trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2673	trn(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt
2674	trb(:,:,:,jp_fabm_o3ba) = trb(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt
2675	END WHERE
2676	#else
2677	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2678	#endif
2679	ENDIF
2680
2681	IF ( ln_po2inc ) THEN
2682	#if defined key_medusa
2683	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2684	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
2685	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
2686	trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
2687	END WHERE
2688	#elif defined key_fabm
2689	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2690	& trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
2691	trn(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt
2692	trb(:,:,:,jp_fabm_o2o) = trb(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt
2693	END WHERE
2694	#else
2695	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2696	#endif
2697	ENDIF
2698
2699	IF ( kt == nitiaufin_r ) THEN
2700	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
2701	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
2702	IF ( ln_ppo4inc ) DEALLOCATE( ppo4_bkginc )
2703	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
2704	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
2705	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
2706	ENDIF
2707
2708	ENDIF
2709
2710	ELSEIF ( ll_asmdin ) THEN
2711
2712	!--------------------------------------------------------------------
2713	! Direct Initialization
2714	!--------------------------------------------------------------------
2715
2716	IF ( kt == nitdin_r ) THEN
2717
2718	neuler = 0 ! Force Euler forward step
2719
2720	! Initialize the now fields with the background + increment
2721	! Background currently is what the model is initialised with
2722	CALL ctl_warn( ' Doing direct initialisation with 3D BGC assimilation', &
2723	& ' Background state is taken from model rather than background file' )
2724
2725	IF ( ln_pno3inc ) THEN
2726	#if defined key_hadocc
2727	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2728	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) > 0.0_wp )
2729	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:)
2730	trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut)
2731	END WHERE
2732	#elif defined key_medusa
2733	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2734	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) > 0.0_wp )
2735	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:)
2736	trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
2737	END WHERE
2738	#elif defined key_fabm
2739	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2740	& trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) > 0.0_wp )
2741	trn(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:)
2742	trb(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n)
2743	END WHERE
2744	#else
2745	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2746	#endif
2747	ENDIF
2748
2749	IF ( ln_psi4inc ) THEN
2750	#if defined key_medusa
2751	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2752	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) > 0.0_wp )
2753	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:)
2754	trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
2755	END WHERE
2756	#elif defined key_fabm
2757	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2758	& trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) > 0.0_wp )
2759	trn(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:)
2760	trb(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s)
2761	END WHERE
2762	#else
2763	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2764	#endif
2765	ENDIF
2766
2767	IF ( ln_ppo4inc ) THEN
2768	#if defined key_fabm
2769	WHERE( ppo4_bkginc(:,:,:) > 0.0_wp .OR. &
2770	& trn(:,:,:,jp_fabm_n1p) + ppo4_bkginc(:,:,:) > 0.0_wp )
2771	trn(:,:,:,jp_fabm_n1p) = trn(:,:,:,jp_fabm_n1p) + ppo4_bkginc(:,:,:)
2772	trb(:,:,:,jp_fabm_n1p) = trn(:,:,:,jp_fabm_n1p)
2773	END WHERE
2774	#else
2775	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2776	#endif
2777	ENDIF
2778
2779	IF ( ln_pdicinc ) THEN
2780	#if defined key_hadocc
2781	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2782	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) > 0.0_wp )
2783	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:)
2784	trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic)
2785	END WHERE
2786	#elif defined key_medusa
2787	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2788	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) > 0.0_wp )
2789	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:)
2790	trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
2791	END WHERE
2792	#elif defined key_fabm
2793	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2794	& trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) > 0.0_wp )
2795	trn(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:)
2796	trb(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c)
2797	END WHERE
2798	#else
2799	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2800	#endif
2801	ENDIF
2802
2803	IF ( ln_palkinc ) THEN
2804	#if defined key_hadocc
2805	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2806	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) > 0.0_wp )
2807	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:)
2808	trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk)
2809	END WHERE
2810	#elif defined key_medusa
2811	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2812	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) > 0.0_wp )
2813	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:)
2814	trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
2815	END WHERE
2816	#elif defined key_fabm
2817	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2818	& trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) > 0.0_wp )
2819	trn(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:)
2820	trb(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba)
2821	END WHERE
2822	#else
2823	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2824	#endif
2825	ENDIF
2826
2827	IF ( ln_po2inc ) THEN
2828	#if defined key_medusa
2829	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2830	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) > 0.0_wp )
2831	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:)
2832	trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
2833	END WHERE
2834	#elif defined key_fabm
2835	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2836	& trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) > 0.0_wp )
2837	trn(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:)
2838	trb(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o)
2839	END WHERE
2840	#else
2841	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2842	#endif
2843	ENDIF
2844
2845	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
2846	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
2847	IF ( ln_ppo4inc ) DEALLOCATE( ppo4_bkginc )
2848	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
2849	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
2850	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
2851	ENDIF
2852	!
2853	ENDIF
2854	!
2855	END SUBROUTINE bgc3d_asm_inc
2856
2857	!!===========================================================================
2858
2859	END MODULE asmbgc

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

Download in other formats: