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

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

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

Download in other formats: