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 @ 9431

Last change on this file since 9431 was 9431, checked in by dford, 6 years ago
Allow assimilation of non-log chlorophyll, and prepare for assimilation of PFTs.
File size: 102.3 KB

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

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

Download in other formats: