New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

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

Context Navigation

source: branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 6650

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

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

Download in other formats: