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 trunk/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago
#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
Property svn:keywords set to `Id`
File size: 48.6 KB

Line
1	MODULE asminc
2	!!======================================================================
3	!! * MODULE asminc *
4	!! Assimilation increment : Apply an increment generated by data
5	!! assimilation
6	!!======================================================================
7	!! History : ! 2007-03 (M. Martin) Met Office version
8	!! ! 2007-04 (A. Weaver) calc_date original code
9	!! ! 2007-04 (A. Weaver) Merge with OPAVAR/NEMOVAR
10	!! NEMO 3.3 ! 2010-05 (D. Lea) Update to work with NEMO v3.2
11	!! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init
12	!! 3.4 ! 2012-10 (A. Weaver and K. Mogensen) Fix for direct initialization
13	!!----------------------------------------------------------------------
14
15	!!----------------------------------------------------------------------
16	!! 'key_asminc' : Switch on the assimilation increment interface
17	!!----------------------------------------------------------------------
18	!! asm_inc_init : Initialize the increment arrays and IAU weights
19	!! calc_date : Compute the calendar date YYYYMMDD on a given step
20	!! tra_asm_inc : Apply the tracer (T and S) increments
21	!! dyn_asm_inc : Apply the dynamic (u and v) increments
22	!! ssh_asm_inc : Apply the SSH increment
23	!! seaice_asm_inc : Apply the seaice increment
24	!!----------------------------------------------------------------------
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 ( ( ln_bkgwri ).AND.( ( ln_asmdin ).OR.( ln_asmiau ) ) ) &
240	& CALL ctl_stop( ' ln_bkgwri and either ln_asmdin or ln_asmiau are set to .true.:', &
241	& ' The background state must be written before applying the increments')
242
243	IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
244	& CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :', &
245	& ' Type IAU weighting function is invalid')
246
247	IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
248	& ) &
249	& CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', &
250	& ' The assimilation increments are not applied')
251
252	IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
253	& CALL ctl_stop( ' nitiaustr = nitiaufin :', &
254	& ' IAU interval is of zero length')
255
256	IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
257	& CALL ctl_stop( ' nitiaustr or nitiaufin :', &
258	& ' IAU starting or final time step is outside the cycle interval', &
259	& ' Valid range nit000 to nitend')
260
261	IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
262	& CALL ctl_stop( ' nitbkg :', &
263	& ' Background time step is outside the cycle interval')
264
265	IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
266	& CALL ctl_stop( ' nitdin :', &
267	& ' Background time step for Direct Initialization is outside', &
268	& ' the cycle interval')
269
270	IF ( nstop > 0 ) RETURN ! if there are any errors then go no further
271
272	!--------------------------------------------------------------------
273	! Initialize the Incremental Analysis Updating weighting function
274	!--------------------------------------------------------------------
275
276	IF ( ln_asmiau ) THEN
277
278	ALLOCATE( wgtiau( icycper ) )
279
280	wgtiau(:) = 0.0
281
282	IF ( niaufn == 0 ) THEN
283
284	!---------------------------------------------------------
285	! Constant IAU forcing
286	!---------------------------------------------------------
287
288	DO jt = 1, iiauper
289	wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
290	END DO
291
292	ELSEIF ( niaufn == 1 ) THEN
293
294	!---------------------------------------------------------
295	! Linear hat-like, centred in middle of IAU interval
296	!---------------------------------------------------------
297
298	! Compute the normalization factor
299	znorm = 0.0
300	IF ( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval
301	imid = iiauper / 2
302	DO jt = 1, imid
303	znorm = znorm + REAL( jt )
304	END DO
305	znorm = 2.0 * znorm
306	ELSE ! Odd number of time steps in IAU interval
307	imid = ( iiauper + 1 ) / 2
308	DO jt = 1, imid - 1
309	znorm = znorm + REAL( jt )
310	END DO
311	znorm = 2.0 * znorm + REAL( imid )
312	ENDIF
313	znorm = 1.0 / znorm
314
315	DO jt = 1, imid - 1
316	wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
317	END DO
318	DO jt = imid, iiauper
319	wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
320	END DO
321
322	ENDIF
323
324	! Test that the integral of the weights over the weighting interval equals 1
325	IF(lwp) THEN
326	WRITE(numout,*)
327	WRITE(numout,*) 'asm_inc_init : IAU weights'
328	WRITE(numout,*) '~~~~~~~~~~~~'
329	WRITE(numout,*) ' time step IAU weight'
330	WRITE(numout,*) ' ========= ====================='
331	ztotwgt = 0.0
332	DO jt = 1, icycper
333	ztotwgt = ztotwgt + wgtiau(jt)
334	WRITE(numout,*) ' ', jt, ' ', wgtiau(jt)
335	END DO
336	WRITE(numout,*) ' ==================================='
337	WRITE(numout,*) ' Time-integrated weight = ', ztotwgt
338	WRITE(numout,*) ' ==================================='
339	ENDIF
340
341	ENDIF
342
343	!--------------------------------------------------------------------
344	! Allocate and initialize the increment arrays
345	!--------------------------------------------------------------------
346
347	ALLOCATE( t_bkginc(jpi,jpj,jpk) )
348	ALLOCATE( s_bkginc(jpi,jpj,jpk) )
349	ALLOCATE( u_bkginc(jpi,jpj,jpk) )
350	ALLOCATE( v_bkginc(jpi,jpj,jpk) )
351	ALLOCATE( ssh_bkginc(jpi,jpj) )
352	ALLOCATE( seaice_bkginc(jpi,jpj))
353	#if defined key_asminc
354	ALLOCATE( ssh_iau(jpi,jpj) )
355	#endif
356	t_bkginc(:,:,:) = 0.0
357	s_bkginc(:,:,:) = 0.0
358	u_bkginc(:,:,:) = 0.0
359	v_bkginc(:,:,:) = 0.0
360	ssh_bkginc(:,:) = 0.0
361	seaice_bkginc(:,:) = 0.0
362	#if defined key_asminc
363	ssh_iau(:,:) = 0.0
364	#endif
365	IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN
366
367	!--------------------------------------------------------------------
368	! Read the increments from file
369	!--------------------------------------------------------------------
370
371	CALL iom_open( c_asminc, inum )
372
373	CALL iom_get( inum, 'time', zdate_inc )
374
375	CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
376	CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
377	z_inc_dateb = zdate_inc
378	z_inc_datef = zdate_inc
379
380	IF(lwp) THEN
381	WRITE(numout,*)
382	WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', &
383	& ' between dates ', NINT( z_inc_dateb ),' and ', &
384	& NINT( z_inc_datef )
385	WRITE(numout,*) '~~~~~~~~~~~~'
386	ENDIF
387
388	IF ( ( NINT( z_inc_dateb ) < ndastp ) &
389	& .OR.( NINT( z_inc_datef ) > iitend_date ) ) &
390	& CALL ctl_warn( ' Validity time of assimilation increments is ', &
391	& ' outside the assimilation interval' )
392
393	IF ( ( ln_asmdin ).AND.( NINT( zdate_inc ) /= iitdin_date ) ) &
394	& CALL ctl_warn( ' Validity time of assimilation increments does ', &
395	& ' not agree with Direct Initialization time' )
396
397	IF ( ln_trainc ) THEN
398	CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
399	CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
400	! Apply the masks
401	t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
402	s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
403	! Set missing increments to 0.0 rather than 1e+20
404	! to allow for differences in masks
405	WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
406	WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
407	ENDIF
408
409	IF ( ln_dyninc ) THEN
410	CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )
411	CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )
412	! Apply the masks
413	u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
414	v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
415	! Set missing increments to 0.0 rather than 1e+20
416	! to allow for differences in masks
417	WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
418	WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
419	ENDIF
420
421	IF ( ln_sshinc ) THEN
422	CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
423	! Apply the masks
424	ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
425	! Set missing increments to 0.0 rather than 1e+20
426	! to allow for differences in masks
427	WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
428	ENDIF
429
430	IF ( ln_seaiceinc ) THEN
431	CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 )
432	! Apply the masks
433	seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1)
434	! Set missing increments to 0.0 rather than 1e+20
435	! to allow for differences in masks
436	WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
437	ENDIF
438
439	CALL iom_close( inum )
440
441	ENDIF
442
443	!-----------------------------------------------------------------------
444	! Apply divergence damping filter
445	!-----------------------------------------------------------------------
446
447	IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN
448
449	CALL wrk_alloc(jpi,jpj,hdiv)
450
451	DO jt = 1, nn_divdmp
452
453	DO jk = 1, jpkm1
454
455	hdiv(:,:) = 0._wp
456
457	DO jj = 2, jpjm1
458	DO ji = fs_2, fs_jpim1 ! vector opt.
459	hdiv(ji,jj) = &
460	( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * u_bkginc(ji ,jj ,jk) &
461	- e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * u_bkginc(ji-1,jj ,jk) &
462	+ e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * v_bkginc(ji ,jj ,jk) &
463	- e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * v_bkginc(ji ,jj-1,jk) ) &
464	/ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
465	END DO
466	END DO
467
468	CALL lbc_lnk( hdiv, 'T', 1. ) ! lateral boundary cond. (no sign change)
469
470	DO jj = 2, jpjm1
471	DO ji = fs_2, fs_jpim1 ! vector opt.
472	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) &
473	- e1t(ji ,jj)e2t(ji ,jj) hdiv(ji ,jj) ) &
474	/ e1u(ji,jj) * umask(ji,jj,jk)
475	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) &
476	- e1t(ji,jj )e2t(ji,jj ) hdiv(ji,jj ) ) &
477	/ e2v(ji,jj) * vmask(ji,jj,jk)
478	END DO
479	END DO
480
481	END DO
482
483	END DO
484
485	CALL wrk_dealloc(jpi,jpj,hdiv)
486
487	ENDIF
488
489
490
491	!-----------------------------------------------------------------------
492	! Allocate and initialize the background state arrays
493	!-----------------------------------------------------------------------
494
495	IF ( ln_asmdin ) THEN
496
497	ALLOCATE( t_bkg(jpi,jpj,jpk) )
498	ALLOCATE( s_bkg(jpi,jpj,jpk) )
499	ALLOCATE( u_bkg(jpi,jpj,jpk) )
500	ALLOCATE( v_bkg(jpi,jpj,jpk) )
501	ALLOCATE( ssh_bkg(jpi,jpj) )
502
503	t_bkg(:,:,:) = 0.0
504	s_bkg(:,:,:) = 0.0
505	u_bkg(:,:,:) = 0.0
506	v_bkg(:,:,:) = 0.0
507	ssh_bkg(:,:) = 0.0
508
509	!--------------------------------------------------------------------
510	! Read from file the background state at analysis time
511	!--------------------------------------------------------------------
512
513	CALL iom_open( c_asmdin, inum )
514
515	CALL iom_get( inum, 'rdastp', zdate_bkg )
516
517	IF(lwp) THEN
518	WRITE(numout,*)
519	WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', &
520	& NINT( zdate_bkg )
521	WRITE(numout,*) '~~~~~~~~~~~~'
522	ENDIF
523
524	IF ( NINT( zdate_bkg ) /= iitdin_date ) &
525	& CALL ctl_warn( ' Validity time of assimilation background state does', &
526	& ' not agree with Direct Initialization time' )
527
528	IF ( ln_trainc ) THEN
529	CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
530	CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
531	t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
532	s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
533	ENDIF
534
535	IF ( ln_dyninc ) THEN
536	CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
537	CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
538	u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
539	v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
540	ENDIF
541
542	IF ( ln_sshinc ) THEN
543	CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
544	ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
545	ENDIF
546
547	CALL iom_close( inum )
548
549	ENDIF
550	!
551	END SUBROUTINE asm_inc_init
552
553
554	SUBROUTINE calc_date( kit000, kt, kdate0, kdate )
555	!!----------------------------------------------------------------------
556	!! * ROUTINE calc_date *
557	!!
558	!! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.
559	!!
560	!! ** Method : Compute the calendar date YYYYMMDD at a given time step.
561	!!
562	!! ** Action :
563	!!----------------------------------------------------------------------
564	INTEGER, INTENT(IN) :: kit000 ! Initial time step
565	INTEGER, INTENT(IN) :: kt ! Current time step referenced to kit000
566	INTEGER, INTENT(IN) :: kdate0 ! Initial date
567	INTEGER, INTENT(OUT) :: kdate ! Current date reference to kdate0
568	!
569	INTEGER :: iyea0 ! Initial year
570	INTEGER :: imon0 ! Initial month
571	INTEGER :: iday0 ! Initial day
572	INTEGER :: iyea ! Current year
573	INTEGER :: imon ! Current month
574	INTEGER :: iday ! Current day
575	INTEGER :: idaystp ! Number of days between initial and current date
576	INTEGER :: idaycnt ! Day counter
577
578	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
579
580	!-----------------------------------------------------------------------
581	! Compute the calendar date YYYYMMDD
582	!-----------------------------------------------------------------------
583
584	! Initial date
585	iyea0 = kdate0 / 10000
586	imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100
587	iday0 = kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 )
588
589	! Check that kt >= kit000 - 1
590	IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')
591
592	! If kt = kit000 - 1 then set the date to the restart date
593	IF ( kt == kit000 - 1 ) THEN
594
595	kdate = ndastp
596	RETURN
597
598	ENDIF
599
600	! Compute the number of days from the initial date
601	idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )
602
603	iday = iday0
604	imon = imon0
605	iyea = iyea0
606	idaycnt = 0
607
608	CALL calc_month_len( iyea, imonth_len )
609
610	DO WHILE ( idaycnt < idaystp )
611	iday = iday + 1
612	IF ( iday > imonth_len(imon) ) THEN
613	iday = 1
614	imon = imon + 1
615	ENDIF
616	IF ( imon > 12 ) THEN
617	imon = 1
618	iyea = iyea + 1
619	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
620	ENDIF
621	idaycnt = idaycnt + 1
622	END DO
623	!
624	kdate = iyea * 10000 + imon * 100 + iday
625	!
626	END SUBROUTINE
627
628
629	SUBROUTINE calc_month_len( iyear, imonth_len )
630	!!----------------------------------------------------------------------
631	!! * ROUTINE calc_month_len *
632	!!
633	!! ** Purpose : Compute the number of days in a months given a year.
634	!!
635	!! ** Method :
636	!!----------------------------------------------------------------------
637	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
638	INTEGER :: iyear !: year
639	!!----------------------------------------------------------------------
640	!
641	! length of the month of the current year (from nleapy, read in namelist)
642	IF ( nleapy < 2 ) THEN
643	imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
644	IF ( nleapy == 1 ) THEN ! we are using calendar with leap years
645	IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
646	imonth_len(2) = 29
647	ENDIF
648	ENDIF
649	ELSE
650	imonth_len(:) = nleapy ! all months with nleapy days per year
651	ENDIF
652	!
653	END SUBROUTINE
654
655
656	SUBROUTINE tra_asm_inc( kt )
657	!!----------------------------------------------------------------------
658	!! * ROUTINE tra_asm_inc *
659	!!
660	!! ** Purpose : Apply the tracer (T and S) assimilation increments
661	!!
662	!! ** Method : Direct initialization or Incremental Analysis Updating
663	!!
664	!! ** Action :
665	!!----------------------------------------------------------------------
666	INTEGER, INTENT(IN) :: kt ! Current time step
667	!
668	INTEGER :: ji,jj,jk
669	INTEGER :: it
670	REAL(wp) :: zincwgt ! IAU weight for current time step
671	REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
672	!!----------------------------------------------------------------------
673
674	! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
675	! used to prevent the applied increments taking the temperature below the local freezing point
676
677	DO jk=1, jpkm1
678	fzptnz (:,:,jk) = tfreez( tsn(:,:,jk,jp_sal), fsdept(:,:,jk) )
679	ENDDO
680
681	IF ( ln_asmiau ) THEN
682
683	!--------------------------------------------------------------------
684	! Incremental Analysis Updating
685	!--------------------------------------------------------------------
686
687	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
688
689	it = kt - nit000 + 1
690	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
691
692	IF(lwp) THEN
693	WRITE(numout,*)
694	WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', &
695	& kt,' with IAU weight = ', wgtiau(it)
696	WRITE(numout,*) '~~~~~~~~~~~~'
697	ENDIF
698
699	! Update the tracer tendencies
700	DO jk = 1, jpkm1
701	IF (ln_temnofreeze) THEN
702	! Do not apply negative increments if the temperature will fall below freezing
703	WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
704	& tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )
705	tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt
706	END WHERE
707	ELSE
708	tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt
709	ENDIF
710	IF (ln_salfix) THEN
711	! Do not apply negative increments if the salinity will fall below a specified
712	! minimum value salfixmin
713	WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
714	& tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )
715	tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
716	END WHERE
717	ELSE
718	tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
719	ENDIF
720	END DO
721
722	ENDIF
723
724	IF ( kt == nitiaufin_r + 1 ) THEN ! For bias crcn to work
725	DEALLOCATE( t_bkginc )
726	DEALLOCATE( s_bkginc )
727	ENDIF
728
729
730	ELSEIF ( ln_asmdin ) THEN
731
732	!--------------------------------------------------------------------
733	! Direct Initialization
734	!--------------------------------------------------------------------
735
736	IF ( kt == nitdin_r ) THEN
737
738	neuler = 0 ! Force Euler forward step
739
740	! Initialize the now fields with the background + increment
741	IF (ln_temnofreeze) THEN
742	! Do not apply negative increments if the temperature will fall below freezing
743	WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. &
744	& tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )
745	tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)
746	END WHERE
747	ELSE
748	tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)
749	ENDIF
750	IF (ln_salfix) THEN
751	! Do not apply negative increments if the salinity will fall below a specified
752	! minimum value salfixmin
753	WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. &
754	& tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )
755	tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)
756	END WHERE
757	ELSE
758	tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)
759	ENDIF
760
761	tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields
762
763	CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities
764
765	IF( ln_zps .AND. .NOT. lk_c1d ) &
766	& CALL zps_hde( nit000, jpts, tsb, & ! Partial steps: before horizontal derivative
767	& gtsu, gtsv, rhd, & ! of T, S, rd at the bottom ocean level
768	& gru , grv )
769
770	#if defined key_zdfkpp
771	CALL eos( tsn, rhd, fsdept_n(:,:,:) ) ! Compute rhd
772	#endif
773
774	DEALLOCATE( t_bkginc )
775	DEALLOCATE( s_bkginc )
776	DEALLOCATE( t_bkg )
777	DEALLOCATE( s_bkg )
778	ENDIF
779	!
780	ENDIF
781	! Perhaps the following call should be in step
782	IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment
783	!
784	END SUBROUTINE tra_asm_inc
785
786
787	SUBROUTINE dyn_asm_inc( kt )
788	!!----------------------------------------------------------------------
789	!! * ROUTINE dyn_asm_inc *
790	!!
791	!! ** Purpose : Apply the dynamics (u and v) assimilation increments.
792	!!
793	!! ** Method : Direct initialization or Incremental Analysis Updating.
794	!!
795	!! ** Action :
796	!!----------------------------------------------------------------------
797	INTEGER, INTENT(IN) :: kt ! Current time step
798	!
799	INTEGER :: jk
800	INTEGER :: it
801	REAL(wp) :: zincwgt ! IAU weight for current time step
802	!!----------------------------------------------------------------------
803
804	IF ( ln_asmiau ) THEN
805
806	!--------------------------------------------------------------------
807	! Incremental Analysis Updating
808	!--------------------------------------------------------------------
809
810	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
811
812	it = kt - nit000 + 1
813	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
814
815	IF(lwp) THEN
816	WRITE(numout,*)
817	WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', &
818	& kt,' with IAU weight = ', wgtiau(it)
819	WRITE(numout,*) '~~~~~~~~~~~~'
820	ENDIF
821
822	! Update the dynamic tendencies
823	DO jk = 1, jpkm1
824	ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
825	va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
826	END DO
827
828	IF ( kt == nitiaufin_r ) THEN
829	DEALLOCATE( u_bkginc )
830	DEALLOCATE( v_bkginc )
831	ENDIF
832
833	ENDIF
834
835	ELSEIF ( ln_asmdin ) THEN
836
837	!--------------------------------------------------------------------
838	! Direct Initialization
839	!--------------------------------------------------------------------
840
841	IF ( kt == nitdin_r ) THEN
842
843	neuler = 0 ! Force Euler forward step
844
845	! Initialize the now fields with the background + increment
846	un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
847	vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)
848
849	ub(:,:,:) = un(:,:,:) ! Update before fields
850	vb(:,:,:) = vn(:,:,:)
851
852	DEALLOCATE( u_bkg )
853	DEALLOCATE( v_bkg )
854	DEALLOCATE( u_bkginc )
855	DEALLOCATE( v_bkginc )
856	ENDIF
857	!
858	ENDIF
859	!
860	END SUBROUTINE dyn_asm_inc
861
862
863	SUBROUTINE ssh_asm_inc( kt )
864	!!----------------------------------------------------------------------
865	!! * ROUTINE ssh_asm_inc *
866	!!
867	!! ** Purpose : Apply the sea surface height assimilation increment.
868	!!
869	!! ** Method : Direct initialization or Incremental Analysis Updating.
870	!!
871	!! ** Action :
872	!!----------------------------------------------------------------------
873	INTEGER, INTENT(IN) :: kt ! Current time step
874	!
875	INTEGER :: it
876	INTEGER :: jk
877	REAL(wp) :: zincwgt ! IAU weight for current time step
878	!!----------------------------------------------------------------------
879
880	IF ( ln_asmiau ) THEN
881
882	!--------------------------------------------------------------------
883	! Incremental Analysis Updating
884	!--------------------------------------------------------------------
885
886	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
887
888	it = kt - nit000 + 1
889	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
890
891	IF(lwp) THEN
892	WRITE(numout,*)
893	WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
894	& kt,' with IAU weight = ', wgtiau(it)
895	WRITE(numout,*) '~~~~~~~~~~~~'
896	ENDIF
897
898	! Save the tendency associated with the IAU weighted SSH increment
899	! (applied in dynspg.*)
900	#if defined key_asminc
901	ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
902	#endif
903	IF ( kt == nitiaufin_r ) THEN
904	DEALLOCATE( ssh_bkginc )
905	ENDIF
906
907	ENDIF
908
909	ELSEIF ( ln_asmdin ) THEN
910
911	!--------------------------------------------------------------------
912	! Direct Initialization
913	!--------------------------------------------------------------------
914
915	IF ( kt == nitdin_r ) THEN
916
917	neuler = 0 ! Force Euler forward step
918
919	! Initialize the now fields the background + increment
920	sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)
921
922	! Update before fields
923	sshb(:,:) = sshn(:,:)
924
925	IF( lk_vvl ) THEN
926	DO jk = 1, jpk
927	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
928	END DO
929	ENDIF
930
931	DEALLOCATE( ssh_bkg )
932	DEALLOCATE( ssh_bkginc )
933
934	ENDIF
935	!
936	ENDIF
937	!
938	END SUBROUTINE ssh_asm_inc
939
940
941	SUBROUTINE seaice_asm_inc( kt, kindic )
942	!!----------------------------------------------------------------------
943	!! * ROUTINE seaice_asm_inc *
944	!!
945	!! ** Purpose : Apply the sea ice assimilation increment.
946	!!
947	!! ** Method : Direct initialization or Incremental Analysis Updating.
948	!!
949	!! ** Action :
950	!!
951	!!----------------------------------------------------------------------
952	IMPLICIT NONE
953	!
954	INTEGER, INTENT(in) :: kt ! Current time step
955	INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation
956	!
957	INTEGER :: it
958	REAL(wp) :: zincwgt ! IAU weight for current time step
959	#if defined key_lim2
960	REAL(wp), DIMENSION(jpi,jpj) :: zofrld, zohicif, zseaicendg, zhicifinc ! LIM
961	REAL(wp) :: zhicifmin = 0.5_wp ! ice minimum depth in metres
962	#endif
963	!!----------------------------------------------------------------------
964
965	IF ( ln_asmiau ) THEN
966
967	!--------------------------------------------------------------------
968	! Incremental Analysis Updating
969	!--------------------------------------------------------------------
970
971	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
972
973	it = kt - nit000 + 1
974	zincwgt = wgtiau(it) ! IAU weight for the current time step
975	! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
976
977	IF(lwp) THEN
978	WRITE(numout,*)
979	WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', &
980	& kt,' with IAU weight = ', wgtiau(it)
981	WRITE(numout,*) '~~~~~~~~~~~~'
982	ENDIF
983
984	! Sea-ice : LIM-3 case (to add)
985
986	#if defined key_lim2
987	! Sea-ice : LIM-2 case
988	zofrld (:,:) = frld(:,:)
989	zohicif(:,:) = hicif(:,:)
990	!
991	frld = MIN( MAX( frld (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
992	pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
993	fr_i(:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction
994	!
995	zseaicendg(:,:) = zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied
996	!
997	! Nudge sea ice depth to bring it up to a required minimum depth
998	WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )
999	zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt
1000	ELSEWHERE
1001	zhicifinc(:,:) = 0.0_wp
1002	END WHERE
1003	!
1004	! nudge ice depth
1005	hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1006	phicif(:,:) = phicif(:,:) + zhicifinc(:,:)
1007	!
1008	! seaice salinity balancing (to add)
1009	#endif
1010
1011	#if defined key_cice && defined key_asminc
1012	! Sea-ice : CICE case. Pass ice increment tendency into CICE
1013	ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
1014	#endif
1015
1016	IF ( kt == nitiaufin_r ) THEN
1017	DEALLOCATE( seaice_bkginc )
1018	ENDIF
1019
1020	ELSE
1021
1022	#if defined key_cice && defined key_asminc
1023	! Sea-ice : CICE case. Zero ice increment tendency into CICE
1024	ndaice_da(:,:) = 0.0_wp
1025	#endif
1026
1027	ENDIF
1028
1029	ELSEIF ( ln_asmdin ) THEN
1030
1031	!--------------------------------------------------------------------
1032	! Direct Initialization
1033	!--------------------------------------------------------------------
1034
1035	IF ( kt == nitdin_r ) THEN
1036
1037	neuler = 0 ! Force Euler forward step
1038
1039	! Sea-ice : LIM-3 case (to add)
1040
1041	#if defined key_lim2
1042	! Sea-ice : LIM-2 case.
1043	zofrld(:,:)=frld(:,:)
1044	zohicif(:,:)=hicif(:,:)
1045	!
1046	! Initialize the now fields the background + increment
1047	frld (:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
1048	pfrld(:,:) = frld(:,:)
1049	fr_i (:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction
1050	zseaicendg(:,:) = zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied
1051	!
1052	! Nudge sea ice depth to bring it up to a required minimum depth
1053	WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )
1054	zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt
1055	ELSEWHERE
1056	zhicifinc(:,:) = 0.0_wp
1057	END WHERE
1058	!
1059	! nudge ice depth
1060	hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1061	phicif(:,:) = phicif(:,:)
1062	!
1063	! seaice salinity balancing (to add)
1064	#endif
1065
1066	#if defined key_cice && defined key_asminc
1067	! Sea-ice : CICE case. Pass ice increment tendency into CICE
1068	ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
1069	#endif
1070	IF ( .NOT. PRESENT(kindic) ) THEN
1071	DEALLOCATE( seaice_bkginc )
1072	END IF
1073
1074	ELSE
1075
1076	#if defined key_cice && defined key_asminc
1077	! Sea-ice : CICE case. Zero ice increment tendency into CICE
1078	ndaice_da(:,:) = 0.0_wp
1079	#endif
1080
1081	ENDIF
1082
1083	!#if defined defined key_lim2 \|\| defined key_cice
1084	!
1085	! IF (ln_seaicebal ) THEN
1086	! !! balancing salinity increments
1087	! !! simple case from limflx.F90 (doesn't include a mass flux)
1088	! !! assumption is that as ice concentration is reduced or increased
1089	! !! the snow and ice depths remain constant
1090	! !! note that snow is being created where ice concentration is being increased
1091	! !! - could be more sophisticated and
1092	! !! not do this (but would need to alter h_snow)
1093	!
1094	! usave(:,:,:)=sb(:,:,:) ! use array as a temporary store
1095	!
1096	! DO jj = 1, jpj
1097	! DO ji = 1, jpi
1098	! ! calculate change in ice and snow mass per unit area
1099	! ! positive values imply adding salt to the ocean (results from ice formation)
1100	! ! fwf : ice formation and melting
1101	!
1102	! zfons = ( -nfresh_da(ji,jj)soce + nfsalt_da(ji,jj) )rdt
1103	!
1104	! ! change salinity down to mixed layer depth
1105	! mld=hmld_kara(ji,jj)
1106	!
1107	! ! prevent small mld
1108	! ! less than 10m can cause salinity instability
1109	! IF (mld < 10) mld=10
1110	!
1111	! ! set to bottom of a level
1112	! DO jk = jpk-1, 2, -1
1113	! IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN
1114	! mld=gdepw(ji,jj,jk+1)
1115	! jkmax=jk
1116	! ENDIF
1117	! ENDDO
1118	!
1119	! ! avoid applying salinity balancing in shallow water or on land
1120	! !
1121	!
1122	! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
1123	!
1124	! dsal_ocn=0.0_wp
1125	! sal_thresh=5.0_wp ! minimum salinity threshold for salinity balancing
1126	!
1127	! if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
1128	! dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
1129	!
1130	! ! put increments in for levels in the mixed layer
1131	! ! but prevent salinity below a threshold value
1132	!
1133	! DO jk = 1, jkmax
1134	!
1135	! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
1136	! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
1137	! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
1138	! ENDIF
1139	!
1140	! ENDDO
1141	!
1142	! ! ! salt exchanges at the ice/ocean interface
1143	! ! zpmess = zfons / rdt_ice ! rdt_ice is ice timestep
1144	! !
1145	! !! Adjust fsalt. A +ve fsalt means adding salt to ocean
1146	! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt
1147	! !!
1148	! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d)
1149	! !! ! E-P (kg m-2 s-2)
1150	! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2)
1151	! ENDDO !ji
1152	! ENDDO !jj!
1153	!
1154	! ENDIF !ln_seaicebal
1155	!
1156	!#endif
1157
1158	ENDIF
1159
1160	END SUBROUTINE seaice_asm_inc
1161
1162	!!======================================================================
1163	END MODULE asminc

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

Download in other formats: