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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 2392

Last change on this file since 2392 was 2392, checked in by gm, 13 years ago
v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections
Property svn:keywords set to `Id`
File size: 35.2 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	!!----------------------------------------------------------------------
13
14	!!----------------------------------------------------------------------
15	!! 'key_asminc' : Switch on the assimilation increment interface
16	!!----------------------------------------------------------------------
17	!! asm_inc_init : Initialize the increment arrays and IAU weights
18	!! calc_date : Compute the calendar date YYYYMMDD on a given step
19	!! tra_asm_inc : Apply the tracer (T and S) increments
20	!! dyn_asm_inc : Apply the dynamic (u and v) increments
21	!! ssh_asm_inc : Apply the SSH increment
22	!!----------------------------------------------------------------------
23	USE in_out_manager ! I/O manager
24	USE par_oce ! Ocean space and time domain variables
25	USE dom_oce ! Ocean space and time domain
26	USE oce ! Dynamics and active tracers defined in memory
27	USE divcur ! Horizontal divergence and relative vorticity
28	USE eosbn2 ! Equation of state - in situ and potential density
29	USE zpshde ! Partial step : Horizontal Derivative
30	USE iom ! Library to read input files
31	USE asmpar ! Parameters for the assmilation interface
32	#if defined key_c1d
33	USE c1d, ONLY : lk_c1d ! 1D initialization
34	#endif
35
36	IMPLICIT NONE
37	PRIVATE
38
39	PUBLIC asm_inc_init !: Initialize the increment arrays and IAU weights
40	PUBLIC calc_date !: Compute the calendar date YYYYMMDD on a given step
41	PUBLIC tra_asm_inc !: Apply the tracer (T and S) increments
42	PUBLIC dyn_asm_inc !: Apply the dynamic (u and v) increments
43	PUBLIC ssh_asm_inc !: Apply the SSH increment
44
45	#if defined key_asminc
46	LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE. !: Logical switch for assimilation increment interface
47	#else
48	LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE. !: No assimilation increments
49	#endif
50	LOGICAL, PUBLIC :: ln_bkgwri = .FALSE. !: No output of the background state fields
51	LOGICAL, PUBLIC :: ln_trjwri = .FALSE. !: No output of the state trajectory fields
52	LOGICAL, PUBLIC :: ln_asmiau = .FALSE. !: No applying forcing with an assimilation increment
53	LOGICAL, PUBLIC :: ln_asmdin = .FALSE. !: No direct initialization
54	LOGICAL, PUBLIC :: ln_trainc = .FALSE. !: No tracer (T and S) assimilation increments
55	LOGICAL, PUBLIC :: ln_dyninc = .FALSE. !: No dynamics (u and v) assimilation increments
56	LOGICAL, PUBLIC :: ln_sshinc = .FALSE. !: No sea surface height assimilation increment
57	LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check
58
59	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkg , s_bkg !: Background temperature and salinity
60	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkg , v_bkg !: Background u- & v- velocity components
61	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkginc, s_bkginc !: Increment to the background T & S
62	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components
63	REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step
64	#if defined key_asminc
65	REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssh_iau !: IAU-weighted sea surface height increment
66	#endif
67	! !!! time steps relative to the cycle interval [0,nitend-nit000-1]
68	INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term
69	INTEGER , PUBLIC :: nitdin !: Time step of the background state for direct initialization
70	INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval
71	INTEGER , PUBLIC :: nitiaufin !: Time step of the end of the IAU interval
72	!
73	INTEGER , PUBLIC :: niaufn !: Type of IAU weighing function: = 0 Constant weighting
74	! !: = 1 Linear hat-like, centred in middle of IAU interval
75	REAL(wp), PUBLIC :: salfixmin !: Ensure that the salinity is larger than this value if (ln_salfix)
76
77	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc ! Background sea surface height and its increment
78
79	!!----------------------------------------------------------------------
80	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
81	!! $Id$
82	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
83	!!----------------------------------------------------------------------
84	CONTAINS
85
86	SUBROUTINE asm_inc_init
87	!!----------------------------------------------------------------------
88	!! * ROUTINE asm_inc_init *
89	!!
90	!! ** Purpose : Initialize the assimilation increment and IAU weights.
91	!!
92	!! ** Method : Initialize the assimilation increment and IAU weights.
93	!!
94	!! ** Action :
95	!!----------------------------------------------------------------------
96	INTEGER :: jt
97	INTEGER :: imid
98	INTEGER :: inum
99	INTEGER :: iiauper ! Number of time steps in the IAU period
100	INTEGER :: icycper ! Number of time steps in the cycle
101	INTEGER :: iitend_date ! Date YYYYMMDD of final time step
102	INTEGER :: iitbkg_date ! Date YYYYMMDD of background time step for Jb term
103	INTEGER :: iitdin_date ! Date YYYYMMDD of background time step for DI
104	INTEGER :: iitiaustr_date ! Date YYYYMMDD of IAU interval start time step
105	INTEGER :: iitiaufin_date ! Date YYYYMMDD of IAU interval final time step
106
107	REAL(wp) :: znorm ! Normalization factor for IAU weights
108	REAL(wp) :: ztotwgt ! Value of time-integrated IAU weights
109	! (should be equal to one)
110	REAL(wp) :: z_inc_dateb ! Start date of interval on which increment is valid
111	REAL(wp) :: z_inc_datef ! End date of interval on which increment is valid
112	REAL(wp) :: zdate_bkg ! Date in background state file for DI
113	REAL(wp) :: zdate_inc ! Time axis in increments file
114	!!
115	NAMELIST/nam_asminc/ ln_bkgwri, ln_trjwri, &
116	& ln_trainc, ln_dyninc, ln_sshinc, &
117	& ln_asmdin, ln_asmiau, &
118	& nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, &
119	& nittrjfrq, ln_salfix, salfixmin
120	!!----------------------------------------------------------------------
121
122	!-----------------------------------------------------------------------
123	! Read Namelist nam_asminc : assimilation increment interface
124	!-----------------------------------------------------------------------
125
126	! Set default values
127	ln_bkgwri = .FALSE.
128	ln_trjwri = .FALSE.
129	ln_trainc = .FALSE.
130	ln_dyninc = .FALSE.
131	ln_sshinc = .FALSE.
132	ln_asmdin = .FALSE.
133	ln_asmiau = .TRUE.
134	ln_salfix = .FALSE.
135	salfixmin = -9999
136	nitbkg = 0
137	nitdin = 0
138	nitiaustr = 1
139	nitiaufin = 150 ! = 10 days with ORCA2
140	niaufn = 0
141	nittrjfrq = 1
142
143	REWIND ( numnam )
144	READ ( numnam, nam_asminc )
145
146	! Control print
147	IF(lwp) THEN
148	WRITE(numout,*)
149	WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
150	WRITE(numout,*) '~~~~~~~~~~~~'
151	WRITE(numout,*) ' Namelist namasm : set assimilation increment parameters'
152	WRITE(numout,*) ' Logical switch for writing out background state ln_bkgwri = ', ln_bkgwri
153	WRITE(numout,*) ' Logical switch for writing out state trajectory ln_trjwri = ', ln_trjwri
154	WRITE(numout,*) ' Logical switch for applying tracer increments ln_trainc = ', ln_trainc
155	WRITE(numout,*) ' Logical switch for applying velocity increments ln_dyninc = ', ln_dyninc
156	WRITE(numout,*) ' Logical switch for applying SSH increments ln_sshinc = ', ln_sshinc
157	WRITE(numout,*) ' Logical switch for Direct Initialization (DI) ln_asmdin = ', ln_asmdin
158	WRITE(numout,*) ' Logical switch for Incremental Analysis Updating (IAU) ln_asmiau = ', ln_asmiau
159	WRITE(numout,*) ' Timestep of background in [0,nitend-nit000-1] nitbkg = ', nitbkg
160	WRITE(numout,*) ' Timestep of background for DI in [0,nitend-nit000-1] nitdin = ', nitdin
161	WRITE(numout,*) ' Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
162	WRITE(numout,*) ' Timestep of end of IAU interval in [0,nitend-nit000-1] nitiaufin = ', nitiaufin
163	WRITE(numout,*) ' Type of IAU weighting function niaufn = ', niaufn
164	WRITE(numout,*) ' Frequency of trajectory output for 4D-VAR nittrjfrq = ', nittrjfrq
165	WRITE(numout,*) ' Logical switch for ensuring that the sa > salfixmin ln_salfix = ', ln_salfix
166	WRITE(numout,*) ' Minimum salinity after applying the increments salfixmin = ', salfixmin
167	ENDIF
168
169	nitbkg_r = nitbkg + nit000 - 1 ! Background time referenced to nit000
170	nitdin_r = nitdin + nit000 - 1 ! Background time for DI referenced to nit000
171	nitiaustr_r = nitiaustr + nit000 - 1 ! Start of IAU interval referenced to nit000
172	nitiaufin_r = nitiaufin + nit000 - 1 ! End of IAU interval referenced to nit000
173
174	iiauper = nitiaufin_r - nitiaustr_r + 1 ! IAU interval length
175	icycper = nitend - nit000 + 1 ! Cycle interval length
176
177	! Date of final time step
178	CALL calc_date( nit000, nitend, ndate0, iitend_date )
179
180	! Background time for Jb referenced to ndate0
181	CALL calc_date( nit000, nitbkg_r, ndate0, iitbkg_date )
182
183	! Background time for DI referenced to ndate0
184	CALL calc_date( nit000, nitdin_r, ndate0, iitdin_date )
185
186	! IAU start time referenced to ndate0
187	CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date )
188
189	! IAU end time referenced to ndate0
190	CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date )
191
192	IF(lwp) THEN
193	WRITE(numout,*)
194	WRITE(numout,*) ' Time steps referenced to current cycle:'
195	WRITE(numout,*) ' iitrst = ', nit000 - 1
196	WRITE(numout,*) ' nit000 = ', nit000
197	WRITE(numout,*) ' nitend = ', nitend
198	WRITE(numout,*) ' nitbkg_r = ', nitbkg_r
199	WRITE(numout,*) ' nitdin_r = ', nitdin_r
200	WRITE(numout,*) ' nitiaustr_r = ', nitiaustr_r
201	WRITE(numout,*) ' nitiaufin_r = ', nitiaufin_r
202	WRITE(numout,*) ' nittrjfrq = ', nittrjfrq
203	WRITE(numout,*)
204	WRITE(numout,*) ' Dates referenced to current cycle:'
205	WRITE(numout,*) ' ndastp = ', ndastp
206	WRITE(numout,*) ' ndate0 = ', ndate0
207	WRITE(numout,*) ' iitend_date = ', iitend_date
208	WRITE(numout,*) ' iitbkg_date = ', iitbkg_date
209	WRITE(numout,*) ' iitdin_date = ', iitdin_date
210	WRITE(numout,*) ' iitiaustr_date = ', iitiaustr_date
211	WRITE(numout,*) ' iitiaufin_date = ', iitiaufin_date
212	ENDIF
213
214	IF ( nacc /= 0 ) &
215	& CALL ctl_stop( ' nacc /= 0 and key_asminc :', &
216	& ' Assimilation increments have only been implemented', &
217	& ' for synchronous time stepping' )
218
219	IF ( ( ln_asmdin ).AND.( ln_asmiau ) ) &
220	& CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
221	& ' Choose Direct Initialization OR Incremental Analysis Updating')
222
223	IF ( ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
224	.AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) ) ) &
225	& CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc and ln_sshinc is set to .true.', &
226	& ' but ln_asmdin and ln_asmiau are both set to .false. :', &
227	& ' Inconsistent options')
228
229	IF ( ( ln_bkgwri ).AND.( ( ln_asmdin ).OR.( ln_asmiau ) ) ) &
230	& CALL ctl_stop( ' ln_bkgwri and either ln_asmdin or ln_asmiau are set to .true.:', &
231	& ' The background state must be written before applying the increments')
232
233	IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
234	& CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :', &
235	& ' Type IAU weighting function is invalid')
236
237	IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ) &
238	& ) &
239	& CALL ctl_warn( ' ln_trainc, ln_dyninc and ln_sshinc are set to .false. :', &
240	& ' The assimilation increments are not applied')
241
242	IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
243	& CALL ctl_stop( ' nitiaustr = nitiaufin :', &
244	& ' IAU interval is of zero length')
245
246	IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
247	& CALL ctl_stop( ' nitiaustr or nitiaufin :', &
248	& ' IAU starting or final time step is outside the cycle interval', &
249	& ' Valid range nit000 to nitend')
250
251	IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
252	& CALL ctl_stop( ' nitbkg :', &
253	& ' Background time step is outside the cycle interval')
254
255	IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
256	& CALL ctl_stop( ' nitdin :', &
257	& ' Background time step for Direct Initialization is outside', &
258	& ' the cycle interval')
259
260	IF ( nstop > 0 ) RETURN ! if there are any errors then go no further
261
262	!--------------------------------------------------------------------
263	! Initialize the Incremental Analysis Updating weighting function
264	!--------------------------------------------------------------------
265
266	IF ( ln_asmiau ) THEN
267
268	ALLOCATE( wgtiau( icycper ) )
269
270	wgtiau(:) = 0.0
271
272	IF ( niaufn == 0 ) THEN
273
274	!---------------------------------------------------------
275	! Constant IAU forcing
276	!---------------------------------------------------------
277
278	DO jt = 1, iiauper
279	wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
280	END DO
281
282	ELSEIF ( niaufn == 1 ) THEN
283
284	!---------------------------------------------------------
285	! Linear hat-like, centred in middle of IAU interval
286	!---------------------------------------------------------
287
288	! Compute the normalization factor
289	znorm = 0.0
290	IF ( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval
291	imid = iiauper / 2
292	DO jt = 1, imid
293	znorm = znorm + REAL( jt )
294	END DO
295	znorm = 2.0 * znorm
296	ELSE ! Odd number of time steps in IAU interval
297	imid = ( iiauper + 1 ) / 2
298	DO jt = 1, imid - 1
299	znorm = znorm + REAL( jt )
300	END DO
301	znorm = 2.0 * znorm + REAL( imid )
302	ENDIF
303	znorm = 1.0 / znorm
304
305	DO jt = 1, imid - 1
306	wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
307	END DO
308	DO jt = imid, iiauper
309	wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
310	END DO
311
312	ENDIF
313
314	! Test that the integral of the weights over the weighting interval equals 1
315	IF(lwp) THEN
316	WRITE(numout,*)
317	WRITE(numout,*) 'asm_inc_init : IAU weights'
318	WRITE(numout,*) '~~~~~~~~~~~~'
319	WRITE(numout,*) ' time step IAU weight'
320	WRITE(numout,*) ' ========= ====================='
321	ztotwgt = 0.0
322	DO jt = 1, icycper
323	ztotwgt = ztotwgt + wgtiau(jt)
324	WRITE(numout,*) ' ', jt, ' ', wgtiau(jt)
325	END DO
326	WRITE(numout,*) ' ==================================='
327	WRITE(numout,*) ' Time-integrated weight = ', ztotwgt
328	WRITE(numout,*) ' ==================================='
329	ENDIF
330
331	ENDIF
332
333	!--------------------------------------------------------------------
334	! Allocate and initialize the increment arrays
335	!--------------------------------------------------------------------
336
337	ALLOCATE( t_bkginc(jpi,jpj,jpk) )
338	ALLOCATE( s_bkginc(jpi,jpj,jpk) )
339	ALLOCATE( u_bkginc(jpi,jpj,jpk) )
340	ALLOCATE( v_bkginc(jpi,jpj,jpk) )
341	ALLOCATE( ssh_bkginc(jpi,jpj) )
342	#if defined key_asminc
343	ALLOCATE( ssh_iau(jpi,jpj) )
344	#endif
345	t_bkginc(:,:,:) = 0.0
346	s_bkginc(:,:,:) = 0.0
347	u_bkginc(:,:,:) = 0.0
348	v_bkginc(:,:,:) = 0.0
349	ssh_bkginc(:,:) = 0.0
350	#if defined key_asminc
351	ssh_iau(:,:) = 0.0
352	#endif
353	IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) ) THEN
354
355	!--------------------------------------------------------------------
356	! Read the increments from file
357	!--------------------------------------------------------------------
358
359	CALL iom_open( c_asminc, inum )
360
361	CALL iom_get( inum, 'time', zdate_inc )
362
363	CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
364	CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
365	z_inc_dateb = zdate_inc
366	z_inc_datef = zdate_inc
367
368	IF(lwp) THEN
369	WRITE(numout,*)
370	WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', &
371	& ' between dates ', NINT( z_inc_dateb ),' and ', &
372	& NINT( z_inc_datef )
373	WRITE(numout,*) '~~~~~~~~~~~~'
374	ENDIF
375
376	IF ( ( NINT( z_inc_dateb ) < ndastp ) &
377	& .OR.( NINT( z_inc_datef ) > iitend_date ) ) &
378	& CALL ctl_warn( ' Validity time of assimilation increments is ', &
379	& ' outside the assimilation interval' )
380
381	IF ( ( ln_asmdin ).AND.( NINT( zdate_inc ) /= iitdin_date ) ) &
382	& CALL ctl_warn( ' Validity time of assimilation increments does ', &
383	& ' not agree with Direct Initialization time' )
384
385	IF ( ln_trainc ) THEN
386	CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
387	CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
388	! Apply the masks
389	t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
390	s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
391	! Set missing increments to 0.0 rather than 1e+20
392	! to allow for differences in masks
393	WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
394	WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
395	ENDIF
396
397	IF ( ln_dyninc ) THEN
398	CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )
399	CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )
400	! Apply the masks
401	u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
402	v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
403	! Set missing increments to 0.0 rather than 1e+20
404	! to allow for differences in masks
405	WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
406	WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
407	ENDIF
408
409	IF ( ln_sshinc ) THEN
410	CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
411	! Apply the masks
412	ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
413	! Set missing increments to 0.0 rather than 1e+20
414	! to allow for differences in masks
415	WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
416	ENDIF
417
418	CALL iom_close( inum )
419
420	ENDIF
421
422	!-----------------------------------------------------------------------
423	! Allocate and initialize the background state arrays
424	!-----------------------------------------------------------------------
425
426	IF ( ln_asmdin ) THEN
427
428	ALLOCATE( t_bkg(jpi,jpj,jpk) )
429	ALLOCATE( s_bkg(jpi,jpj,jpk) )
430	ALLOCATE( u_bkg(jpi,jpj,jpk) )
431	ALLOCATE( v_bkg(jpi,jpj,jpk) )
432	ALLOCATE( ssh_bkg(jpi,jpj) )
433
434	t_bkg(:,:,:) = 0.0
435	s_bkg(:,:,:) = 0.0
436	u_bkg(:,:,:) = 0.0
437	v_bkg(:,:,:) = 0.0
438	ssh_bkg(:,:) = 0.0
439
440	!--------------------------------------------------------------------
441	! Read from file the background state at analysis time
442	!--------------------------------------------------------------------
443
444	CALL iom_open( c_asmdin, inum )
445
446	CALL iom_get( inum, 'zdate', zdate_bkg )
447
448	IF(lwp) THEN
449	WRITE(numout,*)
450	WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', &
451	& NINT( zdate_bkg )
452	WRITE(numout,*) '~~~~~~~~~~~~'
453	ENDIF
454
455	IF ( NINT( zdate_bkg ) /= iitdin_date ) &
456	& CALL ctl_warn( ' Validity time of assimilation background state does', &
457	& ' not agree with Direct Initialization time' )
458
459	IF ( ln_trainc ) THEN
460	CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
461	CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
462	t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
463	s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
464	ENDIF
465
466	IF ( ln_dyninc ) THEN
467	CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
468	CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
469	u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
470	v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
471	ENDIF
472
473	IF ( ln_sshinc ) THEN
474	CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
475	ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
476	ENDIF
477
478	CALL iom_close( inum )
479
480	ENDIF
481	!
482	END SUBROUTINE asm_inc_init
483
484
485	SUBROUTINE calc_date( kit000, kt, kdate0, kdate )
486	!!----------------------------------------------------------------------
487	!! * ROUTINE calc_date *
488	!!
489	!! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.
490	!!
491	!! ** Method : Compute the calendar date YYYYMMDD at a given time step.
492	!!
493	!! ** Action :
494	!!----------------------------------------------------------------------
495	INTEGER, INTENT(IN) :: kit000 ! Initial time step
496	INTEGER, INTENT(IN) :: kt ! Current time step referenced to kit000
497	INTEGER, INTENT(IN) :: kdate0 ! Initial date
498	INTEGER, INTENT(OUT) :: kdate ! Current date reference to kdate0
499	!
500	INTEGER :: iyea0 ! Initial year
501	INTEGER :: imon0 ! Initial month
502	INTEGER :: iday0 ! Initial day
503	INTEGER :: iyea ! Current year
504	INTEGER :: imon ! Current month
505	INTEGER :: iday ! Current day
506	INTEGER :: idaystp ! Number of days between initial and current date
507	INTEGER :: idaycnt ! Day counter
508
509	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
510
511	!-----------------------------------------------------------------------
512	! Compute the calendar date YYYYMMDD
513	!-----------------------------------------------------------------------
514
515	! Initial date
516	iyea0 = kdate0 / 10000
517	imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100
518	iday0 = kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 )
519
520	! Check that kt >= kit000 - 1
521	IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')
522
523	! If kt = kit000 - 1 then set the date to the restart date
524	IF ( kt == kit000 - 1 ) THEN
525
526	kdate = ndastp
527	RETURN
528
529	ENDIF
530
531	! Compute the number of days from the initial date
532	idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )
533
534	iday = iday0
535	imon = imon0
536	iyea = iyea0
537	idaycnt = 0
538
539	CALL calc_month_len( iyea, imonth_len )
540
541	DO WHILE ( idaycnt < idaystp )
542	iday = iday + 1
543	IF ( iday > imonth_len(imon) ) THEN
544	iday = 1
545	imon = imon + 1
546	ENDIF
547	IF ( imon > 12 ) THEN
548	imon = 1
549	iyea = iyea + 1
550	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
551	ENDIF
552	idaycnt = idaycnt + 1
553	END DO
554	!
555	kdate = iyea * 10000 + imon * 100 + iday
556	!
557	END SUBROUTINE
558
559
560	SUBROUTINE calc_month_len( iyear, imonth_len )
561	!!----------------------------------------------------------------------
562	!! * ROUTINE calc_month_len *
563	!!
564	!! ** Purpose : Compute the number of days in a months given a year.
565	!!
566	!! ** Method :
567	!!----------------------------------------------------------------------
568	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
569	INTEGER :: iyear !: year
570	!!----------------------------------------------------------------------
571	!
572	! length of the month of the current year (from nleapy, read in namelist)
573	IF ( nleapy < 2 ) THEN
574	imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
575	IF ( nleapy == 1 ) THEN ! we are using calendar with leap years
576	IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
577	imonth_len(2) = 29
578	ENDIF
579	ENDIF
580	ELSE
581	imonth_len(:) = nleapy ! all months with nleapy days per year
582	ENDIF
583	!
584	END SUBROUTINE
585
586
587	SUBROUTINE tra_asm_inc( kt )
588	!!----------------------------------------------------------------------
589	!! * ROUTINE tra_asm_inc *
590	!!
591	!! ** Purpose : Apply the tracer (T and S) assimilation increments
592	!!
593	!! ** Method : Direct initialization or Incremental Analysis Updating
594	!!
595	!! ** Action :
596	!!----------------------------------------------------------------------
597	INTEGER, INTENT(IN) :: kt ! Current time step
598	!
599	INTEGER :: ji,jj,jk
600	INTEGER :: it
601	REAL(wp) :: zincwgt ! IAU weight for current time step
602	!!----------------------------------------------------------------------
603
604	IF ( ln_asmiau ) THEN
605
606	!--------------------------------------------------------------------
607	! Incremental Analysis Updating
608	!--------------------------------------------------------------------
609
610	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
611
612	it = kt - nit000 + 1
613	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
614
615	IF(lwp) THEN
616	WRITE(numout,*)
617	WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', &
618	& kt,' with IAU weight = ', wgtiau(it)
619	WRITE(numout,*) '~~~~~~~~~~~~'
620	ENDIF
621
622	! Update the tracer tendencies
623	DO jk = 1, jpkm1
624	tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt
625	tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
626	END DO
627
628	! Salinity fix
629	IF (ln_salfix) THEN
630	DO jk = 1, jpkm1
631	DO jj = 1, jpj
632	DO ji= 1, jpi
633	tsa(ji,jj,jk,jp_sal) = MAX( tsa(ji,jj,jk,jp_sal), salfixmin )
634	END DO
635	END DO
636	END DO
637	ENDIF
638
639	ENDIF
640
641	IF ( kt == nitiaufin_r + 1 ) THEN ! For bias crcn to work
642	DEALLOCATE( t_bkginc )
643	DEALLOCATE( s_bkginc )
644	ENDIF
645
646
647	ELSEIF ( ln_asmdin ) THEN
648
649	!--------------------------------------------------------------------
650	! Direct Initialization
651	!--------------------------------------------------------------------
652
653	IF ( kt == nitdin_r ) THEN
654
655	neuler = 0 ! Force Euler forward step
656
657	! Initialize the now fields with the background + increment
658	tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)
659	tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)
660
661	! Optional salinity fix
662	IF (ln_salfix) THEN
663	DO jk = 1, jpkm1
664	DO jj = 1, jpj
665	DO ji= 1, jpi
666	tsn(ji,jj,jk,jp_sal) = MAX( tsn(ji,jj,jk,jp_sal), salfixmin )
667	END DO
668	END DO
669	END DO
670	ENDIF
671
672	tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields
673
674	CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities
675
676	#if ! defined key_c1d
677	IF( ln_zps ) &
678	& CALL zps_hde( nit000, jpts, tsb, & ! Partial steps: before horizontal derivative
679	& gtsu, gtsv, rhd, & ! of T, S, rd at the bottom ocean level
680	& gru , grv )
681	#endif
682
683	DEALLOCATE( t_bkginc )
684	DEALLOCATE( s_bkginc )
685	DEALLOCATE( t_bkg )
686	DEALLOCATE( s_bkg )
687	ENDIF
688	!
689	ENDIF
690	!
691	END SUBROUTINE tra_asm_inc
692
693
694	SUBROUTINE dyn_asm_inc( kt )
695	!!----------------------------------------------------------------------
696	!! * ROUTINE dyn_asm_inc *
697	!!
698	!! ** Purpose : Apply the dynamics (u and v) assimilation increments.
699	!!
700	!! ** Method : Direct initialization or Incremental Analysis Updating.
701	!!
702	!! ** Action :
703	!!----------------------------------------------------------------------
704	INTEGER, INTENT(IN) :: kt ! Current time step
705	!
706	INTEGER :: jk
707	INTEGER :: it
708	REAL(wp) :: zincwgt ! IAU weight for current time step
709	!!----------------------------------------------------------------------
710
711	IF ( ln_asmiau ) THEN
712
713	!--------------------------------------------------------------------
714	! Incremental Analysis Updating
715	!--------------------------------------------------------------------
716
717	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
718
719	it = kt - nit000 + 1
720	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
721
722	IF(lwp) THEN
723	WRITE(numout,*)
724	WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', &
725	& kt,' with IAU weight = ', wgtiau(it)
726	WRITE(numout,*) '~~~~~~~~~~~~'
727	ENDIF
728
729	! Update the dynamic tendencies
730	DO jk = 1, jpkm1
731	ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
732	va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
733	END DO
734
735	IF ( kt == nitiaufin_r ) THEN
736	DEALLOCATE( u_bkginc )
737	DEALLOCATE( v_bkginc )
738	ENDIF
739
740	ENDIF
741
742	ELSEIF ( ln_asmdin ) THEN
743
744	!--------------------------------------------------------------------
745	! Direct Initialization
746	!--------------------------------------------------------------------
747
748	IF ( kt == nitdin_r ) THEN
749
750	neuler = 0 ! Force Euler forward step
751
752	! Initialize the now fields with the background + increment
753	un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
754	vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)
755
756	ub(:,:,:) = un(:,:,:) ! Update before fields
757	vb(:,:,:) = vn(:,:,:)
758
759	CALL div_cur( kt ) ! Compute divergence and curl for now fields
760
761	rotb (:,:,:) = rotn (:,:,:) ! Update before fields
762	hdivb(:,:,:) = hdivn(:,:,:)
763
764	DEALLOCATE( u_bkg )
765	DEALLOCATE( v_bkg )
766	DEALLOCATE( u_bkginc )
767	DEALLOCATE( v_bkginc )
768	ENDIF
769	!
770	ENDIF
771	!
772	END SUBROUTINE dyn_asm_inc
773
774
775	SUBROUTINE ssh_asm_inc( kt )
776	!!----------------------------------------------------------------------
777	!! * ROUTINE ssh_asm_inc *
778	!!
779	!! ** Purpose : Apply the sea surface height assimilation increment.
780	!!
781	!! ** Method : Direct initialization or Incremental Analysis Updating.
782	!!
783	!! ** Action :
784	!!----------------------------------------------------------------------
785	INTEGER, INTENT(IN) :: kt ! Current time step
786	!
787	INTEGER :: it
788	REAL(wp) :: zincwgt ! IAU weight for current time step
789	!!----------------------------------------------------------------------
790
791	IF ( ln_asmiau ) THEN
792
793	!--------------------------------------------------------------------
794	! Incremental Analysis Updating
795	!--------------------------------------------------------------------
796
797	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
798
799	it = kt - nit000 + 1
800	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
801
802	IF(lwp) THEN
803	WRITE(numout,*)
804	WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
805	& kt,' with IAU weight = ', wgtiau(it)
806	WRITE(numout,*) '~~~~~~~~~~~~'
807	ENDIF
808
809	! Save the tendency associated with the IAU weighted SSH increment
810	! (applied in dynspg.*)
811	#if defined key_asminc
812	ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
813	#endif
814	IF ( kt == nitiaufin_r ) THEN
815	DEALLOCATE( ssh_bkginc )
816	ENDIF
817
818	ENDIF
819
820	ELSEIF ( ln_asmdin ) THEN
821
822	!--------------------------------------------------------------------
823	! Direct Initialization
824	!--------------------------------------------------------------------
825
826	IF ( kt == nitdin_r ) THEN
827
828	neuler = 0 ! Force Euler forward step
829
830	! Initialize the now fields the background + increment
831	sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)
832
833	sshb(:,:) = sshn(:,:) ! Update before fields
834
835	DEALLOCATE( ssh_bkg )
836	DEALLOCATE( ssh_bkginc )
837
838	ENDIF
839	!
840	ENDIF
841	!
842	END SUBROUTINE ssh_asm_inc
843
844	!!======================================================================
845	END MODULE asminc

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

Download in other formats: