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

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 4400

Last change on this file since 4400 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

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

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

Download in other formats: