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

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

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

Download in other formats: