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

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

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

Download in other formats: