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/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 10302

Last change on this file since 10302 was 10302, checked in by dford, 5 years ago
Merge in revisions 8447:10159 of dev_r5518_GO6_package.
File size: 102.3 KB

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

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

Download in other formats: