Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

asminc.F90

Last change on this file was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

Download in other formats: