Context Navigation

asminc.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

Property svn:keywords set to Id

File size: 48.6 KB

Line
1	MODULE asminc
2	!!======================================================================
3	!! * MODULE asminc *
4	!! Assimilation increment : Apply an increment generated by data
5	!! assimilation
6	!!======================================================================
7	!! History : ! 2007-03 (M. Martin) Met Office version
8	!! ! 2007-04 (A. Weaver) calc_date original code
9	!! ! 2007-04 (A. Weaver) Merge with OPAVAR/NEMOVAR
10	!! NEMO 3.3 ! 2010-05 (D. Lea) Update to work with NEMO v3.2
11	!! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init
12	!! 3.4 ! 2012-10 (A. Weaver and K. Mogensen) Fix for direct initialization
13	!! ! 2014-09 (D. Lea) Local calc_date removed use routine from OBS
14	!! ! 2015-11 (D. Lea) Handle non-zero initial time of day
15	!!----------------------------------------------------------------------
16
17	!!----------------------------------------------------------------------
18	!! asm_inc_init : Initialize the increment arrays and IAU weights
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	!! ssh_asm_div : Apply divergence associated with SSH increment
23	!! seaice_asm_inc : Apply the seaice increment
24	!!----------------------------------------------------------------------
25	USE oce ! Dynamics and active tracers defined in memory
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 ldfdyn ! lateral diffusion: eddy viscosity coefficients
30	USE eosbn2 ! Equation of state - in situ and potential density
31	USE zpshde ! Partial step : Horizontal Derivative
32	USE asmpar ! Parameters for the assmilation interface
33	USE asmbkg !
34	USE c1d ! 1D initialization
35	USE sbc_oce ! Surface boundary condition variables.
36	USE diaobs , ONLY : calc_date ! Compute the calendar date on a given step
37	#if defined key_si3
38	USE ice , ONLY : hm_i, at_i, at_i_b
39	#endif
40	!
41	USE in_out_manager ! I/O manager
42	USE iom ! Library to read input files
43	USE lib_mpp ! MPP library
44
45	IMPLICIT NONE
46	PRIVATE
47
48	PUBLIC asm_inc_init !: Initialize the increment arrays and IAU weights
49	PUBLIC tra_asm_inc !: Apply the tracer (T and S) increments
50	PUBLIC dyn_asm_inc !: Apply the dynamic (u and v) increments
51	PUBLIC ssh_asm_inc !: Apply the SSH increment
52	PUBLIC ssh_asm_div !: Apply the SSH divergence
53	PUBLIC seaice_asm_inc !: Apply the seaice increment
54
55	#if defined key_asminc
56	LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE. !: Logical switch for assimilation increment interface
57	#else
58	LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE. !: No assimilation increments
59	#endif
60	LOGICAL, PUBLIC :: ln_bkgwri !: No output of the background state fields
61	LOGICAL, PUBLIC :: ln_asmiau !: No applying forcing with an assimilation increment
62	LOGICAL, PUBLIC :: ln_asmdin !: No direct initialization
63	LOGICAL, PUBLIC :: ln_trainc !: No tracer (T and S) assimilation increments
64	LOGICAL, PUBLIC :: ln_dyninc !: No dynamics (u and v) assimilation increments
65	LOGICAL, PUBLIC :: ln_sshinc !: No sea surface height assimilation increment
66	LOGICAL, PUBLIC :: ln_seaiceinc !: No sea ice concentration increment
67	LOGICAL, PUBLIC :: ln_salfix !: Apply minimum salinity check
68	LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
69	INTEGER, PUBLIC :: nn_divdmp !: Apply divergence damping filter nn_divdmp times
70
71	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkg , s_bkg !: Background temperature and salinity
72	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkg , v_bkg !: Background u- & v- velocity components
73	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkginc, s_bkginc !: Increment to the background T & S
74	REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components
75	REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step
76	#if defined key_asminc
77	REAL(wp), PUBLIC, DIMENSION(:,:) , ALLOCATABLE :: ssh_iau !: IAU-weighted sea surface height increment
78	#endif
79	! !!! time steps relative to the cycle interval [0,nitend-nit000-1]
80	INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term
81	INTEGER , PUBLIC :: nitdin !: Time step of the background state for direct initialization
82	INTEGER , PUBLIC :: nitiaustr !: Time step of the start of the IAU interval
83	INTEGER , PUBLIC :: nitiaufin !: Time step of the end of the IAU interval
84	!
85	INTEGER , PUBLIC :: niaufn !: Type of IAU weighing function: = 0 Constant weighting
86	! !: = 1 Linear hat-like, centred in middle of IAU interval
87	REAL(wp), PUBLIC :: salfixmin !: Ensure that the salinity is larger than this value if (ln_salfix)
88
89	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc ! Background sea surface height and its increment
90	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: seaice_bkginc ! Increment to the background sea ice conc
91	#if defined key_cice && defined key_asminc
92	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ndaice_da ! ice increment tendency into CICE
93	#endif
94
95	!! * Substitutions
96	# include "vectopt_loop_substitute.h90"
97	# include "do_loop_substitute.h90"
98	!!----------------------------------------------------------------------
99	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
100	!! $Id$
101	!! Software governed by the CeCILL license (see ./LICENSE)
102	!!----------------------------------------------------------------------
103	CONTAINS
104
105	SUBROUTINE asm_inc_init( Kbb, Kmm, Krhs )
106	!!----------------------------------------------------------------------
107	!! * ROUTINE asm_inc_init *
108	!!
109	!! ** Purpose : Initialize the assimilation increment and IAU weights.
110	!!
111	!! ** Method : Initialize the assimilation increment and IAU weights.
112	!!
113	!! ** Action :
114	!!----------------------------------------------------------------------
115	INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices
116	!
117	INTEGER :: ji, jj, jk, jt ! dummy loop indices
118	INTEGER :: imid, inum ! local integers
119	INTEGER :: ios ! Local integer output status for namelist read
120	INTEGER :: iiauper ! Number of time steps in the IAU period
121	INTEGER :: icycper ! Number of time steps in the cycle
122	REAL(KIND=dp) :: ditend_date ! Date YYYYMMDD.HHMMSS of final time step
123	REAL(KIND=dp) :: ditbkg_date ! Date YYYYMMDD.HHMMSS of background time step for Jb term
124	REAL(KIND=dp) :: ditdin_date ! Date YYYYMMDD.HHMMSS of background time step for DI
125	REAL(KIND=dp) :: ditiaustr_date ! Date YYYYMMDD.HHMMSS of IAU interval start time step
126	REAL(KIND=dp) :: ditiaufin_date ! Date YYYYMMDD.HHMMSS of IAU interval final time step
127
128	REAL(wp) :: znorm ! Normalization factor for IAU weights
129	REAL(wp) :: ztotwgt ! Value of time-integrated IAU weights (should be equal to one)
130	REAL(wp) :: z_inc_dateb ! Start date of interval on which increment is valid
131	REAL(wp) :: z_inc_datef ! End date of interval on which increment is valid
132	REAL(wp) :: zdate_bkg ! Date in background state file for DI
133	REAL(wp) :: zdate_inc ! Time axis in increments file
134	!
135	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zhdiv ! 2D workspace
136	!!
137	NAMELIST/nam_asminc/ ln_bkgwri, &
138	& ln_trainc, ln_dyninc, ln_sshinc, &
139	& ln_asmdin, ln_asmiau, &
140	& nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, &
141	& ln_salfix, salfixmin, nn_divdmp
142	!!----------------------------------------------------------------------
143
144	!-----------------------------------------------------------------------
145	! Read Namelist nam_asminc : assimilation increment interface
146	!-----------------------------------------------------------------------
147	ln_seaiceinc = .FALSE.
148	ln_temnofreeze = .FALSE.
149
150	READ ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
151	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in reference namelist' )
152	READ ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
153	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_asminc in configuration namelist' )
154	IF(lwm) WRITE ( numond, nam_asminc )
155
156	! Control print
157	IF(lwp) THEN
158	WRITE(numout,*)
159	WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
160	WRITE(numout,*) '~~~~~~~~~~~~'
161	WRITE(numout,*) ' Namelist namasm : set assimilation increment parameters'
162	WRITE(numout,*) ' Logical switch for writing out background state ln_bkgwri = ', ln_bkgwri
163	WRITE(numout,*) ' Logical switch for applying tracer increments ln_trainc = ', ln_trainc
164	WRITE(numout,*) ' Logical switch for applying velocity increments ln_dyninc = ', ln_dyninc
165	WRITE(numout,*) ' Logical switch for applying SSH increments ln_sshinc = ', ln_sshinc
166	WRITE(numout,*) ' Logical switch for Direct Initialization (DI) ln_asmdin = ', ln_asmdin
167	WRITE(numout,*) ' Logical switch for applying sea ice increments ln_seaiceinc = ', ln_seaiceinc
168	WRITE(numout,*) ' Logical switch for Incremental Analysis Updating (IAU) ln_asmiau = ', ln_asmiau
169	WRITE(numout,*) ' Timestep of background in [0,nitend-nit000-1] nitbkg = ', nitbkg
170	WRITE(numout,*) ' Timestep of background for DI in [0,nitend-nit000-1] nitdin = ', nitdin
171	WRITE(numout,*) ' Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
172	WRITE(numout,*) ' Timestep of end of IAU interval in [0,nitend-nit000-1] nitiaufin = ', nitiaufin
173	WRITE(numout,*) ' Type of IAU weighting function niaufn = ', niaufn
174	WRITE(numout,*) ' Logical switch for ensuring that the sa > salfixmin ln_salfix = ', ln_salfix
175	WRITE(numout,*) ' Minimum salinity after applying the increments salfixmin = ', salfixmin
176	ENDIF
177
178	nitbkg_r = nitbkg + nit000 - 1 ! Background time referenced to nit000
179	nitdin_r = nitdin + nit000 - 1 ! Background time for DI referenced to nit000
180	nitiaustr_r = nitiaustr + nit000 - 1 ! Start of IAU interval referenced to nit000
181	nitiaufin_r = nitiaufin + nit000 - 1 ! End of IAU interval referenced to nit000
182
183	iiauper = nitiaufin_r - nitiaustr_r + 1 ! IAU interval length
184	icycper = nitend - nit000 + 1 ! Cycle interval length
185
186	CALL calc_date( nitend , ditend_date ) ! Date of final time step
187	CALL calc_date( nitbkg_r , ditbkg_date ) ! Background time for Jb referenced to ndate0
188	CALL calc_date( nitdin_r , ditdin_date ) ! Background time for DI referenced to ndate0
189	CALL calc_date( nitiaustr_r, ditiaustr_date ) ! IAU start time referenced to ndate0
190	CALL calc_date( nitiaufin_r, ditiaufin_date ) ! IAU end time referenced to ndate0
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,*)
203	WRITE(numout,*) ' Dates referenced to current cycle:'
204	WRITE(numout,*) ' ndastp = ', ndastp
205	WRITE(numout,*) ' ndate0 = ', ndate0
206	WRITE(numout,*) ' nn_time0 = ', nn_time0
207	WRITE(numout,*) ' ditend_date = ', ditend_date
208	WRITE(numout,*) ' ditbkg_date = ', ditbkg_date
209	WRITE(numout,*) ' ditdin_date = ', ditdin_date
210	WRITE(numout,*) ' ditiaustr_date = ', ditiaustr_date
211	WRITE(numout,*) ' ditiaufin_date = ', ditiaufin_date
212	ENDIF
213
214
215	IF ( ( ln_asmdin ).AND.( ln_asmiau ) ) &
216	& CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
217	& ' Choose Direct Initialization OR Incremental Analysis Updating')
218
219	IF ( ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
220	.AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) &
221	& CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', &
222	& ' but ln_asmdin and ln_asmiau are both set to .false. :', &
223	& ' Inconsistent options')
224
225	IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
226	& CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :', &
227	& ' Type IAU weighting function is invalid')
228
229	IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
230	& ) &
231	& CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', &
232	& ' The assimilation increments are not applied')
233
234	IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
235	& CALL ctl_stop( ' nitiaustr = nitiaufin :', &
236	& ' IAU interval is of zero length')
237
238	IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
239	& CALL ctl_stop( ' nitiaustr or nitiaufin :', &
240	& ' IAU starting or final time step is outside the cycle interval', &
241	& ' Valid range nit000 to nitend')
242
243	IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
244	& CALL ctl_stop( ' nitbkg :', &
245	& ' Background time step is outside the cycle interval')
246
247	IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
248	& CALL ctl_stop( ' nitdin :', &
249	& ' Background time step for Direct Initialization is outside', &
250	& ' the cycle interval')
251
252	IF ( nstop > 0 ) RETURN ! if there are any errors then go no further
253
254	!--------------------------------------------------------------------
255	! Initialize the Incremental Analysis Updating weighting function
256	!--------------------------------------------------------------------
257
258	IF( ln_asmiau ) THEN
259	!
260	ALLOCATE( wgtiau( icycper ) )
261	!
262	wgtiau(:) = 0._wp
263	!
264	! !---------------------------------------------------------
265	IF( niaufn == 0 ) THEN ! Constant IAU forcing
266	! !---------------------------------------------------------
267	DO jt = 1, iiauper
268	wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
269	END DO
270	! !---------------------------------------------------------
271	ELSEIF ( niaufn == 1 ) THEN ! Linear hat-like, centred in middle of IAU interval
272	! !---------------------------------------------------------
273	! Compute the normalization factor
274	znorm = 0._wp
275	IF( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval
276	imid = iiauper / 2
277	DO jt = 1, imid
278	znorm = znorm + REAL( jt )
279	END DO
280	znorm = 2.0 * znorm
281	ELSE ! Odd number of time steps in IAU interval
282	imid = ( iiauper + 1 ) / 2
283	DO jt = 1, imid - 1
284	znorm = znorm + REAL( jt )
285	END DO
286	znorm = 2.0 * znorm + REAL( imid )
287	ENDIF
288	znorm = 1.0 / znorm
289	!
290	DO jt = 1, imid - 1
291	wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
292	END DO
293	DO jt = imid, iiauper
294	wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
295	END DO
296	!
297	ENDIF
298
299	! Test that the integral of the weights over the weighting interval equals 1
300	IF(lwp) THEN
301	WRITE(numout,*)
302	WRITE(numout,*) 'asm_inc_init : IAU weights'
303	WRITE(numout,*) '~~~~~~~~~~~~'
304	WRITE(numout,*) ' time step IAU weight'
305	WRITE(numout,*) ' ========= ====================='
306	ztotwgt = 0.0
307	DO jt = 1, icycper
308	ztotwgt = ztotwgt + wgtiau(jt)
309	WRITE(numout,*) ' ', jt, ' ', wgtiau(jt)
310	END DO
311	WRITE(numout,*) ' ==================================='
312	WRITE(numout,*) ' Time-integrated weight = ', ztotwgt
313	WRITE(numout,*) ' ==================================='
314	ENDIF
315
316	ENDIF
317
318	!--------------------------------------------------------------------
319	! Allocate and initialize the increment arrays
320	!--------------------------------------------------------------------
321
322	ALLOCATE( t_bkginc (jpi,jpj,jpk) ) ; t_bkginc (:,:,:) = 0._wp
323	ALLOCATE( s_bkginc (jpi,jpj,jpk) ) ; s_bkginc (:,:,:) = 0._wp
324	ALLOCATE( u_bkginc (jpi,jpj,jpk) ) ; u_bkginc (:,:,:) = 0._wp
325	ALLOCATE( v_bkginc (jpi,jpj,jpk) ) ; v_bkginc (:,:,:) = 0._wp
326	ALLOCATE( ssh_bkginc (jpi,jpj) ) ; ssh_bkginc (:,:) = 0._wp
327	ALLOCATE( seaice_bkginc(jpi,jpj) ) ; seaice_bkginc(:,:) = 0._wp
328	#if defined key_asminc
329	ALLOCATE( ssh_iau (jpi,jpj) ) ; ssh_iau (:,:) = 0._wp
330	#endif
331	#if defined key_cice && defined key_asminc
332	ALLOCATE( ndaice_da (jpi,jpj) ) ; ndaice_da (:,:) = 0._wp
333	#endif
334	!
335	IF ( ln_trainc .OR. ln_dyninc .OR. & !--------------------------------------
336	& ln_sshinc .OR. ln_seaiceinc ) THEN ! Read the increments from file
337	! !--------------------------------------
338	CALL iom_open( c_asminc, inum )
339	!
340	CALL iom_get( inum, 'time' , zdate_inc )
341	CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
342	CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
343	z_inc_dateb = zdate_inc
344	z_inc_datef = zdate_inc
345	!
346	IF(lwp) THEN
347	WRITE(numout,*)
348	WRITE(numout,*) 'asm_inc_init : Assimilation increments valid between dates ', z_inc_dateb,' and ', z_inc_datef
349	WRITE(numout,*) '~~~~~~~~~~~~'
350	ENDIF
351	!
352	IF ( ( z_inc_dateb < ndastp + nn_time0*0.0001_wp ) .OR. &
353	& ( z_inc_datef > ditend_date ) ) &
354	& CALL ctl_warn( ' Validity time of assimilation increments is ', &
355	& ' outside the assimilation interval' )
356
357	IF ( ( ln_asmdin ).AND.( zdate_inc /= ditdin_date ) ) &
358	& CALL ctl_warn( ' Validity time of assimilation increments does ', &
359	& ' not agree with Direct Initialization time' )
360
361	IF ( ln_trainc ) THEN
362	CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
363	CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
364	! Apply the masks
365	t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
366	s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
367	! Set missing increments to 0.0 rather than 1e+20
368	! to allow for differences in masks
369	WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
370	WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
371	ENDIF
372
373	IF ( ln_dyninc ) THEN
374	CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )
375	CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )
376	! Apply the masks
377	u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
378	v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
379	! Set missing increments to 0.0 rather than 1e+20
380	! to allow for differences in masks
381	WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
382	WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
383	ENDIF
384
385	IF ( ln_sshinc ) THEN
386	CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
387	! Apply the masks
388	ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
389	! Set missing increments to 0.0 rather than 1e+20
390	! to allow for differences in masks
391	WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
392	ENDIF
393
394	IF ( ln_seaiceinc ) THEN
395	CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 )
396	! Apply the masks
397	seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1)
398	! Set missing increments to 0.0 rather than 1e+20
399	! to allow for differences in masks
400	WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
401	ENDIF
402	!
403	CALL iom_close( inum )
404	!
405	ENDIF
406	!
407	! !--------------------------------------
408	IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN ! Apply divergence damping filter
409	! !--------------------------------------
410	ALLOCATE( zhdiv(jpi,jpj) )
411	!
412	DO jt = 1, nn_divdmp
413	!
414	DO jk = 1, jpkm1 ! zhdiv = e1e1 * div
415	zhdiv(:,:) = 0._wp
416	DO_2D_00_00
417	zhdiv(ji,jj) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * u_bkginc(ji ,jj,jk) &
418	& - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk) &
419	& + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * v_bkginc(ji,jj ,jk) &
420	& - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk) ) / e3t(ji,jj,jk,Kmm)
421	END_2D
422	CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. ) ! lateral boundary cond. (no sign change)
423	!
424	DO_2D_00_00
425	u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) &
426	& + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
427	v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) &
428	& + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
429	END_2D
430	END DO
431	!
432	END DO
433	!
434	DEALLOCATE( zhdiv )
435	!
436	ENDIF
437	!
438	! !-----------------------------------------------------
439	IF ( ln_asmdin ) THEN ! Allocate and initialize the background state arrays
440	! !-----------------------------------------------------
441	!
442	ALLOCATE( t_bkg (jpi,jpj,jpk) ) ; t_bkg (:,:,:) = 0._wp
443	ALLOCATE( s_bkg (jpi,jpj,jpk) ) ; s_bkg (:,:,:) = 0._wp
444	ALLOCATE( u_bkg (jpi,jpj,jpk) ) ; u_bkg (:,:,:) = 0._wp
445	ALLOCATE( v_bkg (jpi,jpj,jpk) ) ; v_bkg (:,:,:) = 0._wp
446	ALLOCATE( ssh_bkg(jpi,jpj) ) ; ssh_bkg(:,:) = 0._wp
447	!
448	!
449	!--------------------------------------------------------------------
450	! Read from file the background state at analysis time
451	!--------------------------------------------------------------------
452	!
453	CALL iom_open( c_asmdin, inum )
454	!
455	CALL iom_get( inum, 'rdastp', zdate_bkg )
456	!
457	IF(lwp) THEN
458	WRITE(numout,*)
459	WRITE(numout,*) ' ==>>> Assimilation background state valid at : ', zdate_bkg
460	WRITE(numout,*)
461	ENDIF
462	!
463	IF ( zdate_bkg /= ditdin_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	ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
484	ENDIF
485	!
486	CALL iom_close( inum )
487	!
488	ENDIF
489	!
490	IF(lwp) WRITE(numout,*) ' ==>>> Euler time step switch is ', neuler
491	!
492	IF( lk_asminc ) THEN !== data assimilation ==!
493	IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1, Kmm ) ! Output background fields
494	IF( ln_asmdin ) THEN ! Direct initialization
495	IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1, Kbb, Kmm, ts , Krhs ) ! Tracers
496	IF( ln_dyninc ) CALL dyn_asm_inc( nit000 - 1, Kbb, Kmm, uu, vv, Krhs ) ! Dynamics
497	IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1, Kbb, Kmm ) ! SSH
498	ENDIF
499	ENDIF
500	!
501	END SUBROUTINE asm_inc_init
502
503
504	SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs )
505	!!----------------------------------------------------------------------
506	!! * ROUTINE tra_asm_inc *
507	!!
508	!! ** Purpose : Apply the tracer (T and S) assimilation increments
509	!!
510	!! ** Method : Direct initialization or Incremental Analysis Updating
511	!!
512	!! ** Action :
513	!!----------------------------------------------------------------------
514	INTEGER , INTENT(in ) :: kt ! Current time step
515	INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! Time level indices
516	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
517	!
518	INTEGER :: ji, jj, jk
519	INTEGER :: it
520	REAL(wp) :: zincwgt ! IAU weight for current time step
521	REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
522	!!----------------------------------------------------------------------
523	!
524	! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
525	! used to prevent the applied increments taking the temperature below the local freezing point
526	DO jk = 1, jpkm1
527	CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) )
528	END DO
529	!
530	! !--------------------------------------
531	IF ( ln_asmiau ) THEN ! Incremental Analysis Updating
532	! !--------------------------------------
533	!
534	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
535	!
536	it = kt - nit000 + 1
537	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
538	!
539	IF(lwp) THEN
540	WRITE(numout,*)
541	WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
542	WRITE(numout,*) '~~~~~~~~~~~~'
543	ENDIF
544	!
545	! Update the tracer tendencies
546	DO jk = 1, jpkm1
547	IF (ln_temnofreeze) THEN
548	! Do not apply negative increments if the temperature will fall below freezing
549	WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
550	& pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )
551	pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt
552	END WHERE
553	ELSE
554	pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt
555	ENDIF
556	IF (ln_salfix) THEN
557	! Do not apply negative increments if the salinity will fall below a specified
558	! minimum value salfixmin
559	WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
560	& pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )
561	pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt
562	END WHERE
563	ELSE
564	pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt
565	ENDIF
566	END DO
567	!
568	ENDIF
569	!
570	IF ( kt == nitiaufin_r + 1 ) THEN ! For bias crcn to work
571	DEALLOCATE( t_bkginc )
572	DEALLOCATE( s_bkginc )
573	ENDIF
574	! !--------------------------------------
575	ELSEIF ( ln_asmdin ) THEN ! Direct Initialization
576	! !--------------------------------------
577	!
578	IF ( kt == nitdin_r ) THEN
579	!
580	neuler = 0 ! Force Euler forward step
581	!
582	! Initialize the now fields with the background + increment
583	IF (ln_temnofreeze) THEN
584	! Do not apply negative increments if the temperature will fall below freezing
585	WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) )
586	pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)
587	END WHERE
588	ELSE
589	pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)
590	ENDIF
591	IF (ln_salfix) THEN
592	! Do not apply negative increments if the salinity will fall below a specified
593	! minimum value salfixmin
594	WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin )
595	pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)
596	END WHERE
597	ELSE
598	pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)
599	ENDIF
600
601	pts(:,:,:,:,Kbb) = pts(:,:,:,:,Kmm) ! Update before fields
602
603	CALL eos( pts(:,:,:,:,Kbb), rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities
604	!!gm fabien
605	! CALL eos( pts(:,:,:,:,Kbb), rhd, rhop ) ! Before potential and in situ densities
606	!!gm
607
608	IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) &
609	& CALL zps_hde ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient
610	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
611	IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) &
612	& CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
613	& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the last ocean level
614
615	DEALLOCATE( t_bkginc )
616	DEALLOCATE( s_bkginc )
617	DEALLOCATE( t_bkg )
618	DEALLOCATE( s_bkg )
619	ENDIF
620	!
621	ENDIF
622	! Perhaps the following call should be in step
623	IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment
624	!
625	END SUBROUTINE tra_asm_inc
626
627
628	SUBROUTINE dyn_asm_inc( kt, Kbb, Kmm, puu, pvv, Krhs )
629	!!----------------------------------------------------------------------
630	!! * ROUTINE dyn_asm_inc *
631	!!
632	!! ** Purpose : Apply the dynamics (u and v) assimilation increments.
633	!!
634	!! ** Method : Direct initialization or Incremental Analysis Updating.
635	!!
636	!! ** Action :
637	!!----------------------------------------------------------------------
638	INTEGER , INTENT( in ) :: kt ! ocean time-step index
639	INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices
640	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
641	!
642	INTEGER :: jk
643	INTEGER :: it
644	REAL(wp) :: zincwgt ! IAU weight for current time step
645	!!----------------------------------------------------------------------
646	!
647	! !--------------------------------------------
648	IF ( ln_asmiau ) THEN ! Incremental Analysis Updating
649	! !--------------------------------------------
650	!
651	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
652	!
653	it = kt - nit000 + 1
654	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
655	!
656	IF(lwp) THEN
657	WRITE(numout,*)
658	WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
659	WRITE(numout,*) '~~~~~~~~~~~~'
660	ENDIF
661	!
662	! Update the dynamic tendencies
663	DO jk = 1, jpkm1
664	puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + u_bkginc(:,:,jk) * zincwgt
665	pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + v_bkginc(:,:,jk) * zincwgt
666	END DO
667	!
668	IF ( kt == nitiaufin_r ) THEN
669	DEALLOCATE( u_bkginc )
670	DEALLOCATE( v_bkginc )
671	ENDIF
672	!
673	ENDIF
674	! !-----------------------------------------
675	ELSEIF ( ln_asmdin ) THEN ! Direct Initialization
676	! !-----------------------------------------
677	!
678	IF ( kt == nitdin_r ) THEN
679	!
680	neuler = 0 ! Force Euler forward step
681	!
682	! Initialize the now fields with the background + increment
683	puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:)
684	pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:)
685	!
686	puu(:,:,:,Kbb) = puu(:,:,:,Kmm) ! Update before fields
687	pvv(:,:,:,Kbb) = pvv(:,:,:,Kmm)
688	!
689	DEALLOCATE( u_bkg )
690	DEALLOCATE( v_bkg )
691	DEALLOCATE( u_bkginc )
692	DEALLOCATE( v_bkginc )
693	ENDIF
694	!
695	ENDIF
696	!
697	END SUBROUTINE dyn_asm_inc
698
699
700	SUBROUTINE ssh_asm_inc( kt, Kbb, Kmm )
701	!!----------------------------------------------------------------------
702	!! * ROUTINE ssh_asm_inc *
703	!!
704	!! ** Purpose : Apply the sea surface height assimilation increment.
705	!!
706	!! ** Method : Direct initialization or Incremental Analysis Updating.
707	!!
708	!! ** Action :
709	!!----------------------------------------------------------------------
710	INTEGER, INTENT(IN) :: kt ! Current time step
711	INTEGER, INTENT(IN) :: Kbb, Kmm ! Current time step
712	!
713	INTEGER :: it
714	INTEGER :: jk
715	REAL(wp) :: zincwgt ! IAU weight for current time step
716	!!----------------------------------------------------------------------
717	!
718	! !-----------------------------------------
719	IF ( ln_asmiau ) THEN ! Incremental Analysis Updating
720	! !-----------------------------------------
721	!
722	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
723	!
724	it = kt - nit000 + 1
725	zincwgt = wgtiau(it) / rdt ! IAU weight for the current time step
726	!
727	IF(lwp) THEN
728	WRITE(numout,*)
729	WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
730	& kt,' with IAU weight = ', wgtiau(it)
731	WRITE(numout,*) '~~~~~~~~~~~~'
732	ENDIF
733	!
734	! Save the tendency associated with the IAU weighted SSH increment
735	! (applied in dynspg.*)
736	#if defined key_asminc
737	ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
738	#endif
739	!
740	ELSE IF( kt == nitiaufin_r+1 ) THEN
741	!
742	! test on ssh_bkginc needed as ssh_asm_inc is called twice by time step
743	IF (ALLOCATED(ssh_bkginc)) DEALLOCATE( ssh_bkginc )
744	!
745	#if defined key_asminc
746	ssh_iau(:,:) = 0._wp
747	#endif
748	!
749	ENDIF
750	! !-----------------------------------------
751	ELSEIF ( ln_asmdin ) THEN ! Direct Initialization
752	! !-----------------------------------------
753	!
754	IF ( kt == nitdin_r ) THEN
755	!
756	neuler = 0 ! Force Euler forward step
757	!
758	ssh(:,:,Kmm) = ssh_bkg(:,:) + ssh_bkginc(:,:) ! Initialize the now fields the background + increment
759	!
760	ssh(:,:,Kbb) = ssh(:,:,Kmm) ! Update before fields
761	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
762	!!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ????
763	!
764	DEALLOCATE( ssh_bkg )
765	DEALLOCATE( ssh_bkginc )
766	!
767	ENDIF
768	!
769	ENDIF
770	!
771	END SUBROUTINE ssh_asm_inc
772
773
774	SUBROUTINE ssh_asm_div( kt, Kbb, Kmm, phdivn )
775	!!----------------------------------------------------------------------
776	!! * ROUTINE ssh_asm_div *
777	!!
778	!! ** Purpose : ssh increment with z* is incorporated via a correction of the local divergence
779	!! across all the water column
780	!!
781	!! ** Method :
782	!! CAUTION : sshiau is positive (inflow) decreasing the
783	!! divergence and expressed in m/s
784	!!
785	!! ** Action : phdivn decreased by the ssh increment
786	!!----------------------------------------------------------------------
787	INTEGER, INTENT(IN) :: kt ! ocean time-step index
788	INTEGER, INTENT(IN) :: Kbb, Kmm ! time level indices
789	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence
790	!!
791	INTEGER :: jk ! dummy loop index
792	REAL(wp), DIMENSION(:,:) , POINTER :: ztim ! local array
793	!!----------------------------------------------------------------------
794	!
795	#if defined key_asminc
796	CALL ssh_asm_inc( kt, Kbb, Kmm ) !== (calculate increments)
797	!
798	IF( ln_linssh ) THEN
799	phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1)
800	ELSE
801	ALLOCATE( ztim(jpi,jpj) )
802	ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) )
803	DO jk = 1, jpkm1
804	phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk)
805	END DO
806	!
807	DEALLOCATE(ztim)
808	ENDIF
809	#endif
810	!
811	END SUBROUTINE ssh_asm_div
812
813
814	SUBROUTINE seaice_asm_inc( kt, kindic )
815	!!----------------------------------------------------------------------
816	!! * ROUTINE seaice_asm_inc *
817	!!
818	!! ** Purpose : Apply the sea ice assimilation increment.
819	!!
820	!! ** Method : Direct initialization or Incremental Analysis Updating.
821	!!
822	!! ** Action :
823	!!
824	!!----------------------------------------------------------------------
825	INTEGER, INTENT(in) :: kt ! Current time step
826	INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation
827	!
828	INTEGER :: it
829	REAL(wp) :: zincwgt ! IAU weight for current time step
830	#if defined key_si3
831	REAL(wp), DIMENSION(jpi,jpj) :: zofrld, zohicif, zseaicendg, zhicifinc
832	REAL(wp) :: zhicifmin = 0.5_wp ! ice minimum depth in metres
833	#endif
834	!!----------------------------------------------------------------------
835	!
836	! !-----------------------------------------
837	IF ( ln_asmiau ) THEN ! Incremental Analysis Updating
838	! !-----------------------------------------
839	!
840	IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
841	!
842	it = kt - nit000 + 1
843	zincwgt = wgtiau(it) ! IAU weight for the current time step
844	! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
845	!
846	IF(lwp) THEN
847	WRITE(numout,*)
848	WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
849	WRITE(numout,*) '~~~~~~~~~~~~'
850	ENDIF
851	!
852	! Sea-ice : SI3 case
853	!
854	#if defined key_si3
855	zofrld (:,:) = 1._wp - at_i(:,:)
856	zohicif(:,:) = hm_i(:,:)
857	!
858	at_i (:,:) = 1. - MIN( MAX( 1.-at_i (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
859	at_i_b(:,:) = 1. - MIN( MAX( 1.-at_i_b(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
860	fr_i(:,:) = at_i(:,:) ! adjust ice fraction
861	!
862	zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:)) ! find out actual sea ice nudge applied
863	!
864	! Nudge sea ice depth to bring it up to a required minimum depth
865	WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin )
866	zhicifinc(:,:) = (zhicifmin - hm_i(:,:)) * zincwgt
867	ELSEWHERE
868	zhicifinc(:,:) = 0.0_wp
869	END WHERE
870	!
871	! nudge ice depth
872	hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:)
873	!
874	! seaice salinity balancing (to add)
875	#endif
876	!
877	#if defined key_cice && defined key_asminc
878	! Sea-ice : CICE case. Pass ice increment tendency into CICE
879	ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
880	#endif
881	!
882	IF ( kt == nitiaufin_r ) THEN
883	DEALLOCATE( seaice_bkginc )
884	ENDIF
885	!
886	ELSE
887	!
888	#if defined key_cice && defined key_asminc
889	ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE
890	#endif
891	!
892	ENDIF
893	! !-----------------------------------------
894	ELSEIF ( ln_asmdin ) THEN ! Direct Initialization
895	! !-----------------------------------------
896	!
897	IF ( kt == nitdin_r ) THEN
898	!
899	neuler = 0 ! Force Euler forward step
900	!
901	! Sea-ice : SI3 case
902	!
903	#if defined key_si3
904	zofrld (:,:) = 1._wp - at_i(:,:)
905	zohicif(:,:) = hm_i(:,:)
906	!
907	! Initialize the now fields the background + increment
908	at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
909	at_i_b(:,:) = at_i(:,:)
910	fr_i(:,:) = at_i(:,:) ! adjust ice fraction
911	!
912	zseaicendg(:,:) = zofrld(:,:) - (1. - at_i(:,:)) ! find out actual sea ice nudge applied
913	!
914	! Nudge sea ice depth to bring it up to a required minimum depth
915	WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(:,:) < zhicifmin )
916	zhicifinc(:,:) = zhicifmin - hm_i(:,:)
917	ELSEWHERE
918	zhicifinc(:,:) = 0.0_wp
919	END WHERE
920	!
921	! nudge ice depth
922	hm_i (:,:) = hm_i (:,:) + zhicifinc(:,:)
923	!
924	! seaice salinity balancing (to add)
925	#endif
926	!
927	#if defined key_cice && defined key_asminc
928	! Sea-ice : CICE case. Pass ice increment tendency into CICE
929	ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
930	#endif
931	IF ( .NOT. PRESENT(kindic) ) THEN
932	DEALLOCATE( seaice_bkginc )
933	END IF
934	!
935	ELSE
936	!
937	#if defined key_cice && defined key_asminc
938	ndaice_da(:,:) = 0._wp ! Sea-ice : CICE case. Zero ice increment tendency into CICE
939	#endif
940	!
941	ENDIF
942
943	!#if defined defined key_si3 \|\| defined key_cice
944	!
945	! IF (ln_seaicebal ) THEN
946	! !! balancing salinity increments
947	! !! simple case from limflx.F90 (doesn't include a mass flux)
948	! !! assumption is that as ice concentration is reduced or increased
949	! !! the snow and ice depths remain constant
950	! !! note that snow is being created where ice concentration is being increased
951	! !! - could be more sophisticated and
952	! !! not do this (but would need to alter h_snow)
953	!
954	! usave(:,:,:)=sb(:,:,:) ! use array as a temporary store
955	!
956	! DO jj = 1, jpj
957	! DO ji = 1, jpi
958	! ! calculate change in ice and snow mass per unit area
959	! ! positive values imply adding salt to the ocean (results from ice formation)
960	! ! fwf : ice formation and melting
961	!
962	! zfons = ( -nfresh_da(ji,jj)soce + nfsalt_da(ji,jj) )rdt
963	!
964	! ! change salinity down to mixed layer depth
965	! mld=hmld_kara(ji,jj)
966	!
967	! ! prevent small mld
968	! ! less than 10m can cause salinity instability
969	! IF (mld < 10) mld=10
970	!
971	! ! set to bottom of a level
972	! DO jk = jpk-1, 2, -1
973	! IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN
974	! mld=gdepw(ji,jj,jk+1)
975	! jkmax=jk
976	! ENDIF
977	! ENDDO
978	!
979	! ! avoid applying salinity balancing in shallow water or on land
980	! !
981	!
982	! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
983	!
984	! dsal_ocn=0.0_wp
985	! sal_thresh=5.0_wp ! minimum salinity threshold for salinity balancing
986	!
987	! if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
988	! dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
989	!
990	! ! put increments in for levels in the mixed layer
991	! ! but prevent salinity below a threshold value
992	!
993	! DO jk = 1, jkmax
994	!
995	! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
996	! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
997	! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
998	! ENDIF
999	!
1000	! ENDDO
1001	!
1002	! ! ! salt exchanges at the ice/ocean interface
1003	! ! zpmess = zfons / rdt_ice ! rdt_ice is ice timestep
1004	! !
1005	! !! Adjust fsalt. A +ve fsalt means adding salt to ocean
1006	! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt
1007	! !!
1008	! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d)
1009	! !! ! E-P (kg m-2 s-2)
1010	! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2)
1011	! ENDDO !ji
1012	! ENDDO !jj!
1013	!
1014	! ENDIF !ln_seaicebal
1015	!
1016	!#endif
1017	!
1018	ENDIF
1019	!
1020	END SUBROUTINE seaice_asm_inc
1021
1022	!!======================================================================
1023	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: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ASM/asminc.F90 @ 12340

Download in other formats: