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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_phytobal3d/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 13097

Last change on this file since 13097 was 13097, checked in by dford, 4 years ago
Allow nitrogen balancing for both 2D and 3D chlorophyll increments.
File size: 108.2 KB

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

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

Download in other formats: