Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

asmbgc.F90

Last change on this file was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 102.4 KB

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90

Download in other formats: