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/2014/dev_r4650_UKMO14.11_SETTE_OBSASM/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.11_SETTE_OBSASM/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 4781

Last change on this file since 4781 was 4781, checked in by djlea, 10 years ago
Remove redundant error check preventing background write and increments being applied in the same run. Switch to use ctl_stop in diaobs. Fix to jul2greg to deal properly with negative julian dates. Tidying obs_write. Restore agrif sette test.
Property svn:keywords set to `Id`
File size: 48.4 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	!!----------------------------------------------------------------------
25	USE wrk_nemo ! Memory Allocation
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 oce ! Dynamics and active tracers defined in memory
30	USE ldfdyn_oce ! ocean dynamics: lateral physics
31	USE eosbn2 ! Equation of state - in situ and potential density
32	USE zpshde ! Partial step : Horizontal Derivative
33	USE iom ! Library to read input files
34	USE asmpar ! Parameters for the assmilation interface
35	USE c1d ! 1D initialization
36	USE in_out_manager ! I/O manager
37	USE lib_mpp ! MPP library
38	#if defined key_lim2
39	USE ice_2 ! LIM2
40	#endif
41	USE sbc_oce ! Surface boundary condition variables.
42
43	IMPLICIT NONE
44	PRIVATE
45
46	PUBLIC asm_inc_init !: Initialize the increment arrays and IAU weights
47	PUBLIC calc_date !: Compute the calendar date YYYYMMDD on a given step
48	PUBLIC tra_asm_inc !: Apply the tracer (T and S) increments
49	PUBLIC dyn_asm_inc !: Apply the dynamic (u and v) increments
50	PUBLIC ssh_asm_inc !: Apply the SSH increment
51	PUBLIC seaice_asm_inc !: Apply the seaice increment
52
53	#if defined key_asminc
54	LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE. !: Logical switch for assimilation increment interface
55	#else
56	LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE. !: No assimilation increments
57	#endif
58	LOGICAL, PUBLIC :: ln_bkgwri !: No output of the background state fields
59	LOGICAL, PUBLIC :: ln_asmiau !: No applying forcing with an assimilation increment
60	LOGICAL, PUBLIC :: ln_asmdin !: No direct initialization
61	LOGICAL, PUBLIC :: ln_trainc !: No tracer (T and S) assimilation increments
62	LOGICAL, PUBLIC :: ln_dyninc !: No dynamics (u and v) assimilation increments
63	LOGICAL, PUBLIC :: ln_sshinc !: No sea surface height assimilation increment
64	LOGICAL, PUBLIC :: ln_seaiceinc !: No sea ice concentration increment
65	LOGICAL, PUBLIC :: ln_salfix !: Apply minimum salinity check
66	LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
67	INTEGER, PUBLIC :: nn_divdmp !: Apply divergence damping filter nn_divdmp times
68
69	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkg , s_bkg !: Background temperature and salinity
70	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkg , v_bkg !: Background u- & v- velocity components
71	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkginc, s_bkginc !: Increment to the background T & S
72	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components
73	REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step
74	#if defined key_asminc
75	REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssh_iau !: IAU-weighted sea surface height increment
76	#endif
77	! !!! time steps relative to the cycle interval [0,nitend-nit000-1]
78	INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term
79	INTEGER , PUBLIC :: nitdin !: Time step of the background state for direct initialization
80	INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval
81	INTEGER , PUBLIC :: nitiaufin !: Time step of the end of the IAU interval
82	!
83	INTEGER , PUBLIC :: niaufn !: Type of IAU weighing function: = 0 Constant weighting
84	! !: = 1 Linear hat-like, centred in middle of IAU interval
85	REAL(wp), PUBLIC :: salfixmin !: Ensure that the salinity is larger than this value if (ln_salfix)
86
87	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc ! Background sea surface height and its increment
88	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: seaice_bkginc ! Increment to the background sea ice conc
89
90	!! * Substitutions
91	# include "domzgr_substitute.h90"
92	# include "ldfdyn_substitute.h90"
93	# include "vectopt_loop_substitute.h90"
94	!!----------------------------------------------------------------------
95	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
96	!! $Id$
97	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
98	!!----------------------------------------------------------------------
99	CONTAINS
100
101	SUBROUTINE asm_inc_init
102	!!----------------------------------------------------------------------
103	!! * ROUTINE asm_inc_init *
104	!!
105	!! ** Purpose : Initialize the assimilation increment and IAU weights.
106	!!
107	!! ** Method : Initialize the assimilation increment and IAU weights.
108	!!
109	!! ** Action :
110	!!----------------------------------------------------------------------
111	INTEGER :: ji, jj, jk
112	INTEGER :: jt
113	INTEGER :: imid
114	INTEGER :: inum
115	INTEGER :: iiauper ! Number of time steps in the IAU period
116	INTEGER :: icycper ! Number of time steps in the cycle
117	INTEGER :: iitend_date ! Date YYYYMMDD of final time step
118	INTEGER :: iitbkg_date ! Date YYYYMMDD of background time step for Jb term
119	INTEGER :: iitdin_date ! Date YYYYMMDD of background time step for DI
120	INTEGER :: iitiaustr_date ! Date YYYYMMDD of IAU interval start time step
121	INTEGER :: iitiaufin_date ! Date YYYYMMDD of IAU interval final time step
122	INTEGER :: ios ! Local integer output status for namelist read
123
124	REAL(wp) :: znorm ! Normalization factor for IAU weights
125	REAL(wp) :: ztotwgt ! Value of time-integrated IAU weights
126	! (should be equal to one)
127	REAL(wp) :: z_inc_dateb ! Start date of interval on which increment is valid
128	REAL(wp) :: z_inc_datef ! End date of interval on which increment is valid
129	REAL(wp) :: zdate_bkg ! Date in background state file for DI
130	REAL(wp) :: zdate_inc ! Time axis in increments file
131
132	REAL(wp), POINTER, DIMENSION(:,:) :: hdiv
133	!!
134	NAMELIST/nam_asminc/ ln_bkgwri, &
135	& ln_trainc, ln_dyninc, ln_sshinc, &
136	& ln_asmdin, ln_asmiau, &
137	& nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, &
138	& ln_salfix, salfixmin, &
139	& nn_divdmp
140	!!----------------------------------------------------------------------
141
142	!-----------------------------------------------------------------------
143	! Read Namelist nam_asminc : assimilation increment interface
144	!-----------------------------------------------------------------------
145
146	ln_seaiceinc = .FALSE.
147	ln_temnofreeze = .FALSE.
148
149	REWIND( numnam_ref ) ! Namelist nam_asminc in reference namelist : Assimilation increment
150	READ ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
151	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in reference namelist', lwp )
152
153	REWIND( numnam_cfg ) ! Namelist nam_asminc in configuration namelist : Assimilation increment
154	READ ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
155	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in configuration namelist', lwp )
156	IF(lwm) WRITE ( numond, nam_asminc )
157
158	! Control print
159	IF(lwp) THEN
160	WRITE(numout,*)
161	WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
162	WRITE(numout,*) '~~~~~~~~~~~~'
163	WRITE(numout,*) ' Namelist namasm : set assimilation increment parameters'
164	WRITE(numout,*) ' Logical switch for writing out background state ln_bkgwri = ', ln_bkgwri
165	WRITE(numout,*) ' Logical switch for applying tracer increments ln_trainc = ', ln_trainc
166	WRITE(numout,*) ' Logical switch for applying velocity increments ln_dyninc = ', ln_dyninc
167	WRITE(numout,*) ' Logical switch for applying SSH increments ln_sshinc = ', ln_sshinc
168	WRITE(numout,*) ' Logical switch for Direct Initialization (DI) ln_asmdin = ', ln_asmdin
169	WRITE(numout,*) ' Logical switch for applying sea ice increments ln_seaiceinc = ', ln_seaiceinc
170	WRITE(numout,*) ' Logical switch for Incremental Analysis Updating (IAU) ln_asmiau = ', ln_asmiau
171	WRITE(numout,*) ' Timestep of background in [0,nitend-nit000-1] nitbkg = ', nitbkg
172	WRITE(numout,*) ' Timestep of background for DI in [0,nitend-nit000-1] nitdin = ', nitdin
173	WRITE(numout,*) ' Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
174	WRITE(numout,*) ' Timestep of end of IAU interval in [0,nitend-nit000-1] nitiaufin = ', nitiaufin
175	WRITE(numout,*) ' Type of IAU weighting function niaufn = ', niaufn
176	WRITE(numout,*) ' Logical switch for ensuring that the sa > salfixmin ln_salfix = ', ln_salfix
177	WRITE(numout,*) ' Minimum salinity after applying the increments salfixmin = ', salfixmin
178	ENDIF
179
180	nitbkg_r = nitbkg + nit000 - 1 ! Background time referenced to nit000
181	nitdin_r = nitdin + nit000 - 1 ! Background time for DI referenced to nit000
182	nitiaustr_r = nitiaustr + nit000 - 1 ! Start of IAU interval referenced to nit000
183	nitiaufin_r = nitiaufin + nit000 - 1 ! End of IAU interval referenced to nit000
184
185	iiauper = nitiaufin_r - nitiaustr_r + 1 ! IAU interval length
186	icycper = nitend - nit000 + 1 ! Cycle interval length
187
188	! Date of final time step
189	CALL calc_date( nit000, nitend, ndate0, iitend_date )
190
191	! Background time for Jb referenced to ndate0
192	CALL calc_date( nit000, nitbkg_r, ndate0, iitbkg_date )
193
194	! Background time for DI referenced to ndate0
195	CALL calc_date( nit000, nitdin_r, ndate0, iitdin_date )
196
197	! IAU start time referenced to ndate0
198	CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date )
199
200	! IAU end time referenced to ndate0
201	CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date )
202
203	IF(lwp) THEN
204	WRITE(numout,*)
205	WRITE(numout,*) ' Time steps referenced to current cycle:'
206	WRITE(numout,*) ' iitrst = ', nit000 - 1
207	WRITE(numout,*) ' nit000 = ', nit000
208	WRITE(numout,*) ' nitend = ', nitend
209	WRITE(numout,*) ' nitbkg_r = ', nitbkg_r
210	WRITE(numout,*) ' nitdin_r = ', nitdin_r
211	WRITE(numout,*) ' nitiaustr_r = ', nitiaustr_r
212	WRITE(numout,*) ' nitiaufin_r = ', nitiaufin_r
213	WRITE(numout,*)
214	WRITE(numout,*) ' Dates referenced to current cycle:'
215	WRITE(numout,*) ' ndastp = ', ndastp
216	WRITE(numout,*) ' ndate0 = ', ndate0
217	WRITE(numout,*) ' iitend_date = ', iitend_date
218	WRITE(numout,*) ' iitbkg_date = ', iitbkg_date
219	WRITE(numout,*) ' iitdin_date = ', iitdin_date
220	WRITE(numout,*) ' iitiaustr_date = ', iitiaustr_date
221	WRITE(numout,*) ' iitiaufin_date = ', iitiaufin_date
222	ENDIF
223
224	IF ( nacc /= 0 ) &
225	& CALL ctl_stop( ' nacc /= 0 and key_asminc :', &
226	& ' Assimilation increments have only been implemented', &
227	& ' for synchronous time stepping' )
228
229	IF ( ( ln_asmdin ).AND.( ln_asmiau ) ) &
230	& CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
231	& ' Choose Direct Initialization OR Incremental Analysis Updating')
232
233	IF ( ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
234	.AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) &
235	& CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', &
236	& ' but ln_asmdin and ln_asmiau are both set to .false. :', &
237	& ' Inconsistent options')
238
239	IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
240	& CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :', &
241	& ' Type IAU weighting function is invalid')
242
243	IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
244	& ) &
245	& CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', &
246	& ' The assimilation increments are not applied')
247
248	IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
249	& CALL ctl_stop( ' nitiaustr = nitiaufin :', &
250	& ' IAU interval is of zero length')
251
252	IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
253	& CALL ctl_stop( ' nitiaustr or nitiaufin :', &
254	& ' IAU starting or final time step is outside the cycle interval', &
255	& ' Valid range nit000 to nitend')
256
257	IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
258	& CALL ctl_stop( ' nitbkg :', &
259	& ' Background time step is outside the cycle interval')
260
261	IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
262	& CALL ctl_stop( ' nitdin :', &
263	& ' Background time step for Direct Initialization is outside', &
264	& ' the cycle interval')
265
266	IF ( nstop > 0 ) RETURN ! if there are any errors then go no further
267
268	!--------------------------------------------------------------------
269	! Initialize the Incremental Analysis Updating weighting function
270	!--------------------------------------------------------------------
271
272	IF ( ln_asmiau ) THEN
273
274	ALLOCATE( wgtiau( icycper ) )
275
276	wgtiau(:) = 0.0
277
278	IF ( niaufn == 0 ) THEN
279
280	!---------------------------------------------------------
281	! Constant IAU forcing
282	!---------------------------------------------------------
283
284	DO jt = 1, iiauper
285	wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
286	END DO
287
288	ELSEIF ( niaufn == 1 ) THEN
289
290	!---------------------------------------------------------
291	! Linear hat-like, centred in middle of IAU interval
292	!---------------------------------------------------------
293
294	! Compute the normalization factor
295	znorm = 0.0
296	IF ( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval
297	imid = iiauper / 2
298	DO jt = 1, imid
299	znorm = znorm + REAL( jt )
300	END DO
301	znorm = 2.0 * znorm
302	ELSE ! Odd number of time steps in IAU interval
303	imid = ( iiauper + 1 ) / 2
304	DO jt = 1, imid - 1
305	znorm = znorm + REAL( jt )
306	END DO
307	znorm = 2.0 * znorm + REAL( imid )
308	ENDIF
309	znorm = 1.0 / znorm
310
311	DO jt = 1, imid - 1
312	wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
313	END DO
314	DO jt = imid, iiauper
315	wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
316	END DO
317
318	ENDIF
319
320	! Test that the integral of the weights over the weighting interval equals 1
321	IF(lwp) THEN
322	WRITE(numout,*)
323	WRITE(numout,*) 'asm_inc_init : IAU weights'
324	WRITE(numout,*) '~~~~~~~~~~~~'
325	WRITE(numout,*) ' time step IAU weight'
326	WRITE(numout,*) ' ========= ====================='
327	ztotwgt = 0.0
328	DO jt = 1, icycper
329	ztotwgt = ztotwgt + wgtiau(jt)
330	WRITE(numout,*) ' ', jt, ' ', wgtiau(jt)
331	END DO
332	WRITE(numout,*) ' ==================================='
333	WRITE(numout,*) ' Time-integrated weight = ', ztotwgt
334	WRITE(numout,*) ' ==================================='
335	ENDIF
336
337	ENDIF
338
339	!--------------------------------------------------------------------
340	! Allocate and initialize the increment arrays
341	!--------------------------------------------------------------------
342
343	ALLOCATE( t_bkginc(jpi,jpj,jpk) )
344	ALLOCATE( s_bkginc(jpi,jpj,jpk) )
345	ALLOCATE( u_bkginc(jpi,jpj,jpk) )
346	ALLOCATE( v_bkginc(jpi,jpj,jpk) )
347	ALLOCATE( ssh_bkginc(jpi,jpj) )
348	ALLOCATE( seaice_bkginc(jpi,jpj))
349	#if defined key_asminc
350	ALLOCATE( ssh_iau(jpi,jpj) )
351	#endif
352	t_bkginc(:,:,:) = 0.0
353	s_bkginc(:,:,:) = 0.0
354	u_bkginc(:,:,:) = 0.0
355	v_bkginc(:,:,:) = 0.0
356	ssh_bkginc(:,:) = 0.0
357	seaice_bkginc(:,:) = 0.0
358	#if defined key_asminc
359	ssh_iau(:,:) = 0.0
360	#endif
361	IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN
362
363	!--------------------------------------------------------------------
364	! Read the increments from file
365	!--------------------------------------------------------------------
366
367	CALL iom_open( c_asminc, inum )
368
369	CALL iom_get( inum, 'time', zdate_inc )
370
371	CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
372	CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
373	z_inc_dateb = zdate_inc
374	z_inc_datef = zdate_inc
375
376	IF(lwp) THEN
377	WRITE(numout,*)
378	WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', &
379	& ' between dates ', NINT( z_inc_dateb ),' and ', &
380	& NINT( z_inc_datef )
381	WRITE(numout,*) '~~~~~~~~~~~~'
382	ENDIF
383
384	IF ( ( NINT( z_inc_dateb ) < ndastp ) &
385	& .OR.( NINT( z_inc_datef ) > iitend_date ) ) &
386	& CALL ctl_warn( ' Validity time of assimilation increments is ', &
387	& ' outside the assimilation interval' )
388
389	IF ( ( ln_asmdin ).AND.( NINT( zdate_inc ) /= iitdin_date ) ) &
390	& CALL ctl_warn( ' Validity time of assimilation increments does ', &
391	& ' not agree with Direct Initialization time' )
392
393	IF ( ln_trainc ) THEN
394	CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
395	CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
396	! Apply the masks
397	t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
398	s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
399	! Set missing increments to 0.0 rather than 1e+20
400	! to allow for differences in masks
401	WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
402	WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
403	ENDIF
404
405	IF ( ln_dyninc ) THEN
406	CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )
407	CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )
408	! Apply the masks
409	u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
410	v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
411	! Set missing increments to 0.0 rather than 1e+20
412	! to allow for differences in masks
413	WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
414	WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
415	ENDIF
416
417	IF ( ln_sshinc ) THEN
418	CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
419	! Apply the masks
420	ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
421	! Set missing increments to 0.0 rather than 1e+20
422	! to allow for differences in masks
423	WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
424	ENDIF
425
426	IF ( ln_seaiceinc ) THEN
427	CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 )
428	! Apply the masks
429	seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1)
430	! Set missing increments to 0.0 rather than 1e+20
431	! to allow for differences in masks
432	WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
433	ENDIF
434
435	CALL iom_close( inum )
436
437	ENDIF
438
439	!-----------------------------------------------------------------------
440	! Apply divergence damping filter
441	!-----------------------------------------------------------------------
442
443	IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN
444
445	CALL wrk_alloc(jpi,jpj,hdiv)
446
447	DO jt = 1, nn_divdmp
448
449	DO jk = 1, jpkm1
450
451	hdiv(:,:) = 0._wp
452
453	DO jj = 2, jpjm1
454	DO ji = fs_2, fs_jpim1 ! vector opt.
455	hdiv(ji,jj) = &
456	( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * u_bkginc(ji ,jj ,jk) &
457	- e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * u_bkginc(ji-1,jj ,jk) &
458	+ e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * v_bkginc(ji ,jj ,jk) &
459	- e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * v_bkginc(ji ,jj-1,jk) ) &
460	/ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
461	END DO
462	END DO
463
464	CALL lbc_lnk( hdiv, 'T', 1. ) ! lateral boundary cond. (no sign change)
465
466	DO jj = 2, jpjm1
467	DO ji = fs_2, fs_jpim1 ! vector opt.
468	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) &
469	- e1t(ji ,jj)e2t(ji ,jj) hdiv(ji ,jj) ) &
470	/ e1u(ji,jj) * umask(ji,jj,jk)
471	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) &
472	- e1t(ji,jj )e2t(ji,jj ) hdiv(ji,jj ) ) &
473	/ e2v(ji,jj) * vmask(ji,jj,jk)
474	END DO
475	END DO
476
477	END DO
478
479	END DO
480
481	CALL wrk_dealloc(jpi,jpj,hdiv)
482
483	ENDIF
484
485
486
487	!-----------------------------------------------------------------------
488	! Allocate and initialize the background state arrays
489	!-----------------------------------------------------------------------
490
491	IF ( ln_asmdin ) THEN
492
493	ALLOCATE( t_bkg(jpi,jpj,jpk) )
494	ALLOCATE( s_bkg(jpi,jpj,jpk) )
495	ALLOCATE( u_bkg(jpi,jpj,jpk) )
496	ALLOCATE( v_bkg(jpi,jpj,jpk) )
497	ALLOCATE( ssh_bkg(jpi,jpj) )
498
499	t_bkg(:,:,:) = 0.0
500	s_bkg(:,:,:) = 0.0
501	u_bkg(:,:,:) = 0.0
502	v_bkg(:,:,:) = 0.0
503	ssh_bkg(:,:) = 0.0
504
505	!--------------------------------------------------------------------
506	! Read from file the background state at analysis time
507	!--------------------------------------------------------------------
508
509	CALL iom_open( c_asmdin, inum )
510
511	CALL iom_get( inum, 'rdastp', zdate_bkg )
512
513	IF(lwp) THEN
514	WRITE(numout,*)
515	WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', &
516	& NINT( zdate_bkg )
517	WRITE(numout,*) '~~~~~~~~~~~~'
518	ENDIF
519
520	IF ( NINT( zdate_bkg ) /= iitdin_date ) &
521	& CALL ctl_warn( ' Validity time of assimilation background state does', &
522	& ' not agree with Direct Initialization time' )
523
524	IF ( ln_trainc ) THEN
525	CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
526	CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
527	t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
528	s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
529	ENDIF
530
531	IF ( ln_dyninc ) THEN
532	CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
533	CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
534	u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
535	v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
536	ENDIF
537
538	IF ( ln_sshinc ) THEN
539	CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
540	ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
541	ENDIF
542
543	CALL iom_close( inum )
544
545	ENDIF
546	!
547	END SUBROUTINE asm_inc_init
548
549
550	SUBROUTINE calc_date( kit000, kt, kdate0, kdate )
551	!!----------------------------------------------------------------------
552	!! * ROUTINE calc_date *
553	!!
554	!! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.
555	!!
556	!! ** Method : Compute the calendar date YYYYMMDD at a given time step.
557	!!
558	!! ** Action :
559	!!----------------------------------------------------------------------
560	INTEGER, INTENT(IN) :: kit000 ! Initial time step
561	INTEGER, INTENT(IN) :: kt ! Current time step referenced to kit000
562	INTEGER, INTENT(IN) :: kdate0 ! Initial date
563	INTEGER, INTENT(OUT) :: kdate ! Current date reference to kdate0
564	!
565	INTEGER :: iyea0 ! Initial year
566	INTEGER :: imon0 ! Initial month
567	INTEGER :: iday0 ! Initial day
568	INTEGER :: iyea ! Current year
569	INTEGER :: imon ! Current month
570	INTEGER :: iday ! Current day
571	INTEGER :: idaystp ! Number of days between initial and current date
572	INTEGER :: idaycnt ! Day counter
573
574	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
575
576	!-----------------------------------------------------------------------
577	! Compute the calendar date YYYYMMDD
578	!-----------------------------------------------------------------------
579
580	! Initial date
581	iyea0 = kdate0 / 10000
582	imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100
583	iday0 = kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 )
584
585	! Check that kt >= kit000 - 1
586	IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')
587
588	! If kt = kit000 - 1 then set the date to the restart date
589	IF ( kt == kit000 - 1 ) THEN
590
591	kdate = ndastp
592	RETURN
593
594	ENDIF
595
596	! Compute the number of days from the initial date
597	idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )
598
599	iday = iday0
600	imon = imon0
601	iyea = iyea0
602	idaycnt = 0
603
604	CALL calc_month_len( iyea, imonth_len )
605
606	DO WHILE ( idaycnt < idaystp )
607	iday = iday + 1
608	IF ( iday > imonth_len(imon) ) THEN
609	iday = 1
610	imon = imon + 1
611	ENDIF
612	IF ( imon > 12 ) THEN
613	imon = 1
614	iyea = iyea + 1
615	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
616	ENDIF
617	idaycnt = idaycnt + 1
618	END DO
619	!
620	kdate = iyea * 10000 + imon * 100 + iday
621	!
622	END SUBROUTINE
623
624
625	SUBROUTINE calc_month_len( iyear, imonth_len )
626	!!----------------------------------------------------------------------
627	!! * ROUTINE calc_month_len *
628	!!
629	!! ** Purpose : Compute the number of days in a months given a year.
630	!!
631	!! ** Method :
632	!!----------------------------------------------------------------------
633	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
634	INTEGER :: iyear !: year
635	!!----------------------------------------------------------------------
636	!
637	! length of the month of the current year (from nleapy, read in namelist)
638	IF ( nleapy < 2 ) THEN
639	imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
640	IF ( nleapy == 1 ) THEN ! we are using calendar with leap years
641	IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
642	imonth_len(2) = 29
643	ENDIF
644	ENDIF
645	ELSE
646	imonth_len(:) = nleapy ! all months with nleapy days per year
647	ENDIF
648	!
649	END SUBROUTINE
650
651
652	SUBROUTINE tra_asm_inc( kt )
653	!!----------------------------------------------------------------------
654	!! * ROUTINE tra_asm_inc *
655	!!
656	!! ** Purpose : Apply the tracer (T and S) assimilation increments
657	!!
658	!! ** Method : Direct initialization or Incremental Analysis Updating
659	!!
660	!! ** Action :
661	!!----------------------------------------------------------------------
662	INTEGER, INTENT(IN) :: kt ! Current time step
663	!
664	INTEGER :: ji,jj,jk
665	INTEGER :: it
666	REAL(wp) :: zincwgt ! IAU weight for current time step
667	REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
668	!!----------------------------------------------------------------------
669
670	! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
671	! used to prevent the applied increments taking the temperature below the local freezing point
672
673	DO jk=1, jpkm1
674	fzptnz (:,:,jk) = tfreez( tsn(:,:,jk,jp_sal), fsdept(:,:,jk) )
675	ENDDO
676
677	IF ( ln_asmiau ) THEN
678
679	!--------------------------------------------------------------------
680	! Incremental Analysis Updating
681	!--------------------------------------------------------------------
682
683	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
684
685	it = kt - nit000 + 1
686	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
687
688	IF(lwp) THEN
689	WRITE(numout,*)
690	WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', &
691	& kt,' with IAU weight = ', wgtiau(it)
692	WRITE(numout,*) '~~~~~~~~~~~~'
693	ENDIF
694
695	! Update the tracer tendencies
696	DO jk = 1, jpkm1
697	IF (ln_temnofreeze) THEN
698	! Do not apply negative increments if the temperature will fall below freezing
699	WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
700	& tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )
701	tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt
702	END WHERE
703	ELSE
704	tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt
705	ENDIF
706	IF (ln_salfix) THEN
707	! Do not apply negative increments if the salinity will fall below a specified
708	! minimum value salfixmin
709	WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
710	& tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )
711	tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
712	END WHERE
713	ELSE
714	tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
715	ENDIF
716	END DO
717
718	ENDIF
719
720	IF ( kt == nitiaufin_r + 1 ) THEN ! For bias crcn to work
721	DEALLOCATE( t_bkginc )
722	DEALLOCATE( s_bkginc )
723	ENDIF
724
725
726	ELSEIF ( ln_asmdin ) THEN
727
728	!--------------------------------------------------------------------
729	! Direct Initialization
730	!--------------------------------------------------------------------
731
732	IF ( kt == nitdin_r ) THEN
733
734	neuler = 0 ! Force Euler forward step
735
736	! Initialize the now fields with the background + increment
737	IF (ln_temnofreeze) THEN
738	! Do not apply negative increments if the temperature will fall below freezing
739	WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. &
740	& tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )
741	tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)
742	END WHERE
743	ELSE
744	tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)
745	ENDIF
746	IF (ln_salfix) THEN
747	! Do not apply negative increments if the salinity will fall below a specified
748	! minimum value salfixmin
749	WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. &
750	& tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )
751	tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)
752	END WHERE
753	ELSE
754	tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)
755	ENDIF
756
757	tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields
758
759	CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities
760
761	IF( ln_zps .AND. .NOT. lk_c1d ) &
762	& CALL zps_hde( nit000, jpts, tsb, & ! Partial steps: before horizontal derivative
763	& gtsu, gtsv, rhd, & ! of T, S, rd at the bottom ocean level
764	& gru , grv )
765
766	#if defined key_zdfkpp
767	CALL eos( tsn, rhd, fsdept_n(:,:,:) ) ! Compute rhd
768	#endif
769
770	DEALLOCATE( t_bkginc )
771	DEALLOCATE( s_bkginc )
772	DEALLOCATE( t_bkg )
773	DEALLOCATE( s_bkg )
774	ENDIF
775	!
776	ENDIF
777	! Perhaps the following call should be in step
778	IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment
779	!
780	END SUBROUTINE tra_asm_inc
781
782
783	SUBROUTINE dyn_asm_inc( kt )
784	!!----------------------------------------------------------------------
785	!! * ROUTINE dyn_asm_inc *
786	!!
787	!! ** Purpose : Apply the dynamics (u and v) assimilation increments.
788	!!
789	!! ** Method : Direct initialization or Incremental Analysis Updating.
790	!!
791	!! ** Action :
792	!!----------------------------------------------------------------------
793	INTEGER, INTENT(IN) :: kt ! Current time step
794	!
795	INTEGER :: jk
796	INTEGER :: it
797	REAL(wp) :: zincwgt ! IAU weight for current time step
798	!!----------------------------------------------------------------------
799
800	IF ( ln_asmiau ) THEN
801
802	!--------------------------------------------------------------------
803	! Incremental Analysis Updating
804	!--------------------------------------------------------------------
805
806	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
807
808	it = kt - nit000 + 1
809	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
810
811	IF(lwp) THEN
812	WRITE(numout,*)
813	WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', &
814	& kt,' with IAU weight = ', wgtiau(it)
815	WRITE(numout,*) '~~~~~~~~~~~~'
816	ENDIF
817
818	! Update the dynamic tendencies
819	DO jk = 1, jpkm1
820	ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
821	va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
822	END DO
823
824	IF ( kt == nitiaufin_r ) THEN
825	DEALLOCATE( u_bkginc )
826	DEALLOCATE( v_bkginc )
827	ENDIF
828
829	ENDIF
830
831	ELSEIF ( ln_asmdin ) THEN
832
833	!--------------------------------------------------------------------
834	! Direct Initialization
835	!--------------------------------------------------------------------
836
837	IF ( kt == nitdin_r ) THEN
838
839	neuler = 0 ! Force Euler forward step
840
841	! Initialize the now fields with the background + increment
842	un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
843	vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)
844
845	ub(:,:,:) = un(:,:,:) ! Update before fields
846	vb(:,:,:) = vn(:,:,:)
847
848	DEALLOCATE( u_bkg )
849	DEALLOCATE( v_bkg )
850	DEALLOCATE( u_bkginc )
851	DEALLOCATE( v_bkginc )
852	ENDIF
853	!
854	ENDIF
855	!
856	END SUBROUTINE dyn_asm_inc
857
858
859	SUBROUTINE ssh_asm_inc( kt )
860	!!----------------------------------------------------------------------
861	!! * ROUTINE ssh_asm_inc *
862	!!
863	!! ** Purpose : Apply the sea surface height assimilation increment.
864	!!
865	!! ** Method : Direct initialization or Incremental Analysis Updating.
866	!!
867	!! ** Action :
868	!!----------------------------------------------------------------------
869	INTEGER, INTENT(IN) :: kt ! Current time step
870	!
871	INTEGER :: it
872	INTEGER :: jk
873	REAL(wp) :: zincwgt ! IAU weight for current time step
874	!!----------------------------------------------------------------------
875
876	IF ( ln_asmiau ) THEN
877
878	!--------------------------------------------------------------------
879	! Incremental Analysis Updating
880	!--------------------------------------------------------------------
881
882	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
883
884	it = kt - nit000 + 1
885	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
886
887	IF(lwp) THEN
888	WRITE(numout,*)
889	WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
890	& kt,' with IAU weight = ', wgtiau(it)
891	WRITE(numout,*) '~~~~~~~~~~~~'
892	ENDIF
893
894	! Save the tendency associated with the IAU weighted SSH increment
895	! (applied in dynspg.*)
896	#if defined key_asminc
897	ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
898	#endif
899	IF ( kt == nitiaufin_r ) THEN
900	DEALLOCATE( ssh_bkginc )
901	ENDIF
902
903	ENDIF
904
905	ELSEIF ( ln_asmdin ) THEN
906
907	!--------------------------------------------------------------------
908	! Direct Initialization
909	!--------------------------------------------------------------------
910
911	IF ( kt == nitdin_r ) THEN
912
913	neuler = 0 ! Force Euler forward step
914
915	! Initialize the now fields the background + increment
916	sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)
917
918	! Update before fields
919	sshb(:,:) = sshn(:,:)
920
921	IF( lk_vvl ) THEN
922	DO jk = 1, jpk
923	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
924	END DO
925	ENDIF
926
927	DEALLOCATE( ssh_bkg )
928	DEALLOCATE( ssh_bkginc )
929
930	ENDIF
931	!
932	ENDIF
933	!
934	END SUBROUTINE ssh_asm_inc
935
936
937	SUBROUTINE seaice_asm_inc( kt, kindic )
938	!!----------------------------------------------------------------------
939	!! * ROUTINE seaice_asm_inc *
940	!!
941	!! ** Purpose : Apply the sea ice assimilation increment.
942	!!
943	!! ** Method : Direct initialization or Incremental Analysis Updating.
944	!!
945	!! ** Action :
946	!!
947	!!----------------------------------------------------------------------
948	IMPLICIT NONE
949	!
950	INTEGER, INTENT(in) :: kt ! Current time step
951	INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation
952	!
953	INTEGER :: it
954	REAL(wp) :: zincwgt ! IAU weight for current time step
955	#if defined key_lim2
956	REAL(wp), DIMENSION(jpi,jpj) :: zofrld, zohicif, zseaicendg, zhicifinc ! LIM
957	REAL(wp) :: zhicifmin = 0.5_wp ! ice minimum depth in metres
958	#endif
959	!!----------------------------------------------------------------------
960
961	IF ( ln_asmiau ) THEN
962
963	!--------------------------------------------------------------------
964	! Incremental Analysis Updating
965	!--------------------------------------------------------------------
966
967	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
968
969	it = kt - nit000 + 1
970	zincwgt = wgtiau(it) ! IAU weight for the current time step
971	! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
972
973	IF(lwp) THEN
974	WRITE(numout,*)
975	WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', &
976	& kt,' with IAU weight = ', wgtiau(it)
977	WRITE(numout,*) '~~~~~~~~~~~~'
978	ENDIF
979
980	! Sea-ice : LIM-3 case (to add)
981
982	#if defined key_lim2
983	! Sea-ice : LIM-2 case
984	zofrld (:,:) = frld(:,:)
985	zohicif(:,:) = hicif(:,:)
986	!
987	frld = MIN( MAX( frld (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
988	pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
989	fr_i(:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction
990	!
991	zseaicendg(:,:) = zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied
992	!
993	! Nudge sea ice depth to bring it up to a required minimum depth
994	WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )
995	zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt
996	ELSEWHERE
997	zhicifinc(:,:) = 0.0_wp
998	END WHERE
999	!
1000	! nudge ice depth
1001	hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1002	phicif(:,:) = phicif(:,:) + zhicifinc(:,:)
1003	!
1004	! seaice salinity balancing (to add)
1005	#endif
1006
1007	#if defined key_cice && defined key_asminc
1008	! Sea-ice : CICE case. Pass ice increment tendency into CICE
1009	ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
1010	#endif
1011
1012	IF ( kt == nitiaufin_r ) THEN
1013	DEALLOCATE( seaice_bkginc )
1014	ENDIF
1015
1016	ELSE
1017
1018	#if defined key_cice && defined key_asminc
1019	! Sea-ice : CICE case. Zero ice increment tendency into CICE
1020	ndaice_da(:,:) = 0.0_wp
1021	#endif
1022
1023	ENDIF
1024
1025	ELSEIF ( ln_asmdin ) THEN
1026
1027	!--------------------------------------------------------------------
1028	! Direct Initialization
1029	!--------------------------------------------------------------------
1030
1031	IF ( kt == nitdin_r ) THEN
1032
1033	neuler = 0 ! Force Euler forward step
1034
1035	! Sea-ice : LIM-3 case (to add)
1036
1037	#if defined key_lim2
1038	! Sea-ice : LIM-2 case.
1039	zofrld(:,:)=frld(:,:)
1040	zohicif(:,:)=hicif(:,:)
1041	!
1042	! Initialize the now fields the background + increment
1043	frld (:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
1044	pfrld(:,:) = frld(:,:)
1045	fr_i (:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction
1046	zseaicendg(:,:) = zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied
1047	!
1048	! Nudge sea ice depth to bring it up to a required minimum depth
1049	WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )
1050	zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt
1051	ELSEWHERE
1052	zhicifinc(:,:) = 0.0_wp
1053	END WHERE
1054	!
1055	! nudge ice depth
1056	hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1057	phicif(:,:) = phicif(:,:)
1058	!
1059	! seaice salinity balancing (to add)
1060	#endif
1061
1062	#if defined key_cice && defined key_asminc
1063	! Sea-ice : CICE case. Pass ice increment tendency into CICE
1064	ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
1065	#endif
1066	IF ( .NOT. PRESENT(kindic) ) THEN
1067	DEALLOCATE( seaice_bkginc )
1068	END IF
1069
1070	ELSE
1071
1072	#if defined key_cice && defined key_asminc
1073	! Sea-ice : CICE case. Zero ice increment tendency into CICE
1074	ndaice_da(:,:) = 0.0_wp
1075	#endif
1076
1077	ENDIF
1078
1079	!#if defined defined key_lim2 \|\| defined key_cice
1080	!
1081	! IF (ln_seaicebal ) THEN
1082	! !! balancing salinity increments
1083	! !! simple case from limflx.F90 (doesn't include a mass flux)
1084	! !! assumption is that as ice concentration is reduced or increased
1085	! !! the snow and ice depths remain constant
1086	! !! note that snow is being created where ice concentration is being increased
1087	! !! - could be more sophisticated and
1088	! !! not do this (but would need to alter h_snow)
1089	!
1090	! usave(:,:,:)=sb(:,:,:) ! use array as a temporary store
1091	!
1092	! DO jj = 1, jpj
1093	! DO ji = 1, jpi
1094	! ! calculate change in ice and snow mass per unit area
1095	! ! positive values imply adding salt to the ocean (results from ice formation)
1096	! ! fwf : ice formation and melting
1097	!
1098	! zfons = ( -nfresh_da(ji,jj)soce + nfsalt_da(ji,jj) )rdt
1099	!
1100	! ! change salinity down to mixed layer depth
1101	! mld=hmld_kara(ji,jj)
1102	!
1103	! ! prevent small mld
1104	! ! less than 10m can cause salinity instability
1105	! IF (mld < 10) mld=10
1106	!
1107	! ! set to bottom of a level
1108	! DO jk = jpk-1, 2, -1
1109	! IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN
1110	! mld=gdepw(ji,jj,jk+1)
1111	! jkmax=jk
1112	! ENDIF
1113	! ENDDO
1114	!
1115	! ! avoid applying salinity balancing in shallow water or on land
1116	! !
1117	!
1118	! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
1119	!
1120	! dsal_ocn=0.0_wp
1121	! sal_thresh=5.0_wp ! minimum salinity threshold for salinity balancing
1122	!
1123	! if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
1124	! dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
1125	!
1126	! ! put increments in for levels in the mixed layer
1127	! ! but prevent salinity below a threshold value
1128	!
1129	! DO jk = 1, jkmax
1130	!
1131	! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
1132	! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
1133	! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
1134	! ENDIF
1135	!
1136	! ENDDO
1137	!
1138	! ! ! salt exchanges at the ice/ocean interface
1139	! ! zpmess = zfons / rdt_ice ! rdt_ice is ice timestep
1140	! !
1141	! !! Adjust fsalt. A +ve fsalt means adding salt to ocean
1142	! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt
1143	! !!
1144	! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d)
1145	! !! ! E-P (kg m-2 s-2)
1146	! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2)
1147	! ENDDO !ji
1148	! ENDDO !jj!
1149	!
1150	! ENDIF !ln_seaicebal
1151	!
1152	!#endif
1153
1154	ENDIF
1155
1156	END SUBROUTINE seaice_asm_inc
1157
1158	!!======================================================================
1159	END MODULE asminc

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

Download in other formats: