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.

asminc.F90 in branches/UKMO/dev_r5518_v3.4_asm_nemovar_community_bgc_ersem/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_v3.4_asm_nemovar_community_bgc_ersem/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 6983

Last change on this file since 6983 was 6983, checked in by dford, 8 years ago
An initial version of code to apply chlorophyll increments to FABM-ERSEM.
File size: 99.6 KB

Line
1	MODULE asminc
2	!!======================================================================
3	!! * MODULE asminc *
4	!! Assimilation increment : Apply an increment generated by data
5	!! assimilation
6	!!======================================================================
7	!! History : ! 2007-03 (M. Martin) Met Office version
8	!! ! 2007-04 (A. Weaver) calc_date original code
9	!! ! 2007-04 (A. Weaver) Merge with OPAVAR/NEMOVAR
10	!! NEMO 3.3 ! 2010-05 (D. Lea) Update to work with NEMO v3.2
11	!! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init
12	!! 3.4 ! 2012-10 (A. Weaver and K. Mogensen) Fix for direct initialization
13	!!----------------------------------------------------------------------
14
15	!!----------------------------------------------------------------------
16	!! 'key_asminc' : Switch on the assimilation increment interface
17	!!----------------------------------------------------------------------
18	!! asm_inc_init : Initialize the increment arrays and IAU weights
19	!! calc_date : Compute the calendar date YYYYMMDD on a given step
20	!! tra_asm_inc : Apply the tracer (T and S) increments
21	!! dyn_asm_inc : Apply the dynamic (u and v) increments
22	!! ssh_asm_inc : Apply the SSH increment
23	!! seaice_asm_inc : Apply the seaice increment
24	!! logchl_asm_inc : Apply the logchl increment
25	!!----------------------------------------------------------------------
26	USE wrk_nemo ! Memory Allocation
27	USE par_oce ! Ocean space and time domain variables
28	USE dom_oce ! Ocean space and time domain
29	USE domvvl ! domain: variable volume level
30	USE oce ! Dynamics and active tracers defined in memory
31	USE ldfdyn_oce ! ocean dynamics: lateral physics
32	USE eosbn2 ! Equation of state - in situ and potential density
33	USE zpshde ! Partial step : Horizontal Derivative
34	USE iom ! Library to read input files
35	USE asmpar ! Parameters for the assmilation interface
36	USE c1d ! 1D initialization
37	USE in_out_manager ! I/O manager
38	USE lib_mpp ! MPP library
39	#if defined key_lim2
40	USE ice_2 ! LIM2
41	#endif
42	USE sbc_oce ! Surface boundary condition variables.
43	USE zdfmxl, ONLY : &
44	& hmld_tref, &
45	#if defined key_karaml
46	& hmld_kara, &
47	& ln_kara, &
48	#endif
49	& hmld, &
50	& hmlp, &
51	& hmlpt
52	#if defined key_bdy
53	USE bdy_oce, ONLY: bdytmask
54	#endif
55	#if defined key_fabm
56	USE asmlogchlbal_ersem, ONLY: &
57	& asm_logchl_bal_ersem
58	USE trc, ONLY: &
59	& trn, &
60	& trb
61	USE par_fabm
62	#elif defined key_medusa && defined key_foam_medusa
63	USE asmlogchlbal_medusa, ONLY: &
64	& asm_logchl_bal_medusa
65	USE trc, ONLY: &
66	& trn, &
67	& trb
68	USE par_medusa
69	#elif defined key_hadocc
70	USE asmlogchlbal_hadocc, ONLY: &
71	& asm_logchl_bal_hadocc
72	USE trc, ONLY: &
73	& trn, &
74	& trb
75	USE par_hadocc
76	#endif
77
78	IMPLICIT NONE
79	PRIVATE
80
81	PUBLIC asm_inc_init !: Initialize the increment arrays and IAU weights
82	PUBLIC calc_date !: Compute the calendar date YYYYMMDD on a given step
83	PUBLIC tra_asm_inc !: Apply the tracer (T and S) increments
84	PUBLIC dyn_asm_inc !: Apply the dynamic (u and v) increments
85	PUBLIC ssh_asm_inc !: Apply the SSH increment
86	PUBLIC seaice_asm_inc !: Apply the seaice increment
87	PUBLIC logchl_asm_inc !: Apply the seaice increment
88
89	#if defined key_asminc
90	LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE. !: Logical switch for assimilation increment interface
91	#else
92	LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE. !: No assimilation increments
93	#endif
94	LOGICAL, PUBLIC :: ln_bkgwri = .FALSE. !: No output of the background state fields
95	LOGICAL, PUBLIC :: ln_balwri = .FALSE. !: No output of the assimilation balancing increments
96	LOGICAL, PUBLIC :: ln_asmiau = .FALSE. !: No applying forcing with an assimilation increment
97	LOGICAL, PUBLIC :: ln_asmdin = .FALSE. !: No direct initialization
98	LOGICAL, PUBLIC :: ln_trainc = .FALSE. !: No tracer (T and S) assimilation increments
99	LOGICAL, PUBLIC :: ln_dyninc = .FALSE. !: No dynamics (u and v) assimilation increments
100	LOGICAL, PUBLIC :: ln_sshinc = .FALSE. !: No sea surface height assimilation increment
101	LOGICAL, PUBLIC :: ln_seaiceinc = .FALSE. !: No sea ice concentration increment
102	LOGICAL, PUBLIC :: ln_logchltotinc = .FALSE. !: No total log10(chlorophyll) increment
103	LOGICAL, PUBLIC :: ln_logchlpftinc = .FALSE. !: No PFT log10(chlorophyll) increment
104	LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check
105	LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
106	INTEGER, PUBLIC :: nn_divdmp !: Apply divergence damping filter nn_divdmp times
107
108	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkg , s_bkg !: Background temperature and salinity
109	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkg , v_bkg !: Background u- & v- velocity components
110	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkginc, s_bkginc !: Increment to the background T & S
111	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components
112	REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step
113	#if defined key_asminc
114	REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssh_iau !: IAU-weighted sea surface height increment
115	#endif
116	! !!! time steps relative to the cycle interval [0,nitend-nit000-1]
117	INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term
118	INTEGER , PUBLIC :: nitdin !: Time step of the background state for direct initialization
119	INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval
120	INTEGER , PUBLIC :: nitiaufin !: Time step of the end of the IAU interval
121	!
122	INTEGER , PUBLIC :: niaufn !: Type of IAU weighing function: = 0 Constant weighting
123	! !: = 1 Linear hat-like, centred in middle of IAU interval
124	REAL(wp), PUBLIC :: salfixmin !: Ensure that the salinity is larger than this value if (ln_salfix)
125
126	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc ! Background sea surface height and its increment
127	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: seaice_bkginc ! Increment to the background sea ice conc
128
129	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_bkginc !: Increment to background logchl
130	#if defined key_fabm
131	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_chl1 !: Increment to ERSEM diatom chl from logchl
132	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_chl2 !: Increment to ERSEM nanophy chl from logchl
133	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_chl3 !: Increment to ERSEM picophy chl from logchl
134	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_chl4 !: Increment to ERSEM microphy chl from logchl
135	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p1c !: Increment to ERSEM diatom c from logchl
136	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p1n !: Increment to ERSEM diatom n from logchl
137	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p1p !: Increment to ERSEM diatom p from logchl
138	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p1s !: Increment to ERSEM diatom s from logchl
139	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p2c !: Increment to ERSEM nanophy c from logchl
140	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p2n !: Increment to ERSEM nanophy n from logchl
141	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p2p !: Increment to ERSEM nanophy p from logchl
142	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p3c !: Increment to ERSEM picophy c from logchl
143	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p3n !: Increment to ERSEM picophy n from logchl
144	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p3p !: Increment to ERSEM picophy p from logchl
145	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p4c !: Increment to ERSEM microphy c from logchl
146	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p4n !: Increment to ERSEM microphy n from logchl
147	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_p4p !: Increment to ERSEM microphy p from logchl
148	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_z4c !: Increment to ERSEM mesozoo c from logchl
149	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_z5c !: Increment to ERSEM microzoo c from logchl
150	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_z5n !: Increment to ERSEM microzoo n from logchl
151	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_z5p !: Increment to ERSEM microzoo p from logchl
152	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_z6c !: Increment to ERSEM het flag c from logchl
153	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_z6n !: Increment to ERSEM het flag n from logchl
154	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_z6p !: Increment to ERSEM het flag p from logchl
155	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_n1p !: Increment to ERSEM phosphate from logchl
156	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_n3n !: Increment to ERSEM nitrate from logchl
157	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_n4n !: Increment to ERSEM ammonium from logchl
158	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_ersem_n5s !: Increment to ERSEM silicate from logchl
159	#elif defined key_medusa && defined key_foam_medusa
160	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_chn !: Increment to MEDUSA nondiatom chl from logchl
161	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_chd !: Increment to MEDUSA diatom chl from logchl
162	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_phn !: Increment to MEDUSA nondiatom n from logchl
163	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_phd !: Increment to MEDUSA diatom n from logchl
164	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_pds !: Increment to MEDUSA diatom s from logchl
165	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_zmi !: Increment to MEDUSA microzoop n from logchl
166	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_zme !: Increment to MEDUSA mesozoop n from logchl
167	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_din !: Increment to MEDUSA nitrate from logchl
168	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_sil !: Increment to MEDUSA silicate from logchl
169	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_fer !: Increment to MEDUSA iron from logchl
170	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_det !: Increment to MEDUSA detritus n from logchl
171	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_dtc !: Increment to MEDUSA detritus c from logchl
172	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_dic !: Increment to MEDUSA dic from logchl
173	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_alk !: Increment to MEDUSA alkalinity from logchl
174	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_medusa_oxy !: Increment to MEDUSA oxygen from logchl
175	#elif defined key_hadocc
176	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_hadocc_nut !: Increment to HadOCC nutrient from logchl
177	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_hadocc_phy !: Increment to HadOCC phytoplankton from logchl
178	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_hadocc_zoo !: Increment to HadOCC zooplankton from logchl
179	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_hadocc_det !: Increment to HadOCC detritus from logchl
180	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_hadocc_dic !: Increment to HadOCC DIC from logchl
181	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: logchl_balinc_hadocc_alk !: Increment to HadOCC alkalinity from logchl
182	#endif
183
184	INTEGER :: mld_choice = 4 !: choice of mld criteria to use for physics assimilation
185	!: 1) hmld - Turbocline/mixing depth [W points]
186	!: 2) hmlp - Density criterion (0.01 kg/m^3 change from 10m) [W points]
187	!: 3) hmld_kara - Kara MLD [Interpolated]
188	!: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points]
189
190	INTEGER :: mld_choice_bgc = 5 !: choice of mld criteria to use for physics assimilation
191	!: 1) hmld - Turbocline/mixing depth [W points]
192	!: 2) hmlp - Density criterion (0.01 kg/m^3 change from 10m) [W points]
193	!: 3) hmld_kara - Kara MLD [Interpolated]
194	!: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points]
195	!: 5) hmlpt - Density criterion (0.01 kg/m^3 change from 10m) [T points]
196
197	INTEGER :: nn_asmpfts = 0 !: number of logchl PFTs assimilated
198
199	!! * Substitutions
200	# include "domzgr_substitute.h90"
201	# include "ldfdyn_substitute.h90"
202	# include "vectopt_loop_substitute.h90"
203	!!----------------------------------------------------------------------
204	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
205	!! $Id$
206	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
207	!!----------------------------------------------------------------------
208	CONTAINS
209
210	SUBROUTINE asm_inc_init
211	!!----------------------------------------------------------------------
212	!! * ROUTINE asm_inc_init *
213	!!
214	!! ** Purpose : Initialize the assimilation increment and IAU weights.
215	!!
216	!! ** Method : Initialize the assimilation increment and IAU weights.
217	!!
218	!! ** Action :
219	!!----------------------------------------------------------------------
220	INTEGER :: ji, jj, jk, jt ! dummy loop indices
221	INTEGER :: imid, inum ! local integers
222	INTEGER :: ios ! Local integer output status for namelist read
223	INTEGER :: iiauper ! Number of time steps in the IAU period
224	INTEGER :: icycper ! Number of time steps in the cycle
225	INTEGER :: iitend_date ! Date YYYYMMDD of final time step
226	INTEGER :: iitbkg_date ! Date YYYYMMDD of background time step for Jb term
227	INTEGER :: iitdin_date ! Date YYYYMMDD of background time step for DI
228	INTEGER :: iitiaustr_date ! Date YYYYMMDD of IAU interval start time step
229	INTEGER :: iitiaufin_date ! Date YYYYMMDD of IAU interval final time step
230	INTEGER :: isurfstat ! Local integer for status of reading surft variable
231	!
232	REAL(wp) :: znorm ! Normalization factor for IAU weights
233	REAL(wp) :: ztotwgt ! Value of time-integrated IAU weights (should be equal to one)
234	REAL(wp) :: z_inc_dateb ! Start date of interval on which increment is valid
235	REAL(wp) :: z_inc_datef ! End date of interval on which increment is valid
236	REAL(wp) :: zdate_bkg ! Date in background state file for DI
237	REAL(wp) :: zdate_inc ! Time axis in increments file
238	!
239	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: &
240	& t_bkginc_2d ! file for reading in 2D
241	! ! temperature increments
242	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: &
243	& z_mld ! Mixed layer depth
244
245	REAL(wp), POINTER, DIMENSION(:,:) :: hdiv ! 2D workspace
246	!
247	LOGICAL :: lk_surft ! Logical: T => Increments file contains surft variable
248	! so only apply surft increments.
249	!
250	CHARACTER(LEN=2) :: cl_pftstr
251	!!
252	NAMELIST/nam_asminc/ ln_bkgwri, ln_balwri, &
253	& ln_trainc, ln_dyninc, ln_sshinc, &
254	& ln_logchltotinc, ln_logchlpftinc, &
255	& ln_asmdin, ln_asmiau, &
256	& nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, &
257	& ln_salfix, salfixmin, nn_divdmp, mld_choice, &
258	& mld_choice_bgc
259	!!----------------------------------------------------------------------
260
261	!-----------------------------------------------------------------------
262	! Read Namelist nam_asminc : assimilation increment interface
263	!-----------------------------------------------------------------------
264
265	! Set default values
266	ln_bkgwri = .FALSE.
267	ln_balwri = .FALSE.
268	ln_trainc = .FALSE.
269	ln_dyninc = .FALSE.
270	ln_sshinc = .FALSE.
271	ln_seaiceinc = .FALSE.
272	ln_logchltotinc = .FALSE.
273	ln_logchlpftinc = .FALSE.
274	ln_asmdin = .FALSE.
275	ln_asmiau = .TRUE.
276	ln_salfix = .FALSE.
277	ln_temnofreeze = .FALSE.
278	salfixmin = -9999
279	nitbkg = 0
280	nitdin = 0
281	nitiaustr = 1
282	nitiaufin = 150
283	niaufn = 0
284	#if defined key_fabm
285	nn_asmpfts = 4
286	#elif defined key_medusa && defined key_foam_medusa
287	nn_asmpfts = 2
288	#elif defined key_hadocc
289	nn_asmpfts = 1
290	#else
291	nn_asmpfts = 0
292	#endif
293
294	REWIND( numnam_ref ) ! Namelist nam_asminc in reference namelist : Assimilation increment
295	READ ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
296	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in reference namelist', lwp )
297
298	REWIND( numnam_cfg ) ! Namelist nam_asminc in configuration namelist : Assimilation increment
299	READ ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
300	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in configuration namelist', lwp )
301	IF(lwm) WRITE ( numond, nam_asminc )
302
303	IF ( ln_logchltotinc ) THEN
304	nn_asmpfts = 1
305	ELSE IF ( .NOT.( ln_logchlpftinc ) ) THEN
306	nn_asmpfts = 0
307	ENDIF
308
309	! Control print
310	IF(lwp) THEN
311	WRITE(numout,*)
312	WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
313	WRITE(numout,*) '~~~~~~~~~~~~'
314	WRITE(numout,*) ' Namelist namasm : set assimilation increment parameters'
315	WRITE(numout,*) ' Logical switch for writing out background state ln_bkgwri = ', ln_bkgwri
316	WRITE(numout,*) ' Logical switch for writing out balancing increments ln_balwri = ', ln_balwri
317	WRITE(numout,*) ' Logical switch for applying tracer increments ln_trainc = ', ln_trainc
318	WRITE(numout,*) ' Logical switch for applying velocity increments ln_dyninc = ', ln_dyninc
319	WRITE(numout,*) ' Logical switch for applying SSH increments ln_sshinc = ', ln_sshinc
320	WRITE(numout,*) ' Logical switch for Direct Initialization (DI) ln_asmdin = ', ln_asmdin
321	WRITE(numout,*) ' Logical switch for applying sea ice increments ln_seaiceinc = ', ln_seaiceinc
322	WRITE(numout,*) ' Logical switch for applying total logchl incs ln_logchltotinc = ', ln_logchltotinc
323	WRITE(numout,*) ' Logical switch for applying PFT logchl incs ln_logchlpftinc = ', ln_logchlpftinc
324	WRITE(numout,*) ' Number of logchl PFTs assimilated nn_asmpfts = ', nn_asmpfts
325	WRITE(numout,*) ' Logical switch for Incremental Analysis Updating (IAU) ln_asmiau = ', ln_asmiau
326	WRITE(numout,*) ' Timestep of background in [0,nitend-nit000-1] nitbkg = ', nitbkg
327	WRITE(numout,*) ' Timestep of background for DI in [0,nitend-nit000-1] nitdin = ', nitdin
328	WRITE(numout,*) ' Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
329	WRITE(numout,*) ' Timestep of end of IAU interval in [0,nitend-nit000-1] nitiaufin = ', nitiaufin
330	WRITE(numout,*) ' Type of IAU weighting function niaufn = ', niaufn
331	WRITE(numout,*) ' Logical switch for ensuring that the sa > salfixmin ln_salfix = ', ln_salfix
332	WRITE(numout,*) ' Minimum salinity after applying the increments salfixmin = ', salfixmin
333	WRITE(numout,*) ' Choice of MLD for physics assimilation mld_choice = ', mld_choice
334	WRITE(numout,*) ' Choice of MLD for BGC assimilation mld_choice_bgc = ', mld_choice_bgc
335	ENDIF
336
337	nitbkg_r = nitbkg + nit000 - 1 ! Background time referenced to nit000
338	nitdin_r = nitdin + nit000 - 1 ! Background time for DI referenced to nit000
339	nitiaustr_r = nitiaustr + nit000 - 1 ! Start of IAU interval referenced to nit000
340	nitiaufin_r = nitiaufin + nit000 - 1 ! End of IAU interval referenced to nit000
341
342	iiauper = nitiaufin_r - nitiaustr_r + 1 ! IAU interval length
343	icycper = nitend - nit000 + 1 ! Cycle interval length
344
345	CALL calc_date( nit000, nitend , ndate0, iitend_date ) ! Date of final time step
346	CALL calc_date( nit000, nitbkg_r , ndate0, iitbkg_date ) ! Background time for Jb referenced to ndate0
347	CALL calc_date( nit000, nitdin_r , ndate0, iitdin_date ) ! Background time for DI referenced to ndate0
348	CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date ) ! IAU start time referenced to ndate0
349	CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date ) ! IAU end time referenced to ndate0
350	!
351	IF(lwp) THEN
352	WRITE(numout,*)
353	WRITE(numout,*) ' Time steps referenced to current cycle:'
354	WRITE(numout,*) ' iitrst = ', nit000 - 1
355	WRITE(numout,*) ' nit000 = ', nit000
356	WRITE(numout,*) ' nitend = ', nitend
357	WRITE(numout,*) ' nitbkg_r = ', nitbkg_r
358	WRITE(numout,*) ' nitdin_r = ', nitdin_r
359	WRITE(numout,*) ' nitiaustr_r = ', nitiaustr_r
360	WRITE(numout,*) ' nitiaufin_r = ', nitiaufin_r
361	WRITE(numout,*)
362	WRITE(numout,*) ' Dates referenced to current cycle:'
363	WRITE(numout,*) ' ndastp = ', ndastp
364	WRITE(numout,*) ' ndate0 = ', ndate0
365	WRITE(numout,*) ' iitend_date = ', iitend_date
366	WRITE(numout,*) ' iitbkg_date = ', iitbkg_date
367	WRITE(numout,*) ' iitdin_date = ', iitdin_date
368	WRITE(numout,*) ' iitiaustr_date = ', iitiaustr_date
369	WRITE(numout,*) ' iitiaufin_date = ', iitiaufin_date
370	ENDIF
371
372	IF ( nacc /= 0 ) &
373	& CALL ctl_stop( ' nacc /= 0 and key_asminc :', &
374	& ' Assimilation increments have only been implemented', &
375	& ' for synchronous time stepping' )
376
377	IF ( ( ln_asmdin ).AND.( ln_asmiau ) ) &
378	& CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
379	& ' Choose Direct Initialization OR Incremental Analysis Updating')
380
381	IF ( ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
382	& .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. &
383	& ( ln_logchltotinc ).OR.( ln_logchlpftinc ) )) &
384	& CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
385	& ' ln_logchltotinc, and ln_logchlpftinc is set to .true.', &
386	& ' but ln_asmdin and ln_asmiau are both set to .false. :', &
387	& ' Inconsistent options')
388
389	IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
390	& CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :', &
391	& ' Type IAU weighting function is invalid')
392
393	IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
394	& .AND.( .NOT. ln_logchltotinc ).AND.( .NOT. ln_logchlpftinc ) ) &
395	& CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
396	& ' ln_logchltotinc, and ln_logchlpftinc are set to .false. :', &
397	& ' The assimilation increments are not applied')
398
399	IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
400	& CALL ctl_stop( ' nitiaustr = nitiaufin :', &
401	& ' IAU interval is of zero length')
402
403	IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
404	& CALL ctl_stop( ' nitiaustr or nitiaufin :', &
405	& ' IAU starting or final time step is outside the cycle interval', &
406	& ' Valid range nit000 to nitend')
407
408	IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
409	& CALL ctl_stop( ' nitbkg :', &
410	& ' Background time step is outside the cycle interval')
411
412	IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
413	& CALL ctl_stop( ' nitdin :', &
414	& ' Background time step for Direct Initialization is outside', &
415	& ' the cycle interval')
416
417	IF ( ( ln_logchltotinc ).AND.( ln_logchlpftinc ) ) THEN
418	CALL ctl_stop( ' ln_logchltotinc and ln_logchlpftinc both set:', &
419	& ' These options are not compatible')
420	ENDIF
421
422	IF ( ( ln_balwri ).AND.( .NOT. ( ( ln_logchltotinc ).OR.( ln_logchlpftinc ) ) ) ) THEN
423	CALL ctl_warn( ' Balancing increments are only calculated for logchl', &
424	& ' Not assimilating logchl, so ln_balwri will be set to .false.')
425	ln_balwri = .FALSE.
426	ENDIF
427
428	IF ( nstop > 0 ) RETURN ! if there are any errors then go no further
429
430	!--------------------------------------------------------------------
431	! Initialize the Incremental Analysis Updating weighting function
432	!--------------------------------------------------------------------
433
434	IF ( ln_asmiau ) THEN
435
436	ALLOCATE( wgtiau( icycper ) )
437
438	wgtiau(:) = 0.0
439
440	IF ( niaufn == 0 ) THEN
441
442	!---------------------------------------------------------
443	! Constant IAU forcing
444	!---------------------------------------------------------
445
446	DO jt = 1, iiauper
447	wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
448	END DO
449
450	ELSEIF ( niaufn == 1 ) THEN
451
452	!---------------------------------------------------------
453	! Linear hat-like, centred in middle of IAU interval
454	!---------------------------------------------------------
455
456	! Compute the normalization factor
457	znorm = 0.0
458	IF ( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval
459	imid = iiauper / 2
460	DO jt = 1, imid
461	znorm = znorm + REAL( jt )
462	END DO
463	znorm = 2.0 * znorm
464	ELSE ! Odd number of time steps in IAU interval
465	imid = ( iiauper + 1 ) / 2
466	DO jt = 1, imid - 1
467	znorm = znorm + REAL( jt )
468	END DO
469	znorm = 2.0 * znorm + REAL( imid )
470	ENDIF
471	znorm = 1.0 / znorm
472
473	DO jt = 1, imid - 1
474	wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
475	END DO
476	DO jt = imid, iiauper
477	wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
478	END DO
479
480	ENDIF
481
482	! Test that the integral of the weights over the weighting interval equals 1
483	IF(lwp) THEN
484	WRITE(numout,*)
485	WRITE(numout,*) 'asm_inc_init : IAU weights'
486	WRITE(numout,*) '~~~~~~~~~~~~'
487	WRITE(numout,*) ' time step IAU weight'
488	WRITE(numout,*) ' ========= ====================='
489	ztotwgt = 0.0
490	DO jt = 1, icycper
491	ztotwgt = ztotwgt + wgtiau(jt)
492	WRITE(numout,*) ' ', jt, ' ', wgtiau(jt)
493	END DO
494	WRITE(numout,*) ' ==================================='
495	WRITE(numout,*) ' Time-integrated weight = ', ztotwgt
496	WRITE(numout,*) ' ==================================='
497	ENDIF
498
499	ENDIF
500
501	!--------------------------------------------------------------------
502	! Allocate and initialize the increment arrays
503	!--------------------------------------------------------------------
504
505	IF ( ln_trainc ) THEN
506	ALLOCATE( t_bkginc(jpi,jpj,jpk) )
507	ALLOCATE( s_bkginc(jpi,jpj,jpk) )
508	t_bkginc(:,:,:) = 0.0
509	s_bkginc(:,:,:) = 0.0
510	ENDIF
511	IF ( ln_dyninc ) THEN
512	ALLOCATE( u_bkginc(jpi,jpj,jpk) )
513	ALLOCATE( v_bkginc(jpi,jpj,jpk) )
514	u_bkginc(:,:,:) = 0.0
515	v_bkginc(:,:,:) = 0.0
516	ENDIF
517	IF ( ln_sshinc ) THEN
518	ALLOCATE( ssh_bkginc(jpi,jpj) )
519	ssh_bkginc(:,:) = 0.0
520	ENDIF
521	IF ( ln_seaiceinc ) THEN
522	ALLOCATE( seaice_bkginc(jpi,jpj))
523	seaice_bkginc(:,:) = 0.0
524	ENDIF
525	#if defined key_asminc
526	ALLOCATE( ssh_iau(jpi,jpj) )
527	ssh_iau(:,:) = 0.0
528	#endif
529	IF ( ( ln_logchltotinc ).OR.( ln_logchlpftinc ) ) THEN
530	IF ( ln_logchltotinc ) THEN
531	ALLOCATE( logchl_bkginc(jpi,jpj,1))
532	ELSE IF ( ln_logchlpftinc ) THEN
533	ALLOCATE( logchl_bkginc(jpi,jpj,nn_asmpfts))
534	ENDIF
535	logchl_bkginc(:,:,:) = 0.0
536	#if defined key_fabm
537	ALLOCATE( logchl_balinc_ersem_chl1(jpi,jpj,jpk) )
538	ALLOCATE( logchl_balinc_ersem_chl2(jpi,jpj,jpk) )
539	ALLOCATE( logchl_balinc_ersem_chl3(jpi,jpj,jpk) )
540	ALLOCATE( logchl_balinc_ersem_chl4(jpi,jpj,jpk) )
541	ALLOCATE( logchl_balinc_ersem_p1c(jpi,jpj,jpk) )
542	ALLOCATE( logchl_balinc_ersem_p1n(jpi,jpj,jpk) )
543	ALLOCATE( logchl_balinc_ersem_p1p(jpi,jpj,jpk) )
544	ALLOCATE( logchl_balinc_ersem_p1s(jpi,jpj,jpk) )
545	ALLOCATE( logchl_balinc_ersem_p2c(jpi,jpj,jpk) )
546	ALLOCATE( logchl_balinc_ersem_p2n(jpi,jpj,jpk) )
547	ALLOCATE( logchl_balinc_ersem_p2p(jpi,jpj,jpk) )
548	ALLOCATE( logchl_balinc_ersem_p3c(jpi,jpj,jpk) )
549	ALLOCATE( logchl_balinc_ersem_p3n(jpi,jpj,jpk) )
550	ALLOCATE( logchl_balinc_ersem_p3p(jpi,jpj,jpk) )
551	ALLOCATE( logchl_balinc_ersem_p4c(jpi,jpj,jpk) )
552	ALLOCATE( logchl_balinc_ersem_p4n(jpi,jpj,jpk) )
553	ALLOCATE( logchl_balinc_ersem_p4p(jpi,jpj,jpk) )
554	ALLOCATE( logchl_balinc_ersem_z4c(jpi,jpj,jpk) )
555	ALLOCATE( logchl_balinc_ersem_z5c(jpi,jpj,jpk) )
556	ALLOCATE( logchl_balinc_ersem_z5n(jpi,jpj,jpk) )
557	ALLOCATE( logchl_balinc_ersem_z5p(jpi,jpj,jpk) )
558	ALLOCATE( logchl_balinc_ersem_z6c(jpi,jpj,jpk) )
559	ALLOCATE( logchl_balinc_ersem_z6n(jpi,jpj,jpk) )
560	ALLOCATE( logchl_balinc_ersem_z6p(jpi,jpj,jpk) )
561	ALLOCATE( logchl_balinc_ersem_n1p(jpi,jpj,jpk) )
562	ALLOCATE( logchl_balinc_ersem_n3n(jpi,jpj,jpk) )
563	ALLOCATE( logchl_balinc_ersem_n4n(jpi,jpj,jpk) )
564	ALLOCATE( logchl_balinc_ersem_n5s(jpi,jpj,jpk) )
565	logchl_balinc_ersem_chl1(:,:,:) = 0.0
566	logchl_balinc_ersem_chl2(:,:,:) = 0.0
567	logchl_balinc_ersem_chl3(:,:,:) = 0.0
568	logchl_balinc_ersem_chl4(:,:,:) = 0.0
569	logchl_balinc_ersem_p1c(:,:,:) = 0.0
570	logchl_balinc_ersem_p1n(:,:,:) = 0.0
571	logchl_balinc_ersem_p1p(:,:,:) = 0.0
572	logchl_balinc_ersem_p1s(:,:,:) = 0.0
573	logchl_balinc_ersem_p2c(:,:,:) = 0.0
574	logchl_balinc_ersem_p2n(:,:,:) = 0.0
575	logchl_balinc_ersem_p2p(:,:,:) = 0.0
576	logchl_balinc_ersem_p3c(:,:,:) = 0.0
577	logchl_balinc_ersem_p3n(:,:,:) = 0.0
578	logchl_balinc_ersem_p3p(:,:,:) = 0.0
579	logchl_balinc_ersem_p4c(:,:,:) = 0.0
580	logchl_balinc_ersem_p4n(:,:,:) = 0.0
581	logchl_balinc_ersem_p4p(:,:,:) = 0.0
582	logchl_balinc_ersem_z4c(:,:,:) = 0.0
583	logchl_balinc_ersem_z5c(:,:,:) = 0.0
584	logchl_balinc_ersem_z5n(:,:,:) = 0.0
585	logchl_balinc_ersem_z5p(:,:,:) = 0.0
586	logchl_balinc_ersem_z6c(:,:,:) = 0.0
587	logchl_balinc_ersem_z6n(:,:,:) = 0.0
588	logchl_balinc_ersem_z6p(:,:,:) = 0.0
589	logchl_balinc_ersem_n1p(:,:,:) = 0.0
590	logchl_balinc_ersem_n3n(:,:,:) = 0.0
591	logchl_balinc_ersem_n4n(:,:,:) = 0.0
592	logchl_balinc_ersem_n5s(:,:,:) = 0.0
593	#elif defined key_medusa && defined key_foam_medusa
594	ALLOCATE( logchl_balinc_medusa_chn(jpi,jpj,jpk) )
595	ALLOCATE( logchl_balinc_medusa_chd(jpi,jpj,jpk) )
596	ALLOCATE( logchl_balinc_medusa_phn(jpi,jpj,jpk) )
597	ALLOCATE( logchl_balinc_medusa_phd(jpi,jpj,jpk) )
598	ALLOCATE( logchl_balinc_medusa_pds(jpi,jpj,jpk) )
599	ALLOCATE( logchl_balinc_medusa_zmi(jpi,jpj,jpk) )
600	ALLOCATE( logchl_balinc_medusa_zme(jpi,jpj,jpk) )
601	ALLOCATE( logchl_balinc_medusa_din(jpi,jpj,jpk) )
602	ALLOCATE( logchl_balinc_medusa_sil(jpi,jpj,jpk) )
603	ALLOCATE( logchl_balinc_medusa_fer(jpi,jpj,jpk) )
604	ALLOCATE( logchl_balinc_medusa_det(jpi,jpj,jpk) )
605	ALLOCATE( logchl_balinc_medusa_dtc(jpi,jpj,jpk) )
606	ALLOCATE( logchl_balinc_medusa_dic(jpi,jpj,jpk) )
607	ALLOCATE( logchl_balinc_medusa_alk(jpi,jpj,jpk) )
608	ALLOCATE( logchl_balinc_medusa_oxy(jpi,jpj,jpk) )
609	logchl_balinc_medusa_chn(:,:,:) = 0.0
610	logchl_balinc_medusa_chd(:,:,:) = 0.0
611	logchl_balinc_medusa_phn(:,:,:) = 0.0
612	logchl_balinc_medusa_phd(:,:,:) = 0.0
613	logchl_balinc_medusa_pds(:,:,:) = 0.0
614	logchl_balinc_medusa_zmi(:,:,:) = 0.0
615	logchl_balinc_medusa_zme(:,:,:) = 0.0
616	logchl_balinc_medusa_din(:,:,:) = 0.0
617	logchl_balinc_medusa_sil(:,:,:) = 0.0
618	logchl_balinc_medusa_fer(:,:,:) = 0.0
619	logchl_balinc_medusa_det(:,:,:) = 0.0
620	logchl_balinc_medusa_dtc(:,:,:) = 0.0
621	logchl_balinc_medusa_dic(:,:,:) = 0.0
622	logchl_balinc_medusa_alk(:,:,:) = 0.0
623	logchl_balinc_medusa_oxy(:,:,:) = 0.0
624	#elif defined key_hadocc
625	ALLOCATE( logchl_balinc_hadocc_nut(jpi,jpj,jpk) )
626	ALLOCATE( logchl_balinc_hadocc_phy(jpi,jpj,jpk) )
627	ALLOCATE( logchl_balinc_hadocc_zoo(jpi,jpj,jpk) )
628	ALLOCATE( logchl_balinc_hadocc_det(jpi,jpj,jpk) )
629	ALLOCATE( logchl_balinc_hadocc_dic(jpi,jpj,jpk) )
630	ALLOCATE( logchl_balinc_hadocc_alk(jpi,jpj,jpk) )
631	logchl_balinc_hadocc_nut(:,:,:) = 0.0
632	logchl_balinc_hadocc_phy(:,:,:) = 0.0
633	logchl_balinc_hadocc_zoo(:,:,:) = 0.0
634	logchl_balinc_hadocc_det(:,:,:) = 0.0
635	logchl_balinc_hadocc_dic(:,:,:) = 0.0
636	logchl_balinc_hadocc_alk(:,:,:) = 0.0
637	#endif
638	ENDIF
639	IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) &
640	& .OR.( ln_logchltotinc ).OR.( ln_logchlpftinc ) ) THEN
641
642	!--------------------------------------------------------------------
643	! Read the increments from file
644	!--------------------------------------------------------------------
645
646	CALL iom_open( c_asminc, inum )
647
648	CALL iom_get( inum, 'time', zdate_inc )
649
650	CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
651	CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
652	z_inc_dateb = zdate_inc
653	z_inc_datef = zdate_inc
654
655	IF(lwp) THEN
656	WRITE(numout,*)
657	WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', &
658	& ' between dates ', NINT( z_inc_dateb ),' and ', &
659	& NINT( z_inc_datef )
660	WRITE(numout,*) '~~~~~~~~~~~~'
661	ENDIF
662
663	IF ( ( NINT( z_inc_dateb ) < ndastp ) &
664	& .OR.( NINT( z_inc_datef ) > iitend_date ) ) &
665	& CALL ctl_warn( ' Validity time of assimilation increments is ', &
666	& ' outside the assimilation interval' )
667
668	IF ( ( ln_asmdin ).AND.( NINT( zdate_inc ) /= iitdin_date ) ) &
669	& CALL ctl_warn( ' Validity time of assimilation increments does ', &
670	& ' not agree with Direct Initialization time' )
671
672	IF ( ln_trainc ) THEN
673
674	!Test if the increments file contains the surft variable.
675	isurfstat = iom_varid( inum, 'bckinsurft', ldstop = .FALSE. )
676	IF ( isurfstat == -1 ) THEN
677	lk_surft = .FALSE.
678	ELSE
679	lk_surft = .TRUE.
680	CALL ctl_warn( ' Applying 2D temperature increment to bottom of ML: ', &
681	& ' bckinsurft found in increments file.' )
682	ENDIF
683
684	IF (lk_surft) THEN
685
686	ALLOCATE(z_mld(jpi,jpj))
687	SELECT CASE(mld_choice)
688	CASE(1)
689	z_mld = hmld
690	CASE(2)
691	z_mld = hmlp
692	CASE(3)
693	#if defined key_karaml
694	IF ( ln_kara ) THEN
695	z_mld = hmld_kara
696	ELSE
697	CALL ctl_stop("Kara mixed layer not calculated as ln_kara=.false.")
698	ENDIF
699	#else
700	CALL ctl_stop("Kara mixed layer not defined in current version of NEMO") ! JW: Safety feature, should be removed
701	! once the Kara mixed layer is available
702	#endif
703	CASE(4)
704	z_mld = hmld_tref
705	END SELECT
706
707	ALLOCATE( t_bkginc_2d(jpi,jpj) )
708	CALL iom_get( inum, jpdom_autoglo, 'bckinsurft', t_bkginc_2d, 1)
709	#if defined key_bdy
710	DO jk = 1,jpkm1
711	WHERE( z_mld(:,:) > fsdepw(:,:,jk) )
712	t_bkginc(:,:,jk) = t_bkginc_2d(:,:) * 0.5 * &
713	& ( 1 + cos( (fsdept(:,:,jk)/z_mld(:,:) ) * rpi ) )
714
715	t_bkginc(:,:,jk) = t_bkginc(:,:,jk) * bdytmask(:,:)
716	ELSEWHERE
717	t_bkginc(:,:,jk) = 0.
718	ENDWHERE
719	ENDDO
720	#else
721	t_bkginc(:,:,:) = 0.
722	#endif
723	s_bkginc(:,:,:) = 0.
724
725	DEALLOCATE(z_mld, t_bkginc_2d)
726
727	ELSE
728
729	CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
730	CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
731	! Apply the masks
732	t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
733	s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
734	! Set missing increments to 0.0 rather than 1e+20
735	! to allow for differences in masks
736	WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
737	WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
738
739	ENDIF
740
741	ENDIF
742
743	IF ( ln_dyninc ) THEN
744	CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )
745	CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )
746	! Apply the masks
747	u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
748	v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
749	! Set missing increments to 0.0 rather than 1e+20
750	! to allow for differences in masks
751	WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
752	WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
753	ENDIF
754
755	IF ( ln_sshinc ) THEN
756	CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
757	! Apply the masks
758	ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
759	! Set missing increments to 0.0 rather than 1e+20
760	! to allow for differences in masks
761	WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
762	ENDIF
763
764	IF ( ln_seaiceinc ) THEN
765	CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 )
766	! Apply the masks
767	seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1)
768	! Set missing increments to 0.0 rather than 1e+20
769	! to allow for differences in masks
770	WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
771	ENDIF
772
773	IF ( ln_logchltotinc ) THEN
774	CALL iom_get( inum, jpdom_autoglo, 'bckinlogchl', logchl_bkginc(:,:,1), 1 )
775	! Apply the masks
776	logchl_bkginc(:,:,1) = logchl_bkginc(:,:,1) * tmask(:,:,1)
777	! Set missing increments to 0.0 rather than 1e+20
778	! to allow for differences in masks
779	WHERE( ABS( logchl_bkginc(:,:,1) ) > 1.0e+10 ) logchl_bkginc(:,:,1) = 0.0
780	ENDIF
781
782	IF ( ln_logchlpftinc ) THEN
783	DO jt = 1, nn_asmpfts
784	WRITE(cl_pftstr,'(I2.2)') jt
785	CALL iom_get( inum, jpdom_autoglo, 'bckinlogchl'//cl_pftstr, logchl_bkginc(:,:,jt), 1 )
786	! Apply the masks
787	logchl_bkginc(:,:,jt) = logchl_bkginc(:,:,jt) * tmask(:,:,1)
788	! Set missing increments to 0.0 rather than 1e+20
789	! to allow for differences in masks
790	WHERE( ABS( logchl_bkginc(:,:,jt) ) > 1.0e+10 ) logchl_bkginc(:,:,jt) = 0.0
791	END DO
792	ENDIF
793
794	CALL iom_close( inum )
795
796	ENDIF
797
798	!-----------------------------------------------------------------------
799	! Apply divergence damping filter
800	!-----------------------------------------------------------------------
801
802	IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN
803
804	CALL wrk_alloc(jpi,jpj,hdiv)
805
806	DO jt = 1, nn_divdmp
807
808	DO jk = 1, jpkm1
809
810	hdiv(:,:) = 0._wp
811
812	DO jj = 2, jpjm1
813	DO ji = fs_2, fs_jpim1 ! vector opt.
814	hdiv(ji,jj) = &
815	( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * u_bkginc(ji ,jj ,jk) &
816	- e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * u_bkginc(ji-1,jj ,jk) &
817	+ e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * v_bkginc(ji ,jj ,jk) &
818	- e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * v_bkginc(ji ,jj-1,jk) ) &
819	/ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
820	END DO
821	END DO
822
823	CALL lbc_lnk( hdiv, 'T', 1. ) ! lateral boundary cond. (no sign change)
824
825	DO jj = 2, jpjm1
826	DO ji = fs_2, fs_jpim1 ! vector opt.
827	u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji+1,jj)e2t(ji+1,jj) hdiv(ji+1,jj) &
828	- e1t(ji ,jj)e2t(ji ,jj) hdiv(ji ,jj) ) &
829	/ e1u(ji,jj) * umask(ji,jj,jk)
830	v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji,jj+1)e2t(ji,jj+1) hdiv(ji,jj+1) &
831	- e1t(ji,jj )e2t(ji,jj ) hdiv(ji,jj ) ) &
832	/ e2v(ji,jj) * vmask(ji,jj,jk)
833	END DO
834	END DO
835
836	END DO
837
838	END DO
839
840	CALL wrk_dealloc(jpi,jpj,hdiv)
841
842	ENDIF
843
844
845
846	!-----------------------------------------------------------------------
847	! Allocate and initialize the background state arrays
848	!-----------------------------------------------------------------------
849
850	IF ( ln_asmdin ) THEN
851
852	ALLOCATE( t_bkg(jpi,jpj,jpk) )
853	ALLOCATE( s_bkg(jpi,jpj,jpk) )
854	ALLOCATE( u_bkg(jpi,jpj,jpk) )
855	ALLOCATE( v_bkg(jpi,jpj,jpk) )
856	ALLOCATE( ssh_bkg(jpi,jpj) )
857
858	t_bkg(:,:,:) = 0.0
859	s_bkg(:,:,:) = 0.0
860	u_bkg(:,:,:) = 0.0
861	v_bkg(:,:,:) = 0.0
862	ssh_bkg(:,:) = 0.0
863
864	!--------------------------------------------------------------------
865	! Read from file the background state at analysis time
866	!--------------------------------------------------------------------
867
868	CALL iom_open( c_asmdin, inum )
869
870	CALL iom_get( inum, 'rdastp', zdate_bkg )
871
872	IF(lwp) THEN
873	WRITE(numout,*)
874	WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', &
875	& NINT( zdate_bkg )
876	WRITE(numout,*) '~~~~~~~~~~~~'
877	ENDIF
878
879	IF ( NINT( zdate_bkg ) /= iitdin_date ) &
880	& CALL ctl_warn( ' Validity time of assimilation background state does', &
881	& ' not agree with Direct Initialization time' )
882
883	IF ( ln_trainc ) THEN
884	CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
885	CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
886	t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
887	s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
888	ENDIF
889
890	IF ( ln_dyninc ) THEN
891	CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
892	CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
893	u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
894	v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
895	ENDIF
896
897	IF ( ln_sshinc ) THEN
898	CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
899	ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
900	ENDIF
901
902	CALL iom_close( inum )
903
904	ENDIF
905	!
906	END SUBROUTINE asm_inc_init
907
908
909	SUBROUTINE calc_date( kit000, kt, kdate0, kdate )
910	!!----------------------------------------------------------------------
911	!! * ROUTINE calc_date *
912	!!
913	!! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.
914	!!
915	!! ** Method : Compute the calendar date YYYYMMDD at a given time step.
916	!!
917	!! ** Action :
918	!!----------------------------------------------------------------------
919	INTEGER, INTENT(IN) :: kit000 ! Initial time step
920	INTEGER, INTENT(IN) :: kt ! Current time step referenced to kit000
921	INTEGER, INTENT(IN) :: kdate0 ! Initial date
922	INTEGER, INTENT(OUT) :: kdate ! Current date reference to kdate0
923	!
924	INTEGER :: iyea0 ! Initial year
925	INTEGER :: imon0 ! Initial month
926	INTEGER :: iday0 ! Initial day
927	INTEGER :: iyea ! Current year
928	INTEGER :: imon ! Current month
929	INTEGER :: iday ! Current day
930	INTEGER :: idaystp ! Number of days between initial and current date
931	INTEGER :: idaycnt ! Day counter
932
933	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
934
935	!-----------------------------------------------------------------------
936	! Compute the calendar date YYYYMMDD
937	!-----------------------------------------------------------------------
938
939	! Initial date
940	iyea0 = kdate0 / 10000
941	imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100
942	iday0 = kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 )
943
944	! Check that kt >= kit000 - 1
945	IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')
946
947	! If kt = kit000 - 1 then set the date to the restart date
948	IF ( kt == kit000 - 1 ) THEN
949
950	kdate = ndastp
951	RETURN
952
953	ENDIF
954
955	! Compute the number of days from the initial date
956	idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )
957
958	iday = iday0
959	imon = imon0
960	iyea = iyea0
961	idaycnt = 0
962
963	CALL calc_month_len( iyea, imonth_len )
964
965	DO WHILE ( idaycnt < idaystp )
966	iday = iday + 1
967	IF ( iday > imonth_len(imon) ) THEN
968	iday = 1
969	imon = imon + 1
970	ENDIF
971	IF ( imon > 12 ) THEN
972	imon = 1
973	iyea = iyea + 1
974	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
975	ENDIF
976	idaycnt = idaycnt + 1
977	END DO
978	!
979	kdate = iyea * 10000 + imon * 100 + iday
980	!
981	END SUBROUTINE
982
983
984	SUBROUTINE calc_month_len( iyear, imonth_len )
985	!!----------------------------------------------------------------------
986	!! * ROUTINE calc_month_len *
987	!!
988	!! ** Purpose : Compute the number of days in a months given a year.
989	!!
990	!! ** Method :
991	!!----------------------------------------------------------------------
992	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
993	INTEGER :: iyear !: year
994	!!----------------------------------------------------------------------
995	!
996	! length of the month of the current year (from nleapy, read in namelist)
997	IF ( nleapy < 2 ) THEN
998	imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
999	IF ( nleapy == 1 ) THEN ! we are using calendar with leap years
1000	IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
1001	imonth_len(2) = 29
1002	ENDIF
1003	ENDIF
1004	ELSE
1005	imonth_len(:) = nleapy ! all months with nleapy days per year
1006	ENDIF
1007	!
1008	END SUBROUTINE
1009
1010
1011	SUBROUTINE tra_asm_inc( kt )
1012	!!----------------------------------------------------------------------
1013	!! * ROUTINE tra_asm_inc *
1014	!!
1015	!! ** Purpose : Apply the tracer (T and S) assimilation increments
1016	!!
1017	!! ** Method : Direct initialization or Incremental Analysis Updating
1018	!!
1019	!! ** Action :
1020	!!----------------------------------------------------------------------
1021	INTEGER, INTENT(IN) :: kt ! Current time step
1022	!
1023	INTEGER :: ji,jj,jk
1024	INTEGER :: it
1025	REAL(wp) :: zincwgt ! IAU weight for current time step
1026	REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
1027	!!----------------------------------------------------------------------
1028
1029	! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
1030	! used to prevent the applied increments taking the temperature below the local freezing point
1031
1032	DO jk = 1, jpkm1
1033	fzptnz(:,:,jk) = eos_fzp( tsn(:,:,jk,jp_sal), fsdept(:,:,jk) )
1034	END DO
1035
1036	IF ( ln_asmiau ) THEN
1037
1038	!--------------------------------------------------------------------
1039	! Incremental Analysis Updating
1040	!--------------------------------------------------------------------
1041
1042	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1043
1044	it = kt - nit000 + 1
1045	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
1046
1047	IF(lwp) THEN
1048	WRITE(numout,*)
1049	WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
1050	WRITE(numout,*) '~~~~~~~~~~~~'
1051	ENDIF
1052
1053	! Update the tracer tendencies
1054	DO jk = 1, jpkm1
1055	IF (ln_temnofreeze) THEN
1056	! Do not apply negative increments if the temperature will fall below freezing
1057	WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
1058	& tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )
1059	tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt
1060	END WHERE
1061	ELSE
1062	tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt
1063	ENDIF
1064	IF (ln_salfix) THEN
1065	! Do not apply negative increments if the salinity will fall below a specified
1066	! minimum value salfixmin
1067	WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
1068	& tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )
1069	tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
1070	END WHERE
1071	ELSE
1072	tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
1073	ENDIF
1074	END DO
1075
1076	ENDIF
1077
1078	IF ( kt == nitiaufin_r + 1 ) THEN ! For bias crcn to work
1079	DEALLOCATE( t_bkginc )
1080	DEALLOCATE( s_bkginc )
1081	ENDIF
1082
1083
1084	ELSEIF ( ln_asmdin ) THEN
1085
1086	!--------------------------------------------------------------------
1087	! Direct Initialization
1088	!--------------------------------------------------------------------
1089
1090	IF ( kt == nitdin_r ) THEN
1091
1092	neuler = 0 ! Force Euler forward step
1093
1094	! Initialize the now fields with the background + increment
1095	IF (ln_temnofreeze) THEN
1096	! Do not apply negative increments if the temperature will fall below freezing
1097	WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )
1098	tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)
1099	END WHERE
1100	ELSE
1101	tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)
1102	ENDIF
1103	IF (ln_salfix) THEN
1104	! Do not apply negative increments if the salinity will fall below a specified
1105	! minimum value salfixmin
1106	WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )
1107	tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)
1108	END WHERE
1109	ELSE
1110	tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)
1111	ENDIF
1112
1113	tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields
1114
1115	CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities
1116	!!gm fabien
1117	! CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities
1118	!!gm
1119
1120
1121	IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) &
1122	& CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient
1123	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
1124	IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) &
1125	& CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF)
1126	& rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , &
1127	& gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level
1128
1129	#if defined key_zdfkpp
1130	CALL eos( tsn, rhd, fsdept_n(:,:,:) ) ! Compute rhd
1131	!!gm fabien CALL eos( tsn, rhd ) ! Compute rhd
1132	#endif
1133
1134	DEALLOCATE( t_bkginc )
1135	DEALLOCATE( s_bkginc )
1136	DEALLOCATE( t_bkg )
1137	DEALLOCATE( s_bkg )
1138	ENDIF
1139	!
1140	ENDIF
1141	! Perhaps the following call should be in step
1142	IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment
1143	!
1144	END SUBROUTINE tra_asm_inc
1145
1146
1147	SUBROUTINE dyn_asm_inc( kt )
1148	!!----------------------------------------------------------------------
1149	!! * ROUTINE dyn_asm_inc *
1150	!!
1151	!! ** Purpose : Apply the dynamics (u and v) assimilation increments.
1152	!!
1153	!! ** Method : Direct initialization or Incremental Analysis Updating.
1154	!!
1155	!! ** Action :
1156	!!----------------------------------------------------------------------
1157	INTEGER, INTENT(IN) :: kt ! Current time step
1158	!
1159	INTEGER :: jk
1160	INTEGER :: it
1161	REAL(wp) :: zincwgt ! IAU weight for current time step
1162	!!----------------------------------------------------------------------
1163
1164	IF ( ln_asmiau ) THEN
1165
1166	!--------------------------------------------------------------------
1167	! Incremental Analysis Updating
1168	!--------------------------------------------------------------------
1169
1170	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1171
1172	it = kt - nit000 + 1
1173	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
1174
1175	IF(lwp) THEN
1176	WRITE(numout,*)
1177	WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', &
1178	& kt,' with IAU weight = ', wgtiau(it)
1179	WRITE(numout,*) '~~~~~~~~~~~~'
1180	ENDIF
1181
1182	! Update the dynamic tendencies
1183	DO jk = 1, jpkm1
1184	ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
1185	va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
1186	END DO
1187
1188	IF ( kt == nitiaufin_r ) THEN
1189	DEALLOCATE( u_bkginc )
1190	DEALLOCATE( v_bkginc )
1191	ENDIF
1192
1193	ENDIF
1194
1195	ELSEIF ( ln_asmdin ) THEN
1196
1197	!--------------------------------------------------------------------
1198	! Direct Initialization
1199	!--------------------------------------------------------------------
1200
1201	IF ( kt == nitdin_r ) THEN
1202
1203	neuler = 0 ! Force Euler forward step
1204
1205	! Initialize the now fields with the background + increment
1206	un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
1207	vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)
1208
1209	ub(:,:,:) = un(:,:,:) ! Update before fields
1210	vb(:,:,:) = vn(:,:,:)
1211
1212	DEALLOCATE( u_bkg )
1213	DEALLOCATE( v_bkg )
1214	DEALLOCATE( u_bkginc )
1215	DEALLOCATE( v_bkginc )
1216	ENDIF
1217	!
1218	ENDIF
1219	!
1220	END SUBROUTINE dyn_asm_inc
1221
1222
1223	SUBROUTINE ssh_asm_inc( kt )
1224	!!----------------------------------------------------------------------
1225	!! * ROUTINE ssh_asm_inc *
1226	!!
1227	!! ** Purpose : Apply the sea surface height assimilation increment.
1228	!!
1229	!! ** Method : Direct initialization or Incremental Analysis Updating.
1230	!!
1231	!! ** Action :
1232	!!----------------------------------------------------------------------
1233	INTEGER, INTENT(IN) :: kt ! Current time step
1234	!
1235	INTEGER :: it
1236	INTEGER :: jk
1237	REAL(wp) :: zincwgt ! IAU weight for current time step
1238	!!----------------------------------------------------------------------
1239
1240	IF ( ln_asmiau ) THEN
1241
1242	!--------------------------------------------------------------------
1243	! Incremental Analysis Updating
1244	!--------------------------------------------------------------------
1245
1246	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1247
1248	it = kt - nit000 + 1
1249	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
1250
1251	IF(lwp) THEN
1252	WRITE(numout,*)
1253	WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
1254	& kt,' with IAU weight = ', wgtiau(it)
1255	WRITE(numout,*) '~~~~~~~~~~~~'
1256	ENDIF
1257
1258	! Save the tendency associated with the IAU weighted SSH increment
1259	! (applied in dynspg.*)
1260	#if defined key_asminc
1261	ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
1262	#endif
1263	IF ( kt == nitiaufin_r ) THEN
1264	DEALLOCATE( ssh_bkginc )
1265	ENDIF
1266
1267	ENDIF
1268
1269	ELSEIF ( ln_asmdin ) THEN
1270
1271	!--------------------------------------------------------------------
1272	! Direct Initialization
1273	!--------------------------------------------------------------------
1274
1275	IF ( kt == nitdin_r ) THEN
1276
1277	neuler = 0 ! Force Euler forward step
1278
1279	! Initialize the now fields the background + increment
1280	sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)
1281
1282	! Update before fields
1283	sshb(:,:) = sshn(:,:)
1284
1285	IF( lk_vvl ) THEN
1286	DO jk = 1, jpk
1287	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
1288	END DO
1289	ENDIF
1290
1291	DEALLOCATE( ssh_bkg )
1292	DEALLOCATE( ssh_bkginc )
1293
1294	ENDIF
1295	!
1296	ENDIF
1297	!
1298	END SUBROUTINE ssh_asm_inc
1299
1300
1301	SUBROUTINE seaice_asm_inc( kt, kindic )
1302	!!----------------------------------------------------------------------
1303	!! * ROUTINE seaice_asm_inc *
1304	!!
1305	!! ** Purpose : Apply the sea ice assimilation increment.
1306	!!
1307	!! ** Method : Direct initialization or Incremental Analysis Updating.
1308	!!
1309	!! ** Action :
1310	!!
1311	!!----------------------------------------------------------------------
1312	IMPLICIT NONE
1313	!
1314	INTEGER, INTENT(in) :: kt ! Current time step
1315	INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation
1316	!
1317	INTEGER :: it
1318	REAL(wp) :: zincwgt ! IAU weight for current time step
1319	#if defined key_lim2
1320	REAL(wp), DIMENSION(jpi,jpj) :: zofrld, zohicif, zseaicendg, zhicifinc ! LIM
1321	REAL(wp) :: zhicifmin = 0.5_wp ! ice minimum depth in metres
1322	#endif
1323	!!----------------------------------------------------------------------
1324
1325	IF ( ln_asmiau ) THEN
1326
1327	!--------------------------------------------------------------------
1328	! Incremental Analysis Updating
1329	!--------------------------------------------------------------------
1330
1331	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1332
1333	it = kt - nit000 + 1
1334	zincwgt = wgtiau(it) ! IAU weight for the current time step
1335	! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
1336
1337	IF(lwp) THEN
1338	WRITE(numout,*)
1339	WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', &
1340	& kt,' with IAU weight = ', wgtiau(it)
1341	WRITE(numout,*) '~~~~~~~~~~~~'
1342	ENDIF
1343
1344	! Sea-ice : LIM-3 case (to add)
1345
1346	#if defined key_lim2
1347	! Sea-ice : LIM-2 case
1348	zofrld (:,:) = frld(:,:)
1349	zohicif(:,:) = hicif(:,:)
1350	!
1351	frld = MIN( MAX( frld (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
1352	pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
1353	fr_i(:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction
1354	!
1355	zseaicendg(:,:) = zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied
1356	!
1357	! Nudge sea ice depth to bring it up to a required minimum depth
1358	WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )
1359	zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt
1360	ELSEWHERE
1361	zhicifinc(:,:) = 0.0_wp
1362	END WHERE
1363	!
1364	! nudge ice depth
1365	hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1366	phicif(:,:) = phicif(:,:) + zhicifinc(:,:)
1367	!
1368	! seaice salinity balancing (to add)
1369	#endif
1370
1371	#if defined key_cice && defined key_asminc
1372	! Sea-ice : CICE case. Pass ice increment tendency into CICE
1373	ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
1374	#endif
1375
1376	IF ( kt == nitiaufin_r ) THEN
1377	DEALLOCATE( seaice_bkginc )
1378	ENDIF
1379
1380	ELSE
1381
1382	#if defined key_cice && defined key_asminc
1383	! Sea-ice : CICE case. Zero ice increment tendency into CICE
1384	ndaice_da(:,:) = 0.0_wp
1385	#endif
1386
1387	ENDIF
1388
1389	ELSEIF ( ln_asmdin ) THEN
1390
1391	!--------------------------------------------------------------------
1392	! Direct Initialization
1393	!--------------------------------------------------------------------
1394
1395	IF ( kt == nitdin_r ) THEN
1396
1397	neuler = 0 ! Force Euler forward step
1398
1399	! Sea-ice : LIM-3 case (to add)
1400
1401	#if defined key_lim2
1402	! Sea-ice : LIM-2 case.
1403	zofrld(:,:)=frld(:,:)
1404	zohicif(:,:)=hicif(:,:)
1405	!
1406	! Initialize the now fields the background + increment
1407	frld (:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
1408	pfrld(:,:) = frld(:,:)
1409	fr_i (:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction
1410	zseaicendg(:,:) = zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied
1411	!
1412	! Nudge sea ice depth to bring it up to a required minimum depth
1413	WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )
1414	zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt
1415	ELSEWHERE
1416	zhicifinc(:,:) = 0.0_wp
1417	END WHERE
1418	!
1419	! nudge ice depth
1420	hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1421	phicif(:,:) = phicif(:,:)
1422	!
1423	! seaice salinity balancing (to add)
1424	#endif
1425
1426	#if defined key_cice && defined key_asminc
1427	! Sea-ice : CICE case. Pass ice increment tendency into CICE
1428	ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
1429	#endif
1430	IF ( .NOT. PRESENT(kindic) ) THEN
1431	DEALLOCATE( seaice_bkginc )
1432	END IF
1433
1434	ELSE
1435
1436	#if defined key_cice && defined key_asminc
1437	! Sea-ice : CICE case. Zero ice increment tendency into CICE
1438	ndaice_da(:,:) = 0.0_wp
1439	#endif
1440
1441	ENDIF
1442
1443	!#if defined defined key_lim2 \|\| defined key_cice
1444	!
1445	! IF (ln_seaicebal ) THEN
1446	! !! balancing salinity increments
1447	! !! simple case from limflx.F90 (doesn't include a mass flux)
1448	! !! assumption is that as ice concentration is reduced or increased
1449	! !! the snow and ice depths remain constant
1450	! !! note that snow is being created where ice concentration is being increased
1451	! !! - could be more sophisticated and
1452	! !! not do this (but would need to alter h_snow)
1453	!
1454	! usave(:,:,:)=sb(:,:,:) ! use array as a temporary store
1455	!
1456	! DO jj = 1, jpj
1457	! DO ji = 1, jpi
1458	! ! calculate change in ice and snow mass per unit area
1459	! ! positive values imply adding salt to the ocean (results from ice formation)
1460	! ! fwf : ice formation and melting
1461	!
1462	! zfons = ( -nfresh_da(ji,jj)soce + nfsalt_da(ji,jj) )rdt
1463	!
1464	! ! change salinity down to mixed layer depth
1465	! mld=hmld_kara(ji,jj)
1466	!
1467	! ! prevent small mld
1468	! ! less than 10m can cause salinity instability
1469	! IF (mld < 10) mld=10
1470	!
1471	! ! set to bottom of a level
1472	! DO jk = jpk-1, 2, -1
1473	! IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN
1474	! mld=gdepw(ji,jj,jk+1)
1475	! jkmax=jk
1476	! ENDIF
1477	! ENDDO
1478	!
1479	! ! avoid applying salinity balancing in shallow water or on land
1480	! !
1481	!
1482	! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
1483	!
1484	! dsal_ocn=0.0_wp
1485	! sal_thresh=5.0_wp ! minimum salinity threshold for salinity balancing
1486	!
1487	! if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
1488	! dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
1489	!
1490	! ! put increments in for levels in the mixed layer
1491	! ! but prevent salinity below a threshold value
1492	!
1493	! DO jk = 1, jkmax
1494	!
1495	! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
1496	! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
1497	! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
1498	! ENDIF
1499	!
1500	! ENDDO
1501	!
1502	! ! ! salt exchanges at the ice/ocean interface
1503	! ! zpmess = zfons / rdt_ice ! rdt_ice is ice timestep
1504	! !
1505	! !! Adjust fsalt. A +ve fsalt means adding salt to ocean
1506	! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt
1507	! !!
1508	! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d)
1509	! !! ! E-P (kg m-2 s-2)
1510	! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2)
1511	! ENDDO !ji
1512	! ENDDO !jj!
1513	!
1514	! ENDIF !ln_seaicebal
1515	!
1516	!#endif
1517
1518	ENDIF
1519
1520	END SUBROUTINE seaice_asm_inc
1521
1522	SUBROUTINE logchl_asm_inc( kt )
1523	!!----------------------------------------------------------------------
1524	!! * ROUTINE logchl_asm_inc *
1525	!!
1526	!! ** Purpose : Apply the chlorophyll assimilation increments.
1527	!!
1528	!! ** Method : Calculate increments to state variables using nitrogen
1529	!! balancing.
1530	!! Direct initialization or Incremental Analysis Updating.
1531	!!
1532	!! ** Action :
1533	!!----------------------------------------------------------------------
1534	INTEGER, INTENT(IN) :: kt ! Current time step
1535	!
1536	INTEGER :: jk ! Loop counter
1537	INTEGER :: it ! Index
1538	REAL(wp) :: zincwgt ! IAU weight for current time step
1539	REAL(wp) :: zincper ! IAU interval in seconds
1540	!!----------------------------------------------------------------------
1541
1542	IF ( kt <= nit000 ) THEN
1543
1544	zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
1545
1546	#if defined key_fabm
1547	CALL asm_logchl_bal_ersem( ln_logchlpftinc, nn_asmpfts, &
1548	& mld_choice_bgc, logchl_bkginc, &
1549	& logchl_balinc_ersem_chl1, logchl_balinc_ersem_chl2, &
1550	& logchl_balinc_ersem_chl3, logchl_balinc_ersem_chl4, &
1551	& logchl_balinc_ersem_p1c, logchl_balinc_ersem_p1n, &
1552	& logchl_balinc_ersem_p1p, logchl_balinc_ersem_p1s, &
1553	& logchl_balinc_ersem_p2c, logchl_balinc_ersem_p2n, &
1554	& logchl_balinc_ersem_p2p, logchl_balinc_ersem_p3c, &
1555	& logchl_balinc_ersem_p3n, logchl_balinc_ersem_p3p, &
1556	& logchl_balinc_ersem_p4c, logchl_balinc_ersem_p4n, &
1557	& logchl_balinc_ersem_p4p, logchl_balinc_ersem_z4c, &
1558	& logchl_balinc_ersem_z5c, logchl_balinc_ersem_z5n, &
1559	& logchl_balinc_ersem_z5p, logchl_balinc_ersem_z6c, &
1560	& logchl_balinc_ersem_z6n, logchl_balinc_ersem_z6p, &
1561	& logchl_balinc_ersem_n1p, logchl_balinc_ersem_n3n, &
1562	& logchl_balinc_ersem_n4n, logchl_balinc_ersem_n5s )
1563	#elif defined key_medusa && defined key_foam_medusa
1564	!CALL asm_logchl_bal_medusa()
1565	CALL ctl_stop( 'Attempting to assimilate logchl into MEDUSA, ', &
1566	& 'but not fully implemented yet' )
1567	#elif defined key_hadocc
1568	!CALL asm_logchl_bal_hadocc()
1569	CALL ctl_stop( 'Attempting to assimilate logchl into HadOCC, ', &
1570	& 'but not fully implemented yet' )
1571	#else
1572	CALL ctl_stop( 'Attempting to assimilate logchl, ', &
1573	& 'but not defined a biogeochemical model' )
1574	#endif
1575
1576	ENDIF
1577
1578	IF ( ln_asmiau ) THEN
1579
1580	!--------------------------------------------------------------------
1581	! Incremental Analysis Updating
1582	!--------------------------------------------------------------------
1583
1584	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1585
1586	it = kt - nit000 + 1
1587	zincwgt = wgtiau(it) ! IAU weight for the current time step
1588	! note this is not a tendency so should not be divided by rdt
1589
1590	IF(lwp) THEN
1591	WRITE(numout,*)
1592	WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', &
1593	& kt,' with IAU weight = ', wgtiau(it)
1594	WRITE(numout,*) '~~~~~~~~~~~~'
1595	ENDIF
1596
1597	! Update the biogeochemical tendencies
1598	#if defined key_fabm
1599	DO jk = 1, jpkm1
1600	! Add directly to trn and trb, rather than to tra, as not a tendency
1601	trn(:,:,jk,jp_fabm_m1+jp_fabm_chl1) = trn(:,:,jk,jp_fabm_m1+jp_fabm_chl1) + &
1602	& logchl_balinc_ersem_chl1(:,:,jk) * zincwgt
1603	trn(:,:,jk,jp_fabm_m1+jp_fabm_chl2) = trn(:,:,jk,jp_fabm_m1+jp_fabm_chl2) + &
1604	& logchl_balinc_ersem_chl2(:,:,jk) * zincwgt
1605	trn(:,:,jk,jp_fabm_m1+jp_fabm_chl3) = trn(:,:,jk,jp_fabm_m1+jp_fabm_chl3) + &
1606	& logchl_balinc_ersem_chl3(:,:,jk) * zincwgt
1607	trn(:,:,jk,jp_fabm_m1+jp_fabm_chl4) = trn(:,:,jk,jp_fabm_m1+jp_fabm_chl4) + &
1608	& logchl_balinc_ersem_chl4(:,:,jk) * zincwgt
1609	trn(:,:,jk,jp_fabm_m1+jp_fabm_p1c) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p1c) + &
1610	& logchl_balinc_ersem_p1c(:,:,jk) * zincwgt
1611	trn(:,:,jk,jp_fabm_m1+jp_fabm_p1n) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p1n) + &
1612	& logchl_balinc_ersem_p1n(:,:,jk) * zincwgt
1613	trn(:,:,jk,jp_fabm_m1+jp_fabm_p1p) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p1p) + &
1614	& logchl_balinc_ersem_p1p(:,:,jk) * zincwgt
1615	trn(:,:,jk,jp_fabm_m1+jp_fabm_p1s) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p1s) + &
1616	& logchl_balinc_ersem_p1s(:,:,jk) * zincwgt
1617	trn(:,:,jk,jp_fabm_m1+jp_fabm_p2c) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p2c) + &
1618	& logchl_balinc_ersem_p2c(:,:,jk) * zincwgt
1619	trn(:,:,jk,jp_fabm_m1+jp_fabm_p2n) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p2n) + &
1620	& logchl_balinc_ersem_p2n(:,:,jk) * zincwgt
1621	trn(:,:,jk,jp_fabm_m1+jp_fabm_p2p) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p2p) + &
1622	& logchl_balinc_ersem_p2p(:,:,jk) * zincwgt
1623	trn(:,:,jk,jp_fabm_m1+jp_fabm_p3c) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p3c) + &
1624	& logchl_balinc_ersem_p3c(:,:,jk) * zincwgt
1625	trn(:,:,jk,jp_fabm_m1+jp_fabm_p3n) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p3n) + &
1626	& logchl_balinc_ersem_p3n(:,:,jk) * zincwgt
1627	trn(:,:,jk,jp_fabm_m1+jp_fabm_p3p) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p3p) + &
1628	& logchl_balinc_ersem_p3p(:,:,jk) * zincwgt
1629	trn(:,:,jk,jp_fabm_m1+jp_fabm_p4c) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p4c) + &
1630	& logchl_balinc_ersem_p4c(:,:,jk) * zincwgt
1631	trn(:,:,jk,jp_fabm_m1+jp_fabm_p4n) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p4n) + &
1632	& logchl_balinc_ersem_p4n(:,:,jk) * zincwgt
1633	trn(:,:,jk,jp_fabm_m1+jp_fabm_p4p) = trn(:,:,jk,jp_fabm_m1+jp_fabm_p4p) + &
1634	& logchl_balinc_ersem_p4p(:,:,jk) * zincwgt
1635	trn(:,:,jk,jp_fabm_m1+jp_fabm_z4c) = trn(:,:,jk,jp_fabm_m1+jp_fabm_z4c) + &
1636	& logchl_balinc_ersem_z4c(:,:,jk) * zincwgt
1637	trn(:,:,jk,jp_fabm_m1+jp_fabm_z5c) = trn(:,:,jk,jp_fabm_m1+jp_fabm_z5c) + &
1638	& logchl_balinc_ersem_z5c(:,:,jk) * zincwgt
1639	trn(:,:,jk,jp_fabm_m1+jp_fabm_z5n) = trn(:,:,jk,jp_fabm_m1+jp_fabm_z5n) + &
1640	& logchl_balinc_ersem_z5n(:,:,jk) * zincwgt
1641	trn(:,:,jk,jp_fabm_m1+jp_fabm_z5p) = trn(:,:,jk,jp_fabm_m1+jp_fabm_z5p) + &
1642	& logchl_balinc_ersem_z5p(:,:,jk) * zincwgt
1643	trn(:,:,jk,jp_fabm_m1+jp_fabm_z6c) = trn(:,:,jk,jp_fabm_m1+jp_fabm_z6c) + &
1644	& logchl_balinc_ersem_z6c(:,:,jk) * zincwgt
1645	trn(:,:,jk,jp_fabm_m1+jp_fabm_z6n) = trn(:,:,jk,jp_fabm_m1+jp_fabm_z6n) + &
1646	& logchl_balinc_ersem_z6n(:,:,jk) * zincwgt
1647	trn(:,:,jk,jp_fabm_m1+jp_fabm_z6p) = trn(:,:,jk,jp_fabm_m1+jp_fabm_z6p) + &
1648	& logchl_balinc_ersem_z6p(:,:,jk) * zincwgt
1649	trn(:,:,jk,jp_fabm_m1+jp_fabm_n1p) = trn(:,:,jk,jp_fabm_m1+jp_fabm_n1p) + &
1650	& logchl_balinc_ersem_n1p(:,:,jk) * zincwgt
1651	trn(:,:,jk,jp_fabm_m1+jp_fabm_n3n) = trn(:,:,jk,jp_fabm_m1+jp_fabm_n3n) + &
1652	& logchl_balinc_ersem_n3n(:,:,jk) * zincwgt
1653	trn(:,:,jk,jp_fabm_m1+jp_fabm_n4n) = trn(:,:,jk,jp_fabm_m1+jp_fabm_n4n) + &
1654	& logchl_balinc_ersem_n4n(:,:,jk) * zincwgt
1655	trn(:,:,jk,jp_fabm_m1+jp_fabm_n5s) = trn(:,:,jk,jp_fabm_m1+jp_fabm_n5s) + &
1656	& logchl_balinc_ersem_n5s(:,:,jk) * zincwgt
1657
1658	trb(:,:,jk,jp_fabm_m1+jp_fabm_chl1) = trb(:,:,jk,jp_fabm_m1+jp_fabm_chl1) + &
1659	& logchl_balinc_ersem_chl1(:,:,jk) * zincwgt
1660	trb(:,:,jk,jp_fabm_m1+jp_fabm_chl2) = trb(:,:,jk,jp_fabm_m1+jp_fabm_chl2) + &
1661	& logchl_balinc_ersem_chl2(:,:,jk) * zincwgt
1662	trb(:,:,jk,jp_fabm_m1+jp_fabm_chl3) = trb(:,:,jk,jp_fabm_m1+jp_fabm_chl3) + &
1663	& logchl_balinc_ersem_chl3(:,:,jk) * zincwgt
1664	trb(:,:,jk,jp_fabm_m1+jp_fabm_chl4) = trb(:,:,jk,jp_fabm_m1+jp_fabm_chl4) + &
1665	& logchl_balinc_ersem_chl4(:,:,jk) * zincwgt
1666	trb(:,:,jk,jp_fabm_m1+jp_fabm_p1c) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p1c) + &
1667	& logchl_balinc_ersem_p1c(:,:,jk) * zincwgt
1668	trb(:,:,jk,jp_fabm_m1+jp_fabm_p1n) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p1n) + &
1669	& logchl_balinc_ersem_p1n(:,:,jk) * zincwgt
1670	trb(:,:,jk,jp_fabm_m1+jp_fabm_p1p) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p1p) + &
1671	& logchl_balinc_ersem_p1p(:,:,jk) * zincwgt
1672	trb(:,:,jk,jp_fabm_m1+jp_fabm_p1s) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p1s) + &
1673	& logchl_balinc_ersem_p1s(:,:,jk) * zincwgt
1674	trb(:,:,jk,jp_fabm_m1+jp_fabm_p2c) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p2c) + &
1675	& logchl_balinc_ersem_p2c(:,:,jk) * zincwgt
1676	trb(:,:,jk,jp_fabm_m1+jp_fabm_p2n) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p2n) + &
1677	& logchl_balinc_ersem_p2n(:,:,jk) * zincwgt
1678	trb(:,:,jk,jp_fabm_m1+jp_fabm_p2p) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p2p) + &
1679	& logchl_balinc_ersem_p2p(:,:,jk) * zincwgt
1680	trb(:,:,jk,jp_fabm_m1+jp_fabm_p3c) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p3c) + &
1681	& logchl_balinc_ersem_p3c(:,:,jk) * zincwgt
1682	trb(:,:,jk,jp_fabm_m1+jp_fabm_p3n) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p3n) + &
1683	& logchl_balinc_ersem_p3n(:,:,jk) * zincwgt
1684	trb(:,:,jk,jp_fabm_m1+jp_fabm_p3p) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p3p) + &
1685	& logchl_balinc_ersem_p3p(:,:,jk) * zincwgt
1686	trb(:,:,jk,jp_fabm_m1+jp_fabm_p4c) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p4c) + &
1687	& logchl_balinc_ersem_p4c(:,:,jk) * zincwgt
1688	trb(:,:,jk,jp_fabm_m1+jp_fabm_p4n) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p4n) + &
1689	& logchl_balinc_ersem_p4n(:,:,jk) * zincwgt
1690	trb(:,:,jk,jp_fabm_m1+jp_fabm_p4p) = trb(:,:,jk,jp_fabm_m1+jp_fabm_p4p) + &
1691	& logchl_balinc_ersem_p4p(:,:,jk) * zincwgt
1692	trb(:,:,jk,jp_fabm_m1+jp_fabm_z4c) = trb(:,:,jk,jp_fabm_m1+jp_fabm_z4c) + &
1693	& logchl_balinc_ersem_z4c(:,:,jk) * zincwgt
1694	trb(:,:,jk,jp_fabm_m1+jp_fabm_z5c) = trb(:,:,jk,jp_fabm_m1+jp_fabm_z5c) + &
1695	& logchl_balinc_ersem_z5c(:,:,jk) * zincwgt
1696	trb(:,:,jk,jp_fabm_m1+jp_fabm_z5n) = trb(:,:,jk,jp_fabm_m1+jp_fabm_z5n) + &
1697	& logchl_balinc_ersem_z5n(:,:,jk) * zincwgt
1698	trb(:,:,jk,jp_fabm_m1+jp_fabm_z5p) = trb(:,:,jk,jp_fabm_m1+jp_fabm_z5p) + &
1699	& logchl_balinc_ersem_z5p(:,:,jk) * zincwgt
1700	trb(:,:,jk,jp_fabm_m1+jp_fabm_z6c) = trb(:,:,jk,jp_fabm_m1+jp_fabm_z6c) + &
1701	& logchl_balinc_ersem_z6c(:,:,jk) * zincwgt
1702	trb(:,:,jk,jp_fabm_m1+jp_fabm_z6n) = trb(:,:,jk,jp_fabm_m1+jp_fabm_z6n) + &
1703	& logchl_balinc_ersem_z6n(:,:,jk) * zincwgt
1704	trb(:,:,jk,jp_fabm_m1+jp_fabm_z6p) = trb(:,:,jk,jp_fabm_m1+jp_fabm_z6p) + &
1705	& logchl_balinc_ersem_z6p(:,:,jk) * zincwgt
1706	trb(:,:,jk,jp_fabm_m1+jp_fabm_n1p) = trb(:,:,jk,jp_fabm_m1+jp_fabm_n1p) + &
1707	& logchl_balinc_ersem_n1p(:,:,jk) * zincwgt
1708	trb(:,:,jk,jp_fabm_m1+jp_fabm_n3n) = trb(:,:,jk,jp_fabm_m1+jp_fabm_n3n) + &
1709	& logchl_balinc_ersem_n3n(:,:,jk) * zincwgt
1710	trb(:,:,jk,jp_fabm_m1+jp_fabm_n4n) = trb(:,:,jk,jp_fabm_m1+jp_fabm_n4n) + &
1711	& logchl_balinc_ersem_n4n(:,:,jk) * zincwgt
1712	trb(:,:,jk,jp_fabm_m1+jp_fabm_n5s) = trb(:,:,jk,jp_fabm_m1+jp_fabm_n5s) + &
1713	& logchl_balinc_ersem_n5s(:,:,jk) * zincwgt
1714	END DO
1715
1716	#elif defined key_medusa && defined key_foam_medusa
1717	DO jk = 1, jpkm1
1718	! Add directly to trn and trb, rather than to tra, as not a tendency
1719	trn(:,:,jk,jpchn) = trn(:,:,jk,jpchn) + logchl_balinc_medusa_chn(:,:,jk) * zincwgt
1720	trn(:,:,jk,jpchd) = trn(:,:,jk,jpchd) + logchl_balinc_medusa_chd(:,:,jk) * zincwgt
1721	trn(:,:,jk,jpphn) = trn(:,:,jk,jpphn) + logchl_balinc_medusa_phn(:,:,jk) * zincwgt
1722	trn(:,:,jk,jpphd) = trn(:,:,jk,jpphd) + logchl_balinc_medusa_phd(:,:,jk) * zincwgt
1723	trn(:,:,jk,jpzmi) = trn(:,:,jk,jpzmi) + logchl_balinc_medusa_zmi(:,:,jk) * zincwgt
1724	trn(:,:,jk,jpzme) = trn(:,:,jk,jpzme) + logchl_balinc_medusa_zme(:,:,jk) * zincwgt
1725	trn(:,:,jk,jpdin) = trn(:,:,jk,jpdin) + logchl_balinc_medusa_din(:,:,jk) * zincwgt
1726	trn(:,:,jk,jpsil) = trn(:,:,jk,jpsil) + logchl_balinc_medusa_sil(:,:,jk) * zincwgt
1727	trn(:,:,jk,jpfer) = trn(:,:,jk,jpfer) + logchl_balinc_medusa_fer(:,:,jk) * zincwgt
1728	trn(:,:,jk,jpdet) = trn(:,:,jk,jpdet) + logchl_balinc_medusa_det(:,:,jk) * zincwgt
1729	trn(:,:,jk,jppds) = trn(:,:,jk,jppds) + logchl_balinc_medusa_pds(:,:,jk) * zincwgt
1730	trn(:,:,jk,jpdtc) = trn(:,:,jk,jpdtc) + logchl_balinc_medusa_dtc(:,:,jk) * zincwgt
1731	trn(:,:,jk,jpdic) = trn(:,:,jk,jpdic) + logchl_balinc_medusa_dic(:,:,jk) * zincwgt
1732	trn(:,:,jk,jpalk) = trn(:,:,jk,jpalk) + logchl_balinc_medusa_alk(:,:,jk) * zincwgt
1733	trn(:,:,jk,jpoxy) = trn(:,:,jk,jpoxy) + logchl_balinc_medusa_oxy(:,:,jk) * zincwgt
1734
1735	trb(:,:,jk,jpchn) = trb(:,:,jk,jpchn) + logchl_balinc_medusa_chn(:,:,jk) * zincwgt
1736	trb(:,:,jk,jpchd) = trb(:,:,jk,jpchd) + logchl_balinc_medusa_chd(:,:,jk) * zincwgt
1737	trb(:,:,jk,jpphn) = trb(:,:,jk,jpphn) + logchl_balinc_medusa_phn(:,:,jk) * zincwgt
1738	trb(:,:,jk,jpphd) = trb(:,:,jk,jpphd) + logchl_balinc_medusa_phd(:,:,jk) * zincwgt
1739	trb(:,:,jk,jpzmi) = trb(:,:,jk,jpzmi) + logchl_balinc_medusa_zmi(:,:,jk) * zincwgt
1740	trb(:,:,jk,jpzme) = trb(:,:,jk,jpzme) + logchl_balinc_medusa_zme(:,:,jk) * zincwgt
1741	trb(:,:,jk,jpdin) = trb(:,:,jk,jpdin) + logchl_balinc_medusa_din(:,:,jk) * zincwgt
1742	trb(:,:,jk,jpsil) = trb(:,:,jk,jpsil) + logchl_balinc_medusa_sil(:,:,jk) * zincwgt
1743	trb(:,:,jk,jpfer) = trb(:,:,jk,jpfer) + logchl_balinc_medusa_fer(:,:,jk) * zincwgt
1744	trb(:,:,jk,jpdet) = trb(:,:,jk,jpdet) + logchl_balinc_medusa_det(:,:,jk) * zincwgt
1745	trb(:,:,jk,jppds) = trb(:,:,jk,jppds) + logchl_balinc_medusa_pds(:,:,jk) * zincwgt
1746	trb(:,:,jk,jpdtc) = trb(:,:,jk,jpdtc) + logchl_balinc_medusa_dtc(:,:,jk) * zincwgt
1747	trb(:,:,jk,jpdic) = trb(:,:,jk,jpdic) + logchl_balinc_medusa_dic(:,:,jk) * zincwgt
1748	trb(:,:,jk,jpalk) = trb(:,:,jk,jpalk) + logchl_balinc_medusa_alk(:,:,jk) * zincwgt
1749	trb(:,:,jk,jpoxy) = trb(:,:,jk,jpoxy) + logchl_balinc_medusa_oxy(:,:,jk) * zincwgt
1750	END DO
1751
1752	#elif defined key_hadocc
1753	DO jk = 1, jpkm1
1754	! Add directly to trn and trb, rather than to tra, as not a tendency
1755	trn(:,:,jk,jp_had_nut) = trn(:,:,jk,jp_had_nut) + logchl_balinc_hadocc_nut(:,:,jk) * zincwgt
1756	trn(:,:,jk,jp_had_phy) = trn(:,:,jk,jp_had_phy) + logchl_balinc_hadocc_phy(:,:,jk) * zincwgt
1757	trn(:,:,jk,jp_had_zoo) = trn(:,:,jk,jp_had_zoo) + logchl_balinc_hadocc_zoo(:,:,jk) * zincwgt
1758	trn(:,:,jk,jp_had_pdn) = trn(:,:,jk,jp_had_pdn) + logchl_balinc_hadocc_det(:,:,jk) * zincwgt
1759	trn(:,:,jk,jp_had_dic) = trn(:,:,jk,jp_had_dic) + logchl_balinc_hadocc_dic(:,:,jk) * zincwgt
1760	trn(:,:,jk,jp_had_alk) = trn(:,:,jk,jp_had_alk) + logchl_balinc_hadocc_alk(:,:,jk) * zincwgt
1761
1762	trb(:,:,jk,jp_had_nut) = trb(:,:,jk,jp_had_nut) + logchl_balinc_hadocc_nut(:,:,jk) * zincwgt
1763	trb(:,:,jk,jp_had_phy) = trb(:,:,jk,jp_had_phy) + logchl_balinc_hadocc_phy(:,:,jk) * zincwgt
1764	trb(:,:,jk,jp_had_zoo) = trb(:,:,jk,jp_had_zoo) + logchl_balinc_hadocc_zoo(:,:,jk) * zincwgt
1765	trb(:,:,jk,jp_had_pdn) = trb(:,:,jk,jp_had_pdn) + logchl_balinc_hadocc_det(:,:,jk) * zincwgt
1766	trb(:,:,jk,jp_had_dic) = trb(:,:,jk,jp_had_dic) + logchl_balinc_hadocc_dic(:,:,jk) * zincwgt
1767	trb(:,:,jk,jp_had_alk) = trb(:,:,jk,jp_had_alk) + logchl_balinc_hadocc_alk(:,:,jk) * zincwgt
1768	END DO
1769	#endif
1770
1771	! Do not deallocate arrays - needed by asm_bal_wri
1772	! which is called at end of model run
1773
1774	ENDIF
1775
1776	ELSEIF ( ln_asmdin ) THEN
1777
1778	!--------------------------------------------------------------------
1779	! Direct Initialization
1780	!--------------------------------------------------------------------
1781
1782	IF ( kt == nitdin_r ) THEN
1783
1784	neuler = 0 ! Force Euler forward step
1785
1786	#if defined key_fabm
1787	! Initialize the now fields with the background + increment
1788	! Background currently is what the model is initialised with
1789	CALL ctl_warn( ' Doing direct initialisation of ERSEM with chlorophyll assimilation', &
1790	& ' Background state is taken from model rather than background file' )
1791	trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + &
1792	& logchl_balinc_ersem_chl1(:,:,:)
1793	trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &
1794	& logchl_balinc_ersem_chl2(:,:,:)
1795	trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + &
1796	& logchl_balinc_ersem_chl3(:,:,:)
1797	trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) + &
1798	& logchl_balinc_ersem_chl4(:,:,:)
1799	trn(:,:,:,jp_fabm_m1+jp_fabm_p1c) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1c) + &
1800	& logchl_balinc_ersem_p1c(:,:,:)
1801	trn(:,:,:,jp_fabm_m1+jp_fabm_p1n) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1n) + &
1802	& logchl_balinc_ersem_p1n(:,:,:)
1803	trn(:,:,:,jp_fabm_m1+jp_fabm_p1p) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1p) + &
1804	& logchl_balinc_ersem_p1p(:,:,:)
1805	trn(:,:,:,jp_fabm_m1+jp_fabm_p1s) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1s) + &
1806	& logchl_balinc_ersem_p1s(:,:,:)
1807	trn(:,:,:,jp_fabm_m1+jp_fabm_p2c) = trn(:,:,:,jp_fabm_m1+jp_fabm_p2c) + &
1808	& logchl_balinc_ersem_p2c(:,:,:)
1809	trn(:,:,:,jp_fabm_m1+jp_fabm_p2n) = trn(:,:,:,jp_fabm_m1+jp_fabm_p2n) + &
1810	& logchl_balinc_ersem_p2n(:,:,:)
1811	trn(:,:,:,jp_fabm_m1+jp_fabm_p2p) = trn(:,:,:,jp_fabm_m1+jp_fabm_p2p) + &
1812	& logchl_balinc_ersem_p2p(:,:,:)
1813	trn(:,:,:,jp_fabm_m1+jp_fabm_p3c) = trn(:,:,:,jp_fabm_m1+jp_fabm_p3c) + &
1814	& logchl_balinc_ersem_p3c(:,:,:)
1815	trn(:,:,:,jp_fabm_m1+jp_fabm_p3n) = trn(:,:,:,jp_fabm_m1+jp_fabm_p3n) + &
1816	& logchl_balinc_ersem_p3n(:,:,:)
1817	trn(:,:,:,jp_fabm_m1+jp_fabm_p3p) = trn(:,:,:,jp_fabm_m1+jp_fabm_p3p) + &
1818	& logchl_balinc_ersem_p3p(:,:,:)
1819	trn(:,:,:,jp_fabm_m1+jp_fabm_p4c) = trn(:,:,:,jp_fabm_m1+jp_fabm_p4c) + &
1820	& logchl_balinc_ersem_p4c(:,:,:)
1821	trn(:,:,:,jp_fabm_m1+jp_fabm_p4n) = trn(:,:,:,jp_fabm_m1+jp_fabm_p4n) + &
1822	& logchl_balinc_ersem_p4n(:,:,:)
1823	trn(:,:,:,jp_fabm_m1+jp_fabm_p4p) = trn(:,:,:,jp_fabm_m1+jp_fabm_p4p) + &
1824	& logchl_balinc_ersem_p4p(:,:,:)
1825	trn(:,:,:,jp_fabm_m1+jp_fabm_z4c) = trn(:,:,:,jp_fabm_m1+jp_fabm_z4c) + &
1826	& logchl_balinc_ersem_z4c(:,:,:)
1827	trn(:,:,:,jp_fabm_m1+jp_fabm_z5c) = trn(:,:,:,jp_fabm_m1+jp_fabm_z5c) + &
1828	& logchl_balinc_ersem_z5c(:,:,:)
1829	trn(:,:,:,jp_fabm_m1+jp_fabm_z5n) = trn(:,:,:,jp_fabm_m1+jp_fabm_z5n) + &
1830	& logchl_balinc_ersem_z5n(:,:,:)
1831	trn(:,:,:,jp_fabm_m1+jp_fabm_z5p) = trn(:,:,:,jp_fabm_m1+jp_fabm_z5p) + &
1832	& logchl_balinc_ersem_z5p(:,:,:)
1833	trn(:,:,:,jp_fabm_m1+jp_fabm_z6c) = trn(:,:,:,jp_fabm_m1+jp_fabm_z6c) + &
1834	& logchl_balinc_ersem_z6c(:,:,:)
1835	trn(:,:,:,jp_fabm_m1+jp_fabm_z6n) = trn(:,:,:,jp_fabm_m1+jp_fabm_z6n) + &
1836	& logchl_balinc_ersem_z6n(:,:,:)
1837	trn(:,:,:,jp_fabm_m1+jp_fabm_z6p) = trn(:,:,:,jp_fabm_m1+jp_fabm_z6p) + &
1838	& logchl_balinc_ersem_z6p(:,:,:)
1839	trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) + &
1840	& logchl_balinc_ersem_n1p(:,:,:)
1841	trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) + &
1842	& logchl_balinc_ersem_n3n(:,:,:)
1843	trn(:,:,:,jp_fabm_m1+jp_fabm_n4n) = trn(:,:,:,jp_fabm_m1+jp_fabm_n4n) + &
1844	& logchl_balinc_ersem_n4n(:,:,:)
1845	trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) + &
1846	& logchl_balinc_ersem_n5s(:,:,:)
1847
1848	trb(:,:,:,jp_fabm_m1+jp_fabm_chl1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1)
1849	trb(:,:,:,jp_fabm_m1+jp_fabm_chl2) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl2)
1850	trb(:,:,:,jp_fabm_m1+jp_fabm_chl3) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl3)
1851	trb(:,:,:,jp_fabm_m1+jp_fabm_chl4) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
1852	trb(:,:,:,jp_fabm_m1+jp_fabm_p1c) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1c)
1853	trb(:,:,:,jp_fabm_m1+jp_fabm_p1n) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1n)
1854	trb(:,:,:,jp_fabm_m1+jp_fabm_p1p) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1p)
1855	trb(:,:,:,jp_fabm_m1+jp_fabm_p1s) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1s)
1856	trb(:,:,:,jp_fabm_m1+jp_fabm_p2c) = trn(:,:,:,jp_fabm_m1+jp_fabm_p2c)
1857	trb(:,:,:,jp_fabm_m1+jp_fabm_p2n) = trn(:,:,:,jp_fabm_m1+jp_fabm_p2n)
1858	trb(:,:,:,jp_fabm_m1+jp_fabm_p2p) = trn(:,:,:,jp_fabm_m1+jp_fabm_p2p)
1859	trb(:,:,:,jp_fabm_m1+jp_fabm_p3c) = trn(:,:,:,jp_fabm_m1+jp_fabm_p3c)
1860	trb(:,:,:,jp_fabm_m1+jp_fabm_p3n) = trn(:,:,:,jp_fabm_m1+jp_fabm_p3n)
1861	trb(:,:,:,jp_fabm_m1+jp_fabm_p3p) = trn(:,:,:,jp_fabm_m1+jp_fabm_p3p)
1862	trb(:,:,:,jp_fabm_m1+jp_fabm_p4c) = trn(:,:,:,jp_fabm_m1+jp_fabm_p4c)
1863	trb(:,:,:,jp_fabm_m1+jp_fabm_p4n) = trn(:,:,:,jp_fabm_m1+jp_fabm_p4n)
1864	trb(:,:,:,jp_fabm_m1+jp_fabm_p4p) = trn(:,:,:,jp_fabm_m1+jp_fabm_p4p)
1865	trb(:,:,:,jp_fabm_m1+jp_fabm_z4c) = trn(:,:,:,jp_fabm_m1+jp_fabm_z4c)
1866	trb(:,:,:,jp_fabm_m1+jp_fabm_z5c) = trn(:,:,:,jp_fabm_m1+jp_fabm_z5c)
1867	trb(:,:,:,jp_fabm_m1+jp_fabm_z5n) = trn(:,:,:,jp_fabm_m1+jp_fabm_z5n)
1868	trb(:,:,:,jp_fabm_m1+jp_fabm_z5p) = trn(:,:,:,jp_fabm_m1+jp_fabm_z5p)
1869	trb(:,:,:,jp_fabm_m1+jp_fabm_z6c) = trn(:,:,:,jp_fabm_m1+jp_fabm_z6c)
1870	trb(:,:,:,jp_fabm_m1+jp_fabm_z6n) = trn(:,:,:,jp_fabm_m1+jp_fabm_z6n)
1871	trb(:,:,:,jp_fabm_m1+jp_fabm_z6p) = trn(:,:,:,jp_fabm_m1+jp_fabm_z6p)
1872	trb(:,:,:,jp_fabm_m1+jp_fabm_n1p) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p)
1873	trb(:,:,:,jp_fabm_m1+jp_fabm_n3n) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n)
1874	trb(:,:,:,jp_fabm_m1+jp_fabm_n4n) = trn(:,:,:,jp_fabm_m1+jp_fabm_n4n)
1875	trb(:,:,:,jp_fabm_m1+jp_fabm_n5s) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s)
1876
1877	#elif defined key_medusa && defined key_foam_medusa
1878	! Initialize the now fields with the background + increment
1879	! Background currently is what the model is initialised with
1880	CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', &
1881	& ' Background state is taken from model rather than background file' )
1882	trn(:,:,:,jpchn) = trn(:,:,:,jpchn) + logchl_balinc_medusa_chn(:,:,:)
1883	trn(:,:,:,jpchd) = trn(:,:,:,jpchd) + logchl_balinc_medusa_chd(:,:,:)
1884	trn(:,:,:,jpphn) = trn(:,:,:,jpphn) + logchl_balinc_medusa_phn(:,:,:)
1885	trn(:,:,:,jpphd) = trn(:,:,:,jpphd) + logchl_balinc_medusa_phd(:,:,:)
1886	trn(:,:,:,jpzmi) = trn(:,:,:,jpzmi) + logchl_balinc_medusa_zmi(:,:,:)
1887	trn(:,:,:,jpzme) = trn(:,:,:,jpzme) + logchl_balinc_medusa_zme(:,:,:)
1888	trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + logchl_balinc_medusa_din(:,:,:)
1889	trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + logchl_balinc_medusa_sil(:,:,:)
1890	trn(:,:,:,jpfer) = trn(:,:,:,jpfer) + logchl_balinc_medusa_fer(:,:,:)
1891	trn(:,:,:,jpdet) = trn(:,:,:,jpdet) + logchl_balinc_medusa_det(:,:,:)
1892	trn(:,:,:,jppds) = trn(:,:,:,jppds) + logchl_balinc_medusa_pds(:,:,:)
1893	trn(:,:,:,jpdtc) = trn(:,:,:,jpdtc) + logchl_balinc_medusa_dtc(:,:,:)
1894	trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + logchl_balinc_medusa_dic(:,:,:)
1895	trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + logchl_balinc_medusa_alk(:,:,:)
1896	trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + logchl_balinc_medusa_oxy(:,:,:)
1897
1898	trb(:,:,:,jpchn) = trn(:,:,:,jpchn)
1899	trb(:,:,:,jpchd) = trn(:,:,:,jpchd)
1900	trb(:,:,:,jpphn) = trn(:,:,:,jpphn)
1901	trb(:,:,:,jpphd) = trn(:,:,:,jpphd)
1902	trb(:,:,:,jpzmi) = trn(:,:,:,jpzmi)
1903	trb(:,:,:,jpzme) = trn(:,:,:,jpzme)
1904	trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
1905	trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
1906	trb(:,:,:,jpfer) = trn(:,:,:,jpfer)
1907	trb(:,:,:,jpdet) = trn(:,:,:,jpdet)
1908	trb(:,:,:,jppds) = trn(:,:,:,jppds)
1909	trb(:,:,:,jpdtc) = trn(:,:,:,jpdtc)
1910	trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
1911	trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
1912	trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
1913
1914	#elif defined key_hadocc
1915	! Initialize the now fields with the background + increment
1916	! Background currently is what the model is initialised with
1917	CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', &
1918	& ' Background state is taken from model rather than background file' )
1919	trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + logchl_balinc_hadocc_nut(:,:,:)
1920	trn(:,:,:,jp_had_phy) = trn(:,:,:,jp_had_phy) + logchl_balinc_hadocc_phy(:,:,:)
1921	trn(:,:,:,jp_had_zoo) = trn(:,:,:,jp_had_zoo) + logchl_balinc_hadocc_zoo(:,:,:)
1922	trn(:,:,:,jp_had_pdn) = trn(:,:,:,jp_had_pdn) + logchl_balinc_hadocc_det(:,:,:)
1923	trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + logchl_balinc_hadocc_dic(:,:,:)
1924	trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + logchl_balinc_hadocc_alk(:,:,:)
1925
1926	trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut)
1927	trb(:,:,:,jp_had_phy) = trn(:,:,:,jp_had_phy)
1928	trb(:,:,:,jp_had_zoo) = trn(:,:,:,jp_had_zoo)
1929	trb(:,:,:,jp_had_pdn) = trn(:,:,:,jp_had_pdn)
1930	trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic)
1931	trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk)
1932	#endif
1933
1934	! Do not deallocate arrays - needed by asm_bal_wri
1935	! which is called at end of model run
1936	ENDIF
1937	!
1938	ENDIF
1939	!
1940	END SUBROUTINE logchl_asm_inc
1941
1942	!!======================================================================
1943	END MODULE asminc

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

Download in other formats: