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

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

bias.F90 in branches/UKMO/dev_r5518_pcbias_ipc/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_pcbias_ipc/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90 @ 7082

Last change on this file since 7082 was 7082, checked in by isabella, 8 years ago
changed back to match sum
File size: 44.6 KB

Line
1	MODULE bias
2	!! Is used by OPA and STEP
3	!!======================================================================
4	!! * Module bias *
5	!! Code to estimate and apply bias correction.
6	!! The bias is in T/S and Pressure. It follows the generalized
7	!! bias algorithm presented in Balmaseda et al 2007.
8	!!
9	!! It can be read from a file offline, estimated online from relaxation
10	!! terms or from assimilation increments (this latter estimtd in inner loop)
11	!!
12	!! The parameters for the time evolution of the bias can be different
13	!! from the relaxation terms and for the assim increments. Only the
14	!! parameter for the relaxtion terms are considered here.
15	!!
16	!! The offline bias term can contain the seasonal cycle.
17	!!
18	!! The time evolution of the bias for relaxtion is estimated as followed
19	!! bias_rlx(t+1)=t_rlx_membias_rlx(t)+t_rlx_upd(t/s)trdmp.
20	!!
21	!! The total bias in T/S is partion between the correction to T/S only
22	!! (fb_t) and the correction applied to the pressure gradient (fb_p).
23	!! We impose that (fb_p=1.-fb_t). These factors can be different for the
24	!! different terms(fb_t_rxl,fb_t_asm,fb_t_ofl)
25	!!
26	!! (t/s)bias = fb_t_ofl * (t/s)bias_ofl +
27	!! fb_t_rlx * (t/s)bias_rlx +
28	!! fb_t_asm * (t/s)bias_asm
29	!! (t/s)bias_p =fb_p_ofl * (t/s)bias_ofl+
30	!! fb_p_rlx * (t/s)bias_rlx_p +
31	!! fb_p_asm * (t/s)bias_asm_p
32	!! (t/s)bias is applied directely to correct T and S
33	!! (t/s)bias_p is used to compute rhd_pc and gru/v_pc
34	!!
35	!! Note: the above is an adhoc /simple way of making the partition
36	!! between bias in pressure and bias in T/S. It would be better
37	!! if the partition is done at the time of estimating the time
38	!! evolution of the bias. That would mean doubling the number of
39	!! 3D arrays.
40	!!
41	!! New addtion: (not in Balmaseda et al 2007):
42	!! A more physical alternative to the partition of the bias can be
43	!! done in terms of the inertial frequency: when the time scales of
44	!! adjustment are shorter that >1/f (Equator), the bias correction should
45	!! be in the the pressure term. Otherwise, it can act directly in T/S.
46	!! NOTE that the bias correction in the pressure term here (following
47	!! (Bell et al 2007) is just the "opposite" as the semi-prognostic method
48	!! in Greatbatch et al 2004.
49	!! The use of this partition is controlled by ln_inertial=.true.
50	!!
51	!!
52	!! 2009-03 (M.A. Balmaseda ECMWF)
53	!!======================================================================
54
55	!!----------------------------------------------------------------------
56	!! bias_init : Read in the bias namelist and the bias arrays
57	!! tra_bias : Apply the bias fields on T/S directly
58	!! dyn_bias : Compute density correction for dynamic hpg
59	!! bias_opn : open bias files for restart capabilities
60	!! bias_wrt : write bias fies " " "
61	!!----------------------------------------------------------------------
62	!! * Modules used
63	USE par_kind, ONLY: &
64	& wp
65	USE par_oce, ONLY: &
66	& jpi, &
67	& jpj, &
68	& jpk
69	USE dom_oce, ONLY: &
70	& rdt, &
71	& ln_zps, &
72	& gphit
73	USE phycst, ONLY: &
74	& rday, &
75	& rad
76	USE oce, ONLY: &
77	& tsb, tsn, tsa, &
78	& rhop, &
79	& gtsu, gtsv
80	USE dynhpg, ONLY: &
81	& ln_dynhpg_imp
82	USE tradmp
83	USE dtatsd, ONLY: &
84	& ln_tsd_tradmp
85	USE in_out_manager, ONLY: &
86	& lwp, &
87	& numnam_ref, &
88	& numnam_cfg, &
89	& numond, &
90	& numout, &
91	& lrst_oce, &
92	& nit000
93	USE iom
94	USE eosbn2
95	USE zpshde ! partial step: hor. derivative (zps_hde routine)
96	USE biaspar
97	USE fldread ! read input fields
98	USE lbclnk ! lateral boundary conditions (or mpp link)
99	USE asmpar
100	USE asminc
101	USE lib_mpp, ONLY: &
102	& ctl_stop, &
103	& ctl_nam
104
105	IMPLICIT NONE
106
107	!! * Routine accessibility
108	PRIVATE
109	PUBLIC bias_init, & !: Read in the bias namelist and the bias arrays
110	& tra_bias, & !: Estimate/apply bias on traces
111	& dyn_bias, & !: " density correction for pressure gradient.
112	& bias_opn, &
113	& bias_wrt
114
115
116	!! * Shared variables
117	!! * Private module variables
118	REAL(wp), PRIVATE :: &
119	& bias_time_unit_asm, & !: bias_asm units in s ( per day = 86400 s)
120	& bias_time_unit_rlx, & !: bias_rlx units in s ( 1 second)
121	& bias_time_unit_ofl, & !: bias_ofl units in s ( 1 second)
122	& t_rlx_mem, & !: time param for mem in bias_rlx model
123	& t_rlx_upd, & !: time param for update in bias_rlx model
124	!: (pct of incr for computation of bias)
125	& t_asm_mem, & !: time param for mem in bias_asm model
126	& t_asm_upd, & !: time param for update in bias_asm model
127	!: (pct of incr for computation of bias)
128	& fb_t_rlx, & !: parition of bias in T for rlx bias term
129	& fb_t_asm, & !: parition of bias in T for asm bias term
130	& fb_t_ofl, & !: parition of bias in T for ofl bias term
131	& fb_p_rlx, & !: parition of bias in P for rlx bias term
132	& fb_p_asm, & !: parition of bias in P for asm bias term
133	& fb_p_ofl, & !: parition of bias in P for ofl bias term
134	& fctamp, & !: amplification factor for T if inertial
135	& rn_maxlat_bias, & !: Max lat for latitudinal ramp
136	& rn_minlat_bias !: Min lat for latitudinal ramp
137
138	LOGICAL, PRIVATE :: lalloc
139	REAL(wp), PRIVATE, DIMENSION(:,:,:), ALLOCATABLE :: &
140	& tbias_asm, & !: Temperature bias field
141	& sbias_asm, & !: Salinity bias field
142	& tbias_rlx, & !: Temperature bias field
143	& sbias_rlx, & !: Salinity bias field
144	& tbias_asm_out, & !: Output temperature bias field
145	& sbias_asm_out, & !: Output salinity bias field
146	& tbias_rlx_out, & !: Output temperature bias field
147	& sbias_rlx_out, & !: Output salinity bias field
148	& tbias_p_out, & !: Output temperature bias field for P correction
149	& sbias_p_out !: Output salinity bias field for P correction
150
151	INTEGER, PRIVATE :: nn_lat_ramp ! choice of latitude dependent ramp
152	! for the pressure correction.
153	! 1:inertial ramp, 2:linear ramp, else:no ramp
154	LOGICAL, PRIVATE :: ln_bsyncro ! syncronous or assincrous bias correction
155	LOGICAL, PRIVATE :: ln_itdecay ! evolve bias correction at every time step.
156	LOGICAL, PRIVATE :: ln_incpc ! incremental pressure correction
157	LOGICAL, PRIVATE :: lat_ramp_incpc ! ramp for incremental pressure correction
158
159	REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef
160	REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef_stscale
161
162	INTEGER, PRIVATE :: &
163	& numbias_asm, & ! netcdf id of bias file from assim
164	& numbias_tot, & ! netcdf id of bias file with total bias
165	& nn_bias_itwrt ! time step for outputting bias pressure corr
166
167	CHARACTER(LEN=128), PRIVATE :: &
168	& cn_bias_asm, & ! name of bias file from assim
169	& cn_bias_tot ! name of bias with total/rlx bias
170
171	! Structure of input T and S bias offline (file informations, fields read)
172	TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) :: sf_tbias_ofl
173	TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) :: sf_sbias_ofl
174
175	TYPE(FLD_N), PRIVATE ::& ! information about the fields to be read
176	& sn_tbias_ofl, sn_sbias_ofl
177
178	CONTAINS
179
180	SUBROUTINE bias_init
181	!!----------------------------------------------------------------------
182	!! * ROUTINE bias_init *
183	!!
184	!! ** Purpose : Read in the bias namelist and read in the bias fields.
185	!!
186	!! ** Method : Read in the bias namelist and read in the bias fields.
187	!!
188	!! ** Action :
189	!!
190	!! History :
191	!! ! 08-05 (D. Lea) Initial version
192	!! ! 08-10 (M. Martin) Tidy
193	!! ! 09-03 (M. Balmaseda). Generalize to estimate the bias
194	!! from relax and offline bias term.
195	!! Introduce parameters to control the
196	!! model for the bias
197	!! (variables and time evolution)
198	!!----------------------------------------------------------------------
199
200	IMPLICIT NONE
201
202	!! * Local declarations
203
204	CHARACTER(len=100) :: cn_dir ! dir for location ofline bias
205	INTEGER :: ierror
206	INTEGER :: ios ! Local integer output status for namelist read
207	REAL(wp) :: eft_rlx, & ! efolding time (bias memory) in days
208	& eft_asm, & ! " "
209	& log2, &
210	& lenscl_bias !lengthscale of the pressure bias decay between minlat and maxlat.
211
212	NAMELIST/nambias/ ln_bias, ln_bias_asm, ln_bias_rlx, ln_bias_ofl, &
213	& ln_bias_ts_app, ln_bias_pc_app, &
214	& fb_t_asm, fb_t_rlx, fb_t_ofl, fb_p_asm, fb_p_rlx, fb_p_ofl, &
215	& eft_rlx, t_rlx_upd, eft_asm, t_asm_upd, nn_lat_ramp, &
216	& bias_time_unit_asm, bias_time_unit_rlx, bias_time_unit_ofl, &
217	& cn_bias_tot, cn_bias_asm, cn_dir, sn_tbias_ofl, sn_sbias_ofl, &
218	& ln_bsyncro, fctamp, rn_maxlat_bias, rn_minlat_bias, &
219	& nn_bias_itwrt, ln_itdecay, ln_incpc, lat_ramp_incpc
220
221
222	!-----------------------------------------------------------------------
223	! Read Namelist : bias interface
224	!-----------------------------------------------------------------------
225
226
227	REWIND( numnam_ref ) ! Namelist nambias in reference namelist : Bias pressure correction
228	READ ( numnam_ref, nambias, IOSTAT = ios, ERR = 901)
229	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in reference namelist', lwp )
230
231
232	! Set additional default values (note that most values are set in the reference namelist)
233
234	IF ( ln_asmiau ) nn_bias_itwrt = nitiaufin
235
236	! ... default values (NB: frequency positive => hours, negative => months)
237	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or !
238	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' !
239	sn_tbias_ofl = FLD_N( 'tbias_ofl' , -1. , 'tbias' , .TRUE. , .FALSE. , 'yearly', '', '', '' )
240	sn_sbias_ofl = FLD_N( 'sbias_ofl' , -1. , 'sbias' , .TRUE. , .FALSE. , 'yearly', '', '', '' )
241
242
243	REWIND( numnam_cfg ) ! Namelist nambias in configuration namelist : Bias pressure correction
244	READ ( numnam_cfg, nambias, IOSTAT = ios, ERR = 902 )
245	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in configuration namelist', lwp )
246	IF(lwm) WRITE ( numond, nambias )
247
248
249	IF ( ( .NOT. ln_bias_asm ) .AND. ( .NOT. ln_bias_ofl ) .AND. ( .NOT. ln_bias_rlx ) ) THEN
250	ln_bias_ts_app = .FALSE.
251	ln_bias_pc_app = .FALSE.
252	ln_bias = .FALSE.
253	ENDIF
254
255	! set up decay scales
256	log2 = LOG( 2.0_wp )
257	t_rlx_mem = EXP( - log2 * rdt / ( eft_rlx * rday ) )
258	t_asm_mem = EXP( - log2 * bias_time_unit_asm/ ( eft_asm * rday ) )
259
260	! Control print
261	IF(lwp) THEN
262	WRITE(numout,*)
263	WRITE(numout,*) 'bias_init : '
264	WRITE(numout,*) '~~~~~~~~~~~ '
265	WRITE(numout,*) ' Namelist nambias : '
266
267	WRITE(numout,*) ' Bias switches/options/variables '
268	WRITE(numout,*) ' bias main switch ln_bias = ',ln_bias
269	WRITE(numout,*) ' bias from assim ln_bias_asm = ',ln_bias_asm
270	WRITE(numout,*) ' bias from relax ln_bias_rlx = ',ln_bias_rlx
271	WRITE(numout,*) ' bias from offln ln_bias_ofl = ',ln_bias_ofl
272	WRITE(numout,*) ' bias T and S apply ln_bias_ts_app = ',ln_bias_ts_app
273	WRITE(numout,*) ' bias pressure correctn apply ln_bias_pc_app = ',ln_bias_pc_app
274	WRITE(numout,*) ' bias pressure correctn apply ln_bias_pc_app = ',ln_bias_pc_app
275	WRITE(numout,*) ' lat ramp for bias correction nn_lat_ramp = ',nn_lat_ramp
276	WRITE(numout,*) ' time step for writing bias fld nn_bias_itwrt = ',nn_bias_itwrt
277	WRITE(numout,*) ' evolve pcbias at each timestep ln_itdecay = ',ln_itdecay
278	WRITE(numout,*) ' incremental press. correction ln_incpc = ',ln_incpc
279	WRITE(numout,*) ' lat ramp for inc.press.corr lat_ramp_incpc = ',lat_ramp_incpc
280	WRITE(numout,*) ' Parameters for parition of bias term '
281	WRITE(numout,*) ' fb_t_rlx = ',fb_t_rlx
282	WRITE(numout,*) ' fb_t_asm = ',fb_t_asm
283	WRITE(numout,*) ' fb_t_ofl = ',fb_t_ofl
284	WRITE(numout,*) ' fb_p_rlx = ',fb_p_rlx
285	WRITE(numout,*) ' fb_p_asm = ',fb_p_asm
286	WRITE(numout,*) ' fb_p_ofl = ',fb_p_ofl
287	WRITE(numout,*) ' Parameters for time evolution of bias '
288	WRITE(numout,) ' Rlx efolding time (mem) eft_rlx,t_rlx_mem = ', eft_rlx, t_rlx_mem, 1. - log2 rdt / (eft_rlx * rday)
289	WRITE(numout,*) ' uptdate factor t_rlx_upd = ',t_rlx_upd
290	WRITE(numout,) ' Asm efolding time (mem) eft_asm,t_asm_mem = ', eft_asm, t_asm_mem, 1. - log2 rdt / (eft_asm * rday)
291	WRITE(numout,*) ' uptdate factor t_asm_upd = ',t_asm_upd
292	WRITE(numout,*) ' Filenames and input structures'
293	WRITE(numout,*) ' bias_tot filename cn_bias_to = ',cn_bias_tot
294	WRITE(numout,*) ' bias_asm filename cn_bias_asm = ',cn_bias_asm
295	WRITE(numout,*) ' bias_asm time unit (secs) bias_time_unit_asm = ',bias_time_unit_asm
296	WRITE(numout,*) ' structure Tem bias ofl sn_tbias_ofl = ',sn_tbias_ofl
297	WRITE(numout,*) ' structure Sal bias ofl sn_sbias_ofl = ',sn_sbias_ofl
298
299	IF ( ( (.NOT. ln_tsd_tradmp) .OR. (.NOT. ln_tradmp) ) .AND. ln_bias_rlx ) &
300	& CALL ctl_stop (' lk_dtatem, lk_dtasal and lk_tradmp need to be true with ln_bias_rlx' )
301
302	ENDIF
303	IF( .NOT. ln_bias ) RETURN
304
305	IF( .NOT. lalloc ) THEN
306
307	ALLOCATE( tbias(jpi,jpj,jpk) , &
308	& sbias(jpi,jpj,jpk) , &
309	& tbias_p(jpi,jpj,jpk), &
310	& sbias_p(jpi,jpj,jpk), &
311	& tbias_i(jpi,jpj,jpk), &
312	& sbias_i(jpi,jpj,jpk), &
313	& rhd_pc(jpi,jpj,jpk) , &
314	& gru_pc(jpi,jpj) , &
315	& grv_pc(jpi,jpj) )
316
317	ALLOCATE( fbcoef(jpi,jpj), fbcoef_stscale(jpi,jpj) )
318
319	IF( ln_bias_asm ) ALLOCATE( tbias_asm(jpi,jpj,jpk), &
320	& sbias_asm(jpi,jpj,jpk), &
321	tbias_asm_out(jpi,jpj,jpk), &
322	sbias_asm_out(jpi,jpj,jpk), &
323	tbias_p_out(jpi,jpj,jpk), &
324	sbias_p_out(jpi,jpj,jpk) )
325
326	IF( ln_bias_rlx ) ALLOCATE( tbias_rlx(jpi,jpj,jpk), &
327	& sbias_rlx(jpi,jpj,jpk), &
328	tbias_rlx_out(jpi,jpj,jpk), &
329	sbias_rlx_out(jpi,jpj,jpk) )
330
331	IF( ln_incpc ) ALLOCATE( tbias_asm_stscale(jpi,jpj,jpk), &
332	& sbias_asm_stscale(jpi,jpj,jpk))
333
334	lalloc = .TRUE.
335
336	ENDIF
337
338	IF( ln_bias_ofl ) THEN ! set sf_tbias_ofl and sf_sbias_ofl strctrs
339	!
340	! tbias
341	!
342	ALLOCATE( sf_tbias_ofl(1), STAT=ierror )
343	IF( ierror > 0 ) THEN
344	CALL ctl_stop( 'bias_init: unable to allocate sf_tbias_ofl structure' ) ; RETURN
345	ENDIF
346	ALLOCATE( sf_tbias_ofl(1)%fnow(jpi,jpj,jpk) )
347	ALLOCATE( sf_tbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
348
349	! fill structure with values and control print
350	CALL fld_fill( sf_tbias_ofl, (/ sn_tbias_ofl /), cn_dir, 'bias_init', 'Offline T bias term ', 'nam_tbias_ofl' )
351	!
352	! salinity bias
353	!
354	ALLOCATE( sf_sbias_ofl(1), STAT=ierror )
355	IF( ierror > 0 ) THEN
356	CALL ctl_stop( 'bias_init: unable to allocate sf_sbias_ofl structure' ) ; RETURN
357	ENDIF
358	ALLOCATE( sf_sbias_ofl(1)%fnow(jpi,jpj,jpk) )
359	ALLOCATE( sf_sbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
360
361	! fill structure with values and control print
362	CALL fld_fill( sf_sbias_ofl, (/ sn_sbias_ofl /), cn_dir, 'bias_init', 'Offline S bias term ', 'nam_sbias_ofl' )
363
364	ENDIF
365
366	! Read total bias
367	IF ( ln_bias ) THEN
368	tbias(:,:,:) = 0.0_wp
369	sbias(:,:,:) = 0.0_wp
370	tbias_p(:,:,:) = 0.0_wp
371	sbias_p(:,:,:) = 0.0_wp
372	gru_pc(:,:) = 0.0_wp
373	grv_pc(:,:) = 0.0_wp
374	IF ( ln_bias_rlx ) THEN
375	tbias_rlx(:,:,:) = 0.0_wp
376	sbias_rlx(:,:,:) = 0.0_wp
377	ENDIF
378	IF ( ln_bias_asm ) THEN !now rlx and asm bias in same file
379	tbias_asm(:,:,:) = 0.0_wp
380	sbias_asm(:,:,:) = 0.0_wp
381	ENDIF
382
383	IF ( ln_incpc ) THEN !incr pressure correction
384	tbias_asm_stscale(:,:,:) = 0.0_wp
385	sbias_asm_stscale(:,:,:) = 0.0_wp
386	ENDIF
387
388
389	numbias_tot = 0
390	! Get bias from file and prevent fail if the file does not exist
391	IF(lwp) WRITE(numout,*) 'Opening ',TRIM( cn_bias_tot )
392	CALL iom_open( cn_bias_tot, numbias_tot, ldstop=.FALSE. )
393
394	IF ( numbias_tot > 0 ) THEN
395	! Could check validity time of bias fields here...
396	! Get the T and S bias data
397	IF(lwp) WRITE(numout,*) 'Reading bias fields from tot...'
398
399	!Search for bias from relaxation term if needed. Use same file
400	IF ( ln_bias_rlx ) THEN
401	IF(lwp) WRITE(numout,*) 'Reading bias fields for bias rlx from file ',cn_bias_tot
402	IF( iom_varid( numbias_tot, 'tbias_rlx' ) > 0 ) THEN
403	! Get the T and S bias data
404	CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_rlx', tbias_rlx )
405	CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_rlx', sbias_rlx )
406	ELSE
407	CALL ctl_stop( 'Bias relaxation variables not found in ',cn_bias_tot )
408	ENDIF
409	ENDIF
410
411
412	!Search for bias from assim term if needed. Use same file
413	IF ( ln_bias_asm .and. .not.ln_bsyncro ) THEN
414	IF(lwp) WRITE(numout,*) 'Reading a-syncro bias fields for bias asm from file ',cn_bias_tot
415	IF( iom_varid( numbias_tot, 'tbias_asm' ) > 0 ) THEN
416	! Get the T and S bias data
417	CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm', tbias_asm )
418	CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm', sbias_asm )
419	ELSE
420	CALL ctl_stop( 'Bias assim variables not found in ',cn_bias_tot )
421	ENDIF
422
423
424	IF ( ln_incpc ) THEN
425	IF(lwp) WRITE(numout,*) 'Reading short time scale bias correction fields for bias asm from file ',cn_bias_tot
426	IF( iom_varid( numbias_tot, 'tbias_asm_stscale' ) > 0 ) THEN
427	! Get the T and S bias data
428	CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm_stscale', tbias_asm_stscale )
429	CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm_stscale', sbias_asm_stscale )
430	ELSE
431	CALL ctl_stop( 'Short time scale bias assim variables not found in ',cn_bias_tot )
432	ENDIF
433	ENDIF
434
435	ENDIF
436
437
438
439	! Close the file
440	CALL iom_close(numbias_tot)
441
442	ELSE
443	IF(lwp) WRITE(numout,*) 'No bias file found so T and S bias fields are set to zero'
444	ENDIF
445
446	ENDIF
447
448	! for the time being, the bias_asm is read in the same file as
449	! bias_rlx
450	! Implications: bias_asm is estimated/evolved in time in the second outer
451	! loop only, when the assimilation increments are ready.
452	! bias_asm is kept cte during the first outer loop.
453	! => Assyncronous bias correction.
454	! Alternative: Syncronous bias correction:
455	! bias_asm estimated/evolved in the first outer loop
456	! with the asm incrments of the previous cycle.
457	! bias_asm kept cte during the second outer loop.
458	! Implication: bias_asm should be estimated really in the
459	! inner loop.
460	IF ( ln_bsyncro ) THEN
461	! Read bias from assimilation from a separate file
462	IF ( ln_bias_asm ) THEN
463	tbias_asm(:,:,:) = 0.0_wp
464	sbias_asm(:,:,:) = 0.0_wp
465	numbias_asm = 0
466	! Get bias from file and prevent fail if the file does not exist
467	IF(lwp) WRITE(numout,*) 'Opening file for syncro assim bias ',TRIM( cn_bias_asm )
468	CALL iom_open( cn_bias_asm, numbias_asm, ldstop=.FALSE. )
469
470	IF ( numbias_asm > 0 ) THEN
471	! Could check validity time of bias fields here...
472
473	! Get the T and S bias data
474	IF(lwp) WRITE(numout,*) 'Reading syncro bias fields from asm from file ',cn_bias_asm
475	CALL iom_get( numbias_asm, jpdom_autoglo, 'tbias_asm', tbias_asm )
476	CALL iom_get( numbias_asm, jpdom_autoglo, 'sbias_asm', sbias_asm )
477
478	! this is only applicable if tbias_asm were to be calculated in the inner loop
479	tbias_asm(:,:,:) = tbias_asm(:,:,:) * rdt / bias_time_unit_asm
480	sbias_asm(:,:,:) = sbias_asm(:,:,:) * rdt / bias_time_unit_asm
481
482	! Close the file
483	CALL iom_close(numbias_asm)
484
485	ELSE
486	IF(lwp) WRITE(numout,*) 'No bias file found from asm so T and S bias fields are set to zero'
487	ENDIF
488
489	ENDIF
490
491	ENDIF
492
493	!latitudinal dependence of partition coeficients. Adhoc
494	IF ( nn_lat_ramp == 1 ) THEN
495	! Use the inertial ramp.
496	lenscl_bias = ( rn_maxlat_bias - rn_minlat_bias )*2._wp
497	WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )
498	fbcoef(:,:) = 0._wp
499	ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias )
500	fbcoef(:,:) = 1._wp
501	ELSEWHERE
502	fbcoef(:,:) = 1._wp - exp( -( abs( gphit(:,:) ) - rn_minlat_bias ) &
503	* ( abs( gphit(:,:) ) - rn_minlat_bias ) / lenscl_bias )
504	ENDWHERE
505	ELSEIF ( nn_lat_ramp == 2 ) THEN
506	! Use a linear ramp consist with the geostrophic velocity balance ramp in NEMOVAR
507
508	WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )
509	fbcoef(:,:) = 0._wp
510	ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias )
511	fbcoef(:,:) = 1._wp
512	ELSEWHERE
513	fbcoef(:,:) = 1._wp - ((rn_maxlat_bias - abs( gphit(:,:)))/(rn_maxlat_bias - rn_minlat_bias))
514	ENDWHERE
515	ELSE
516	fbcoef(:,:) = 0.0_wp
517	fctamp = 0.0_wp
518	fbcoef_stscale(:,:) = 0.0_wp
519	ENDIF
520
521
522	IF ( lat_ramp_incpc) THEN
523	! not sure if this should be a special case of nn_lat_ramp == 2
524	minlat_bias = 3.0_wp
525	maxlat_bias = 8.0_wp
526	WHERE ( abs( gphit(:,:) ) <= minlat_bias )
527	fbcoef_stscale(:,:)=0._wp
528	ELSEWHERE ( abs( gphit(:,:) ) >= maxlat_bias )
529	fbcoef_stscale(:,:)=1._wp
530	ELSEWHERE
531	fbcoef_stscale(:,:)=1._wp - ((maxlat_bias - abs( gphit(:,:)))/(maxlat_bias-minlat_bias))
532	ENDWHERE
533	ENDIF
534
535
536	END SUBROUTINE bias_init
537
538	SUBROUTINE tra_bias ( kt )
539	!!----------------------------------------------------------------------
540	!! * ROUTINE tra_bias *
541	!!
542	!! ** Purpose : Update bias field every time step
543	!!
544	!! ** Method : add contributions to bias from 3 terms
545	!!
546	!! ** Action : Bias from assimilation (read in bias_init)
547	!! Bias from relaxation term is estimated according to
548	!! the prescribed time evolution of the bias
549	!! Bias from ofl term is read from external file
550	!! The difference contributions are added and the partition
551	!! into direct bias in T/S and pressure perfomed.
552	!!
553	!! History : 09-03 (M. Balmaseda)
554	!!----------------------------------------------------------------------
555	!! called every timestep after dta_sst if ln_bias true.
556
557	IMPLICIT NONE
558
559	!! * Arguments
560	INTEGER, INTENT(IN) :: kt ! ocean time-step index
561	!! * Local variables
562	INTEGER :: ji,jj,jk, it ! local loop index
563	REAL(wp) :: tsclf ! time saling factor
564	REAL(wp) :: fb_t_asm_max, fb_t_rlx_max, fb_t_ofl_max
565	REAL(wp) :: ztfrac, ztsday
566	REAL(wp) :: zfrac, zfrac1 ! temporal weights for inst pcbias (names could be changed)
567	REAL(wp) :: ztscale ! reduce the inst pcbias by this amount per 24 hours
568	REAL(wp) :: zwgt ! Weight for the inst pcorr term
569	REAL(wp) :: zdecay ! used in inst pcorr
570	REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2
571
572	IF ( .NOT. ln_bias ) RETURN
573	fb_t_rlx_max = MIN(fb_t_rlx*fctamp,1.0_wp)
574	fb_t_asm_max = MIN(fb_t_asm*fctamp,1.0_wp)
575	fb_t_ofl_max = MIN(fb_t_ofl*fctamp,1.0_wp)
576
577	tbias(:,:,:) = 0.0_wp
578	sbias(:,:,:) = 0.0_wp
579	tbias_p(:,:,:) = 0.0_wp
580	sbias_p(:,:,:) = 0.0_wp
581
582	ztscale = 0.1_wp
583	zwgt = 1.0_wp
584	tbias_i(:,:,:) = 0.0_wp
585	sbias_i(:,:,:) = 0.0_wp
586
587	IF ( ln_bias_asm ) THEN
588
589	tsclf = 1
590	IF ( .NOT.ln_bsyncro ) tsclf = rdt / bias_time_unit_asm
591	zcof1(:,:) = tsclf * ( ( 1.0_wp - fbcoef(:,:) ) * fb_t_asm + &
592	& fbcoef(:,:) * fb_t_asm_max )
593	zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_asm
594
595	IF ( ln_itdecay ) THEN !decay the pressure correction at each time step
596
597	ztsday = rday / real(rdt)
598
599	zdecay = (1-ztscale)**(1/real(ztsday)) ! used in ipc
600	zfrac1 = max(0.0_wp, zdecay**real(kt)) ! used in ipc
601
602
603	IF( ln_asmiau .and. ln_trainc ) THEN ! nudge in increments and decay historical contribution
604
605
606	IF ( kt <= nitiaufin ) THEN ! During IAU calculate the fraction of increments to apply at each time step
607
608	ztfrac = real(kt) / real(nitiaufin) ! nudging factor during the IAU
609
610
611	IF (lwp) THEN
612	WRITE(numout,*) 'tra_bias : bias weights'
613	WRITE(numout,*) '~~~~~~~~~~~~'
614	WRITE(numout,* ) "proportion of increment applied in pcbias ",ztfrac
615	WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
616	ENDIF
617
618	DO jk = 1, jpkm1
619	tbias(:,:,jk) = tbias(:,:,jk) + &
620	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
621	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
622	sbias(:,:,jk) = sbias(:,:,jk) + &
623	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
624	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
625	tbias_p(:,:,jk) = tbias_p(:,:,jk) + &
626	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
627	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
628	sbias_p(:,:,jk) = sbias_p(:,:,jk) + &
629	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
630	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
631	ENDDO
632
633
634	IF (ln_incpc) THEN
635
636	DO jk = 1, jpkm1
637
638	tbias_i(:,:,jk) = ( t_bkginc(:,:,jk) * zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) ) &
639	& - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
640	sbias_i(:,:,jk) = ( s_bkginc(:,:,jk) * zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) ) &
641	& - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
642
643	ENDDO
644
645	ENDIF
646
647	IF ( .not.ln_bsyncro ) THEN
648	IF ( kt == nn_bias_itwrt ) THEN
649	DO jk = 1, jpkm1
650	tbias_asm_out(:,:,jk) = t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
651	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
652	sbias_asm_out(:,:,jk) = t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
653	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
654	END DO
655	ENDIF
656	ENDIF
657
658	! update the historical component with the increments at the end of the IAU
659	IF ( kt == nitiaufin ) THEN
660	DO jk = 1, jpkm1
661	tbias_asm(:,:,jk) = t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
662	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
663	sbias_asm(:,:,jk) = t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
664	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
665	END DO
666	ENDIF
667
668	ELSE ! decay pressure correction from combined historical component and increments after IAU
669
670	ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin)) ! decay from end of IAU
671
672
673	DO jk = 1, jpkm1
674	tbias(:,:,jk) = tbias(:,:,jk) + &
675	& ( ztfrac * tbias_asm(:,:,jk) ) * zcof1(:,:)
676	sbias(:,:,jk) = sbias(:,:,jk) + &
677	& ( ztfrac * sbias_asm(:,:,jk) ) * zcof1(:,:)
678	tbias_p(:,:,jk) = tbias_p(:,:,jk) + &
679	& ( ztfrac * tbias_asm(:,:,jk) ) * zcof2(:,:)
680	sbias_p(:,:,jk) = sbias_p(:,:,jk) + &
681	& ( ztfrac * sbias_asm(:,:,jk) ) * zcof2(:,:)
682
683	ENDDO
684
685	IF (ln_incpc) THEN
686
687	zfrac = max(0.0_wp, zdecay**real(kt-nitiaufin) )
688
689	DO jk = 1, jpkm1
690	tbias_i(:,:,jk) = ( t_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) ) &
691	& - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
692	sbias_i(:,:,jk) = ( s_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) ) &
693	& - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
694	ENDDO
695
696	ENDIF
697
698	IF (lwp) THEN
699	WRITE(numout,*) 'tra_bias : bias weights'
700	WRITE(numout,*) '~~~~~~~~~~~~'
701	WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac
702	ENDIF
703
704	IF ( ln_trainc .and. .not.ln_bsyncro ) THEN
705	IF ( kt == nn_bias_itwrt ) THEN
706	DO jk = 1, jpkm1
707	tbias_asm_out(:,:,jk) = ztfrac * tbias_asm(:,:,jk)
708	sbias_asm_out(:,:,jk) = ztfrac * sbias_asm(:,:,jk)
709	END DO
710	ENDIF
711	ENDIF
712
713	ENDIF
714
715
716	ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts)
717
718	DO jk = 1, jpkm1
719	tbias(:,:,jk) = tbias(:,:,jk) + &
720	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) ) * zcof1(:,:)
721	sbias(:,:,jk) = sbias(:,:,jk) + &
722	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) ) * zcof1(:,:)
723	tbias_p(:,:,jk) = tbias_p(:,:,jk) + &
724	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) ) * zcof2(:,:)
725	sbias_p(:,:,jk) = sbias_p(:,:,jk) + &
726	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) ) * zcof2(:,:)
727	ENDDO
728
729	IF (lwp) THEN
730	WRITE(numout,*) 'tra_bias : bias weights'
731	WRITE(numout,*) '~~~~~~~~~~~~'
732	WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
733	ENDIF
734
735	ENDIF
736
737	ELSE ! maintain a constant pressure correction
738
739	DO jk = 1, jpkm1
740	tbias(:,:,jk) = tbias(:,:,jk) + tbias_asm(:,:,jk) * zcof1(:,:)
741	sbias(:,:,jk) = sbias(:,:,jk) + sbias_asm(:,:,jk) * zcof1(:,:)
742	tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_asm(:,:,jk) * zcof2(:,:)
743	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_asm(:,:,jk) * zcof2(:,:)
744	END DO
745
746	IF( ln_asmiau .and. ln_trainc .and. .not.ln_bsyncro ) THEN
747	! if last outer loop (ln_asmiau=true and ln_trainc=true). t/sbias_asm
748	! is updated, only once (end of run) taking into account units.
749	IF ( kt == nn_bias_itwrt ) THEN
750	IF(lwp) WRITE(numout,*)' estimating asm bias at timestep: ',kt
751	DO jk = 1, jpkm1
752	tbias_asm_out(:,:,jk) = t_asm_mem * tbias_asm(:,:,jk) + &
753	& t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
754	sbias_asm_out(:,:,jk) = t_asm_mem * sbias_asm(:,:,jk) + &
755	& t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
756	END DO
757	ENDIF
758	ENDIF
759
760	ENDIF
761
762	ENDIF
763
764
765	#if defined key_tradmp
766	! Time evolution of bias from relaxation
767	IF ( ln_bias_rlx ) THEN
768	tbias_rlx(:,:,:) = t_rlx_mem * tbias_rlx(:,:,:) + &
769	& t_rlx_upd * ttrdmp(:,:,:) * rdt
770	sbias_rlx(:,:,:) = t_rlx_mem * sbias_rlx(:,:,:) + &
771	& t_rlx_upd * strdmp(:,:,:) * rdt
772	zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_rlx + &
773	& fbcoef(:,:) * fb_t_rlx_max
774	zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_rlx
775	DO jk = 1, jpkm1
776	tbias(:,:,jk) = tbias(:,:,jk) + tbias_rlx(:,:,jk) * zcof1(:,:)
777	sbias(:,:,jk) = sbias(:,:,jk) + sbias_rlx(:,:,jk) * zcof1(:,:)
778	tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_rlx(:,:,jk) * zcof2(:,:)
779	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_rlx(:,:,jk) * zcof2(:,:)
780	ENDDO
781	IF ( kt == nn_bias_itwrt ) THEN
782	tbias_rlx_out(:,:,:) = tbias_rlx(:,:,:)
783	sbias_rlx_out(:,:,:) = sbias_rlx(:,:,:)
784	ENDIF
785	ENDIF
786	#endif
787	! ofline bias
788	IF ( kt == nit000 ) THEN
789	IF(lwp) WRITE(numout,*) ' tra_bias: ln_bias_ofl = ',ln_bias_ofl
790	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~'
791	ENDIF
792	IF ( ln_bias_ofl ) THEN
793	IF(lwp) WRITE(numout,*) 'reading offline bias'
794	CALL fld_read( kt, 1, sf_tbias_ofl )
795	CALL fld_read( kt, 1, sf_sbias_ofl )
796
797	zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_ofl + &
798	& fbcoef(:,:) * fb_t_ofl_max
799	zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_ofl
800	DO jk = 1, jpkm1
801	tbias(:,:,jk) = tbias(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
802	sbias(:,:,jk) = sbias(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
803	tbias_p(:,:,jk) = tbias_p(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
804	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
805	ENDDO
806	ENDIF
807
808
809	!apply bias on tracers if needed
810	IF ( ln_bias_ts_app ) THEN
811
812	! Add the bias directely to T/s
813	tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + tmask(:,:,:) * tbias(:,:,:) / rdt
814	tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + tmask(:,:,:) * sbias(:,:,:) / rdt
815
816	! lateral boundary conditions (is this needed?)
817	CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1.0_wp )
818	CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1.0_wp )
819
820	ENDIF
821
822	IF ( kt == nn_bias_itwrt ) THEN
823	tbias_p_out(:,:,:) = tbias_p(:,:,:)
824	sbias_p_out(:,:,:) = sbias_p(:,:,:)
825	ENDIF
826
827	END SUBROUTINE tra_bias
828
829	SUBROUTINE dyn_bias( kt )
830	!!----------------------------------------------------------------------
831	!! * ROUTINE dyn_bias *
832	!!
833	!! ** Purpose : Computes rhd_pc, gru/v_pc bias corrected
834	!! for hydrostatic pressure gradient
835	!! depending on time step (semi-implicit/centered)
836	!! If partial steps computes bottom pressure gradient.
837	!! These correction terms will affect only the dynamical
838	!! component (pressure gradient calculation), but not
839	!! the isopycnal calculation for the lateral mixing.
840	!!
841	!! ** Method : At this stage of the computation, ta and sa are the
842	!! after temperature and salinity. If semi-implicit, these
843	!! are used to compute rho and bottom pressure gradient.
844	!! If centered, tb,sb are used instead.
845	!! If bias key is activated, the temperature,salinity are
846	!! bias corrected in the calculation of the density fields
847	!! used in the pressure gradient calculation.
848	!!
849	!!
850	!! ** Action : - rhd_pc ready. rhop will be overwriten later
851	!! - if ln_zps, bottom density gradients gru/v_pc ready.
852	!!----------------------------------------------------------------------
853	!!
854	!! * Arguments
855	INTEGER, INTENT(IN) :: kt ! ocean time-step index
856	!! * Local variables
857	REAL(wp) :: tsw(jpi,jpj,jpk,jpts)
858	!!
859	!!----------------------------------------------------------------------
860	!
861	! gtu,gsu,gtv,gsv rhop will be overwritten later in step.
862	!
863	IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg
864	tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:)
865	tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:)
866
867	IF (ln_incpc) THEN
868	tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) -tbias_i(:,:,:) - tbias_p(:,:,:)
869	tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) -sbias_i(:,:,:) - sbias_p(:,:,:)
870	ENDIF
871
872	ELSE
873	tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:)
874	tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:)
875
876	IF (ln_incpc) THEN
877	tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_i(:,:,:) - tbias_p(:,:,:)
878	tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_i(:,:,:) - sbias_p(:,:,:)
879	ENDIF
880
881	ENDIF
882
883	CALL eos( tsw, rhd_pc, rhop )
884
885	! is this needed?
886	CALL lbc_lnk( rhd_pc, 'T', 1.0_wp )
887
888	! Partial steps: now horizontal gradient of t,s,rd
889	! at the bottom ocean level
890	IF( ln_zps ) THEN
891	CALL zps_hde( kt, jpts, tsw, gtsu, gtsv, &
892	& rhd_pc, gru_pc , grv_pc )
893	ENDIF
894
895	END SUBROUTINE dyn_bias
896
897	SUBROUTINE bias_opn( kt )
898	!!---------------------------------------------------------------------
899	!! * ROUTINE bias_opn *
900	!!
901	!! ** Purpose : open bias restart file following the same logic as the
902	!! standard restarts.
903	!!----------------------------------------------------------------------
904	!! * Arguments
905	INTEGER, INTENT(IN) :: kt ! ocean time-step
906	!! * Local variables
907	CHARACTER(LEN=20) :: clbkt ! ocean time-step deine as a character
908	CHARACTER(LEN=50) :: clbias_tot ! total bias restart file name
909	!!----------------------------------------------------------------------
910	!
911	IF( lrst_oce .AND. .NOT.lrst_bias ) THEN ! create bias file
912	IF( nitend > 999999999 ) THEN ; WRITE(clbkt, * ) nitend
913	ELSE ; WRITE(clbkt, '(i8.8)') nitend
914	ENDIF
915	clbias_tot = TRIM(cexper)//"_"//TRIM(ADJUSTL(clbkt))//"_"//TRIM(cn_bias_tot)
916	IF(lwp) THEN
917	WRITE(numout,*)
918	SELECT CASE ( jprstlib )
919	CASE ( jprstdimg ) ; WRITE(numout,*) ' open tot bias restart binary file: '//clbias_tot
920	CASE DEFAULT ; WRITE(numout,*) ' open tot bias restart NetCDF file: '//clbias_tot
921	END SELECT
922	ENDIF
923	CALL iom_open( clbias_tot, numbias_tot , ldwrt = .TRUE., kiolib = jprstlib )
924	lrst_bias=.TRUE.
925
926	ENDIF
927	!
928	END SUBROUTINE bias_opn
929
930	SUBROUTINE bias_wrt( kt )
931	!!---------------------------------------------------------------------
932	!! * ROUTINE bias_wrt *
933	!!
934	!! ** Purpose : Write bias restart fields in the format corresponding to jprstlib
935	!!
936	!! ** Method : Write in numbias_tot when kt == nitend in output
937	!! file, save fields which are necessary for restart.
938	!!
939	!! Changed the timestep for writing to nitend rather than nitrst as this causes a
940	!! problem if the nstock list is used to determine the restart output times and
941	!! means that the bias is not output at all. M. Martin. 08/2011.
942	!! Need to check with Magdalena that this is ok for ECMWF.
943	!!----------------------------------------------------------------------
944	!! * Arguments
945	INTEGER, INTENT(IN) :: kt ! ocean time-step
946	!!----------------------------------------------------------------------
947	! ! the begining of the run [s]
948
949	IF ( ln_bias_rlx ) THEN
950	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_rlx' , tbias_rlx_out )
951	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_rlx' , sbias_rlx_out )
952	ENDIF
953
954	IF ( ln_bias_asm ) THEN
955	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm' , tbias_asm_out )
956	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm' , sbias_asm_out )
957	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_p' , tbias_p_out )
958	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p' , sbias_p_out )
959	ENDIF
960
961	IF( kt == nitend ) THEN
962	CALL iom_close( numbias_tot ) ! close the restart file (only at last time step)
963	lrst_bias = .FALSE.
964	ENDIF
965	!
966	END SUBROUTINE bias_wrt
967
968	END MODULE bias

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

Download in other formats: