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

Last change on this file since 9333 was 9333, checked in by dford, 6 years ago
Some minor tidying.
File size: 78.8 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	WHERE( phyto2d_balinc(:,:,:,:) > 0.0_wp .OR. &
867	& trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,:) * zincwgt > 0.0_wp )
868	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
869	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
870	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
871	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
872	END WHERE
873	#elif defined key_hadocc
874	WHERE( phyto2d_balinc(:,:,:,:) > 0.0_wp .OR. &
875	& trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,:) * zincwgt > 0.0_wp )
876	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
877	& phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
878	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
879	& phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
880	END WHERE
881	#endif
882
883	! Do not deallocate arrays - needed by asm_bgc_bal_wri
884	! which is called at end of model run
885	ENDIF
886
887	ELSEIF ( ll_asmdin ) THEN
888
889	!--------------------------------------------------------------------
890	! Direct Initialization
891	!--------------------------------------------------------------------
892
893	IF ( kt == nitdin_r ) THEN
894
895	neuler = 0 ! Force Euler forward step
896
897	! Initialize the now fields with the background + increment
898	! Background currently is what the model is initialised with
899	CALL ctl_warn( ' Doing direct initialisation with phyto2d assimilation', &
900	& ' Background state is taken from model rather than background file' )
901	#if defined key_medusa && defined key_foam_medusa
902	WHERE( phyto2d_balinc(:,:,:,:) > 0.0_wp .OR. &
903	& trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,:) * zincwgt > 0.0_wp )
904	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
905	& phyto2d_balinc(:,:,:,jp_msa0:jp_msa1)
906	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
907	END WHERE
908	#elif defined key_hadocc
909	WHERE( phyto2d_balinc(:,:,:,:) > 0.0_wp .OR. &
910	& trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,:) * zincwgt > 0.0_wp )
911	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
912	& phyto2d_balinc(:,:,:,jp_had0:jp_had1)
913	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
914	END WHERE
915	#endif
916
917	! Do not deallocate arrays - needed by asm_bgc_bal_wri
918	! which is called at end of model run
919	ENDIF
920	!
921	ENDIF
922	!
923	END SUBROUTINE phyto2d_asm_inc
924
925	!!===========================================================================
926	!!===========================================================================
927	!!===========================================================================
928
929	SUBROUTINE phyto3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
930	!!------------------------------------------------------------------------
931	!! * ROUTINE phyto3d_asm_inc *
932	!!
933	!! ** Purpose : Apply the profile chlorophyll assimilation increments.
934	!!
935	!! ** Method : Calculate increments to state variables.
936	!! Direct initialization or Incremental Analysis Updating.
937	!!
938	!! ** Action :
939	!!------------------------------------------------------------------------
940	INTEGER, INTENT(IN) :: kt ! Current time step
941	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
942	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
943	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
944	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
945	!
946	INTEGER :: ji, jj, jk ! Loop counters
947	INTEGER :: it ! Index
948	REAL(wp) :: zincwgt ! IAU weight for time step
949	REAL(wp) :: zfrac_chn ! Fraction of jpchn
950	REAL(wp) :: zfrac_chd ! Fraction of jpchd
951	REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc ! Chlorophyll increments
952	REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl ! Chlorophyll background
953	!!------------------------------------------------------------------------
954
955	IF ( kt <= nit000 ) THEN
956
957	IF ( ln_plchltotinc ) THEN
958	! Convert log10(chlorophyll) increment back to a chlorophyll increment
959	! In order to transform logchl incs to chl incs, need to account for model
960	! background, cannot simply do 10^logchl_bkginc. Need to:
961	! 1) Add logchl inc to log10(background) to get log10(analysis)
962	! 2) Take 10^log10(analysis) to get analysis
963	! 3) Subtract background from analysis to get chl incs
964	! If rn_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
965	#if defined key_medusa && defined key_foam_medusa
966	bkg_chl(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd)
967	#elif defined key_hadocc
968	bkg_chl(:,:,:) = chl_bkg(:,:,:)
969	#endif
970	DO jk = 1, jpk
971	DO jj = 1, jpj
972	DO ji = 1, jpi
973	IF ( bkg_chl(ji,jj,jk) > 0.0 ) THEN
974	chl_inc(ji,jj,jk) = 10**( LOG10( bkg_chl(ji,jj,jk) ) + plchltot_bkginc(ji,jj,jk) ) - bkg_chl(ji,jj,jk)
975	IF ( rn_maxchlinc > 0.0 ) THEN
976	chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( chl_inc(ji,jj,jk), rn_maxchlinc ) )
977	ENDIF
978	ELSE
979	chl_inc(ji,jj,jk) = 0.0
980	ENDIF
981	END DO
982	END DO
983	END DO
984	ELSE IF ( ln_pchltotinc ) THEN
985	DO jk = 1, jpk
986	DO jj = 1, jpj
987	DO ji = 1, jpi
988	IF ( rn_maxchlinc > 0.0 ) THEN
989	chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( pchltot_bkginc(ji,jj,jk), rn_maxchlinc ) )
990	ELSE
991	chl_inc(ji,jj,jk) = pchltot_bkginc(ji,jj,jk)
992	ENDIF
993	END DO
994	END DO
995	END DO
996	ENDIF
997
998	#if defined key_medusa && defined key_foam_medusa
999	! Loop over each grid point partioning the increments based on existing ratios
1000	DO jk = 1, jpk
1001	DO jj = 1, jpj
1002	DO ji = 1, jpi
1003	IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN
1004	zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd))
1005	zfrac_chd = 1.0 - zfrac_chn
1006	phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn
1007	phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd
1008	ENDIF
1009	END DO
1010	END DO
1011	END DO
1012	#elif defined key_hadocc
1013	phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:)
1014	#else
1015	CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', &
1016	& 'but not defined a biogeochemical model' )
1017	#endif
1018
1019	ENDIF
1020
1021	IF ( ll_asmiau ) THEN
1022
1023	!--------------------------------------------------------------------
1024	! Incremental Analysis Updating
1025	!--------------------------------------------------------------------
1026
1027	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1028
1029	it = kt - nit000 + 1
1030	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1031	! note this is not a tendency so should not be divided by rdt
1032
1033	IF(lwp) THEN
1034	WRITE(numout,*)
1035	WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', &
1036	& kt,' with IAU weight = ', pwgtiau(it)
1037	WRITE(numout,*) '~~~~~~~~~~~~'
1038	ENDIF
1039
1040	! Update the biogeochemical variables
1041	! Add directly to trn and trb, rather than to tra, because tra gets
1042	! reset to zero at the start of trc_stp, called after this routine
1043	#if defined key_medusa && defined key_foam_medusa
1044	WHERE( phyto3d_balinc(:,:,:,:) > 0.0_wp .OR. &
1045	& trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,:) * zincwgt > 0.0_wp )
1046	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1047	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1048	trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1049	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1050	END WHERE
1051	#elif defined key_hadocc
1052	WHERE( phyto3d_balinc(:,:,:,:) > 0.0_wp .OR. &
1053	& trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,:) * zincwgt > 0.0_wp )
1054	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1055	& phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1056	trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1057	& phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1058	END WHERE
1059	#endif
1060
1061	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1062	! which is called at end of model run
1063	ENDIF
1064
1065	ELSEIF ( ll_asmdin ) THEN
1066
1067	!--------------------------------------------------------------------
1068	! Direct Initialization
1069	!--------------------------------------------------------------------
1070
1071	IF ( kt == nitdin_r ) THEN
1072
1073	neuler = 0 ! Force Euler forward step
1074
1075	! Initialize the now fields with the background + increment
1076	! Background currently is what the model is initialised with
1077	CALL ctl_warn( ' Doing direct initialisation with phyto3d assimilation', &
1078	& ' Background state is taken from model rather than background file' )
1079	#if defined key_medusa && defined key_foam_medusa
1080	WHERE( phyto3d_balinc(:,:,:,:) > 0.0_wp .OR. &
1081	& trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,:) * zincwgt > 0.0_wp )
1082	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1083	& phyto3d_balinc(:,:,:,jp_msa0:jp_msa1)
1084	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1085	END WHERE
1086	#elif defined key_hadocc
1087	WHERE( phyto3d_balinc(:,:,:,:) > 0.0_wp .OR. &
1088	& trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,:) * zincwgt > 0.0_wp )
1089	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1090	& phyto3d_balinc(:,:,:,jp_had0:jp_had1)
1091	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1092	END WHERE
1093	#endif
1094
1095	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1096	! which is called at end of model run
1097	ENDIF
1098	!
1099	ENDIF
1100	!
1101	END SUBROUTINE phyto3d_asm_inc
1102
1103	!!===========================================================================
1104	!!===========================================================================
1105	!!===========================================================================
1106
1107	SUBROUTINE pco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
1108	& ll_trainc, pt_bkginc, ps_bkginc )
1109	!!------------------------------------------------------------------------
1110	!! * ROUTINE pco2_asm_inc *
1111	!!
1112	!! ** Purpose : Apply the pco2/fco2 assimilation increments.
1113	!!
1114	!! ** Method : Calculate increments to state variables using carbon
1115	!! balancing.
1116	!! Direct initialization or Incremental Analysis Updating.
1117	!!
1118	!! ** Action :
1119	!!------------------------------------------------------------------------
1120	INTEGER, INTENT(IN) :: kt ! Current time step
1121	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1122	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1123	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1124	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1125	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
1126	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
1127	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
1128	!
1129	INTEGER :: ji, jj, jk ! Loop counters
1130	INTEGER :: it ! Index
1131	INTEGER :: jkmax ! Index of mixed layer depth
1132	REAL(wp) :: zincwgt ! IAU weight for current time step
1133	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth for balancing
1134	REAL(wp), DIMENSION(jpi,jpj) :: pco2_bkginc_temp ! For fCO2 conversion
1135	REAL(wp), DIMENSION(jpi,jpj) :: dic_bkg_temp ! DIC background
1136	REAL(wp), DIMENSION(jpi,jpj) :: alk_bkg_temp ! Alkalinity background
1137	REAL(wp), DIMENSION(jpi,jpj) :: tem_bkg_temp ! Temperature background
1138	REAL(wp), DIMENSION(jpi,jpj) :: sal_bkg_temp ! Salinity background
1139	REAL(wp), DIMENSION(1) :: patm ! Atmospheric pressure
1140	REAL(wp), DIMENSION(1) :: phyd ! Hydrostatic pressure
1141
1142	! Coefficients for fCO2 to pCO2 conversion
1143	REAL(wp) :: zcoef_fco2_1 = -1636.75
1144	REAL(wp) :: zcoef_fco2_2 = 12.0408
1145	REAL(wp) :: zcoef_fco2_3 = 0.0327957
1146	REAL(wp) :: zcoef_fco2_4 = 0.0000316528
1147	REAL(wp) :: zcoef_fco2_5 = 2.0
1148	REAL(wp) :: zcoef_fco2_6 = 57.7
1149	REAL(wp) :: zcoef_fco2_7 = 0.118
1150	REAL(wp) :: zcoef_fco2_8 = 82.0578
1151	!!------------------------------------------------------------------------
1152
1153	IF ( kt <= nit000 ) THEN
1154
1155	!----------------------------------------------------------------------
1156	! DIC and alkalinity balancing
1157	!----------------------------------------------------------------------
1158
1159	IF ( ln_sfco2inc ) THEN
1160	#if defined key_medusa && defined key_roam
1161	! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
1162	patm(1) = 1.0
1163	phyd(1) = 0.0
1164	DO jj = 1, jpj
1165	DO ji = 1, jpi
1166	CALL f2pCO2( sfco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
1167	END DO
1168	END DO
1169	#else
1170	! If assimilating fCO2, then convert to pCO2 using temperature
1171	! See flux_gas.F90 within HadOCC for details of calculation
1172	pco2_bkginc_temp(:,:) = sfco2_bkginc(:,:) / &
1173	& EXP((zcoef_fco2_1 + &
1174	& zcoef_fco2_2 * (tsn(:,:,1,1)+rt0) - &
1175	& zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
1176	& zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0) + &
1177	& zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0))) / &
1178	& (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
1179	#endif
1180	ELSE
1181	pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
1182	ENDIF
1183
1184	! Account for physics assimilation if required
1185	IF ( ll_trainc ) THEN
1186	tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
1187	sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
1188	ELSE
1189	tem_bkg_temp(:,:) = tsn(:,:,1,1)
1190	sal_bkg_temp(:,:) = tsn(:,:,1,2)
1191	ENDIF
1192
1193	#if defined key_medusa
1194	! Account for logchl balancing if required
1195	IF ( ln_slchltotinc .AND. ln_phytobal ) THEN
1196	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
1197	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk)
1198	ELSE
1199	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
1200	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
1201	ENDIF
1202
1203	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1204	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
1205	& pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
1206
1207	#elif defined key_hadocc
1208	! Account for slchltot balancing if required
1209	IF ( ln_slchltotinc .AND. ln_phytobal ) THEN
1210	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
1211	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk)
1212	ELSE
1213	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
1214	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
1215	ENDIF
1216
1217	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1218	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
1219	& pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
1220
1221	#else
1222	CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
1223	& 'but not defined a biogeochemical model' )
1224	#endif
1225
1226	! Select mixed layer
1227	IF ( ll_asmdin ) THEN
1228	#if defined key_hadocc \|\| (defined key_medusa && defined key_foam_medusa)
1229	CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
1230	& ' Mixed layer depth taken to be background maximum mld_max_bkg' )
1231	zmld(:,:) = mld_max_bkg(:,:)
1232	#else
1233	CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
1234	#endif
1235	ELSE
1236	SELECT CASE( mld_choice_bgc )
1237	CASE ( 1 ) ! Turbocline/mixing depth [W points]
1238	zmld(:,:) = hmld(:,:)
1239	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
1240	zmld(:,:) = hmlp(:,:)
1241	CASE ( 3 ) ! Kara MLD [Interpolated]
1242	#if defined key_karaml
1243	IF ( ln_kara ) THEN
1244	zmld(:,:) = hmld_kara(:,:)
1245	ELSE
1246	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1247	& ' but ln_kara=.false.' )
1248	ENDIF
1249	#else
1250	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1251	& ' but is not defined' )
1252	#endif
1253	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
1254	!zmld(:,:) = hmld_tref(:,:)
1255	CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
1256	& ' but is not available in this version' )
1257	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
1258	zmld(:,:) = hmlpt(:,:)
1259	END SELECT
1260	ENDIF
1261
1262	! Propagate through mixed layer
1263	DO jj = 1, jpj
1264	DO ji = 1, jpi
1265	!
1266	jkmax = jpk-1
1267	DO jk = jpk-1, 1, -1
1268	IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. &
1269	& ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
1270	zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
1271	jkmax = jk
1272	ENDIF
1273	END DO
1274	!
1275	#if defined key_top
1276	DO jk = 2, jkmax
1277	pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
1278	END DO
1279	#endif
1280	!
1281	END DO
1282	END DO
1283
1284	ENDIF
1285
1286	IF ( ll_asmiau ) THEN
1287
1288	!--------------------------------------------------------------------
1289	! Incremental Analysis Updating
1290	!--------------------------------------------------------------------
1291
1292	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1293
1294	it = kt - nit000 + 1
1295	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1296	! note this is not a tendency so should not be divided by rdt
1297
1298	IF(lwp) THEN
1299	WRITE(numout,*)
1300	IF ( ln_spco2inc ) THEN
1301	WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
1302	& kt,' with IAU weight = ', pwgtiau(it)
1303	ELSE IF ( ln_sfco2inc ) THEN
1304	WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
1305	& kt,' with IAU weight = ', pwgtiau(it)
1306	ENDIF
1307	WRITE(numout,*) '~~~~~~~~~~~~'
1308	ENDIF
1309
1310	! Update the biogeochemical variables
1311	! Add directly to trn and trb, rather than to tra, because tra gets
1312	! reset to zero at the start of trc_stp, called after this routine
1313	#if defined key_medusa && defined key_foam_medusa
1314	DO jk = 1, jpkm1
1315	trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
1316	& pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
1317	trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
1318	& pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
1319	END DO
1320	#elif defined key_hadocc
1321	DO jk = 1, jpkm1
1322	trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
1323	& pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
1324	trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
1325	& pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
1326	END DO
1327	#endif
1328
1329	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1330	! which is called at end of model run
1331
1332	ENDIF
1333
1334	ELSEIF ( ll_asmdin ) THEN
1335
1336	!--------------------------------------------------------------------
1337	! Direct Initialization
1338	!--------------------------------------------------------------------
1339
1340	IF ( kt == nitdin_r ) THEN
1341
1342	neuler = 0 ! Force Euler forward step
1343
1344	#if defined key_medusa && defined key_foam_medusa
1345	! Initialize the now fields with the background + increment
1346	! Background currently is what the model is initialised with
1347	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', &
1348	& ' Background state is taken from model rather than background file' )
1349	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1350	& pco2_balinc(:,:,:,jp_msa0:jp_msa1)
1351	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1352	#elif defined key_hadocc
1353	! Initialize the now fields with the background + increment
1354	! Background currently is what the model is initialised with
1355	CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', &
1356	& ' Background state is taken from model rather than background file' )
1357	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1358	& pco2_balinc(:,:,:,jp_had0:jp_had1)
1359	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1360	#endif
1361
1362	! Do not deallocate arrays - needed by asm_bgc_bal_wri
1363	! which is called at end of model run
1364	ENDIF
1365	!
1366	ENDIF
1367	!
1368	END SUBROUTINE pco2_asm_inc
1369
1370	!!===========================================================================
1371	!!===========================================================================
1372	!!===========================================================================
1373
1374	SUBROUTINE ph_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
1375	& ll_trainc, pt_bkginc, ps_bkginc )
1376	!!------------------------------------------------------------------------
1377	!! * ROUTINE ph_asm_inc *
1378	!!
1379	!! ** Purpose : Apply the pH assimilation increments.
1380	!!
1381	!! ** Method : Calculate increments to state variables using carbon
1382	!! balancing.
1383	!! Direct initialization or Incremental Analysis Updating.
1384	!!
1385	!! ** Action :
1386	!!------------------------------------------------------------------------
1387	INTEGER, INTENT(IN) :: kt ! Current time step
1388	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1389	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1390	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1391	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1392	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
1393	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
1394	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
1395	!!------------------------------------------------------------------------
1396
1397	CALL ctl_stop( ' pH balancing not yet implemented' )
1398
1399	! See solve_at_general line 281 of mocsy_phsolvers.F90
1400	!
1401	! Or, call to mocsy_interface line 158 of carb_chem.F90
1402	! Input variables (rest are output) are:
1403	! ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj), &
1404	! zdic(ji,jj),zsil(ji,jj),zpho(ji,jj), &
1405	! f_pp0(ji,jj),fsdept(ji,jj,jk), &
1406	! gphit(ji,jj),f_kw660(ji,jj), &
1407	! f_xco2a(ji,jj),1
1408	!
1409	! ztmp = tsn(:,:,:,jp_tem)
1410	! zsal = tsn(:,:,:,jp_sal)
1411	! zalk = trn(:,:,:,jpalk)
1412	! zdic = trn(:,:,:,jpdic)
1413	! zsil = trn(:,:,:,jpsil)
1414	! zpho = trn(:,:,:,jpdin) / 16.0
1415	! f_pp0 = 1.0
1416	! fsdept = fsdept(:,:,:)
1417	! gphit = gphit(:,:)
1418	! f_kw660 = 1.0
1419	! f_xco2a = f_xco2a(:,:)
1420
1421	!
1422	END SUBROUTINE ph_asm_inc
1423
1424	!!===========================================================================
1425	!!===========================================================================
1426	!!===========================================================================
1427
1428	SUBROUTINE bgc3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1429	!!----------------------------------------------------------------------
1430	!! * ROUTINE dyn_asm_inc *
1431	!!
1432	!! ** Purpose : Apply generic 3D biogeochemistry assimilation increments.
1433	!!
1434	!! ** Method : Direct initialization or Incremental Analysis Updating.
1435	!!
1436	!! ** Action :
1437	!!----------------------------------------------------------------------
1438	INTEGER, INTENT(IN) :: kt ! Current time step
1439	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1440	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1441	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
1442	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
1443	!
1444	INTEGER :: jk ! Loop counter
1445	INTEGER :: it ! Index
1446	REAL(wp) :: zincwgt ! IAU weight for current time step
1447	REAL(wp) :: zincper ! IAU interval in seconds
1448	!!----------------------------------------------------------------------
1449
1450	IF ( ll_asmiau ) THEN
1451
1452	!--------------------------------------------------------------------
1453	! Incremental Analysis Updating
1454	!--------------------------------------------------------------------
1455
1456	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1457
1458	it = kt - nit000 + 1
1459	zincwgt = pwgtiau(it) ! IAU weight for the current time step
1460	! note this is not a tendency so should not be divided by rdt
1461
1462	IF(lwp) THEN
1463	WRITE(numout,*)
1464	WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', &
1465	& kt,' with IAU weight = ', pwgtiau(it)
1466	WRITE(numout,*) '~~~~~~~~~~~~'
1467	ENDIF
1468
1469	! Update the 3D BGC variables
1470	! Add directly to trn and trb, rather than to tra, because tra gets
1471	! reset to zero at the start of trc_stp, called after this routine
1472	! Don't apply increments if they'll take concentrations negative
1473
1474	IF ( ln_pno3inc ) THEN
1475	#if defined key_hadocc
1476	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
1477	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
1478	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
1479	trb(:,:,:,jp_had_nut) = trb(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
1480	END WHERE
1481	#elif defined key_medusa && defined key_foam_medusa
1482	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
1483	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
1484	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
1485	trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
1486	END WHERE
1487	#else
1488	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1489	#endif
1490	ENDIF
1491
1492	IF ( ln_psi4inc ) THEN
1493	#if defined key_medusa && defined key_foam_medusa
1494	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
1495	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
1496	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
1497	trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
1498	END WHERE
1499	#else
1500	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1501	#endif
1502	ENDIF
1503
1504	IF ( ln_pdicinc ) THEN
1505	#if defined key_hadocc
1506	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
1507	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
1508	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
1509	trb(:,:,:,jp_had_dic) = trb(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
1510	END WHERE
1511	#elif defined key_medusa && defined key_foam_medusa
1512	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
1513	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
1514	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
1515	trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
1516	END WHERE
1517	#else
1518	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1519	#endif
1520	ENDIF
1521
1522	IF ( ln_palkinc ) THEN
1523	#if defined key_hadocc
1524	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
1525	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
1526	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
1527	trb(:,:,:,jp_had_alk) = trb(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
1528	END WHERE
1529	#elif defined key_medusa && defined key_foam_medusa
1530	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
1531	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
1532	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
1533	trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
1534	END WHERE
1535	#else
1536	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1537	#endif
1538	ENDIF
1539
1540	IF ( ln_po2inc ) THEN
1541	#if defined key_medusa && defined key_foam_medusa
1542	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
1543	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
1544	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
1545	trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
1546	END WHERE
1547	#else
1548	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1549	#endif
1550	ENDIF
1551
1552	IF ( kt == nitiaufin_r ) THEN
1553	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
1554	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
1555	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
1556	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
1557	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
1558	ENDIF
1559
1560	ENDIF
1561
1562	ELSEIF ( ll_asmdin ) THEN
1563
1564	!--------------------------------------------------------------------
1565	! Direct Initialization
1566	!--------------------------------------------------------------------
1567
1568	IF ( kt == nitdin_r ) THEN
1569
1570	neuler = 0 ! Force Euler forward step
1571
1572	! Initialize the now fields with the background + increment
1573	! Background currently is what the model is initialised with
1574	#if defined key_hadocc
1575	CALL ctl_warn( ' Doing direct initialisation of HadOCC with 3D BGC assimilation', &
1576	& ' Background state is taken from model rather than background file' )
1577	#elif defined key_medusa && defined key_foam_medusa
1578	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with 3D BGC assimilation', &
1579	& ' Background state is taken from model rather than background file' )
1580	#endif
1581
1582	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1583	& pco2_balinc(:,:,:,jp_msa0:jp_msa1)
1584	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1585
1586	IF ( ln_pno3inc ) THEN
1587	#if defined key_hadocc
1588	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
1589	& trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) > 0.0_wp )
1590	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:)
1591	trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut)
1592	END WHERE
1593	#elif defined key_medusa && defined key_foam_medusa
1594	WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
1595	& trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) > 0.0_wp )
1596	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:)
1597	trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
1598	END WHERE
1599	#else
1600	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1601	#endif
1602	ENDIF
1603
1604	IF ( ln_psi4inc ) THEN
1605	#if defined key_medusa && defined key_foam_medusa
1606	WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
1607	& trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) > 0.0_wp )
1608	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:)
1609	trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
1610	END WHERE
1611	#else
1612	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1613	#endif
1614	ENDIF
1615
1616	IF ( ln_pdicinc ) THEN
1617	#if defined key_hadocc
1618	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
1619	& trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) > 0.0_wp )
1620	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:)
1621	trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic)
1622	END WHERE
1623	#elif defined key_medusa && defined key_foam_medusa
1624	WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
1625	& trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) > 0.0_wp )
1626	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:)
1627	trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
1628	END WHERE
1629	#else
1630	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1631	#endif
1632	ENDIF
1633
1634	IF ( ln_palkinc ) THEN
1635	#if defined key_hadocc
1636	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
1637	& trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) > 0.0_wp )
1638	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:)
1639	trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk)
1640	END WHERE
1641	#elif defined key_medusa && defined key_foam_medusa
1642	WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
1643	& trn(:,:,:,jpalk) + palk_bkginc(:,:,:) > 0.0_wp )
1644	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:)
1645	trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
1646	END WHERE
1647	#else
1648	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1649	#endif
1650	ENDIF
1651
1652	IF ( ln_po2inc ) THEN
1653	#if defined key_medusa && defined key_foam_medusa
1654	WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
1655	& trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) > 0.0_wp )
1656	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:)
1657	trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
1658	END WHERE
1659	#else
1660	CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1661	#endif
1662	ENDIF
1663
1664	IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
1665	IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
1666	IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
1667	IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
1668	IF ( ln_po2inc ) DEALLOCATE( po2_bkginc )
1669	ENDIF
1670	!
1671	ENDIF
1672	!
1673	END SUBROUTINE bgc3d_asm_inc
1674
1675	!!===========================================================================
1676
1677	END MODULE asmbgc

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

Download in other formats: