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

Last change on this file since 9297 was 9297, checked in by dford, 6 years ago
Cosmetic changes.
File size: 45.3 KB

Line
1	MODULE asmbgc
2	!!===========================================================================
3	!! * MODULE asmbgc *
4	!! Routines and declarations for biogeochemical assimilation
5	!!===========================================================================
6	!! History : 3.6 ! 2018-02 (D. Ford) Adapted from existing modules
7	!!---------------------------------------------------------------------------
8	!! 'key_asminc' : assimilation increment interface
9	!! 'key_top' : passive tracers
10	!! 'key_hadocc' : HadOCC model
11	!! 'key_medusa' : MEDUSA model
12	!! 'key_foam_medusa' : MEDUSA extras for FOAM OBS and ASM
13	!! 'key_roam' : MEDUSA extras for carbonate cycle
14	!! 'key_karaml' : Kara mixed layer depth
15	!!---------------------------------------------------------------------------
16	!! asm_bgc_check_options : check bgc assimilation options
17	!! asm_bgc_init_incs : initialise bgc increments
18	!! asm_bgc_init_bkg : initialise bgc background
19	!! asm_bgc_bal_wri : write out bgc balancing increments
20	!! asm_bgc_bkg_wri : write out bgc background
21	!! slchltot_asm_inc : apply the slchltot increment
22	!! spco2_asm_inc : apply the pCO2/fCO2 increment
23	!!---------------------------------------------------------------------------
24	USE par_kind, ONLY: & ! kind parameters
25	& wp
26	USE par_oce, ONLY: & ! domain array sizes
27	& jpi, &
28	& jpj, &
29	& jpk
30	USE iom ! i/o
31	USE oce, ONLY: & ! active tracer variables
32	& tsn
33	USE zdfmxl, ONLY : & ! mixed layer depth
34	#if defined key_karaml
35	& hmld_kara, &
36	& ln_kara, &
37	#endif
38	& hmld, &
39	& hmlp, &
40	& hmlpt
41	USE asmpar, ONLY: & ! assimilation parameters
42	& c_asmbkg, &
43	& c_asmbal, &
44	& nitbkg_r, &
45	& nitdin_r, &
46	& nitiaustr_r, &
47	& nitiaufin_r
48	#if defined key_top
49	USE trc, ONLY: & ! passive tracer variables
50	& trn, &
51	& trb
52	USE par_trc, ONLY: & ! passive tracer parameters
53	& jptra
54	#endif
55	#if defined key_medusa && defined key_foam_medusa
56	USE asmlogchlbal_medusa, ONLY: & ! logchl balancing for MEDUSA
57	& asm_logchl_bal_medusa
58	USE asmpco2bal, ONLY: & ! pCO2 balancing for MEDUSA
59	& asm_pco2_bal
60	USE par_medusa ! MEDUSA parameters
61	USE mocsy_mainmod, ONLY: & ! MEDUSA carbonate system
62	& f2pCO2
63	USE sms_medusa, ONLY: & ! MEDUSA diagnostic variables
64	& pgrow_avg, &
65	& ploss_avg, &
66	& phyt_avg, &
67	& mld_max
68	#elif defined key_hadocc
69	USE asmlogchlbal_hadocc, ONLY: & ! logchl balancing for HadOCC
70	& asm_logchl_bal_hadocc
71	USE asmpco2bal, ONLY: & ! pCO2 balancing for HadOCC
72	& asm_pco2_bal
73	USE par_hadocc ! HadOCC parameters
74	USE trc, ONLY: & ! HadOCC diagnostic variables
75	& pgrow_avg, &
76	& ploss_avg, &
77	& phyt_avg, &
78	& mld_max, &
79	& HADOCC_CHL
80	USE had_bgc_const, ONLY: & ! HadOCC C:Chl ratio
81	& cchl_p
82	#endif
83
84	IMPLICIT NONE
85	PRIVATE
86
87	PUBLIC asm_bgc_check_options ! called by asm_inc_init in asminc.F90
88	PUBLIC asm_bgc_init_incs ! called by asm_inc_init in asminc.F90
89	PUBLIC asm_bgc_init_bkg ! called by asm_inc_init in asminc.F90
90	PUBLIC asm_bgc_bal_wri ! called by nemo_gcm in nemogcm.F90
91	PUBLIC asm_bgc_bkg_wri ! called by asm_bkg_wri in asmbkg.F90
92	PUBLIC slchltot_asm_inc ! called by bgc_asm_inc in asminc.F90
93	PUBLIC spco2_asm_inc ! called by bgc_asm_inc in asminc.F90
94
95	LOGICAL, PUBLIC :: ln_balwri = .FALSE. !: No output of balancing incs
96	LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE. !: No slchltot increment
97	LOGICAL, PUBLIC :: ln_slchltotbal = .FALSE. !: No slchltot balancing
98	LOGICAL, PUBLIC :: ln_spco2inc = .FALSE. !: No pCO2 increment
99	LOGICAL, PUBLIC :: ln_sfco2inc = .FALSE. !: No fCO2 increment
100
101	INTEGER, PUBLIC :: mld_choice_bgc = 1 !: choice of mld for bgc assimilation
102	! 1) hmld - Turbocline/mixing depth
103	! [W points]
104	! 2) hmlp - Density criterion
105	! (0.01 kg/m^3 change from 10m)
106	! [W points]
107	! 3) hmld_kara - Kara MLD
108	! [Interpolated]
109	! 4) hmld_tref - Temperature criterion
110	! (0.2 K change from surface)
111	! [T points]
112	! 5) hmlpt - Density criterion
113	! (0.01 kg/m^3 change from 10m)
114	! [T points]
115
116	REAL(wp), PUBLIC :: rn_maxchlinc = -999.0 !: maximum absolute non-log chlorophyll inc
117	! <= 0 no maximum applied (switch turned off)
118	! > 0 absolute chl inc capped at this value
119
120	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: slchltot_bkginc ! slchltot inc
121	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: spco2_bkginc ! sp(/f)co2 inc
122	#if defined key_top
123	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: slchltot_balinc ! Balancing incs from slchltot asm
124	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: spco2_balinc ! Balancing incs from sp(/f)co2 asm
125	#endif
126
127	#if defined key_hadocc \|\| (defined key_medusa && defined key_foam_medusa)
128	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pgrow_avg_bkg ! Background phyto growth
129	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ploss_avg_bkg ! Background phyto loss
130	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: phyt_avg_bkg ! Background phyto conc
131	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mld_max_bkg ! Background max MLD
132	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg ! Background tracer state
133	#endif
134	#if defined key_hadocc
135	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: chl_bkg ! Background surface chl
136	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: cchl_p_bkg ! Background surface c:chl
137	#endif
138
139	CONTAINS
140
141	SUBROUTINE asm_bgc_check_options
142	!!------------------------------------------------------------------------
143	!! * ROUTINE asm_bgc_check_options *
144	!!
145	!! ** Purpose : check validity of biogeochemical assimilation options
146	!!
147	!! ** Method : check settings of logicals
148	!!
149	!! ** Action : call ctl_stop/ctl_warn if required
150	!!
151	!! References : asm_inc_init
152	!!------------------------------------------------------------------------
153
154	IF ( ( ln_balwri ).AND.( .NOT. ln_slchltotinc ).AND. &
155	& ( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ) ) THEN
156	CALL ctl_warn( ' Balancing increments are only calculated for logchl and pCO2/fCO2', &
157	& ' Not assimilating logchl, pCO2 or fCO2, so ln_balwri will be set to .false.')
158	ln_balwri = .FALSE.
159	ENDIF
160
161	IF ( ( ln_slchltotbal ).AND.( .NOT. ln_slchltotinc ) ) THEN
162	CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', &
163	& ' Not assimilating logchl, so ln_slchltotbal will be set to .false.')
164	ln_slchltotbal = .FALSE.
165	ENDIF
166
167	IF ( ( ln_spco2inc ).AND.( ln_sfco2inc ) ) THEN
168	CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
169	ENDIF
170
171	END SUBROUTINE asm_bgc_check_options
172
173	!!===========================================================================
174	!!===========================================================================
175	!!===========================================================================
176
177	SUBROUTINE asm_bgc_init_incs( knum )
178	!!------------------------------------------------------------------------
179	!! * ROUTINE asm_bgc_init_incs *
180	!!
181	!! ** Purpose : initialise bgc increments
182	!!
183	!! ** Method : allocate arrays and read increments from file
184	!!
185	!! ** Action : allocate arrays and read increments from file
186	!!
187	!! References : asm_inc_init
188	!!------------------------------------------------------------------------
189	!!
190	INTEGER, INTENT(in ) :: knum ! i/o unit of increments file
191	!!
192	!!------------------------------------------------------------------------
193
194	IF ( ln_slchltotinc ) THEN
195	ALLOCATE( slchltot_bkginc(jpi,jpj) )
196	slchltot_bkginc(:,:) = 0.0
197	#if defined key_top
198	ALLOCATE( slchltot_balinc(jpi,jpj,jpk,jptra) )
199	slchltot_balinc(:,:,:,:) = 0.0
200	#endif
201	ENDIF
202	IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
203	ALLOCATE( spco2_bkginc(jpi,jpj) )
204	spco2_bkginc(:,:) = 0.0
205	#if defined key_top
206	ALLOCATE( spco2_balinc(jpi,jpj,jpk,jptra) )
207	spco2_balinc(:,:,:,:) = 0.0
208	#endif
209	ENDIF
210
211	IF ( ln_slchltotinc ) THEN
212	CALL iom_get( knum, jpdom_autoglo, 'bckinslchltot', slchltot_bkginc(:,:), 1 )
213	! Apply the masks
214	slchltot_bkginc(:,:) = slchltot_bkginc(:,:) * tmask(:,:,1)
215	! Set missing increments to 0.0 rather than 1e+20
216	! to allow for differences in masks
217	WHERE( ABS( slchltot_bkginc(:,:) ) > 1.0e+10 ) slchltot_bkginc(:,:) = 0.0
218	ENDIF
219
220	IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
221	IF ( ln_spco2inc ) THEN
222	CALL iom_get( knum, jpdom_autoglo, 'bckinspco2', spco2_bkginc(:,:), 1 )
223	ELSE IF ( ln_sfco2inc ) THEN
224	CALL iom_get( knum, jpdom_autoglo, 'bckinsfco2', spco2_bkginc(:,:), 1 )
225	ENDIF
226	! Apply the masks
227	spco2_bkginc(:,:) = spco2_bkginc(:,:) * tmask(:,:,1)
228	! Set missing increments to 0.0 rather than 1e+20
229	! to allow for differences in masks
230	WHERE( ABS( spco2_bkginc(:,:) ) > 1.0e+10 ) spco2_bkginc(:,:) = 0.0
231	ENDIF
232
233	END SUBROUTINE asm_bgc_init_incs
234
235	!!===========================================================================
236	!!===========================================================================
237	!!===========================================================================
238
239	SUBROUTINE asm_bgc_init_bkg
240	!!------------------------------------------------------------------------
241	!! * ROUTINE asm_bgc_init_bkg *
242	!!
243	!! ** Purpose : initialise bgc background
244	!!
245	!! ** Method : allocate arrays and read background from file
246	!!
247	!! ** Action : allocate arrays and read background from file
248	!!
249	!! References : asm_inc_init
250	!!------------------------------------------------------------------------
251	!!
252	INTEGER :: inum ! i/o unit of background file
253	INTEGER :: jt ! loop counter
254	!!
255	!!------------------------------------------------------------------------
256
257	#if defined key_hadocc \|\| (defined key_medusa && defined key_foam_medusa)
258	IF ( ln_slchltotinc ) THEN
259
260	ALLOCATE( pgrow_avg_bkg(jpi,jpj) )
261	ALLOCATE( ploss_avg_bkg(jpi,jpj) )
262	ALLOCATE( phyt_avg_bkg(jpi,jpj) )
263	ALLOCATE( mld_max_bkg(jpi,jpj) )
264	ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
265	pgrow_avg_bkg(:,:) = 0.0
266	ploss_avg_bkg(:,:) = 0.0
267	phyt_avg_bkg(:,:) = 0.0
268	mld_max_bkg(:,:) = 0.0
269	tracer_bkg(:,:,:,:) = 0.0
270
271	#if defined key_hadocc
272	ALLOCATE( chl_bkg(jpi,jpj) )
273	ALLOCATE( cchl_p_bkg(jpi,jpj) )
274	chl_bkg(:,:) = 0.0
275	cchl_p_bkg(:,:) = 0.0
276	#endif
277
278	!--------------------------------------------------------------------
279	! Read background variables for slchltot assimilation
280	! Some only required if performing balancing
281	!--------------------------------------------------------------------
282
283	CALL iom_open( c_asmbkg, inum )
284
285	#if defined key_hadocc
286	CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl', chl_bkg )
287	CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg )
288	chl_bkg(:,:) = chl_bkg(:,:) * tmask(:,:,1)
289	cchl_p_bkg(:,:) = cchl_p_bkg(:,:) * tmask(:,:,1)
290	#elif defined key_medusa
291	CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) )
292	CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) )
293	#endif
294
295	IF ( ln_slchltotbal ) THEN
296
297	CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg )
298	CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg )
299	CALL iom_get( inum, jpdom_autoglo, 'phyt_avg', phyt_avg_bkg )
300	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
301	pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1)
302	ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1)
303	phyt_avg_bkg(:,:) = phyt_avg_bkg(:,:) * tmask(:,:,1)
304	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
305
306	#if defined key_hadocc
307	CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) )
308	CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) )
309	CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) )
310	CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) )
311	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
312	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
313	#elif defined key_medusa
314	CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) )
315	CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
316	CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
317	CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) )
318	CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) )
319	CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) )
320	CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) )
321	CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) )
322	CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) )
323	CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) )
324	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
325	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
326	CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
327	#endif
328	ELSE IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
329	#if defined key_hadocc
330	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
331	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
332	#elif defined key_medusa
333	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
334	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
335	#endif
336	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
337	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
338	ENDIF
339
340	CALL iom_close( inum )
341
342	DO jt = 1, jptra
343	tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
344	END DO
345
346	ELSE IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
347
348	ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
349	ALLOCATE( mld_max_bkg(jpi,jpj) )
350	tracer_bkg(:,:,:,:) = 0.0
351	mld_max_bkg(:,:) = 0.0
352
353	CALL iom_open( c_asmbkg, inum )
354
355	#if defined key_hadocc
356	CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
357	CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
358	#elif defined key_medusa
359	CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
360	CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
361	#endif
362	CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg )
363
364	CALL iom_close( inum )
365
366	DO jt = 1, jptra
367	tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
368	END DO
369	mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
370
371	ENDIF
372	#endif
373
374	END SUBROUTINE asm_bgc_init_bkg
375
376	!!===========================================================================
377	!!===========================================================================
378	!!===========================================================================
379
380	SUBROUTINE asm_bgc_bal_wri( kt )
381	!!------------------------------------------------------------------------
382	!!
383	!! * ROUTINE asm_bgc_bal_wri *
384	!!
385	!! ** Purpose : Write to file the balancing increments
386	!! calculated for biogeochemistry
387	!!
388	!! ** Method : Write to file the balancing increments
389	!! calculated for biogeochemistry
390	!!
391	!! ** Action :
392	!!
393	!! References :
394	!!
395	!! History :
396	!! ! 2014-08 (D. Ford) Initial version, based on asm_bkg_wri
397	!!------------------------------------------------------------------------
398	!! * Arguments
399	INTEGER, INTENT( IN ) :: kt ! Current time-step
400	!! * Local declarations
401	CHARACTER(LEN=50) :: cl_asmbal ! Filename (with extension)
402	LOGICAL :: llok ! Check if file exists
403	INTEGER :: inum ! File unit number
404	REAL(wp) :: zdate ! Date
405	!!------------------------------------------------------------------------
406
407	! Set things up
408	zdate = REAL( ndastp )
409	WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal )
410
411	! Check if file exists
412	INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok )
413	IF( .NOT. llok ) THEN
414	IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// &
415	& TRIM( c_asmbal ) // ' at timestep = ', kt
416
417	! Define the output file
418	CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
419
420	! Write the information
421	CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate )
422
423	IF ( ln_slchltotinc ) THEN
424	#if defined key_medusa
425	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', slchltot_balinc(:,:,:,jpchn) )
426	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', slchltot_balinc(:,:,:,jpchd) )
427	IF ( ln_slchltotbal ) THEN
428	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phn', slchltot_balinc(:,:,:,jpphn) )
429	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phd', slchltot_balinc(:,:,:,jpphd) )
430	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_pds', slchltot_balinc(:,:,:,jppds) )
431	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zmi', slchltot_balinc(:,:,:,jpzmi) )
432	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zme', slchltot_balinc(:,:,:,jpzme) )
433	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_din', slchltot_balinc(:,:,:,jpdin) )
434	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_sil', slchltot_balinc(:,:,:,jpsil) )
435	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_fer', slchltot_balinc(:,:,:,jpfer) )
436	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', slchltot_balinc(:,:,:,jpdet) )
437	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dtc', slchltot_balinc(:,:,:,jpdtc) )
438	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', slchltot_balinc(:,:,:,jpdic) )
439	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', slchltot_balinc(:,:,:,jpalk) )
440	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_oxy', slchltot_balinc(:,:,:,jpoxy) )
441	ENDIF
442	#elif defined key_hadocc
443	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phy', slchltot_balinc(:,:,:,jp_had_phy) )
444	IF ( ln_slchltotbal ) THEN
445	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_nut', slchltot_balinc(:,:,:,jp_had_nut) )
446	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zoo', slchltot_balinc(:,:,:,jp_had_zoo) )
447	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', slchltot_balinc(:,:,:,jp_had_pdn) )
448	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', slchltot_balinc(:,:,:,jp_had_dic) )
449	CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', slchltot_balinc(:,:,:,jp_had_alk) )
450	ENDIF
451	#endif
452	ENDIF
453
454	IF ( ln_spco2inc ) THEN
455	#if defined key_medusa
456	CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', spco2_balinc(:,:,:,jpdic) )
457	CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', spco2_balinc(:,:,:,jpalk) )
458	#elif defined key_hadocc
459	CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', spco2_balinc(:,:,:,jp_had_dic) )
460	CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', spco2_balinc(:,:,:,jp_had_alk) )
461	#endif
462	ELSE IF ( ln_sfco2inc ) THEN
463	#if defined key_medusa
464	CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', spco2_balinc(:,:,:,jpdic) )
465	CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', spco2_balinc(:,:,:,jpalk) )
466	#elif defined key_hadocc
467	CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', spco2_balinc(:,:,:,jp_had_dic) )
468	CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', spco2_balinc(:,:,:,jp_had_alk) )
469	#endif
470	ENDIF
471
472	CALL iom_close( inum )
473	ELSE
474	CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', &
475	& ' Therefore not writing out balancing increments at this timestep', &
476	& ' Something has probably gone wrong somewhere' )
477	IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// &
478	& TRIM( c_asmbal ) // ' at timestep = ', kt
479	ENDIF
480
481	END SUBROUTINE asm_bgc_bal_wri
482
483	!!===========================================================================
484	!!===========================================================================
485	!!===========================================================================
486
487	SUBROUTINE asm_bgc_bkg_wri( kt, knum )
488	!!------------------------------------------------------------------------
489	!! * ROUTINE asm_bgc_bkg_wri *
490	!!
491	!! ** Purpose : write out bgc background
492	!!
493	!! ** Method : write out bgc background
494	!!
495	!! ** Action : write out bgc background
496	!!
497	!! References : asm_bkg_wri
498	!!------------------------------------------------------------------------
499	!!
500	INTEGER, INTENT(in ) :: kt ! Current time-step
501	INTEGER, INTENT(in ) :: knum ! i/o unit of increments file
502	!!
503	!!------------------------------------------------------------------------
504
505	#if defined key_hadocc
506	CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg' , pgrow_avg )
507	CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg' , ploss_avg )
508	CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg' , phyt_avg )
509	CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max' , mld_max )
510	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_nut' , trn(:,:,:,jp_had_nut) )
511	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_phy' , trn(:,:,:,jp_had_phy) )
512	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_zoo' , trn(:,:,:,jp_had_zoo) )
513	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_pdn' , trn(:,:,:,jp_had_pdn) )
514	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_dic' , trn(:,:,:,jp_had_dic) )
515	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_alk' , trn(:,:,:,jp_had_alk) )
516	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_chl' , HADOCC_CHL(:,:,1) )
517	CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,1) )
518	#elif defined key_medusa && defined key_foam_medusa
519	CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg' , pgrow_avg )
520	CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg' , ploss_avg )
521	CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg' , phyt_avg )
522	CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max' , mld_max )
523	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn' , trn(:,:,:,jpchn) )
524	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chd' , trn(:,:,:,jpchd) )
525	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phn' , trn(:,:,:,jpphn) )
526	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phd' , trn(:,:,:,jpphd) )
527	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_pds' , trn(:,:,:,jppds) )
528	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zmi' , trn(:,:,:,jpzmi) )
529	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zme' , trn(:,:,:,jpzme) )
530	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_din' , trn(:,:,:,jpdin) )
531	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_sil' , trn(:,:,:,jpsil) )
532	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_fer' , trn(:,:,:,jpfer) )
533	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_det' , trn(:,:,:,jpdet) )
534	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dtc' , trn(:,:,:,jpdtc) )
535	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dic' , trn(:,:,:,jpdic) )
536	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_alk' , trn(:,:,:,jpalk) )
537	CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_oxy' , trn(:,:,:,jpoxy) )
538	#endif
539
540	END SUBROUTINE asm_bgc_bkg_wri
541
542	!!===========================================================================
543	!!===========================================================================
544	!!===========================================================================
545
546	SUBROUTINE slchltot_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
547	!!------------------------------------------------------------------------
548	!! * ROUTINE logchl_asm_inc *
549	!!
550	!! ** Purpose : Apply the chlorophyll assimilation increments.
551	!!
552	!! ** Method : Calculate increments to state variables using nitrogen
553	!! balancing.
554	!! Direct initialization or Incremental Analysis Updating.
555	!!
556	!! ** Action :
557	!!------------------------------------------------------------------------
558	INTEGER, INTENT(IN) :: kt ! Current time step
559	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
560	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
561	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
562	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
563	!
564	INTEGER :: jk ! Loop counter
565	INTEGER :: it ! Index
566	REAL(wp) :: zincwgt ! IAU weight for current time step
567	REAL(wp) :: zincper ! IAU interval in seconds
568	!!------------------------------------------------------------------------
569
570	IF ( kt <= nit000 ) THEN
571
572	zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
573
574	#if defined key_medusa && defined key_foam_medusa
575	CALL asm_logchl_bal_medusa( slchltot_bkginc, zincper, mld_choice_bgc, &
576	& rn_maxchlinc, ln_slchltotbal, ll_asmdin, &
577	& pgrow_avg_bkg, ploss_avg_bkg, &
578	& phyt_avg_bkg, mld_max_bkg, &
579	& tracer_bkg, slchltot_balinc )
580	#elif defined key_hadocc
581	CALL asm_logchl_bal_hadocc( slchltot_bkginc, zincper, mld_choice_bgc, &
582	& rn_maxchlinc, ln_slchltotbal, ll_asmdin, &
583	& pgrow_avg_bkg, ploss_avg_bkg, &
584	& phyt_avg_bkg, mld_max_bkg, &
585	& chl_bkg, cchl_p_bkg, &
586	& tracer_bkg, slchltot_balinc )
587	#else
588	CALL ctl_stop( 'Attempting to assimilate slchltot, ', &
589	& 'but not defined a biogeochemical model' )
590	#endif
591
592	ENDIF
593
594	IF ( ll_asmiau ) THEN
595
596	!--------------------------------------------------------------------
597	! Incremental Analysis Updating
598	!--------------------------------------------------------------------
599
600	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
601
602	it = kt - nit000 + 1
603	zincwgt = pwgtiau(it) ! IAU weight for the current time step
604	! note this is not a tendency so should not be divided by rdt
605
606	IF(lwp) THEN
607	WRITE(numout,*)
608	WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', &
609	& kt,' with IAU weight = ', pwgtiau(it)
610	WRITE(numout,*) '~~~~~~~~~~~~'
611	ENDIF
612
613	! Update the biogeochemical variables
614	! Add directly to trn and trb, rather than to tra, because tra gets
615	! reset to zero at the start of trc_stp, called after this routine
616	#if defined key_medusa && defined key_foam_medusa
617	DO jk = 1, jpkm1
618	trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
619	& slchltot_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
620	trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
621	& slchltot_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
622	END DO
623	#elif defined key_hadocc
624	DO jk = 1, jpkm1
625	trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
626	& slchltot_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
627	trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
628	& slchltot_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
629	END DO
630	#endif
631
632	! Do not deallocate arrays - needed by asm_bgc_bal_wri
633	! which is called at end of model run
634	ENDIF
635
636	ELSEIF ( ll_asmdin ) THEN
637
638	!--------------------------------------------------------------------
639	! Direct Initialization
640	!--------------------------------------------------------------------
641
642	IF ( kt == nitdin_r ) THEN
643
644	neuler = 0 ! Force Euler forward step
645
646	#if defined key_medusa && defined key_foam_medusa
647	! Initialize the now fields with the background + increment
648	! Background currently is what the model is initialised with
649	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', &
650	& ' Background state is taken from model rather than background file' )
651	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
652	& slchltot_balinc(:,:,:,jp_msa0:jp_msa1)
653	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
654	#elif defined key_hadocc
655	! Initialize the now fields with the background + increment
656	! Background currently is what the model is initialised with
657	CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', &
658	& ' Background state is taken from model rather than background file' )
659	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
660	& slchltot_balinc(:,:,:,jp_had0:jp_had1)
661	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
662	#endif
663
664	! Do not deallocate arrays - needed by asm_bgc_bal_wri
665	! which is called at end of model run
666	ENDIF
667	!
668	ENDIF
669	!
670	END SUBROUTINE slchltot_asm_inc
671
672	!!===========================================================================
673	!!===========================================================================
674	!!===========================================================================
675
676	SUBROUTINE spco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
677	& ll_trainc, pt_bkginc, ps_bkginc )
678	!!------------------------------------------------------------------------
679	!! * ROUTINE pco2_asm_inc *
680	!!
681	!! ** Purpose : Apply the pco2/fco2 assimilation increments.
682	!!
683	!! ** Method : Calculate increments to state variables using carbon
684	!! balancing.
685	!! Direct initialization or Incremental Analysis Updating.
686	!!
687	!! ** Action :
688	!!------------------------------------------------------------------------
689	INTEGER, INTENT(IN) :: kt ! Current time step
690	LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
691	LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
692	INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau
693	REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights
694	LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments
695	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
696	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
697	!
698	INTEGER :: ji, jj, jk ! Loop counters
699	INTEGER :: it ! Index
700	INTEGER :: jkmax ! Index of mixed layer depth
701	REAL(wp) :: zincwgt ! IAU weight for current time step
702	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth for balancing
703	REAL(wp), DIMENSION(jpi,jpj) :: pco2_bkginc_temp ! For fCO2 conversion
704	REAL(wp), DIMENSION(jpi,jpj) :: dic_bkg_temp ! DIC background
705	REAL(wp), DIMENSION(jpi,jpj) :: alk_bkg_temp ! Alkalinity background
706	REAL(wp), DIMENSION(jpi,jpj) :: tem_bkg_temp ! Temperature background
707	REAL(wp), DIMENSION(jpi,jpj) :: sal_bkg_temp ! Salinity background
708	REAL(wp), DIMENSION(1) :: patm ! Atmospheric pressure
709	REAL(wp), DIMENSION(1) :: phyd ! Hydrostatic pressure
710
711	! Coefficients for fCO2 to pCO2 conversion
712	REAL(wp) :: zcoef_fco2_1 = -1636.75
713	REAL(wp) :: zcoef_fco2_2 = 12.0408
714	REAL(wp) :: zcoef_fco2_3 = 0.0327957
715	REAL(wp) :: zcoef_fco2_4 = 0.0000316528
716	REAL(wp) :: zcoef_fco2_5 = 2.0
717	REAL(wp) :: zcoef_fco2_6 = 57.7
718	REAL(wp) :: zcoef_fco2_7 = 0.118
719	REAL(wp) :: zcoef_fco2_8 = 82.0578
720	!!------------------------------------------------------------------------
721
722	IF ( kt <= nit000 ) THEN
723
724	!----------------------------------------------------------------------
725	! DIC and alkalinity balancing
726	!----------------------------------------------------------------------
727
728	IF ( ln_sfco2inc ) THEN
729	#if defined key_medusa && defined key_roam
730	! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
731	patm(1) = 1.0
732	phyd(1) = 0.0
733	DO jj = 1, jpj
734	DO ji = 1, jpi
735	CALL f2pCO2( spco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
736	END DO
737	END DO
738	#else
739	! If assimilating fCO2, then convert to pCO2 using temperature
740	! See flux_gas.F90 within HadOCC for details of calculation
741	pco2_bkginc_temp(:,:) = spco2_bkginc(:,:) / &
742	& EXP((zcoef_fco2_1 + &
743	& zcoef_fco2_2 * (tsn(:,:,1,1)+rt0) - &
744	& zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
745	& zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0)(tsn(:,:,1,1)+rt0) + &
746	& zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0))) / &
747	& (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
748	#endif
749	ELSE
750	pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
751	ENDIF
752
753	! Account for physics assimilation if required
754	IF ( ll_trainc ) THEN
755	tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
756	sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
757	ELSE
758	tem_bkg_temp(:,:) = tsn(:,:,1,1)
759	sal_bkg_temp(:,:) = tsn(:,:,1,2)
760	ENDIF
761
762	#if defined key_medusa
763	! Account for logchl balancing if required
764	IF ( ln_slchltotinc .AND. ln_slchltotbal ) THEN
765	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + slchltot_balinc(:,:,1,jpdic)
766	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + slchltot_balinc(:,:,1,jpalk)
767	ELSE
768	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
769	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
770	ENDIF
771
772	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
773	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
774	& spco2_balinc(:,:,1,jpdic), spco2_balinc(:,:,1,jpalk) )
775
776	#elif defined key_hadocc
777	! Account for slchltot balancing if required
778	IF ( ln_slchltotinc .AND. ln_slchltotbal ) THEN
779	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + slchltot_balinc(:,:,1,jp_had_dic)
780	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + slchltot_balinc(:,:,1,jp_had_alk)
781	ELSE
782	dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
783	alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
784	ENDIF
785
786	CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
787	& tem_bkg_temp(:,:), sal_bkg_temp(:,:), &
788	& spco2_balinc(:,:,1,jp_had_dic), spco2_balinc(:,:,1,jp_had_alk) )
789
790	#else
791	CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
792	& 'but not defined a biogeochemical model' )
793	#endif
794
795	! Select mixed layer
796	IF ( ll_asmdin ) THEN
797	#if defined key_hadocc \|\| (defined key_medusa && defined key_foam_medusa)
798	CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
799	& ' Mixed layer depth taken to be background maximum mld_max_bkg' )
800	zmld(:,:) = mld_max_bkg(:,:)
801	#else
802	CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
803	#endif
804	ELSE
805	SELECT CASE( mld_choice_bgc )
806	CASE ( 1 ) ! Turbocline/mixing depth [W points]
807	zmld(:,:) = hmld(:,:)
808	CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
809	zmld(:,:) = hmlp(:,:)
810	CASE ( 3 ) ! Kara MLD [Interpolated]
811	#if defined key_karaml
812	IF ( ln_kara ) THEN
813	zmld(:,:) = hmld_kara(:,:)
814	ELSE
815	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
816	& ' but ln_kara=.false.' )
817	ENDIF
818	#else
819	CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
820	& ' but is not defined' )
821	#endif
822	CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points]
823	!zmld(:,:) = hmld_tref(:,:)
824	CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
825	& ' but is not available in this version' )
826	CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
827	zmld(:,:) = hmlpt(:,:)
828	END SELECT
829	ENDIF
830
831	! Propagate through mixed layer
832	DO jj = 1, jpj
833	DO ji = 1, jpi
834	!
835	jkmax = jpk-1
836	DO jk = jpk-1, 1, -1
837	IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. &
838	& ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
839	zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
840	jkmax = jk
841	ENDIF
842	END DO
843	!
844	#if defined key_top
845	DO jk = 2, jkmax
846	spco2_balinc(ji,jj,jk,:) = spco2_balinc(ji,jj,1,:)
847	END DO
848	#endif
849	!
850	END DO
851	END DO
852
853	ENDIF
854
855	IF ( ll_asmiau ) THEN
856
857	!--------------------------------------------------------------------
858	! Incremental Analysis Updating
859	!--------------------------------------------------------------------
860
861	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
862
863	it = kt - nit000 + 1
864	zincwgt = pwgtiau(it) ! IAU weight for the current time step
865	! note this is not a tendency so should not be divided by rdt
866
867	IF(lwp) THEN
868	WRITE(numout,*)
869	IF ( ln_spco2inc ) THEN
870	WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
871	& kt,' with IAU weight = ', pwgtiau(it)
872	ELSE IF ( ln_sfco2inc ) THEN
873	WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
874	& kt,' with IAU weight = ', pwgtiau(it)
875	ENDIF
876	WRITE(numout,*) '~~~~~~~~~~~~'
877	ENDIF
878
879	! Update the biogeochemical variables
880	! Add directly to trn and trb, rather than to tra, because tra gets
881	! reset to zero at the start of trc_stp, called after this routine
882	#if defined key_medusa && defined key_foam_medusa
883	DO jk = 1, jpkm1
884	trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
885	& spco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
886	trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
887	& spco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
888	END DO
889	#elif defined key_hadocc
890	DO jk = 1, jpkm1
891	trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
892	& spco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
893	trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
894	& spco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
895	END DO
896	#endif
897
898	! Do not deallocate arrays - needed by asm_bgc_bal_wri
899	! which is called at end of model run
900
901	ENDIF
902
903	ELSEIF ( ll_asmdin ) THEN
904
905	!--------------------------------------------------------------------
906	! Direct Initialization
907	!--------------------------------------------------------------------
908
909	IF ( kt == nitdin_r ) THEN
910
911	neuler = 0 ! Force Euler forward step
912
913	#if defined key_medusa && defined key_foam_medusa
914	! Initialize the now fields with the background + increment
915	! Background currently is what the model is initialised with
916	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', &
917	& ' Background state is taken from model rather than background file' )
918	trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
919	& spco2_balinc(:,:,:,jp_msa0:jp_msa1)
920	trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
921	#elif defined key_hadocc
922	! Initialize the now fields with the background + increment
923	! Background currently is what the model is initialised with
924	CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', &
925	& ' Background state is taken from model rather than background file' )
926	trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
927	& spco2_balinc(:,:,:,jp_had0:jp_had1)
928	trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
929	#endif
930
931	! Do not deallocate arrays - needed by asm_bgc_bal_wri
932	! which is called at end of model run
933	ENDIF
934	!
935	ENDIF
936	!
937	END SUBROUTINE spco2_asm_inc
938
939	!!===========================================================================
940
941	END MODULE asmbgc

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

Download in other formats: