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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 9322

Last change on this file since 9322 was 9322, checked in by dford, 6 years ago
Add increments from various BGC profiles.
File size: 68.7 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_foam_medusa' : MEDUSA extras for FOAM OBS and ASM
13	!! 'key_roam' : MEDUSA extras for carbonate cycle
14	!! 'key_karaml' : Kara mixed layer depth
15	!!---------------------------------------------------------------------------
16	!! asm_bgc_check_options : check bgc assimilation options
17	!! asm_bgc_init_incs : initialise bgc increments
18	!! asm_bgc_init_bkg : initialise bgc background
19	!! asm_bgc_bal_wri : write out bgc balancing increments
20	!! asm_bgc_bkg_wri : write out bgc background
21	!! phyto_asm_inc : apply the ocean colour increments
22	!! pco2_asm_inc : apply the pCO2/fCO2 increments
23	!! ph_asm_inc : apply the pH increments
24	!! bgc3d_asm_inc : apply the generic 3D BGC increments
25	!!---------------------------------------------------------------------------
26	USE par_kind, ONLY: & ! kind parameters
27	& wp
28	USE par_oce, ONLY: & ! domain array sizes
29	& jpi, &
30	& jpj, &
31	& jpk
32	USE iom ! i/o
33	USE oce, ONLY: & ! active tracer variables
34	& tsn
35	USE zdfmxl, ONLY : & ! mixed layer depth
36	#if defined key_karaml
37	& hmld_kara, &
38	& ln_kara, &
39	#endif
40	& hmld, &
41	& hmlp, &
42	& hmlpt
43	USE asmpar, ONLY: & ! assimilation parameters
44	& c_asmbkg, &
45	& c_asmbal, &
46	& nitbkg_r, &
47	& nitdin_r, &
48	& nitiaustr_r, &
49	& nitiaufin_r
50	#if defined key_top
51	USE trc, ONLY: & ! passive tracer variables
52	& trn, &
53	& trb
54	USE par_trc, ONLY: & ! passive tracer parameters
55	& jptra
56	#endif
57	#if defined key_medusa && defined key_foam_medusa
58	USE asmlogchlbal_medusa, ONLY: & ! logchl balancing for MEDUSA
59	& asm_logchl_bal_medusa
60	USE asmpco2bal, ONLY: & ! pCO2 balancing for MEDUSA
61	& asm_pco2_bal
62	USE par_medusa ! MEDUSA parameters
63	USE mocsy_mainmod, ONLY: & ! MEDUSA carbonate system
64	& f2pCO2
65	USE sms_medusa, ONLY: & ! MEDUSA diagnostic variables
66	& pgrow_avg, &
67	& ploss_avg, &
68	& phyt_avg, &
69	& mld_max
70	#elif defined key_hadocc
71	USE asmlogchlbal_hadocc, ONLY: & ! logchl balancing for HadOCC
72	& asm_logchl_bal_hadocc
73	USE asmpco2bal, ONLY: & ! pCO2 balancing for HadOCC
74	& asm_pco2_bal
75	USE par_hadocc ! HadOCC parameters
76	USE trc, ONLY: & ! HadOCC diagnostic variables
77	& pgrow_avg, &
78	& ploss_avg, &
79	& phyt_avg, &
80	& mld_max, &
81	& HADOCC_CHL
82	USE had_bgc_const, ONLY: & ! HadOCC C:Chl ratio
83	& cchl_p
84	#endif
85
86	IMPLICIT NONE
87	PRIVATE
88
89	PUBLIC asm_bgc_check_options ! called by asm_inc_init in asminc.F90
90	PUBLIC asm_bgc_init_incs ! called by asm_inc_init in asminc.F90
91	PRIVATE asm_bgc_read_incs_2d ! called by asm_bgc_init_incs
92	PRIVATE asm_bgc_read_incs_3d ! called by asm_bgc_init_incs
93	PUBLIC asm_bgc_init_bkg ! called by asm_inc_init in asminc.F90
94	PUBLIC asm_bgc_bal_wri ! called by nemo_gcm in nemogcm.F90
95	PUBLIC asm_bgc_bkg_wri ! called by asm_bkg_wri in asmbkg.F90
96	PUBLIC phyto_asm_inc ! called by bgc_asm_inc in asminc.F90
97	PUBLIC pco2_asm_inc ! called by bgc_asm_inc in asminc.F90
98	PUBLIC ph_asm_inc ! called by bgc_asm_inc in asminc.F90
99	PUBLIC bgc3d_asm_inc ! called by bgc_asm_inc in asminc.F90
100
101	LOGICAL, PUBLIC :: ln_balwri = .FALSE. !: No output of balancing incs
102	LOGICAL, PUBLIC :: ln_phytobal = .FALSE. !: No phytoplankton balancing
103	LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE. !: No surface total log10(chlorophyll) increment
104	LOGICAL, PUBLIC :: ln_slchldiainc = .FALSE. !: No surface diatom log10(chlorophyll) increment
105	LOGICAL, PUBLIC :: ln_slchlnoninc = .FALSE. !: No surface non-diatom log10(chlorophyll) increment
106	LOGICAL, PUBLIC :: ln_schltotinc = .FALSE. !: No surface total chlorophyll increment
107	LOGICAL, PUBLIC :: ln_slphytotinc = .FALSE. !: No surface total log10(phyto C) increment
108	LOGICAL, PUBLIC :: ln_slphydiainc = .FALSE. !: No surface diatom log10(phyto C) increment
109	LOGICAL, PUBLIC :: ln_slphynoninc = .FALSE. !: No surface non-diatom log10(phyto C) increment
110	LOGICAL, PUBLIC :: ln_spco2inc = .FALSE. !: No surface pCO2 increment
111	LOGICAL, PUBLIC :: ln_sfco2inc = .FALSE. !: No surface fCO2 increment
112	LOGICAL, PUBLIC :: ln_plchltotinc = .FALSE. !: No profile total log10(chlorophyll) increment
113	LOGICAL, PUBLIC :: ln_pchltotinc = .FALSE. !: No profile total chlorophyll increment
114	LOGICAL, PUBLIC :: ln_pno3inc = .FALSE. !: No profile nitrate increment
115	LOGICAL, PUBLIC :: ln_psi4inc = .FALSE. !: No profile silicate increment
116	LOGICAL, PUBLIC :: ln_pdicinc = .FALSE. !: No profile dissolved inorganic carbon increment
117	LOGICAL, PUBLIC :: ln_palkinc = .FALSE. !: No profile alkalinity increment
118	LOGICAL, PUBLIC :: ln_pphinc = .FALSE. !: No profile pH increment
119	LOGICAL, PUBLIC :: ln_po2inc = .FALSE. !: No profile oxygen increment
120
121	INTEGER, PUBLIC :: mld_choice_bgc = 1 !: choice of mld for bgc assimilation
122	! 1) hmld - Turbocline/mixing depth
123	! [W points]
124	! 2) hmlp - Density criterion
125	! (0.01 kg/m^3 change from 10m)
126	! [W points]
127	! 3) hmld_kara - Kara MLD
128	! [Interpolated]
129	! 4) hmld_tref - Temperature criterion
130	! (0.2 K change from surface)
131	! [T points]
132	! 5) hmlpt - Density criterion
133	! (0.01 kg/m^3 change from 10m)
134	! [T points]
135
136	REAL(wp), PUBLIC :: rn_maxchlinc = -999.0 !: maximum absolute non-log chlorophyll inc
137	! <= 0 no maximum applied (switch turned off)
138	! > 0 absolute chl inc capped at this value
139
140	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchltot_bkginc ! slchltot inc
141	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchldia_bkginc ! slchldia inc
142	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchlnon_bkginc ! slchlnon inc
143	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: schltot_bkginc ! schltot inc
144	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphytot_bkginc ! slphytot inc
145	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphydia_bkginc ! slphydia inc
146	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slphynon_bkginc ! slphynon inc
147	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sfco2_bkginc ! sfco2 inc
148	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: spco2_bkginc ! spco2 inc
149	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: plchltot_bkginc ! plchltot inc
150	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pchltot_bkginc ! pchltot inc
151	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pno3_bkginc ! pno3 inc
152	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: psi4_bkginc ! psi4 inc
153	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pdic_bkginc ! pdic inc
154	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: palk_bkginc ! palk inc
155	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pph_bkginc ! pph inc
156	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: po2_bkginc ! po2 inc
157	#if defined key_top
158	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto_balinc ! Balancing incs from ocean colour
159	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc ! Balancing incs from spco2/sfco2
160	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ph_balinc ! Balancing incs from pph
161	#endif
162
163	#if defined key_hadocc \|\| (defined key_medusa && defined key_foam_medusa)
164	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pgrow_avg_bkg ! Background phyto growth
165	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ploss_avg_bkg ! Background phyto loss
166	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: phyt_avg_bkg ! Background phyto conc
167	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mld_max_bkg ! Background max MLD
168	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg ! Background tracer state
169	#endif
170	#if defined key_hadocc
171	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: chl_bkg ! Background chl
172	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: cchl_p_bkg ! Background c:chl
173	#endif
174
175	CONTAINS
176
177	SUBROUTINE asm_bgc_check_options
178	!!------------------------------------------------------------------------
179	!! * ROUTINE asm_bgc_check_options *
180	!!
181	!! ** Purpose : check validity of biogeochemical assimilation options
182	!!
183	!! ** Method : check settings of logicals
184	!!
185	!! ** Action : call ctl_stop/ctl_warn if required
186	!!
187	!! References : asm_inc_init
188	!!------------------------------------------------------------------------
189
190	#if ! defined key_top \|\| ( ! defined key_hadocc && ( ! defined key_medusa \|\| ! defined key_foam_medusa ) )
191	CALL ctl_stop( ' Attempting to assimilate biogeochemical observations', &
192	& ' but no compatible biogeochemical model is available' )
193	#endif
194
195	#if defined key_hadocc
196	IF ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slphydiainc .OR. &
197	& ln_slphynoninc .OR. ln_psi4inc .OR. ln_pphinc .OR. ln_po2inc ) THEN
198	CALL ctl_stop( ' Cannot assimilate PFTs, Si4, pH or O2 into HadOCC' )
199	ENDIF
200	#endif
201
202	IF ( ( ln_phytobal ).AND. &
203	& ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
204	& ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc ).AND. &
205	& ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
206	& ( .NOT. ln_slphynoninc ) ) THEN
207	CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', &
208	& ' if not assimilating ocean colour,', &
209	& ' so ln_phytobal will be set to .false.')
210	ln_phytobal = .FALSE.
211	ENDIF
212
213	IF ( ( ln_balwri ).AND.( .NOT. ln_phytobal ).AND. &
214	& ( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ).AND.( .NOT. ln_pphinc ) ) THEN
215	CALL ctl_warn( ' Balancing increments are only calculated for ocean colour', &
216	& ' with ln_phytobal, or for pCO2, fCO2, or pH.', &
217	& ' Not the case, so ln_balwri will be set to .false.')
218	ln_balwri = .FALSE.
219	ENDIF
220
221	IF ( ln_spco2inc .AND. ln_sfco2inc ) THEN
222	CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
223	ENDIF
224
225	IF ( ln_slchltotinc .AND. ln_schltotinc ) THEN
226	CALL ctl_stop( ' Can only assimilate surface log10(chlorophyll) or chlorophyll, not both' )
227	ENDIF
228
229	IF ( ln_plchltotinc .AND. ln_pchltotinc ) THEN
230	CALL ctl_stop( ' Can only assimilate profile log10(chlorophyll) or chlorophyll, not both' )
231	ENDIF
232
233	IF ( ( ln_slchltotinc .OR. ln_schltotinc ) .AND. &
234	& ( ln_slchldiainc .OR. ln_slchlnoninc ) ) THEN
235	CALL ctl_stop( ' Can only assimilate total or PFT surface chlorophyll, not both' )
236	ENDIF
237
238	IF ( ln_slphytotinc .AND. ( ln_slphydiainc .OR. ln_slphynoninc ) ) THEN
239	CALL ctl_stop( ' Can only assimilate total or PFT surface phytoplankton carbon, not both' )
240	ENDIF
241
242	END SUBROUTINE asm_bgc_check_options
243
244	!!===========================================================================
245	!!===========================================================================
246	!!===========================================================================
247
248	SUBROUTINE asm_bgc_init_incs( knum )
249	!!------------------------------------------------------------------------
250	!! * ROUTINE asm_bgc_init_incs *
251	!!
252	!! ** Purpose : initialise bgc increments
253	!!
254	!! ** Method : allocate arrays and read increments from file
255	!!
256	!! ** Action : allocate arrays and read increments from file
257	!!
258	!! References : asm_inc_init
259	!!------------------------------------------------------------------------
260	!!
261	INTEGER, INTENT(in ) :: knum ! i/o unit of increments file
262	!!
263	!!------------------------------------------------------------------------
264
265	! Allocate and read increments
266
267	IF ( ln_slchltotinc ) THEN
268	ALLOCATE( slchltot_bkginc(jpi,jpj) )
269	CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc )
270	ENDIF
271
272	IF ( ln_slchldiainc ) THEN
273	ALLOCATE( slchldia_bkginc(jpi,jpj) )
274	CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc )
275	ENDIF
276
277	IF ( ln_slchlnoninc ) THEN
278	ALLOCATE( slchlnon_bkginc(jpi,jpj) )
279	CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc )
280	ENDIF
281
282	IF ( ln_schltotinc ) THEN
283	ALLOCATE( schltot_bkginc(jpi,jpj) )
284	CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc )
285	ENDIF
286
287	IF ( ln_slphytotinc ) THEN
288	ALLOCATE( slphytot_bkginc(jpi,jpj) )
289	CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc )
290	ENDIF
291
292	IF ( ln_slphydiainc ) THEN
293	ALLOCATE( slphydia_bkginc(jpi,jpj) )
294	CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc )
295	ENDIF
296
297	IF ( ln_slphynoninc ) THEN
298	ALLOCATE( slphynon_bkginc(jpi,jpj) )
299	CALL asm_bgc_read_incs_2d( knum, 'bckinslphynon', slphynon_bkginc )
300	ENDIF
301
302	IF ( ln_sfco2inc ) THEN
303	ALLOCATE( sfco2_bkginc(jpi,jpj) )
304	CALL asm_bgc_read_incs_2d( knum, 'bckinsfco2', sfco2_bkginc )
305	ENDIF
306
307	IF ( ln_spco2inc ) THEN
308	ALLOCATE( spco2_bkginc(jpi,jpj) )
309	CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc )
310	ENDIF
311
312	IF ( ln_plchltotinc ) THEN
313	ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) )
314	CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc )
315	ENDIF
316
317	IF ( ln_pchltotinc ) THEN
318	ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) )
319	CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc )
320	ENDIF
321
322	IF ( ln_pno3inc ) THEN
323	ALLOCATE( pno3_bkginc(jpi,jpj,jpk) )
324	CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc )
325	ENDIF
326
327	IF ( ln_psi4inc ) THEN
328	ALLOCATE( psi4_bkginc(jpi,jpj,jpk) )
329	CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc )
330	ENDIF
331
332	IF ( ln_pdicinc ) THEN
333	ALLOCATE( pdic_bkginc(jpi,jpj,jpk) )
334	CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc )
335	ENDIF
336
337	IF ( ln_palkinc ) THEN
338	ALLOCATE( palk_bkginc(jpi,jpj,jpk) )
339	CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc )
340	ENDIF
341
342	IF ( ln_pphinc ) THEN
343	ALLOCATE( pph_bkginc(jpi,jpj,jpk) )
344	CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc )
345	ENDIF
346
347	IF ( ln_po2inc ) THEN
348	ALLOCATE( po2_bkginc(jpi,jpj,jpk) )
349	CALL asm_bgc_read_incs_3d( knum, 'bckinpo2', po2_bkginc )
350	ENDIF
351
352	! Allocate balancing increments
353
354	IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
355	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
356	& ln_slphynoninc ) THEN
357	#if defined key_top
358	ALLOCATE( phyto_balinc(jpi,jpj,jpk,jptra) )
359	phyto_balinc(:,:,:,:) = 0.0
360	#else
361	CALL ctl_stop( ' key_top must be set for balancing increments' )
362	#endif
363	ENDIF
364
365	IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
366	#if defined key_top
367	ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) )
368	pco2_balinc(:,:,:,:) = 0.0
369	#else
370	CALL ctl_stop( ' key_top must be set for balancing increments' )
371	#endif
372	ENDIF
373
374	IF ( ln_pphinc ) THEN
375	#if defined key_top
376	ALLOCATE( ph_balinc(jpi,jpj,jpk,jptra) )
377	ph_balinc(:,:,:,:) = 0.0
378	#else
379	CALL ctl_stop( ' key_top must be set for balancing increments' )
380	#endif
381	ENDIF
382
383	END SUBROUTINE asm_bgc_init_incs
384
385	!!===========================================================================
386	!!===========================================================================
387	!!===========================================================================
388
389	SUBROUTINE asm_bgc_read_incs_2d( knum, cd_bgcname, p_incs )
390	!!------------------------------------------------------------------------
391	!! * ROUTINE asm_bgc_init_incs *
392	!!
393	!! ** Purpose : read 2d bgc increments
394	!!
395	!! ** Method : read increments from file
396	!!
397	!! ** Action : read increments from file
398	!!
399	!! References : asm_inc_init
400	!!------------------------------------------------------------------------
401	!!
402	INTEGER, INTENT(in ) :: knum ! i/o unit
403	CHARACTER(LEN=13), INTENT(in ) :: cd_bgcname ! variable
404	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: p_incs ! increments
405	!!
406	!!------------------------------------------------------------------------
407
408	! Initialise
409	p_incs(:,:) = 0.0
410
411	! read from file
412	CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 )
413
414	! Apply the masks
415	p_incs(:,:) = p_incs(:,:) * tmask(:,:,1)
416
417	! Set missing increments to 0.0 rather than 1e+20
418	! to allow for differences in masks
419	WHERE( ABS( p_incs(:,:) ) > 1.0e+10 ) p_incs(:,:) = 0.0
420
421	END SUBROUTINE asm_bgc_read_incs_2d
422
423	!!===========================================================================
424	!!===========================================================================
425	!!===========================================================================
426
427	SUBROUTINE asm_bgc_read_incs_3d( knum, cd_bgcname, p_incs )
428	!!------------------------------------------------------------------------
429	!! * ROUTINE asm_bgc_init_incs *
430	!!
431	!! ** Purpose : read 3d bgc increments
432	!!
433	!! ** Method : read increments from file
434	!!
435	!! ** Action : read increments from file
436	!!
437	!! References : asm_inc_init
438	!!------------------------------------------------------------------------
439	!!
440	INTEGER, INTENT(in ) :: knum ! i/o unit
441	CHARACTER(LEN=13), INTENT(in ) :: cd_bgcname ! variable
442	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: p_incs ! increments
443	!!
444	!!------------------------------------------------------------------------
445
446	! Initialise
447	p_incs(:,:,:) = 0.0
448
449	! read from file
450	CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 )
451
452	! Apply the masks
453	p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:)
454
455	! Set missing increments to 0.0 rather than 1e+20
456	! to allow for differences in masks
457	WHERE( ABS( p_incs(:,:,:) ) > 1.0e+10 ) p_incs(:,:,:) = 0.0
458
459	END SUBROUTINE asm_bgc_read_incs_3d
460
461	!!===========================================================================
462	!!===========================================================================
463	!!===========================================================================
464
465	SUBROUTINE asm_bgc_init_bkg
466	!!------------------------------------------------------------------------
467	!! * ROUTINE asm_bgc_init_bkg *
468	!!
469	!! ** Purpose : initialise bgc background
470	!!
471	!! ** Method : allocate arrays and read background from file
472	!!
473	!! ** Action : allocate arrays and read background from file
474	!!
475	!! References : asm_inc_init
476	!!------------------------------------------------------------------------
477	!!
478	INTEGER :: inum ! i/o unit of background file
479	INTEGER :: jt ! loop counter
480	!!
481	!!------------------------------------------------------------------------
482
483	#if defined key_top && ( defined key_hadocc \|\| (defined key_medusa && defined key_foam_medusa) )
484	IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
485	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
486	& ln_slphynoninc .OR. ln_plchltotinc .OR. ln_pchltotinc ) THEN
487
488	ALLOCATE( pgrow_avg_bkg(jpi,jpj) )
489	ALLOCATE( ploss_avg_bkg(jpi,jpj) )
490	ALLOCATE( phyt_avg_bkg(jpi,jpj) )
491	ALLOCATE( mld_max_bkg(jpi,jpj) )
492	ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
493	pgrow_avg_bkg(:,:) = 0.0
494	ploss_avg_bkg(:,:) = 0.0
495	phyt_avg_bkg(:,:) = 0.0
496	mld_max_bkg(:,:) = 0.0
497	tracer_bkg(:,:,:,:) = 0.0
498
499	#if defined key_hadocc
500	ALLOCATE( chl_bkg(jpi,jpj,jpk) )
501	ALLOCATE( cchl_p_bkg(jpi,jpj,jpk) )
502	chl_bkg(:,:,:) = 0.0
503	cchl_p_bkg(:,:,:) = 0.0
504	#endif
505
506	!--------------------------------------------------------------------
507	! Read background variables for phytoplankton assimilation
508	! Some only required if performing balancing
509	!--------------------------------------------------------------------
510
511	CALL iom_open( c_asmbkg, inum )
512
513	#if defined key_hadocc
514	CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl', chl_bkg )
515	CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg )
516	chl_bkg(:,:,:) = chl_bkg(:,:,:) * tmask(:,:,:)
517	cchl_p_bkg(:,:,:) = cchl_p_bkg(:,:,:) * tmask(:,:,:)
518	#elif defined key_medusa
519	CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) )
520	CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) )
521	#endif
522
523	IF ( ln_phytobal ) THEN
524
525	CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg )
526	CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg )
527	CALL iom_get( inum, jpdom_autoglo, 'phyt_avg', phyt_avg_bkg )
528	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
529	pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1)
530	ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1)
531	phyt_avg_bkg(:,:) = phyt_avg_bkg(:,:) * tmask(:,:,1)
532	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
533
534	#if defined key_hadocc
535	CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) )
536	CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) )
537	CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) )
538	CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) )
539	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
540	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
541	#elif defined key_medusa
542	CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) )
543	CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
544	CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
545	CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) )
546	CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) )
547	CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) )
548	CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) )
549	CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) )
550	CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) )
551	CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) )
552	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
553	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
554	CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
555	#endif
556	ELSE IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
557	#if defined key_hadocc
558	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
559	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
560	#elif defined key_medusa
561	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
562	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
563	#endif
564	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
565	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
566	ENDIF
567
568	CALL iom_close( inum )
569
570	DO jt = 1, jptra
571	tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
572	END DO
573
574	ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
575
576	ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
577	ALLOCATE( mld_max_bkg(jpi,jpj) )
578	tracer_bkg(:,:,:,:) = 0.0
579	mld_max_bkg(:,:) = 0.0
580
581	CALL iom_open( c_asmbkg, inum )
582
583	#if defined key_hadocc
584	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
585	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
586	#elif defined key_medusa
587	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
588	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
589	#endif
590	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
591
592	CALL iom_close( inum )
593
594	DO jt = 1, jptra
595	tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
596	END DO
597	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
598
599	ENDIF
600	#else
601	CALL ctl_stop( ' asm_bgc_init_bkg: key_top and a compatible biogeochemical model required' )
602	#endif
603
604	END SUBROUTINE asm_bgc_init_bkg
605
606	!!===========================================================================
607	!!===========================================================================
608	!!===========================================================================
609
610	SUBROUTINE asm_bgc_bal_wri( kt )
611	!!------------------------------------------------------------------------
612	!!
613	!! * ROUTINE asm_bgc_bal_wri *
614	!!
615	!! ** Purpose : Write to file the balancing increments
616	!! calculated for biogeochemistry
617	!!
618	!! ** Method : Write to file the balancing increments
619	!! calculated for biogeochemistry
620	!!
621	!! ** Action :
622	!!
623	!! References :
624	!!
625	!! History :
626	!! ! 2014-08 (D. Ford) Initial version, based on asm_bkg_wri
627	!!------------------------------------------------------------------------
628	!! * Arguments
629	INTEGER, INTENT( IN ) :: kt ! Current time-step
630	!! * Local declarations
631	CHARACTER(LEN=50) :: cl_asmbal ! Filename (with extension)
632	LOGICAL :: llok ! Check if file exists
633	INTEGER :: inum ! File unit number
634	REAL(wp) :: zdate ! Date
635	!!------------------------------------------------------------------------
636
637	! Set things up
638	zdate = REAL( ndastp )
639	WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal )
640
641	! Check if file exists
642	INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok )
643	IF( .NOT. llok ) THEN
644	IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// &
645	& TRIM( c_asmbal ) // ' at timestep = ', kt
646
647	! Define the output file
648	CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
649
650	! Write the information
651	CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate )
652
653	IF ( ln_slchltotinc ) THEN
654	#if defined key_medusa
655	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', phyto_balinc(:,:,:,jpchn) )
656	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', phyto_balinc(:,:,:,jpchd) )
657	IF ( ln_phytobal ) THEN
658	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phn', phyto_balinc(:,:,:,jpphn) )
659	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phd', phyto_balinc(:,:,:,jpphd) )
660	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_pds', phyto_balinc(:,:,:,jppds) )
661	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zmi', phyto_balinc(:,:,:,jpzmi) )
662	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zme', phyto_balinc(:,:,:,jpzme) )
663	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_din', phyto_balinc(:,:,:,jpdin) )
664	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_sil', phyto_balinc(:,:,:,jpsil) )
665	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_fer', phyto_balinc(:,:,:,jpfer) )
666	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', phyto_balinc(:,:,:,jpdet) )
667	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dtc', phyto_balinc(:,:,:,jpdtc) )
668	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', phyto_balinc(:,:,:,jpdic) )
669	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', phyto_balinc(:,:,:,jpalk) )
670	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_oxy', phyto_balinc(:,:,:,jpoxy) )
671	ENDIF
672	#elif defined key_hadocc
673	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phy', phyto_balinc(:,:,:,jp_had_phy) )
674	IF ( ln_phytobal ) THEN
675	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_nut', phyto_balinc(:,:,:,jp_had_nut) )
676	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zoo', phyto_balinc(:,:,:,jp_had_zoo) )
677	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', phyto_balinc(:,:,:,jp_had_pdn) )
678	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', phyto_balinc(:,:,:,jp_had_dic) )
679	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', phyto_balinc(:,:,:,jp_had_alk) )
680	ENDIF
681	#endif
682	ENDIF
683
684	IF ( ln_spco2inc ) THEN
685	#if defined key_medusa
686	CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jpdic) )
687	CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jpalk) )
688	#elif defined key_hadocc
689	CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) )
690	CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) )
691	#endif
692	ELSE IF ( ln_sfco2inc ) THEN
693	#if defined key_medusa
694	CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jpdic) )
695	CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jpalk) )
696	#elif defined key_hadocc
697	CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) )
698	CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) )
699	#endif
700	ENDIF
701
702	CALL iom_close( inum )
703	ELSE
704	CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', &
705	& ' Therefore not writing out balancing increments at this timestep', &
706	& ' Something has probably gone wrong somewhere' )
707	IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// &
708	& TRIM( c_asmbal ) // ' at timestep = ', kt
709	ENDIF
710
711	END SUBROUTINE asm_bgc_bal_wri
712
713	!!===========================================================================
714	!!===========================================================================
715	!!===========================================================================
716
717	SUBROUTINE asm_bgc_bkg_wri( kt, knum )
718	!!------------------------------------------------------------------------
719	!! * ROUTINE asm_bgc_bkg_wri *
720	!!
721	!! ** Purpose : write out bgc background
722	!!
723	!! ** Method : write out bgc background
724	!!
725	!! ** Action : write out bgc background
726	!!
727	!! References : asm_bkg_wri
728	!!------------------------------------------------------------------------
729	!!
730	INTEGER, INTENT(in ) :: kt ! Current time-step
731	INTEGER, INTENT(in ) :: knum ! i/o unit of increments file
732	!!
733	!!------------------------------------------------------------------------
734
735	#if defined key_hadocc
736	CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg' , pgrow_avg )
737	CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg' , ploss_avg )
738	CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg' , phyt_avg )
739	CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max' , mld_max )
740	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_nut' , trn(:,:,:,jp_had_nut) )
741	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_phy' , trn(:,:,:,jp_had_phy) )
742	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_zoo' , trn(:,:,:,jp_had_zoo) )
743	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_pdn' , trn(:,:,:,jp_had_pdn) )
744	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_dic' , trn(:,:,:,jp_had_dic) )
745	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_alk' , trn(:,:,:,jp_had_alk) )
746	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_chl' , HADOCC_CHL(:,:,:) )
747	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:) )
748	#elif defined key_medusa && defined key_foam_medusa
749	CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg' , pgrow_avg )
750	CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg' , ploss_avg )
751	CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg' , phyt_avg )
752	CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max' , mld_max )
753	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn' , trn(:,:,:,jpchn) )
754	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chd' , trn(:,:,:,jpchd) )
755	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phn' , trn(:,:,:,jpphn) )
756	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phd' , trn(:,:,:,jpphd) )
757	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_pds' , trn(:,:,:,jppds) )
758	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zmi' , trn(:,:,:,jpzmi) )
759	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zme' , trn(:,:,:,jpzme) )
760	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_din' , trn(:,:,:,jpdin) )
761	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_sil' , trn(:,:,:,jpsil) )
762	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_fer' , trn(:,:,:,jpfer) )
763	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_det' , trn(:,:,:,jpdet) )
764	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dtc' , trn(:,:,:,jpdtc) )
765	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dic' , trn(:,:,:,jpdic) )
766	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_alk' , trn(:,:,:,jpalk) )
767	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_oxy' , trn(:,:,:,jpoxy) )
768	#endif
769
770	END SUBROUTINE asm_bgc_bkg_wri
771
772	!!===========================================================================
773	!!===========================================================================
774	!!===========================================================================
775
776	SUBROUTINE phyto_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
777	!!------------------------------------------------------------------------
778	!! * ROUTINE phyto_asm_inc *
779	!!
780	!! ** Purpose : Apply the chlorophyll assimilation increments.
781	!!
782	!! ** Method : Calculate increments to state variables using nitrogen
783	!! balancing.
784	!! Direct initialization or Incremental Analysis Updating.
785	!!
786	!! ** Action :
787	!!------------------------------------------------------------------------
788	INTEGER, INTENT(IN) :: kt ! Current time step
789	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
790	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
791	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
792	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
793	!
794	INTEGER :: jk ! Loop counter
795	INTEGER :: it ! Index
796	REAL(wp) :: zincwgt ! IAU weight for current time step
797	REAL(wp) :: zincper ! IAU interval in seconds
798	!!------------------------------------------------------------------------
799
800	IF ( ln_slchldiainc .OR. ln_slchlnoninc .OR. &
801	& ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
802	& ln_slphynoninc ) THEN
803	CALL ctl_stop( ' No PFT assimilation quite yet' )
804	ENDIF
805
806	IF ( kt <= nit000 ) THEN
807
808	zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
809
810	#if defined key_medusa && defined key_foam_medusa
811	CALL asm_logchl_bal_medusa( slchltot_bkginc, zincper, mld_choice_bgc, &
812	& rn_maxchlinc, ln_phytobal, ll_asmdin, &
813	& pgrow_avg_bkg, ploss_avg_bkg, &
814	& phyt_avg_bkg, mld_max_bkg, &
815	& tracer_bkg, phyto_balinc )
816	#elif defined key_hadocc
817	CALL asm_logchl_bal_hadocc( slchltot_bkginc, zincper, mld_choice_bgc, &
818	& rn_maxchlinc, ln_phytobal, ll_asmdin, &
819	& pgrow_avg_bkg, ploss_avg_bkg, &
820	& phyt_avg_bkg, mld_max_bkg, &
821	& chl_bkg(:,:,1), cchl_p_bkg(:,:,1), &
822	& tracer_bkg, phyto_balinc )
823	#else
824	CALL ctl_stop( 'Attempting to assimilate slchltot, ', &
825	& 'but not defined a biogeochemical model' )
826	#endif
827
828	ENDIF
829
830	IF ( ll_asmiau ) THEN
831
832	!--------------------------------------------------------------------
833	! Incremental Analysis Updating
834	!--------------------------------------------------------------------
835
836	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
837
838	it = kt - nit000 + 1
839	zincwgt = pwgtiau(it) ! IAU weight for the current time step
840	! note this is not a tendency so should not be divided by rdt
841
842	IF(lwp) THEN
843	WRITE(numout,*)
844	WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', &
845	& kt,' with IAU weight = ', pwgtiau(it)
846	WRITE(numout,*) '~~~~~~~~~~~~'
847	ENDIF
848
849	! Update the biogeochemical variables
850	! Add directly to trn and trb, rather than to tra, because tra gets
851	! reset to zero at the start of trc_stp, called after this routine
852	#if defined key_medusa && defined key_foam_medusa
853	DO jk = 1, jpkm1
854	trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
855	& phyto_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
856	trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
857	& phyto_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
858	END DO
859	#elif defined key_hadocc
860	DO jk = 1, jpkm1
861	trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
862	& phyto_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
863	trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
864	& phyto_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
865	END DO
866	#endif
867
868	! Do not deallocate arrays - needed by asm_bgc_bal_wri
869	! which is called at end of model run
870	ENDIF
871
872	ELSEIF ( ll_asmdin ) THEN
873
874	!--------------------------------------------------------------------
875	! Direct Initialization
876	!--------------------------------------------------------------------
877
878	IF ( kt == nitdin_r ) THEN
879
880	neuler = 0 ! Force Euler forward step
881
882	#if defined key_medusa && defined key_foam_medusa
883	! Initialize the now fields with the background + increment
884	! Background currently is what the model is initialised with
885	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', &
886	& ' Background state is taken from model rather than background file' )
887	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
888	& phyto_balinc(:,:,:,jp_msa0:jp_msa1)
889	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
890	#elif defined key_hadocc
891	! Initialize the now fields with the background + increment
892	! Background currently is what the model is initialised with
893	CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', &
894	& ' Background state is taken from model rather than background file' )
895	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
896	& phyto_balinc(:,:,:,jp_had0:jp_had1)
897	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
898	#endif
899
900	! Do not deallocate arrays - needed by asm_bgc_bal_wri
901	! which is called at end of model run
902	ENDIF
903	!
904	ENDIF
905	!
906	END SUBROUTINE phyto_asm_inc
907
908	!!===========================================================================
909	!!===========================================================================
910	!!===========================================================================
911
912	SUBROUTINE pco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
913	& ll_trainc, pt_bkginc, ps_bkginc )
914	!!------------------------------------------------------------------------
915	!! * ROUTINE pco2_asm_inc *
916	!!
917	!! ** Purpose : Apply the pco2/fco2 assimilation increments.
918	!!
919	!! ** Method : Calculate increments to state variables using carbon
920	!! balancing.
921	!! Direct initialization or Incremental Analysis Updating.
922	!!
923	!! ** Action :
924	!!------------------------------------------------------------------------
925	INTEGER, INTENT(IN) :: kt ! Current time step
926	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
927	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
928	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
929	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
930	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
931	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
932	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
933	!
934	INTEGER :: ji, jj, jk ! Loop counters
935	INTEGER :: it ! Index
936	INTEGER :: jkmax ! Index of mixed layer depth
937	REAL(wp) :: zincwgt ! IAU weight for current time step
938	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth for balancing
939	REAL(wp), DIMENSION(jpi,jpj) :: pco2_bkginc_temp ! For fCO2 conversion
940	REAL(wp), DIMENSION(jpi,jpj) :: dic_bkg_temp ! DIC background
941	REAL(wp), DIMENSION(jpi,jpj) :: alk_bkg_temp ! Alkalinity background
942	REAL(wp), DIMENSION(jpi,jpj) :: tem_bkg_temp ! Temperature background
943	REAL(wp), DIMENSION(jpi,jpj) :: sal_bkg_temp ! Salinity background
944	REAL(wp), DIMENSION(1) :: patm ! Atmospheric pressure
945	REAL(wp), DIMENSION(1) :: phyd ! Hydrostatic pressure
946
947	! Coefficients for fCO2 to pCO2 conversion
948	REAL(wp) :: zcoef_fco2_1 = -1636.75
949	REAL(wp) :: zcoef_fco2_2 = 12.0408
950	REAL(wp) :: zcoef_fco2_3 = 0.0327957
951	REAL(wp) :: zcoef_fco2_4 = 0.0000316528
952	REAL(wp) :: zcoef_fco2_5 = 2.0
953	REAL(wp) :: zcoef_fco2_6 = 57.7
954	REAL(wp) :: zcoef_fco2_7 = 0.118
955	REAL(wp) :: zcoef_fco2_8 = 82.0578
956	!!------------------------------------------------------------------------
957
958	IF ( kt <= nit000 ) THEN
959
960	!----------------------------------------------------------------------
961	! DIC and alkalinity balancing
962	!----------------------------------------------------------------------
963
964	IF ( ln_sfco2inc ) THEN
965	#if defined key_medusa && defined key_roam
966	! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
967	patm(1) = 1.0
968	phyd(1) = 0.0
969	DO jj = 1, jpj
970	DO ji = 1, jpi
971	CALL f2pCO2( sfco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
972	END DO
973	END DO
974	#else
975	! If assimilating fCO2, then convert to pCO2 using temperature
976	! See flux_gas.F90 within HadOCC for details of calculation
977	pco2_bkginc_temp(:,:) = sfco2_bkginc(:,:) / &
978	& EXP((zcoef_fco2_1 + &
979	& zcoef_fco2_2 * (tsn(:,:,1,1)+rt0) - &
980	& zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
981	& zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0) + &
982	& zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0))) / &
983	& (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
984	#endif
985	ELSE
986	pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
987	ENDIF
988
989	! Account for physics assimilation if required
990	IF ( ll_trainc ) THEN
991	tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
992	sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
993	ELSE
994	tem_bkg_temp(:,:) = tsn(:,:,1,1)
995	sal_bkg_temp(:,:) = tsn(:,:,1,2)
996	ENDIF
997
998	#if defined key_medusa
999	! Account for logchl balancing if required
1000	IF ( ln_slchltotinc .AND. ln_phytobal ) THEN
1001	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto_balinc(:,:,1,jpdic)
1002	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto_balinc(:,:,1,jpalk)
1003	ELSE
1004	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
1005	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
1006	ENDIF
1007
1008	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1009	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
1010	& pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
1011
1012	#elif defined key_hadocc
1013	! Account for slchltot balancing if required
1014	IF ( ln_slchltotinc .AND. ln_phytobal ) THEN
1015	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto_balinc(:,:,1,jp_had_dic)
1016	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto_balinc(:,:,1,jp_had_alk)
1017	ELSE
1018	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
1019	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
1020	ENDIF
1021
1022	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1023	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
1024	& pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
1025
1026	#else
1027	CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
1028	& 'but not defined a biogeochemical model' )
1029	#endif
1030
1031	! Select mixed layer
1032	IF ( ll_asmdin ) THEN
1033	#if defined key_hadocc \|\| (defined key_medusa && defined key_foam_medusa)
1034	CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
1035	& ' Mixed layer depth taken to be background maximum mld_max_bkg' )
1036	zmld(:,:) = mld_max_bkg(:,:)
1037	#else
1038	CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
1039	#endif
1040	ELSE
1041	SELECT CASE( mld_choice_bgc )
1042	CASE ( 1 ) ! Turbocline/mixing depth [W points]
1043	zmld(:,:) = hmld(:,:)
1044	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
1045	zmld(:,:) = hmlp(:,:)
1046	CASE ( 3 ) ! Kara MLD [Interpolated]
1047	#if defined key_karaml
1048	IF ( ln_kara ) THEN
1049	zmld(:,:) = hmld_kara(:,:)
1050	ELSE
1051	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1052	& ' but ln_kara=.false.' )
1053	ENDIF
1054	#else
1055	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1056	& ' but is not defined' )
1057	#endif
1058	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
1059	!zmld(:,:) = hmld_tref(:,:)
1060	CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
1061	& ' but is not available in this version' )
1062	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
1063	zmld(:,:) = hmlpt(:,:)
1064	END SELECT
1065	ENDIF
1066
1067	! Propagate through mixed layer
1068	DO jj = 1, jpj
1069	DO ji = 1, jpi
1070	!
1071	jkmax = jpk-1
1072	DO jk = jpk-1, 1, -1
1073	IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. &
1074	& ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
1075	zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
1076	jkmax = jk
1077	ENDIF
1078	END DO
1079	!
1080	#if defined key_top
1081	DO jk = 2, jkmax
1082	pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
1083	END DO
1084	#endif
1085	!
1086	END DO
1087	END DO
1088
1089	ENDIF
1090
1091	IF ( ll_asmiau ) THEN
1092
1093	!--------------------------------------------------------------------
1094	! Incremental Analysis Updating
1095	!--------------------------------------------------------------------
1096
1097	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1098
1099	it = kt - nit000 + 1
1100	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1101	! note this is not a tendency so should not be divided by rdt
1102
1103	IF(lwp) THEN
1104	WRITE(numout,*)
1105	IF ( ln_spco2inc ) THEN
1106	WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
1107	& kt,' with IAU weight = ', pwgtiau(it)
1108	ELSE IF ( ln_sfco2inc ) THEN
1109	WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
1110	& kt,' with IAU weight = ', pwgtiau(it)
1111	ENDIF
1112	WRITE(numout,*) '~~~~~~~~~~~~'
1113	ENDIF
1114
1115	! Update the biogeochemical variables
1116	! Add directly to trn and trb, rather than to tra, because tra gets
1117	! reset to zero at the start of trc_stp, called after this routine
1118	#if defined key_medusa && defined key_foam_medusa
1119	DO jk = 1, jpkm1
1120	trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
1121	& pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
1122	trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
1123	& pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
1124	END DO
1125	#elif defined key_hadocc
1126	DO jk = 1, jpkm1
1127	trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
1128	& pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
1129	trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
1130	& pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
1131	END DO
1132	#endif
1133
1134	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1135	! which is called at end of model run
1136
1137	ENDIF
1138
1139	ELSEIF ( ll_asmdin ) THEN
1140
1141	!--------------------------------------------------------------------
1142	! Direct Initialization
1143	!--------------------------------------------------------------------
1144
1145	IF ( kt == nitdin_r ) THEN
1146
1147	neuler = 0 ! Force Euler forward step
1148
1149	#if defined key_medusa && defined key_foam_medusa
1150	! Initialize the now fields with the background + increment
1151	! Background currently is what the model is initialised with
1152	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', &
1153	& ' Background state is taken from model rather than background file' )
1154	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1155	& pco2_balinc(:,:,:,jp_msa0:jp_msa1)
1156	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1157	#elif defined key_hadocc
1158	! Initialize the now fields with the background + increment
1159	! Background currently is what the model is initialised with
1160	CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', &
1161	& ' Background state is taken from model rather than background file' )
1162	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1163	& pco2_balinc(:,:,:,jp_had0:jp_had1)
1164	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1165	#endif
1166
1167	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1168	! which is called at end of model run
1169	ENDIF
1170	!
1171	ENDIF
1172	!
1173	END SUBROUTINE pco2_asm_inc
1174
1175	!!===========================================================================
1176	!!===========================================================================
1177	!!===========================================================================
1178
1179	SUBROUTINE ph_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
1180	& ll_trainc, pt_bkginc, ps_bkginc )
1181	!!------------------------------------------------------------------------
1182	!! * ROUTINE ph_asm_inc *
1183	!!
1184	!! ** Purpose : Apply the pH assimilation increments.
1185	!!
1186	!! ** Method : Calculate increments to state variables using carbon
1187	!! balancing.
1188	!! Direct initialization or Incremental Analysis Updating.
1189	!!
1190	!! ** Action :
1191	!!------------------------------------------------------------------------
1192	INTEGER, INTENT(IN) :: kt ! Current time step
1193	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1194	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1195	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1196	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1197	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
1198	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
1199	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
1200	!!------------------------------------------------------------------------
1201
1202	CALL ctl_stop( ' pH balancing not yet implemented' )
1203
1204	!
1205	END SUBROUTINE ph_asm_inc
1206
1207	!!===========================================================================
1208	!!===========================================================================
1209	!!===========================================================================
1210
1211	SUBROUTINE bgc3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1212	!!----------------------------------------------------------------------
1213	!! * ROUTINE dyn_asm_inc *
1214	!!
1215	!! ** Purpose : Apply generic 3D biogeochemistry assimilation increments.
1216	!!
1217	!! ** Method : Direct initialization or Incremental Analysis Updating.
1218	!!
1219	!! ** Action :
1220	!!----------------------------------------------------------------------
1221	INTEGER, INTENT(IN) :: kt ! Current time step
1222	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1223	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1224	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1225	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1226	!
1227	INTEGER :: jk ! Loop counter
1228	INTEGER :: it ! Index
1229	REAL(wp) :: zincwgt ! IAU weight for current time step
1230	REAL(wp) :: zincper ! IAU interval in seconds
1231	!!----------------------------------------------------------------------
1232
1233	IF ( ll_asmiau ) THEN
1234
1235	!--------------------------------------------------------------------
1236	! Incremental Analysis Updating
1237	!--------------------------------------------------------------------
1238
1239	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1240
1241	it = kt - nit000 + 1
1242	!zincwgt = pwgtiau(it) / rdt ! IAU weight for the current time step
1243	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1244	! Check which we should use both here and for all others
1245
1246	IF(lwp) THEN
1247	WRITE(numout,*)
1248	WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', &
1249	& kt,' with IAU weight = ', pwgtiau(it)
1250	WRITE(numout,*) '~~~~~~~~~~~~'
1251	ENDIF
1252
1253	! Update the 3D BGC variables
1254	! Add directly to trn and trb, rather than to tra, because tra gets
1255	! reset to zero at the start of trc_stp, called after this routine
1256	! Don't apply increments if they'll take concentrations negative
1257
1258	IF ( ln_pno3inc ) THEN
1259	#if defined key_hadocc
1260	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
1261	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
1262	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
1263	trb(:,:,:,jp_had_nut) = trb(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
1264	END WHERE
1265	#elif defined key_medusa && defined key_foam_medusa
1266	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
1267	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
1268	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
1269	trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
1270	END WHERE
1271	#else
1272	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1273	#endif
1274	ENDIF
1275
1276	IF ( ln_psi4inc ) THEN
1277	#if defined key_medusa && defined key_foam_medusa
1278	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
1279	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
1280	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
1281	trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
1282	END WHERE
1283	#else
1284	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1285	#endif
1286	ENDIF
1287
1288	IF ( ln_pdicinc ) THEN
1289	#if defined key_hadocc
1290	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
1291	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
1292	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
1293	trb(:,:,:,jp_had_dic) = trb(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
1294	END WHERE
1295	#elif defined key_medusa && defined key_foam_medusa
1296	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
1297	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
1298	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
1299	trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
1300	END WHERE
1301	#else
1302	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1303	#endif
1304	ENDIF
1305
1306	IF ( ln_palkinc ) THEN
1307	#if defined key_hadocc
1308	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
1309	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
1310	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
1311	trb(:,:,:,jp_had_alk) = trb(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
1312	END WHERE
1313	#elif defined key_medusa && defined key_foam_medusa
1314	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
1315	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
1316	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
1317	trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
1318	END WHERE
1319	#else
1320	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1321	#endif
1322	ENDIF
1323
1324	IF ( ln_po2inc ) THEN
1325	#if defined key_medusa && defined key_foam_medusa
1326	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
1327	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
1328	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
1329	trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
1330	END WHERE
1331	#else
1332	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1333	#endif
1334	ENDIF
1335
1336	IF ( kt == nitiaufin_r ) THEN
1337	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
1338	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
1339	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
1340	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
1341	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
1342	ENDIF
1343
1344	ENDIF
1345
1346	ELSEIF ( ll_asmdin ) THEN
1347
1348	!--------------------------------------------------------------------
1349	! Direct Initialization
1350	!--------------------------------------------------------------------
1351
1352	IF ( kt == nitdin_r ) THEN
1353
1354	neuler = 0 ! Force Euler forward step
1355
1356	! Initialize the now fields with the background + increment
1357	! Background currently is what the model is initialised with
1358	#if defined key_hadocc
1359	CALL ctl_warn( ' Doing direct initialisation of HadOCC with 3D BGC assimilation', &
1360	& ' Background state is taken from model rather than background file' )
1361	#elif defined key_medusa && defined key_foam_medusa
1362	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with 3D BGC assimilation', &
1363	& ' Background state is taken from model rather than background file' )
1364	#endif
1365
1366	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1367	& pco2_balinc(:,:,:,jp_msa0:jp_msa1)
1368	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1369
1370	IF ( ln_pno3inc ) THEN
1371	#if defined key_hadocc
1372	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
1373	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) > 0.0_wp )
1374	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:)
1375	trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut)
1376	END WHERE
1377	#elif defined key_medusa && defined key_foam_medusa
1378	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
1379	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) > 0.0_wp )
1380	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:)
1381	trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
1382	END WHERE
1383	#else
1384	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1385	#endif
1386	ENDIF
1387
1388	IF ( ln_psi4inc ) THEN
1389	#if defined key_medusa && defined key_foam_medusa
1390	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
1391	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) > 0.0_wp )
1392	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:)
1393	trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
1394	END WHERE
1395	#else
1396	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1397	#endif
1398	ENDIF
1399
1400	IF ( ln_pdicinc ) THEN
1401	#if defined key_hadocc
1402	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
1403	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) > 0.0_wp )
1404	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:)
1405	trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic)
1406	END WHERE
1407	#elif defined key_medusa && defined key_foam_medusa
1408	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
1409	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) > 0.0_wp )
1410	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:)
1411	trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
1412	END WHERE
1413	#else
1414	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1415	#endif
1416	ENDIF
1417
1418	IF ( ln_palkinc ) THEN
1419	#if defined key_hadocc
1420	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
1421	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) > 0.0_wp )
1422	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:)
1423	trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk)
1424	END WHERE
1425	#elif defined key_medusa && defined key_foam_medusa
1426	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
1427	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) > 0.0_wp )
1428	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:)
1429	trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
1430	END WHERE
1431	#else
1432	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1433	#endif
1434	ENDIF
1435
1436	IF ( ln_po2inc ) THEN
1437	#if defined key_medusa && defined key_foam_medusa
1438	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
1439	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) > 0.0_wp )
1440	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:)
1441	trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
1442	END WHERE
1443	#else
1444	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1445	#endif
1446	ENDIF
1447
1448	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
1449	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
1450	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
1451	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
1452	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
1453	ENDIF
1454	!
1455	ENDIF
1456	!
1457	END SUBROUTINE bgc3d_asm_inc
1458
1459	!!===========================================================================
1460
1461	END MODULE asmbgc

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

Download in other formats: