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 NEMO/branches/UKMO/NEMO_4.0.4_FOAM_IAU_SI3_SIC/src/OCE/ASM – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_IAU_SI3_SIC/src/OCE/ASM/asminc.F90 @ 15241

Last change on this file since 15241 was 15241, checked in by dcarneir, 21 months ago
Cleaning the SIC IAU code and keeping only options used
File size: 54.0 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	!! ! 2014-09 (D. Lea) Local calc_date removed use routine from OBS
14	!! ! 2015-11 (D. Lea) Handle non-zero initial time of day
15	!!----------------------------------------------------------------------
16
17	!!----------------------------------------------------------------------
18	!! asm_inc_init : Initialize the increment arrays and IAU weights
19	!! tra_asm_inc : Apply the tracer (T and S) increments
20	!! dyn_asm_inc : Apply the dynamic (u and v) increments
21	!! ssh_asm_inc : Apply the SSH increment
22	!! ssh_asm_div : Apply divergence associated with SSH increment
23	!! seaice_asm_inc : Apply the seaice increment
24	!!----------------------------------------------------------------------
25	USE oce ! Dynamics and active tracers defined in memory
26	USE par_oce ! Ocean space and time domain variables
27	USE dom_oce ! Ocean space and time domain
28	USE domvvl ! domain: variable volume level
29	USE ldfdyn ! lateral diffusion: eddy viscosity coefficients
30	USE eosbn2 ! Equation of state - in situ and potential density
31	USE zpshde ! Partial step : Horizontal Derivative
32	USE asmpar ! Parameters for the assmilation interface
33	USE asmbkg !
34	USE c1d ! 1D initialization
35	USE sbc_oce ! Surface boundary condition variables.
36	USE diaobs , ONLY : calc_date ! Compute the calendar date on a given step
37	#if defined key_si3
38	USE phycst ! physical constants
39	USE ice1D ! sea-ice: thermodynamics variables
40	USE icetab ! sea-ice: 2D <==> 1D
41	USE icethd_do
42	USE ice
43	USE icevar , ONLY : ice_var_zapsmall
44	#endif
45	!
46	USE in_out_manager ! I/O manager
47	USE iom ! Library to read input files
48	USE lib_mpp ! MPP library
49
50	IMPLICIT NONE
51	PRIVATE
52
53	PUBLIC asm_inc_init !: Initialize the increment arrays and IAU weights
54	PUBLIC tra_asm_inc !: Apply the tracer (T and S) increments
55	PUBLIC dyn_asm_inc !: Apply the dynamic (u and v) increments
56	PUBLIC ssh_asm_inc !: Apply the SSH increment
57	PUBLIC ssh_asm_div !: Apply the SSH divergence
58	PUBLIC seaice_asm_inc !: Apply the seaice increment
59
60	#if defined key_asminc
61	LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE. !: Logical switch for assimilation increment interface
62	#else
63	LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE. !: No assimilation increments
64	#endif
65	LOGICAL, PUBLIC :: ln_bkgwri !: No output of the background state fields
66	LOGICAL, PUBLIC :: ln_asmiau !: No applying forcing with an assimilation increment
67	LOGICAL, PUBLIC :: ln_asmdin !: No direct initialization
68	LOGICAL, PUBLIC :: ln_trainc !: No tracer (T and S) assimilation increments
69	LOGICAL, PUBLIC :: ln_dyninc !: No dynamics (u and v) assimilation increments
70	LOGICAL, PUBLIC :: ln_sshinc !: No sea surface height assimilation increment
71	LOGICAL, PUBLIC :: ln_seaiceinc !: No sea ice concentration increment
72	LOGICAL, PUBLIC :: ln_salfix !: Apply minimum salinity check
73	LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
74	INTEGER, PUBLIC :: nn_divdmp !: Apply divergence damping filter nn_divdmp times
75
76	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkg , s_bkg !: Background temperature and salinity
77	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkg , v_bkg !: Background u- & v- velocity components
78	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkginc, s_bkginc !: Increment to the background T & S
79	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components
80	REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step
81	#if defined key_asminc
82	REAL(wp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: ssh_iau !: IAU-weighted sea surface height increment
83	#endif
84	! !!! time steps relative to the cycle interval [0,nitend-nit000-1]
85	INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term
86	INTEGER , PUBLIC :: nitdin !: Time step of the background state for direct initialization
87	INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval
88	INTEGER , PUBLIC :: nitiaufin !: Time step of the end of the IAU interval
89	!
90	INTEGER , PUBLIC :: niaufn !: Type of IAU weighing function: = 0 Constant weighting
91	! !: = 1 Linear hat-like, centred in middle of IAU interval
92	REAL(wp), PUBLIC :: salfixmin !: Ensure that the salinity is larger than this value if (ln_salfix)
93
94	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc ! Background sea surface height and its increment
95	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: seaice_bkginc ! Increment to the background sea ice conc
96	#if defined key_si3
97	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: a_i_bkginc ! Increment to the background sea ice conc categories
98	#endif
99	REAL(wp) :: zhi_damin = 0.5_wp !: ice thickness for new sea ice from da increment
100	INTEGER :: nhi_damin !: thickness category corresponding to zhi_damin
101	LOGICAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: incr_newice !: mask .TRUE.=DA positive ice increment to open water
102	#if defined key_cice && defined key_asminc
103	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ndaice_da ! ice increment tendency into CICE
104	#endif
105
106	!! * Substitutions
107	# include "vectopt_loop_substitute.h90"
108	!!----------------------------------------------------------------------
109	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
110	!! $Id$
111	!! Software governed by the CeCILL license (see ./LICENSE)
112	!!----------------------------------------------------------------------
113	CONTAINS
114
115	SUBROUTINE asm_inc_init
116	!!----------------------------------------------------------------------
117	!! * ROUTINE asm_inc_init *
118	!!
119	!! ** Purpose : Initialize the assimilation increment and IAU weights.
120	!!
121	!! ** Method : Initialize the assimilation increment and IAU weights.
122	!!
123	!! ** Action :
124	!!----------------------------------------------------------------------
125	INTEGER :: ji, jj, jk, jt ! dummy loop indices
126	INTEGER :: imid, inum ! local integers
127	INTEGER :: ios ! Local integer output status for namelist read
128	INTEGER :: iiauper ! Number of time steps in the IAU period
129	INTEGER :: icycper ! Number of time steps in the cycle
130	REAL(KIND=dp) :: ditend_date ! Date YYYYMMDD.HHMMSS of final time step
131	REAL(KIND=dp) :: ditbkg_date ! Date YYYYMMDD.HHMMSS of background time step for Jb term
132	REAL(KIND=dp) :: ditdin_date ! Date YYYYMMDD.HHMMSS of background time step for DI
133	REAL(KIND=dp) :: ditiaustr_date ! Date YYYYMMDD.HHMMSS of IAU interval start time step
134	REAL(KIND=dp) :: ditiaufin_date ! Date YYYYMMDD.HHMMSS of IAU interval final time step
135
136	REAL(wp) :: znorm ! Normalization factor for IAU weights
137	REAL(wp) :: ztotwgt ! Value of time-integrated IAU weights (should be equal to one)
138	REAL(wp) :: z_inc_dateb ! Start date of interval on which increment is valid
139	REAL(wp) :: z_inc_datef ! End date of interval on which increment is valid
140	REAL(wp) :: zdate_bkg ! Date in background state file for DI
141	REAL(wp) :: zdate_inc ! Time axis in increments file
142	!
143	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zhdiv ! 2D workspace
144	REAL(wp) :: zremaining_increment
145
146	#if defined key_si3
147	REAL(wp) :: zopenwater_lim
148	INTEGER :: jl
149	#endif
150	!!
151	NAMELIST/nam_asminc/ ln_bkgwri, &
152	& ln_trainc, ln_dyninc, ln_sshinc, &
153	& ln_seaiceinc, ln_asmdin, ln_asmiau, &
154	& nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, &
155	& ln_salfix, salfixmin, nn_divdmp
156	!!----------------------------------------------------------------------
157
158	!-----------------------------------------------------------------------
159	! Read Namelist nam_asminc : assimilation increment interface
160	!-----------------------------------------------------------------------
161	ln_seaiceinc = .FALSE.
162	ln_temnofreeze = .FALSE.
163	zhi_damin = 0.5_wp
164	zopenwater_lim = 1.0e-2_wp
165
166	REWIND( numnam_ref ) ! Namelist nam_asminc in reference namelist : Assimilation increment
167	READ ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
168	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in reference namelist' )
169	REWIND( numnam_cfg ) ! Namelist nam_asminc in configuration namelist : Assimilation increment
170	READ ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
171	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_asminc in configuration namelist' )
172	IF(lwm) WRITE ( numond, nam_asminc )
173
174	! Control print
175	IF(lwp) THEN
176	WRITE(numout,*)
177	WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
178	WRITE(numout,*) '~~~~~~~~~~~~'
179	WRITE(numout,*) ' Namelist namasm : set assimilation increment parameters'
180	WRITE(numout,*) ' Logical switch for writing out background state ln_bkgwri = ', ln_bkgwri
181	WRITE(numout,*) ' Logical switch for applying tracer increments ln_trainc = ', ln_trainc
182	WRITE(numout,*) ' Logical switch for applying velocity increments ln_dyninc = ', ln_dyninc
183	WRITE(numout,*) ' Logical switch for applying SSH increments ln_sshinc = ', ln_sshinc
184	WRITE(numout,*) ' Logical switch for Direct Initialization (DI) ln_asmdin = ', ln_asmdin
185	WRITE(numout,*) ' Logical switch for applying sea ice increments ln_seaiceinc = ', ln_seaiceinc
186	WRITE(numout,*) ' Logical switch for Incremental Analysis Updating (IAU) ln_asmiau = ', ln_asmiau
187	WRITE(numout,*) ' Timestep of background in [0,nitend-nit000-1] nitbkg = ', nitbkg
188	WRITE(numout,*) ' Timestep of background for DI in [0,nitend-nit000-1] nitdin = ', nitdin
189	WRITE(numout,*) ' Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
190	WRITE(numout,*) ' Timestep of end of IAU interval in [0,nitend-nit000-1] nitiaufin = ', nitiaufin
191	WRITE(numout,*) ' Type of IAU weighting function niaufn = ', niaufn
192	WRITE(numout,*) ' Logical switch for ensuring that the sa > salfixmin ln_salfix = ', ln_salfix
193	WRITE(numout,*) ' Minimum salinity after applying the increments salfixmin = ', salfixmin
194	WRITE(numout,*) ' Minimum ice thickness for new ice from da zhi_damin = ', zhi_damin
195	WRITE(numout,*) ' Limit for open water detection for ice da zopenwater_lim = ', zopenwater_lim
196	ENDIF
197
198	nitbkg_r = nitbkg + nit000 - 1 ! Background time referenced to nit000
199	nitdin_r = nitdin + nit000 - 1 ! Background time for DI referenced to nit000
200	nitiaustr_r = nitiaustr + nit000 - 1 ! Start of IAU interval referenced to nit000
201	nitiaufin_r = nitiaufin + nit000 - 1 ! End of IAU interval referenced to nit000
202
203	iiauper = nitiaufin_r - nitiaustr_r + 1 ! IAU interval length
204	icycper = nitend - nit000 + 1 ! Cycle interval length
205
206	CALL calc_date( nitend , ditend_date ) ! Date of final time step
207	CALL calc_date( nitbkg_r , ditbkg_date ) ! Background time for Jb referenced to ndate0
208	CALL calc_date( nitdin_r , ditdin_date ) ! Background time for DI referenced to ndate0
209	CALL calc_date( nitiaustr_r, ditiaustr_date ) ! IAU start time referenced to ndate0
210	CALL calc_date( nitiaufin_r, ditiaufin_date ) ! IAU end time referenced to ndate0
211
212	IF(lwp) THEN
213	WRITE(numout,*)
214	WRITE(numout,*) ' Time steps referenced to current cycle:'
215	WRITE(numout,*) ' iitrst = ', nit000 - 1
216	WRITE(numout,*) ' nit000 = ', nit000
217	WRITE(numout,*) ' nitend = ', nitend
218	WRITE(numout,*) ' nitbkg_r = ', nitbkg_r
219	WRITE(numout,*) ' nitdin_r = ', nitdin_r
220	WRITE(numout,*) ' nitiaustr_r = ', nitiaustr_r
221	WRITE(numout,*) ' nitiaufin_r = ', nitiaufin_r
222	WRITE(numout,*)
223	WRITE(numout,*) ' Dates referenced to current cycle:'
224	WRITE(numout,*) ' ndastp = ', ndastp
225	WRITE(numout,*) ' ndate0 = ', ndate0
226	WRITE(numout,*) ' nn_time0 = ', nn_time0
227	WRITE(numout,*) ' ditend_date = ', ditend_date
228	WRITE(numout,*) ' ditbkg_date = ', ditbkg_date
229	WRITE(numout,*) ' ditdin_date = ', ditdin_date
230	WRITE(numout,*) ' ditiaustr_date = ', ditiaustr_date
231	WRITE(numout,*) ' ditiaufin_date = ', ditiaufin_date
232	ENDIF
233
234
235	IF ( ( ln_asmdin ).AND.( ln_asmiau ) ) &
236	& CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
237	& ' Choose Direct Initialization OR Incremental Analysis Updating')
238
239	IF ( ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
240	.AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) &
241	& CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', &
242	& ' but ln_asmdin and ln_asmiau are both set to .false. :', &
243	& ' Inconsistent options')
244
245	IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
246	& CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :', &
247	& ' Type IAU weighting function is invalid')
248
249	IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
250	& ) &
251	& CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', &
252	& ' The assimilation increments are not applied')
253
254	IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
255	& CALL ctl_stop( ' nitiaustr = nitiaufin :', &
256	& ' IAU interval is of zero length')
257
258	IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
259	& CALL ctl_stop( ' nitiaustr or nitiaufin :', &
260	& ' IAU starting or final time step is outside the cycle interval', &
261	& ' Valid range nit000 to nitend')
262
263	IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
264	& CALL ctl_stop( ' nitbkg :', &
265	& ' Background time step is outside the cycle interval')
266
267	IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
268	& CALL ctl_stop( ' nitdin :', &
269	& ' Background time step for Direct Initialization is outside', &
270	& ' the cycle interval')
271
272	IF ( nstop > 0 ) RETURN ! if there are any errors then go no further
273
274	!--------------------------------------------------------------------
275	! Initialize the Incremental Analysis Updating weighting function
276	!--------------------------------------------------------------------
277
278	IF( ln_asmiau ) THEN
279	!
280	ALLOCATE( wgtiau( icycper ) )
281	!
282	wgtiau(:) = 0._wp
283	!
284	! !---------------------------------------------------------
285	IF( niaufn == 0 ) THEN ! Constant IAU forcing
286	! !---------------------------------------------------------
287	DO jt = 1, iiauper
288	wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
289	END DO
290	! !---------------------------------------------------------
291	ELSEIF ( niaufn == 1 ) THEN ! Linear hat-like, centred in middle of IAU interval
292	! !---------------------------------------------------------
293	! Compute the normalization factor
294	znorm = 0._wp
295	IF( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval
296	imid = iiauper / 2
297	DO jt = 1, imid
298	znorm = znorm + REAL( jt )
299	END DO
300	znorm = 2.0 * znorm
301	ELSE ! Odd number of time steps in IAU interval
302	imid = ( iiauper + 1 ) / 2
303	DO jt = 1, imid - 1
304	znorm = znorm + REAL( jt )
305	END DO
306	znorm = 2.0 * znorm + REAL( imid )
307	ENDIF
308	znorm = 1.0 / znorm
309	!
310	DO jt = 1, imid - 1
311	wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
312	END DO
313	DO jt = imid, iiauper
314	wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
315	END DO
316	!
317	ENDIF
318
319	! Test that the integral of the weights over the weighting interval equals 1
320	IF(lwp) THEN
321	WRITE(numout,*)
322	WRITE(numout,*) 'asm_inc_init : IAU weights'
323	WRITE(numout,*) '~~~~~~~~~~~~'
324	WRITE(numout,*) ' time step IAU weight'
325	WRITE(numout,*) ' ========= ====================='
326	ztotwgt = 0.0
327	DO jt = 1, icycper
328	ztotwgt = ztotwgt + wgtiau(jt)
329	WRITE(numout,*) ' ', jt, ' ', wgtiau(jt)
330	END DO
331	WRITE(numout,*) ' ==================================='
332	WRITE(numout,*) ' Time-integrated weight = ', ztotwgt
333	WRITE(numout,*) ' ==================================='
334	ENDIF
335
336	ENDIF
337
338	!--------------------------------------------------------------------
339	! Allocate and initialize the increment arrays
340	!--------------------------------------------------------------------
341
342	ALLOCATE( t_bkginc (jpi,jpj,jpk) ) ; t_bkginc (:,:,:) = 0._wp
343	ALLOCATE( s_bkginc (jpi,jpj,jpk) ) ; s_bkginc (:,:,:) = 0._wp
344	ALLOCATE( u_bkginc (jpi,jpj,jpk) ) ; u_bkginc (:,:,:) = 0._wp
345	ALLOCATE( v_bkginc (jpi,jpj,jpk) ) ; v_bkginc (:,:,:) = 0._wp
346	ALLOCATE( ssh_bkginc (jpi,jpj) ) ; ssh_bkginc (:,:) = 0._wp
347	ALLOCATE( seaice_bkginc(jpi,jpj) ) ; seaice_bkginc(:,:) = 0._wp
348	#if defined key_si3
349	ALLOCATE( a_i_bkginc (jpi,jpj,jpl) ) ; a_i_bkginc (:,:,:) = 0._wp
350	ALLOCATE( incr_newice(jpi,jpj,jpl) ) ; incr_newice(:,:,:) = .FALSE.
351	#endif
352	#if defined key_asminc
353	ALLOCATE( ssh_iau (jpi,jpj) ) ; ssh_iau (:,:) = 0._wp
354	#endif
355	#if defined key_cice && defined key_asminc
356	ALLOCATE( ndaice_da (jpi,jpj) ) ; ndaice_da (:,:) = 0._wp
357	#endif
358	!
359	IF ( ln_trainc .OR. ln_dyninc .OR. & !--------------------------------------
360	& ln_sshinc .OR. ln_seaiceinc ) THEN ! Read the increments from file
361	! !--------------------------------------
362	CALL iom_open( c_asminc, inum )
363	!
364	CALL iom_get( inum, 'time' , zdate_inc )
365	CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
366	CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
367	z_inc_dateb = zdate_inc
368	z_inc_datef = zdate_inc
369	!
370	IF(lwp) THEN
371	WRITE(numout,*)
372	WRITE(numout,*) 'asm_inc_init : Assimilation increments valid between dates ', z_inc_dateb,' and ', z_inc_datef
373	WRITE(numout,*) '~~~~~~~~~~~~'
374	ENDIF
375	!
376	IF ( ( z_inc_dateb < ndastp + nn_time0*0.0001_wp ) .OR. &
377	& ( z_inc_datef > ditend_date ) ) &
378	& CALL ctl_warn( ' Validity time of assimilation increments is ', &
379	& ' outside the assimilation interval' )
380
381	IF ( ( ln_asmdin ).AND.( zdate_inc /= ditdin_date ) ) &
382	& CALL ctl_warn( ' Validity time of assimilation increments does ', &
383	& ' not agree with Direct Initialization time' )
384
385	IF ( ln_trainc ) THEN
386	CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
387	CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
388	! Apply the masks
389	t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
390	s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
391	! Set missing increments to 0.0 rather than 1e+20
392	! to allow for differences in masks
393	WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
394	WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
395	ENDIF
396
397	IF ( ln_dyninc ) THEN
398	CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )
399	CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )
400	! Apply the masks
401	u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
402	v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
403	! Set missing increments to 0.0 rather than 1e+20
404	! to allow for differences in masks
405	WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
406	WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
407	ENDIF
408
409	IF ( ln_sshinc ) THEN
410	CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
411	! Apply the masks
412	ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
413	! Set missing increments to 0.0 rather than 1e+20
414	! to allow for differences in masks
415	WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
416	ENDIF
417
418	IF ( ln_seaiceinc ) THEN
419	CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 )
420	! Apply the masks
421	seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1)
422	! Set missing increments to 0.0 rather than 1e+20
423	! to allow for differences in masks
424	WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
425
426	IF (lwp) THEN
427	WRITE(numout,*)
428	WRITE(numout,*) 'asm_inc_init : Converting seaice_bkginc to a_i_bkginc'
429	WRITE(numout,*) '~~~~~~~~~~~~'
430	END IF
431	!single category increment for sea ice conc
432	!convert single category increment to multi category
433	a_i_bkginc = 0.0_wp
434	at_i = SUM(a_i(:,:,:), DIM=3)
435
436	! Calculate which category corresponds to zhi_damin
437	! find which category to fill
438	DO jl = 1, jpl
439	IF( zhi_damin > hi_max(jl-1) .AND. zhi_damin <= hi_max(jl) ) THEN
440	nhi_damin = jl
441	END IF
442	END DO
443
444	IF (lwp) THEN
445	WRITE(numout,*)
446	WRITE(numout,*) 'asm_inc_init : Converting seaice_bkginc to a_i_bkginc using Peterson splitting'
447	WRITE(numout,*) '~~~~~~~~~~~~'
448	END IF
449	DO jj = 1, jpj
450	DO ji = 1, jpi
451	IF ( seaice_bkginc(ji,jj) > 0.0_wp) THEN
452	!Positive ice concentration increments are always
453	!added to the thinnest category of ice
454	a_i_bkginc(ji,jj,nhi_damin) = seaice_bkginc(ji,jj)
455	ELSE
456	!negative increments are first removed from the thinnest
457	!available category until it reaches zero concentration
458	!and then progressively removed from thicker categories
459	zremaining_increment = seaice_bkginc(ji,jj)
460	DO jl = 1, jpl
461	! assign as much increment as possible to current category
462	a_i_bkginc(ji,jj,jl) = -min(a_i(ji,jj,jl), -zremaining_increment)
463	! update remaining amount of increment
464	zremaining_increment = zremaining_increment - a_i_bkginc(ji,jj,jl)
465	END DO
466	END IF
467	END DO
468	END DO
469	DO jl = 1,jpl
470	WHERE (at_i(:,:) < zopenwater_lim .AND. seaice_bkginc(:,:) > 0.0_wp)
471	incr_newice(:,:,jl) = .TRUE.
472	END WHERE
473	END DO
474	ENDIF
475	!
476	CALL iom_close( inum )
477	!
478	ENDIF
479	!
480	! !--------------------------------------
481	IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN ! Apply divergence damping filter
482	! !--------------------------------------
483	ALLOCATE( zhdiv(jpi,jpj) )
484	!
485	DO jt = 1, nn_divdmp
486	!
487	DO jk = 1, jpkm1 ! zhdiv = e1e1 * div
488	zhdiv(:,:) = 0._wp
489	DO jj = 2, jpjm1
490	DO ji = fs_2, fs_jpim1 ! vector opt.
491	zhdiv(ji,jj) = ( e2u(ji ,jj) * e3u_n(ji ,jj,jk) * u_bkginc(ji ,jj,jk) &
492	& - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * u_bkginc(ji-1,jj,jk) &
493	& + e1v(ji,jj ) * e3v_n(ji,jj ,jk) * v_bkginc(ji,jj ,jk) &
494	& - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * v_bkginc(ji,jj-1,jk) ) / e3t_n(ji,jj,jk)
495	END DO
496	END DO
497	CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. ) ! lateral boundary cond. (no sign change)
498	!
499	DO jj = 2, jpjm1
500	DO ji = fs_2, fs_jpim1 ! vector opt.
501	u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) &
502	& + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
503	v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) &
504	& + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
505	END DO
506	END DO
507	END DO
508	!
509	END DO
510	!
511	DEALLOCATE( zhdiv )
512	!
513	ENDIF
514	!
515	! !-----------------------------------------------------
516	IF ( ln_asmdin ) THEN ! Allocate and initialize the background state arrays
517	! !-----------------------------------------------------
518	!
519	ALLOCATE( t_bkg (jpi,jpj,jpk) ) ; t_bkg (:,:,:) = 0._wp
520	ALLOCATE( s_bkg (jpi,jpj,jpk) ) ; s_bkg (:,:,:) = 0._wp
521	ALLOCATE( u_bkg (jpi,jpj,jpk) ) ; u_bkg (:,:,:) = 0._wp
522	ALLOCATE( v_bkg (jpi,jpj,jpk) ) ; v_bkg (:,:,:) = 0._wp
523	ALLOCATE( ssh_bkg(jpi,jpj) ) ; ssh_bkg(:,:) = 0._wp
524	!
525	!
526	!--------------------------------------------------------------------
527	! Read from file the background state at analysis time
528	!--------------------------------------------------------------------
529	!
530	CALL iom_open( c_asmdin, inum )
531	!
532	CALL iom_get( inum, 'rdastp', zdate_bkg )
533	!
534	IF(lwp) THEN
535	WRITE(numout,*)
536	WRITE(numout,*) ' ==>>> Assimilation background state valid at : ', zdate_bkg
537	WRITE(numout,*)
538	ENDIF
539	!
540	IF ( zdate_bkg /= ditdin_date ) &
541	& CALL ctl_warn( ' Validity time of assimilation background state does', &
542	& ' not agree with Direct Initialization time' )
543	!
544	IF ( ln_trainc ) THEN
545	CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
546	CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
547	t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
548	s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
549	ENDIF
550	!
551	IF ( ln_dyninc ) THEN
552	CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
553	CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
554	u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
555	v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
556	ENDIF
557	!
558	IF ( ln_sshinc ) THEN
559	CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
560	ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
561	ENDIF
562	!
563	CALL iom_close( inum )
564	!
565	ENDIF
566	!
567	IF(lwp) WRITE(numout,*) ' ==>>> Euler time step switch is ', neuler
568	!
569	IF( lk_asminc ) THEN !== data assimilation ==!
570	IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 ) ! Output background fields
571	IF( ln_asmdin ) THEN ! Direct initialization
572	IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 ) ! Tracers
573	IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1 ) ! Dynamics
574	IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 ) ! SSH
575	ENDIF
576	ENDIF
577	!
578	END SUBROUTINE asm_inc_init
579
580
581	SUBROUTINE tra_asm_inc( kt )
582	!!----------------------------------------------------------------------
583	!! * ROUTINE tra_asm_inc *
584	!!
585	!! ** Purpose : Apply the tracer (T and S) assimilation increments
586	!!
587	!! ** Method : Direct initialization or Incremental Analysis Updating
588	!!
589	!! ** Action :
590	!!----------------------------------------------------------------------
591	INTEGER, INTENT(IN) :: kt ! Current time step
592	!
593	INTEGER :: ji, jj, jk
594	INTEGER :: it
595	REAL(wp) :: zincwgt ! IAU weight for current time step
596	REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
597	!!----------------------------------------------------------------------
598	!
599	! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
600	! used to prevent the applied increments taking the temperature below the local freezing point
601	DO jk = 1, jpkm1
602	CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), gdept_n(:,:,jk) )
603	END DO
604	!
605	! !--------------------------------------
606	IF ( ln_asmiau ) THEN ! Incremental Analysis Updating
607	! !--------------------------------------
608	!
609	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
610	!
611	it = kt - nit000 + 1
612	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
613	!
614	IF(lwp) THEN
615	WRITE(numout,*)
616	WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
617	WRITE(numout,*) '~~~~~~~~~~~~'
618	ENDIF
619	!
620	! Update the tracer tendencies
621	DO jk = 1, jpkm1
622	IF (ln_temnofreeze) THEN
623	! Do not apply negative increments if the temperature will fall below freezing
624	WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
625	& tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )
626	tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt
627	END WHERE
628	ELSE
629	tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt
630	ENDIF
631	IF (ln_salfix) THEN
632	! Do not apply negative increments if the salinity will fall below a specified
633	! minimum value salfixmin
634	WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
635	& tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )
636	tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
637	END WHERE
638	ELSE
639	tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
640	ENDIF
641	END DO
642	!
643	ENDIF
644	!
645	IF ( kt == nitiaufin_r + 1 ) THEN ! For bias crcn to work
646	DEALLOCATE( t_bkginc )
647	DEALLOCATE( s_bkginc )
648	ENDIF
649	! !--------------------------------------
650	ELSEIF ( ln_asmdin ) THEN ! Direct Initialization
651	! !--------------------------------------
652	!
653	IF ( kt == nitdin_r ) THEN
654	!
655	neuler = 0 ! Force Euler forward step
656	!
657	! Initialize the now fields with the background + increment
658	IF (ln_temnofreeze) THEN
659	! Do not apply negative increments if the temperature will fall below freezing
660	WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )
661	tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)
662	END WHERE
663	ELSE
664	tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)
665	ENDIF
666	IF (ln_salfix) THEN
667	! Do not apply negative increments if the salinity will fall below a specified
668	! minimum value salfixmin
669	WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )
670	tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)
671	END WHERE
672	ELSE
673	tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)
674	ENDIF
675
676	tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields
677
678	CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities
679	!!gm fabien
680	! CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities
681	!!gm
682
683	IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) &
684	& CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient
685	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
686	IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) &
687	& CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
688	& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the last ocean level
689
690	DEALLOCATE( t_bkginc )
691	DEALLOCATE( s_bkginc )
692	DEALLOCATE( t_bkg )
693	DEALLOCATE( s_bkg )
694	ENDIF
695	!
696	ENDIF
697	! Perhaps the following call should be in step
698	IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment
699	!
700	END SUBROUTINE tra_asm_inc
701
702
703	SUBROUTINE dyn_asm_inc( kt )
704	!!----------------------------------------------------------------------
705	!! * ROUTINE dyn_asm_inc *
706	!!
707	!! ** Purpose : Apply the dynamics (u and v) assimilation increments.
708	!!
709	!! ** Method : Direct initialization or Incremental Analysis Updating.
710	!!
711	!! ** Action :
712	!!----------------------------------------------------------------------
713	INTEGER, INTENT(IN) :: kt ! Current time step
714	!
715	INTEGER :: jk
716	INTEGER :: it
717	REAL(wp) :: zincwgt ! IAU weight for current time step
718	!!----------------------------------------------------------------------
719	!
720	! !--------------------------------------------
721	IF ( ln_asmiau ) THEN ! Incremental Analysis Updating
722	! !--------------------------------------------
723	!
724	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
725	!
726	it = kt - nit000 + 1
727	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
728	!
729	IF(lwp) THEN
730	WRITE(numout,*)
731	WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
732	WRITE(numout,*) '~~~~~~~~~~~~'
733	ENDIF
734	!
735	! Update the dynamic tendencies
736	DO jk = 1, jpkm1
737	ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
738	va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
739	END DO
740	!
741	IF ( kt == nitiaufin_r ) THEN
742	DEALLOCATE( u_bkginc )
743	DEALLOCATE( v_bkginc )
744	ENDIF
745	!
746	ENDIF
747	! !-----------------------------------------
748	ELSEIF ( ln_asmdin ) THEN ! Direct Initialization
749	! !-----------------------------------------
750	!
751	IF ( kt == nitdin_r ) THEN
752	!
753	neuler = 0 ! Force Euler forward step
754	!
755	! Initialize the now fields with the background + increment
756	un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
757	vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)
758	!
759	ub(:,:,:) = un(:,:,:) ! Update before fields
760	vb(:,:,:) = vn(:,:,:)
761	!
762	DEALLOCATE( u_bkg )
763	DEALLOCATE( v_bkg )
764	DEALLOCATE( u_bkginc )
765	DEALLOCATE( v_bkginc )
766	ENDIF
767	!
768	ENDIF
769	!
770	END SUBROUTINE dyn_asm_inc
771
772
773	SUBROUTINE ssh_asm_inc( kt )
774	!!----------------------------------------------------------------------
775	!! * ROUTINE ssh_asm_inc *
776	!!
777	!! ** Purpose : Apply the sea surface height assimilation increment.
778	!!
779	!! ** Method : Direct initialization or Incremental Analysis Updating.
780	!!
781	!! ** Action :
782	!!----------------------------------------------------------------------
783	INTEGER, INTENT(IN) :: kt ! Current time step
784	!
785	INTEGER :: it
786	INTEGER :: jk
787	REAL(wp) :: zincwgt ! IAU weight for current time step
788	!!----------------------------------------------------------------------
789	!
790	! !-----------------------------------------
791	IF ( ln_asmiau ) THEN ! Incremental Analysis Updating
792	! !-----------------------------------------
793	!
794	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
795	!
796	it = kt - nit000 + 1
797	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
798	!
799	IF(lwp) THEN
800	WRITE(numout,*)
801	WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
802	& kt,' with IAU weight = ', wgtiau(it)
803	WRITE(numout,*) '~~~~~~~~~~~~'
804	ENDIF
805	!
806	! Save the tendency associated with the IAU weighted SSH increment
807	! (applied in dynspg.*)
808	#if defined key_asminc
809	ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
810	#endif
811	!
812	ELSE IF( kt == nitiaufin_r+1 ) THEN
813	!
814	! test on ssh_bkginc needed as ssh_asm_inc is called twice by time step
815	IF (ALLOCATED(ssh_bkginc)) DEALLOCATE( ssh_bkginc )
816	!
817	#if defined key_asminc
818	ssh_iau(:,:) = 0._wp
819	#endif
820	!
821	ENDIF
822	! !-----------------------------------------
823	ELSEIF ( ln_asmdin ) THEN ! Direct Initialization
824	! !-----------------------------------------
825	!
826	IF ( kt == nitdin_r ) THEN
827	!
828	neuler = 0 ! Force Euler forward step
829	!
830	sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) ! Initialize the now fields the background + increment
831	!
832	sshb(:,:) = sshn(:,:) ! Update before fields
833	e3t_b(:,:,:) = e3t_n(:,:,:)
834	!!gm why not e3u_b, e3v_b, gdept_b ????
835	!
836	DEALLOCATE( ssh_bkg )
837	DEALLOCATE( ssh_bkginc )
838	!
839	ENDIF
840	!
841	ENDIF
842	!
843	END SUBROUTINE ssh_asm_inc
844
845
846	SUBROUTINE ssh_asm_div( kt, phdivn )
847	!!----------------------------------------------------------------------
848	!! * ROUTINE ssh_asm_div *
849	!!
850	!! ** Purpose : ssh increment with z* is incorporated via a correction of the local divergence
851	!! across all the water column
852	!!
853	!! ** Method :
854	!! CAUTION : sshiau is positive (inflow) decreasing the
855	!! divergence and expressed in m/s
856	!!
857	!! ** Action : phdivn decreased by the ssh increment
858	!!----------------------------------------------------------------------
859	INTEGER, INTENT(IN) :: kt ! ocean time-step index
860	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence
861	!!
862	INTEGER :: jk ! dummy loop index
863	REAL(wp), DIMENSION(:,:) , POINTER :: ztim ! local array
864	!!----------------------------------------------------------------------
865	!
866	#if defined key_asminc
867	CALL ssh_asm_inc( kt ) !== (calculate increments)
868	!
869	IF( ln_linssh ) THEN
870	phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t_n(:,:,1) * tmask(:,:,1)
871	ELSE
872	ALLOCATE( ztim(jpi,jpj) )
873	ztim(:,:) = ssh_iau(:,:) / ( ht_n(:,:) + 1.0 - ssmask(:,:) )
874	DO jk = 1, jpkm1
875	phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk)
876	END DO
877	!
878	DEALLOCATE(ztim)
879	ENDIF
880	#endif
881	!
882	END SUBROUTINE ssh_asm_div
883
884
885	SUBROUTINE seaice_asm_inc( kt, kindic )
886	!!----------------------------------------------------------------------
887	!! * ROUTINE seaice_asm_inc *
888	!!
889	!! ** Purpose : Apply the sea ice assimilation increment.
890	!!
891	!! ** Method : Direct initialization or Incremental Analysis Updating.
892	!!
893	!! ** Action :
894	!!
895	!!----------------------------------------------------------------------
896	INTEGER, INTENT(in) :: kt ! Current time step
897	INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation
898	!
899	INTEGER :: it, jl
900	REAL(wp) :: zincwgt ! IAU weight for current time step
901	#if defined key_si3
902	REAL(wp), DIMENSION(jpi,jpj,jpl) :: da_i ! change in ice concentration
903	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i ! inverse of old ice concentration
904	REAL(wp) :: rn_hinew_save
905	LOGICAL, SAVE :: initial_step=.TRUE.
906	REAL(wp), DIMENSION(jpi,jpj) :: at_i_bkginc ! sum of conc increment over categories
907	#endif
908	!!----------------------------------------------------------------------
909	!
910	! !-----------------------------------------
911	IF ( ln_asmiau ) THEN ! Incremental Analysis Updating
912	! !-----------------------------------------
913	!
914	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
915	!
916	it = kt - nit000 + 1
917	zincwgt = wgtiau(it) ! IAU weight for the current time step
918	! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
919	!
920	IF(lwp) THEN
921	WRITE(numout,*)
922	WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
923	WRITE(numout,*) '~~~~~~~~~~~~'
924	ENDIF
925	!
926	! Sea-ice : SI3 case
927	!
928	#if defined key_si3
929	WHERE( a_i(:,:,:) > epsi10 )
930	z1_a_i(:,:,:) = 1.0_wp/a_i(:,:,:)
931	ELSEWHERE
932	z1_a_i(:,:,:) = 0.0_wp
933	END WHERE
934
935	WRITE(numout,*) "zhi_damin = ",zhi_damin
936	! add positive concentration increments with zhi_damin to regions where ice
937	! is already present and bound them to 1
938	! ice volume is added based on the sea ice conc increment and assigned thickness
939	WHERE ( .NOT. incr_newice .AND. a_i_bkginc(:,:,:) > 0.0_wp )
940	a_i(:,:,:) = a_i(:,:,:) + MIN(1.0_wp - a_i(:,:,:), a_i_bkginc(:,:,:) * zincwgt)
941	v_i(:,:,:) = v_i(:,:,:) + MIN(1.0_wp - a_i(:,:,:), a_i_bkginc(:,:,:) * zincwgt) * zhi_damin
942	END WHERE
943
944	! add negative concentration increments to regions where ice
945	! is already present and bound them to 0
946	! in this case ice volume is changed based on the current thickness
947	WHERE ( .NOT. incr_newice .AND. a_i_bkginc(:,:,:) < 0.0_wp )
948	a_i(:,:,:) = MAX(a_i(:,:,:) + a_i_bkginc(:,:,:) * zincwgt, 0.0_wp)
949	v_i(:,:,:) = a_i(:,:,:) * h_i(:,:,:)
950	END WHERE
951
952	! compute change in ice concentration (new / old)
953	WHERE ( incr_newice )
954	da_i(:,:,:) = 1.0_wp
955	ELSEWHERE
956	da_i(:,:,:) = a_i(:,:,:) * z1_a_i(:,:,:)
957	END WHERE
958
959	! compute heat balance that adds specified ice thickness
960	! to open water
961	IF (initial_step) THEN
962	IF(lwp) THEN
963	WRITE(numout,*)
964	WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' compute heat balance qlead'
965	WRITE(numout,*) '~~~~~~~~~~~~'
966	ENDIF
967
968	! compute qlead to ensure thickness of zhi_damin
969	! for the new ice concentration values
970	WHERE (ANY(incr_newice, DIM=3))
971	ht_i_new(:,:) = zhi_damin
972	ELSEWHERE
973	ht_i_new(:,:) = 0.0_wp
974	ENDWHERE
975
976	! ensure all categories of new ice are zero
977	WHERE ( incr_newice )
978	v_i (:,:,:) = 0.0_wp
979	v_s (:,:,:) = 0.0_wp
980	sv_i(:,:,:) = 0.0_wp
981	a_ip(:,:,:) = 0.0_wp
982	v_ip(:,:,:) = 0.0_wp
983	END WHERE
984	DO jl = 1, nlay_i
985	WHERE ( incr_newice )
986	e_i(:,:,jl,:) = 0.0_wp
987	END WHERE
988	END DO
989	DO jl = 1, nlay_s
990	WHERE ( incr_newice )
991	e_s(:,:,jl,:) = 0.0_wp
992	END WHERE
993	END DO
994
995	qlead = 0.0_wp
996	at_i_bkginc(:,:) = SUM( a_i_bkginc(:,:,:)*zincwgt, DIM=3 )
997
998	call seaice_asm_qlead(ht_i_new, at_i_bkginc, qlead)
999
1000	! Call ice_thd_do to create the new ice from open water
1001	rn_hinew_save = rn_hinew
1002	rn_hinew = zhi_damin
1003	call ice_thd_do
1004	rn_hinew = rn_hinew_save
1005
1006	! compute thickness now
1007	WHERE( a_i(:,:,:) > epsi10 )
1008	z1_a_i(:,:,:) = 1.0_wp/a_i(:,:,:)
1009	ELSEWHERE
1010	z1_a_i(:,:,:) = 0.0_wp
1011	END WHERE
1012	h_i(:,:,:) = v_i(:,:,:) * z1_a_i(:,:,:)
1013
1014	incr_newice(:,:,:) = .FALSE.
1015	initial_step = .FALSE.
1016	END IF
1017
1018	! make volume consistent with concentration and thickness
1019	! note where conc is zero this will lead to zero volume
1020	! and soon zero thickness after the call to ice_var_zapsmall
1021	v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
1022
1023	! as ice concentration has changed, maintain equivalent values for prognostic variables
1024	v_s (:,:,:) = v_s (:,:,:) * da_i(:,:,:)
1025	sv_i (:,:,:) = sv_i(:,:,:) * da_i(:,:,:)
1026	a_ip (:,:,:) = a_ip(:,:,:) * da_i(:,:,:)
1027	v_ip (:,:,:) = v_ip(:,:,:) * da_i(:,:,:)
1028	DO jl = 1, nlay_i
1029	e_i(:,:,jl,:) = e_i(:,:,jl,:) * da_i(:,:,:)
1030	END DO
1031	DO jl = 1, nlay_s
1032	e_s(:,:,jl,:) = e_s(:,:,jl,:) * da_i(:,:,:)
1033	END DO
1034
1035	call ice_var_zapsmall()
1036
1037	!
1038	! seaice salinity balancing (to add)
1039	#endif
1040	!
1041	#if defined key_cice && defined key_asminc
1042	! Sea-ice : CICE case. Pass ice increment tendency into CICE
1043	ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
1044	#endif
1045	!
1046	IF ( kt == nitiaufin_r ) THEN
1047	DEALLOCATE( seaice_bkginc )
1048	ENDIF
1049	!
1050	ELSE
1051	!
1052	#if defined key_cice && defined key_asminc
1053	ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1054	#endif
1055	!
1056	ENDIF
1057	! !-----------------------------------------
1058	ELSEIF ( ln_asmdin ) THEN ! Direct Initialization
1059	! !-----------------------------------------
1060	!
1061	IF ( kt == nitdin_r ) THEN
1062	!
1063	neuler = 0 ! Force Euler forward step
1064
1065	IF(lwp) THEN
1066	WRITE(numout,*)
1067	WRITE(numout,*) 'seaice_asm_inc : sea ice direct initialization at time step = ', kt
1068	WRITE(numout,*) '~~~~~~~~~~~~'
1069	ENDIF
1070	!
1071	#if defined key_cice && defined key_asminc
1072	! Sea-ice : CICE case. Pass ice increment tendency into CICE
1073	ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
1074	#endif
1075	IF ( .NOT. PRESENT(kindic) ) THEN
1076	DEALLOCATE( seaice_bkginc )
1077	END IF
1078	!
1079	ELSE
1080	!
1081	#if defined key_cice && defined key_asminc
1082	ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1083	#endif
1084	!
1085	ENDIF
1086	!
1087	ENDIF
1088	!
1089	END SUBROUTINE seaice_asm_inc
1090
1091	SUBROUTINE seaice_asm_qlead(ht_i_new, at_i_new, qlead_new)
1092	!!----------------------------------------------------------------------
1093	!! * ROUTINE seaice_asm_qlead *
1094	!!
1095	!! ** Purpose : Compute the value of qlead that will correspond to new ice of
1096	!! thickness ht_i_new concentration of at_i_new
1097	!!
1098	!! ** Method : Based on symbolic inversion of the icethd_do routine
1099	!!
1100	!! ** Action : return qlead_new
1101	!!----------------------------------------------------------------------
1102	REAL(wp), DIMENSION(:,:), INTENT(in) :: ht_i_new ! total thickness of new ice
1103	REAL(wp), DIMENSION(:,:), INTENT(in) :: at_i_new ! total concentration of new ice
1104	REAL(wp), DIMENSION(:,:), INTENT(inout) :: qlead_new ! output flux? value
1105
1106	INTEGER :: jj, ji
1107
1108	REAL(wp), DIMENSION(jpij) :: at_i_new_1d ! 1d version of at_i_new
1109	REAL(wp), DIMENSION(jpij) :: zs_newice ! salinity of accreted ice
1110	REAL(wp), DIMENSION(jpij) :: zh_newice ! thickness of accreted ice
1111	!!----------------------------------------------------------------------
1112
1113	! Identify grid points where new ice forms
1114	npti = 0 ; nptidx(:) = 0
1115	DO jj = 1, jpj
1116	DO ji = 1, jpi
1117	IF ( ht_i_new(ji,jj) > 0._wp ) THEN
1118	npti = npti + 1
1119	nptidx( npti ) = (jj - 1) * jpi + ji
1120	ENDIF
1121	END DO
1122	END DO
1123
1124	! Move from 2-D to 1-D vectors
1125	IF ( npti > 0 ) THEN
1126	CALL tab_2d_1d( npti, nptidx(1:npti), at_i_new_1d(1:npti), at_i_new )
1127	CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new )
1128	CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d (1:npti) , sss_m )
1129	CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d (1:npti) , qlead_new )
1130	CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d (1:npti) , t_bo )
1131
1132
1133	! --- Salinity of new ice --- !
1134	SELECT CASE ( nn_icesal )
1135	CASE ( 1 ) ! Sice = constant
1136	zs_newice(1:npti) = rn_icesal
1137	CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)]
1138	DO ji = 1, npti
1139	zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
1140	END DO
1141	CASE ( 3 ) ! Sice = F(z) [multiyear ice]
1142	zs_newice(1:npti) = 2.3
1143	END SELECT
1144
1145
1146	DO ji = 1, npti
1147	qlead_1d(ji) = at_i_new_1d(ji)zh_newice(ji)(-r1_rhoirLfusrTmltrhoizs_newice(ji) + &
1148	& r1_rhoirhoi(-rLfus - rTmltrcpzs_newice(ji) + rTmltrcpizs_newice(ji) - &
1149	& rcpirt0 + rcpit_bo_1d(ji))* &
1150	& min(-epsi10, -rt0 + t_bo_1d(ji)) + &
1151	& rcp(rt0 - t_bo_1d(ji))min(-epsi10, -rt0 + t_bo_1d(ji))) / &
1152	& (r1_rhoi*min(-epsi10, -rt0 + t_bo_1d(ji)))
1153	END DO
1154	CALL tab_1d_2d( npti, nptidx(1:npti), qlead_1d(1:npti), qlead_new )
1155	ENDIF ! npti > 0
1156	END SUBROUTINE seaice_asm_qlead
1157
1158	!!======================================================================
1159	END MODULE asminc

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

Download in other formats: