Context Navigation

asmbgc.F90 @ 13355

Last change on this file since 13355 was 13355, checked in by jwhile, 4 years ago
Merged in changes to fix 1d running - documented in UKMO ocean utils ticket 367
File size: 108.0 KB

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

Download in other formats: