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 @ 7092

Last change on this file since 7092 was 7092, checked in by isabella, 8 years ago
declaring variables
File size: 44.1 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	& tbias_asm_stscale, & !: Short time scale temperature bias field
151	& sbias_asm_stscale !: Short time scale salinity bias field
152
153	INTEGER, PRIVATE :: nn_lat_ramp ! choice of latitude dependent ramp
154	! for the pressure correction.
155	! 1:inertial ramp, 2:linear ramp, else:no ramp
156	LOGICAL, PRIVATE :: ln_bsyncro ! syncronous or assincrous bias correction
157	LOGICAL, PRIVATE :: ln_itdecay ! evolve bias correction at every time step.
158	LOGICAL, PRIVATE :: ln_incpc ! incremental pressure correction
159
160	REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef
161	REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef_stscale
162
163	INTEGER, PRIVATE :: &
164	& numbias_asm, & ! netcdf id of bias file from assim
165	& numbias_tot, & ! netcdf id of bias file with total bias
166	& nn_bias_itwrt ! time step for outputting bias pressure corr
167
168	CHARACTER(LEN=128), PRIVATE :: &
169	& cn_bias_asm, & ! name of bias file from assim
170	& cn_bias_tot ! name of bias with total/rlx bias
171
172	! Structure of input T and S bias offline (file informations, fields read)
173	TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) :: sf_tbias_ofl
174	TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) :: sf_sbias_ofl
175
176	TYPE(FLD_N), PRIVATE ::& ! information about the fields to be read
177	& sn_tbias_ofl, sn_sbias_ofl
178
179	CONTAINS
180
181	SUBROUTINE bias_init
182	!!----------------------------------------------------------------------
183	!! * ROUTINE bias_init *
184	!!
185	!! ** Purpose : Read in the bias namelist and read in the bias fields.
186	!!
187	!! ** Method : Read in the bias namelist and read in the bias fields.
188	!!
189	!! ** Action :
190	!!
191	!! History :
192	!! ! 08-05 (D. Lea) Initial version
193	!! ! 08-10 (M. Martin) Tidy
194	!! ! 09-03 (M. Balmaseda). Generalize to estimate the bias
195	!! from relax and offline bias term.
196	!! Introduce parameters to control the
197	!! model for the bias
198	!! (variables and time evolution)
199	!!----------------------------------------------------------------------
200
201	IMPLICIT NONE
202
203	!! * Local declarations
204
205	CHARACTER(len=100) :: cn_dir ! dir for location ofline bias
206	INTEGER :: ierror
207	INTEGER :: ios ! Local integer output status for namelist read
208	REAL(wp) :: eft_rlx, & ! efolding time (bias memory) in days
209	& eft_asm, & ! " "
210	& log2, &
211	& lenscl_bias !lengthscale of the pressure bias decay between minlat and maxlat.
212
213	NAMELIST/nambias/ ln_bias, ln_bias_asm, ln_bias_rlx, ln_bias_ofl, &
214	& ln_bias_ts_app, ln_bias_pc_app, &
215	& fb_t_asm, fb_t_rlx, fb_t_ofl, fb_p_asm, fb_p_rlx, fb_p_ofl, &
216	& eft_rlx, t_rlx_upd, eft_asm, t_asm_upd, nn_lat_ramp, &
217	& bias_time_unit_asm, bias_time_unit_rlx, bias_time_unit_ofl, &
218	& cn_bias_tot, cn_bias_asm, cn_dir, sn_tbias_ofl, sn_sbias_ofl, &
219	& ln_bsyncro, fctamp, rn_maxlat_bias, rn_minlat_bias, &
220	& nn_bias_itwrt, ln_itdecay, ln_incpc
221
222
223	!-----------------------------------------------------------------------
224	! Read Namelist : bias interface
225	!-----------------------------------------------------------------------
226
227
228	REWIND( numnam_ref ) ! Namelist nambias in reference namelist : Bias pressure correction
229	READ ( numnam_ref, nambias, IOSTAT = ios, ERR = 901)
230	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in reference namelist', lwp )
231
232
233	! Set additional default values (note that most values are set in the reference namelist)
234
235	IF ( ln_asmiau ) nn_bias_itwrt = nitiaufin
236
237	! ... default values (NB: frequency positive => hours, negative => months)
238	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or !
239	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' !
240	sn_tbias_ofl = FLD_N( 'tbias_ofl' , -1. , 'tbias' , .TRUE. , .FALSE. , 'yearly', '', '', '' )
241	sn_sbias_ofl = FLD_N( 'sbias_ofl' , -1. , 'sbias' , .TRUE. , .FALSE. , 'yearly', '', '', '' )
242
243
244	REWIND( numnam_cfg ) ! Namelist nambias in configuration namelist : Bias pressure correction
245	READ ( numnam_cfg, nambias, IOSTAT = ios, ERR = 902 )
246	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in configuration namelist', lwp )
247	IF(lwm) WRITE ( numond, nambias )
248
249
250	IF ( ( .NOT. ln_bias_asm ) .AND. ( .NOT. ln_bias_ofl ) .AND. ( .NOT. ln_bias_rlx ) ) THEN
251	ln_bias_ts_app = .FALSE.
252	ln_bias_pc_app = .FALSE.
253	ln_bias = .FALSE.
254	ENDIF
255
256	! set up decay scales
257	log2 = LOG( 2.0_wp )
258	t_rlx_mem = EXP( - log2 * rdt / ( eft_rlx * rday ) )
259	t_asm_mem = EXP( - log2 * bias_time_unit_asm/ ( eft_asm * rday ) )
260
261	! Control print
262	IF(lwp) THEN
263	WRITE(numout,*)
264	WRITE(numout,*) 'bias_init : '
265	WRITE(numout,*) '~~~~~~~~~~~ '
266	WRITE(numout,*) ' Namelist nambias : '
267
268	WRITE(numout,*) ' Bias switches/options/variables '
269	WRITE(numout,*) ' bias main switch ln_bias = ',ln_bias
270	WRITE(numout,*) ' bias from assim ln_bias_asm = ',ln_bias_asm
271	WRITE(numout,*) ' bias from relax ln_bias_rlx = ',ln_bias_rlx
272	WRITE(numout,*) ' bias from offln ln_bias_ofl = ',ln_bias_ofl
273	WRITE(numout,*) ' bias T and S apply ln_bias_ts_app = ',ln_bias_ts_app
274	WRITE(numout,*) ' bias pressure correctn apply ln_bias_pc_app = ',ln_bias_pc_app
275	WRITE(numout,*) ' bias pressure correctn apply ln_bias_pc_app = ',ln_bias_pc_app
276	WRITE(numout,*) ' lat ramp for bias correction nn_lat_ramp = ',nn_lat_ramp
277	WRITE(numout,*) ' time step for writing bias fld nn_bias_itwrt = ',nn_bias_itwrt
278	WRITE(numout,*) ' evolve pcbias at each timestep ln_itdecay = ',ln_itdecay
279	WRITE(numout,*) ' incremental press. correction ln_incpc = ',ln_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	& rhd_pc(jpi,jpj,jpk) , &
312	& gru_pc(jpi,jpj) , &
313	& grv_pc(jpi,jpj) )
314
315	ALLOCATE( fbcoef(jpi,jpj), fbcoef_stscale(jpi,jpj) )
316
317	IF( ln_bias_asm ) ALLOCATE( tbias_asm(jpi,jpj,jpk), &
318	& sbias_asm(jpi,jpj,jpk), &
319	tbias_asm_out(jpi,jpj,jpk), &
320	sbias_asm_out(jpi,jpj,jpk), &
321	tbias_p_out(jpi,jpj,jpk), &
322	sbias_p_out(jpi,jpj,jpk) )
323
324	IF( ln_bias_rlx ) ALLOCATE( tbias_rlx(jpi,jpj,jpk), &
325	& sbias_rlx(jpi,jpj,jpk), &
326	tbias_rlx_out(jpi,jpj,jpk), &
327	sbias_rlx_out(jpi,jpj,jpk) )
328
329	IF( ln_incpc ) ALLOCATE( tbias_asm_stscale(jpi,jpj,jpk), &
330	& sbias_asm_stscale(jpi,jpj,jpk))
331
332	lalloc = .TRUE.
333
334	ENDIF
335
336	IF( ln_bias_ofl ) THEN ! set sf_tbias_ofl and sf_sbias_ofl strctrs
337	!
338	! tbias
339	!
340	ALLOCATE( sf_tbias_ofl(1), STAT=ierror )
341	IF( ierror > 0 ) THEN
342	CALL ctl_stop( 'bias_init: unable to allocate sf_tbias_ofl structure' ) ; RETURN
343	ENDIF
344	ALLOCATE( sf_tbias_ofl(1)%fnow(jpi,jpj,jpk) )
345	ALLOCATE( sf_tbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
346
347	! fill structure with values and control print
348	CALL fld_fill( sf_tbias_ofl, (/ sn_tbias_ofl /), cn_dir, 'bias_init', 'Offline T bias term ', 'nam_tbias_ofl' )
349	!
350	! salinity bias
351	!
352	ALLOCATE( sf_sbias_ofl(1), STAT=ierror )
353	IF( ierror > 0 ) THEN
354	CALL ctl_stop( 'bias_init: unable to allocate sf_sbias_ofl structure' ) ; RETURN
355	ENDIF
356	ALLOCATE( sf_sbias_ofl(1)%fnow(jpi,jpj,jpk) )
357	ALLOCATE( sf_sbias_ofl(1)%fdta(jpi,jpj,jpk,2) )
358
359	! fill structure with values and control print
360	CALL fld_fill( sf_sbias_ofl, (/ sn_sbias_ofl /), cn_dir, 'bias_init', 'Offline S bias term ', 'nam_sbias_ofl' )
361
362	ENDIF
363
364	! Read total bias
365	IF ( ln_bias ) THEN
366	tbias(:,:,:) = 0.0_wp
367	sbias(:,:,:) = 0.0_wp
368	tbias_p(:,:,:) = 0.0_wp
369	sbias_p(:,:,:) = 0.0_wp
370	gru_pc(:,:) = 0.0_wp
371	grv_pc(:,:) = 0.0_wp
372	IF ( ln_bias_rlx ) THEN
373	tbias_rlx(:,:,:) = 0.0_wp
374	sbias_rlx(:,:,:) = 0.0_wp
375	ENDIF
376	IF ( ln_bias_asm ) THEN !now rlx and asm bias in same file
377	tbias_asm(:,:,:) = 0.0_wp
378	sbias_asm(:,:,:) = 0.0_wp
379	ENDIF
380
381	IF ( ln_incpc ) THEN !incr pressure correction
382	tbias_asm_stscale(:,:,:) = 0.0_wp
383	sbias_asm_stscale(:,:,:) = 0.0_wp
384	ENDIF
385
386
387	numbias_tot = 0
388	! Get bias from file and prevent fail if the file does not exist
389	IF(lwp) WRITE(numout,*) 'Opening ',TRIM( cn_bias_tot )
390	CALL iom_open( cn_bias_tot, numbias_tot, ldstop=.FALSE. )
391
392	IF ( numbias_tot > 0 ) THEN
393	! Could check validity time of bias fields here...
394	! Get the T and S bias data
395	IF(lwp) WRITE(numout,*) 'Reading bias fields from tot...'
396
397	!Search for bias from relaxation term if needed. Use same file
398	IF ( ln_bias_rlx ) THEN
399	IF(lwp) WRITE(numout,*) 'Reading bias fields for bias rlx from file ',cn_bias_tot
400	IF( iom_varid( numbias_tot, 'tbias_rlx' ) > 0 ) THEN
401	! Get the T and S bias data
402	CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_rlx', tbias_rlx )
403	CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_rlx', sbias_rlx )
404	ELSE
405	CALL ctl_stop( 'Bias relaxation variables not found in ',cn_bias_tot )
406	ENDIF
407	ENDIF
408
409
410	!Search for bias from assim term if needed. Use same file
411	IF ( ln_bias_asm .and. .not.ln_bsyncro ) THEN
412	IF(lwp) WRITE(numout,*) 'Reading a-syncro bias fields for bias asm from file ',cn_bias_tot
413	IF( iom_varid( numbias_tot, 'tbias_asm' ) > 0 ) THEN
414	! Get the T and S bias data
415	CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm', tbias_asm )
416	CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm', sbias_asm )
417	ELSE
418	CALL ctl_stop( 'Bias assim variables not found in ',cn_bias_tot )
419	ENDIF
420
421
422	IF ( ln_incpc ) THEN
423	IF(lwp) WRITE(numout,*) 'Reading short time scale bias correction fields for bias asm from file ',cn_bias_tot
424	IF( iom_varid( numbias_tot, 'tbias_asm_stscale' ) > 0 ) THEN
425	! Get the T and S bias data
426	CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm_stscale', tbias_asm_stscale )
427	CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm_stscale', sbias_asm_stscale )
428	ELSE
429	CALL ctl_stop( 'Short time scale bias assim variables not found in ',cn_bias_tot )
430	ENDIF
431	ENDIF
432
433	ENDIF
434
435
436
437	! Close the file
438	CALL iom_close(numbias_tot)
439
440	ELSE
441	IF(lwp) WRITE(numout,*) 'No bias file found so T and S bias fields are set to zero'
442	ENDIF
443
444	ENDIF
445
446	! for the time being, the bias_asm is read in the same file as
447	! bias_rlx
448	! Implications: bias_asm is estimated/evolved in time in the second outer
449	! loop only, when the assimilation increments are ready.
450	! bias_asm is kept cte during the first outer loop.
451	! => Assyncronous bias correction.
452	! Alternative: Syncronous bias correction:
453	! bias_asm estimated/evolved in the first outer loop
454	! with the asm incrments of the previous cycle.
455	! bias_asm kept cte during the second outer loop.
456	! Implication: bias_asm should be estimated really in the
457	! inner loop.
458	IF ( ln_bsyncro ) THEN
459	! Read bias from assimilation from a separate file
460	IF ( ln_bias_asm ) THEN
461	tbias_asm(:,:,:) = 0.0_wp
462	sbias_asm(:,:,:) = 0.0_wp
463	numbias_asm = 0
464	! Get bias from file and prevent fail if the file does not exist
465	IF(lwp) WRITE(numout,*) 'Opening file for syncro assim bias ',TRIM( cn_bias_asm )
466	CALL iom_open( cn_bias_asm, numbias_asm, ldstop=.FALSE. )
467
468	IF ( numbias_asm > 0 ) THEN
469	! Could check validity time of bias fields here...
470
471	! Get the T and S bias data
472	IF(lwp) WRITE(numout,*) 'Reading syncro bias fields from asm from file ',cn_bias_asm
473	CALL iom_get( numbias_asm, jpdom_autoglo, 'tbias_asm', tbias_asm )
474	CALL iom_get( numbias_asm, jpdom_autoglo, 'sbias_asm', sbias_asm )
475
476	! this is only applicable if tbias_asm were to be calculated in the inner loop
477	tbias_asm(:,:,:) = tbias_asm(:,:,:) * rdt / bias_time_unit_asm
478	sbias_asm(:,:,:) = sbias_asm(:,:,:) * rdt / bias_time_unit_asm
479
480	! Close the file
481	CALL iom_close(numbias_asm)
482
483	ELSE
484	IF(lwp) WRITE(numout,*) 'No bias file found from asm so T and S bias fields are set to zero'
485	ENDIF
486
487	ENDIF
488
489	ENDIF
490
491	!latitudinal dependence of partition coeficients. Adhoc
492	IF ( nn_lat_ramp == 1 ) THEN
493	! Use the inertial ramp.
494	lenscl_bias = ( rn_maxlat_bias - rn_minlat_bias )*2._wp
495	WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )
496	fbcoef(:,:) = 0._wp
497	ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias )
498	fbcoef(:,:) = 1._wp
499	ELSEWHERE
500	fbcoef(:,:) = 1._wp - exp( -( abs( gphit(:,:) ) - rn_minlat_bias ) &
501	* ( abs( gphit(:,:) ) - rn_minlat_bias ) / lenscl_bias )
502	ENDWHERE
503	ELSEIF ( nn_lat_ramp == 2 ) THEN
504	! Use a linear ramp consist with the geostrophic velocity balance ramp in NEMOVAR
505
506	WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias )
507	fbcoef(:,:) = 0._wp
508	ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias )
509	fbcoef(:,:) = 1._wp
510	ELSEWHERE
511	fbcoef(:,:) = 1._wp - ((rn_maxlat_bias - abs( gphit(:,:)))/(rn_maxlat_bias - rn_minlat_bias))
512	ENDWHERE
513	ELSE
514	fbcoef(:,:) = 0.0_wp
515	fctamp = 0.0_wp
516	fbcoef_stscale(:,:) = 0.0_wp
517	ENDIF
518
519
520	IF ( ln_incpc) THEN
521	! not sure if this should be a special case of nn_lat_ramp == 2
522	minlat_bias = 3.0_wp
523	maxlat_bias = 8.0_wp
524	WHERE ( abs( gphit(:,:) ) <= minlat_bias )
525	fbcoef_stscale(:,:)=0._wp
526	ELSEWHERE ( abs( gphit(:,:) ) >= maxlat_bias )
527	fbcoef_stscale(:,:)=1._wp
528	ELSEWHERE
529	fbcoef_stscale(:,:)=1._wp - ((maxlat_bias - abs( gphit(:,:)))/(maxlat_bias-minlat_bias))
530	ENDWHERE
531	ENDIF
532
533
534	END SUBROUTINE bias_init
535
536	SUBROUTINE tra_bias ( kt )
537	!!----------------------------------------------------------------------
538	!! * ROUTINE tra_bias *
539	!!
540	!! ** Purpose : Update bias field every time step
541	!!
542	!! ** Method : add contributions to bias from 3 terms
543	!!
544	!! ** Action : Bias from assimilation (read in bias_init)
545	!! Bias from relaxation term is estimated according to
546	!! the prescribed time evolution of the bias
547	!! Bias from ofl term is read from external file
548	!! The difference contributions are added and the partition
549	!! into direct bias in T/S and pressure perfomed.
550	!!
551	!! History : 09-03 (M. Balmaseda)
552	!!----------------------------------------------------------------------
553	!! called every timestep after dta_sst if ln_bias true.
554
555	IMPLICIT NONE
556
557	!! * Arguments
558	INTEGER, INTENT(IN) :: kt ! ocean time-step index
559	!! * Local variables
560	INTEGER :: ji,jj,jk, it ! local loop index
561	REAL(wp) :: tsclf ! time saling factor
562	REAL(wp) :: fb_t_asm_max, fb_t_rlx_max, fb_t_ofl_max
563	REAL(wp) :: ztfrac, ztsday
564	REAL(wp) :: zfrac, zfrac1 ! temporal weights for inst pcbias (names could be changed)
565	REAL(wp) :: ztscale ! reduce the inst pcbias by this amount per 24 hours
566	REAL(wp) :: zwgt ! Weight for the inst pcorr term
567	REAL(wp) :: zdecay ! used in inst pcorr
568	REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2
569
570	IF ( .NOT. ln_bias ) RETURN
571	fb_t_rlx_max = MIN(fb_t_rlx*fctamp,1.0_wp)
572	fb_t_asm_max = MIN(fb_t_asm*fctamp,1.0_wp)
573	fb_t_ofl_max = MIN(fb_t_ofl*fctamp,1.0_wp)
574
575	tbias(:,:,:) = 0.0_wp
576	sbias(:,:,:) = 0.0_wp
577	tbias_p(:,:,:) = 0.0_wp
578	sbias_p(:,:,:) = 0.0_wp
579
580	ztscale = 0.1_wp
581	zwgt = 1.0_wp
582
583	IF ( ln_bias_asm ) THEN
584
585	tsclf = 1
586	IF ( .NOT.ln_bsyncro ) tsclf = rdt / bias_time_unit_asm
587	zcof1(:,:) = tsclf * ( ( 1.0_wp - fbcoef(:,:) ) * fb_t_asm + &
588	& fbcoef(:,:) * fb_t_asm_max )
589	zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_asm
590
591	IF ( ln_itdecay ) THEN !decay the pressure correction at each time step
592
593	ztsday = rday / real(rdt)
594
595	zdecay = (1-ztscale)**(1/real(ztsday)) ! used in ipc
596	zfrac1 = max(0.0_wp, zdecay**real(kt)) ! used in ipc
597
598
599	IF( ln_asmiau .and. ln_trainc ) THEN ! nudge in increments and decay historical contribution
600
601
602	IF ( kt <= nitiaufin ) THEN ! During IAU calculate the fraction of increments to apply at each time step
603
604	ztfrac = real(kt) / real(nitiaufin) ! nudging factor during the IAU
605
606
607	IF (lwp) THEN
608	WRITE(numout,*) 'tra_bias : bias weights'
609	WRITE(numout,*) '~~~~~~~~~~~~'
610	WRITE(numout,* ) "proportion of increment applied in pcbias ",ztfrac
611	WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
612	ENDIF
613
614	DO jk = 1, jpkm1
615	tbias(:,:,jk) = tbias(:,:,jk) + &
616	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
617	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
618	sbias(:,:,jk) = sbias(:,:,jk) + &
619	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
620	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:)
621	tbias_p(:,:,jk) = tbias_p(:,:,jk) + &
622	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
623	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
624	sbias_p(:,:,jk) = sbias_p(:,:,jk) + &
625	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
626	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
627	ENDDO
628
629
630	IF (ln_incpc) THEN
631
632	DO jk = 1, jpkm1
633	tbias_p(:,:,jk) = tbias_p(:,:,jk) + ( t_bkginc(:,:,jk) * zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) ) &
634	& - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
635	sbias_p(:,:,jk) = sbias_p(:,:,jk) + ( s_bkginc(:,:,jk) * zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) ) &
636	& - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
637	ENDDO
638
639	ENDIF
640
641	IF ( .not.ln_bsyncro ) THEN
642	IF ( kt == nn_bias_itwrt ) THEN
643	DO jk = 1, jpkm1
644	tbias_asm_out(:,:,jk) = t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
645	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
646	sbias_asm_out(:,:,jk) = t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
647	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
648	END DO
649	ENDIF
650	ENDIF
651
652	! update the historical component with the increments at the end of the IAU
653	IF ( kt == nitiaufin ) THEN
654	DO jk = 1, jpkm1
655	tbias_asm(:,:,jk) = t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
656	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
657	sbias_asm(:,:,jk) = t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
658	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
659	END DO
660	ENDIF
661
662	ELSE ! decay pressure correction from combined historical component and increments after IAU
663
664	ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin)) ! decay from end of IAU
665
666
667	DO jk = 1, jpkm1
668	tbias(:,:,jk) = tbias(:,:,jk) + &
669	& ( ztfrac * tbias_asm(:,:,jk) ) * zcof1(:,:)
670	sbias(:,:,jk) = sbias(:,:,jk) + &
671	& ( ztfrac * sbias_asm(:,:,jk) ) * zcof1(:,:)
672	tbias_p(:,:,jk) = tbias_p(:,:,jk) + &
673	& ( ztfrac * tbias_asm(:,:,jk) ) * zcof2(:,:)
674	sbias_p(:,:,jk) = sbias_p(:,:,jk) + &
675	& ( ztfrac * sbias_asm(:,:,jk) ) * zcof2(:,:)
676
677	ENDDO
678
679	IF (ln_incpc) THEN
680
681	zfrac = max(0.0_wp, zdecay**real(kt-nitiaufin) )
682
683	DO jk = 1, jpkm1
684	tbias_p(:,:,jk) = tbias_p(:,:,jk) + ( t_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) ) &
685	& - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
686	sbias_p(:,:,jk) = sbias_p(:,:,jk) + ( s_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) ) &
687	& - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
688	ENDDO
689
690	ENDIF
691
692	IF (lwp) THEN
693	WRITE(numout,*) 'tra_bias : bias weights'
694	WRITE(numout,*) '~~~~~~~~~~~~'
695	WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac
696	ENDIF
697
698	IF ( ln_trainc .and. .not.ln_bsyncro ) THEN
699	IF ( kt == nn_bias_itwrt ) THEN
700	DO jk = 1, jpkm1
701	tbias_asm_out(:,:,jk) = ztfrac * tbias_asm(:,:,jk)
702	sbias_asm_out(:,:,jk) = ztfrac * sbias_asm(:,:,jk)
703	END DO
704	ENDIF
705	ENDIF
706
707	ENDIF
708
709
710	ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts)
711
712	DO jk = 1, jpkm1
713	tbias(:,:,jk) = tbias(:,:,jk) + &
714	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) ) * zcof1(:,:)
715	sbias(:,:,jk) = sbias(:,:,jk) + &
716	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) ) * zcof1(:,:)
717	tbias_p(:,:,jk) = tbias_p(:,:,jk) + &
718	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) ) * zcof2(:,:)
719	sbias_p(:,:,jk) = sbias_p(:,:,jk) + &
720	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) ) * zcof2(:,:)
721	ENDDO
722
723	IF (lwp) THEN
724	WRITE(numout,*) 'tra_bias : bias weights'
725	WRITE(numout,*) '~~~~~~~~~~~~'
726	WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
727	ENDIF
728
729	ENDIF
730
731	ELSE ! maintain a constant pressure correction
732
733	DO jk = 1, jpkm1
734	tbias(:,:,jk) = tbias(:,:,jk) + tbias_asm(:,:,jk) * zcof1(:,:)
735	sbias(:,:,jk) = sbias(:,:,jk) + sbias_asm(:,:,jk) * zcof1(:,:)
736	tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_asm(:,:,jk) * zcof2(:,:)
737	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_asm(:,:,jk) * zcof2(:,:)
738	END DO
739
740	IF( ln_asmiau .and. ln_trainc .and. .not.ln_bsyncro ) THEN
741	! if last outer loop (ln_asmiau=true and ln_trainc=true). t/sbias_asm
742	! is updated, only once (end of run) taking into account units.
743	IF ( kt == nn_bias_itwrt ) THEN
744	IF(lwp) WRITE(numout,*)' estimating asm bias at timestep: ',kt
745	DO jk = 1, jpkm1
746	tbias_asm_out(:,:,jk) = t_asm_mem * tbias_asm(:,:,jk) + &
747	& t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
748	sbias_asm_out(:,:,jk) = t_asm_mem * sbias_asm(:,:,jk) + &
749	& t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
750	END DO
751	ENDIF
752	ENDIF
753
754	ENDIF
755
756	ENDIF
757
758
759	#if defined key_tradmp
760	! Time evolution of bias from relaxation
761	IF ( ln_bias_rlx ) THEN
762	tbias_rlx(:,:,:) = t_rlx_mem * tbias_rlx(:,:,:) + &
763	& t_rlx_upd * ttrdmp(:,:,:) * rdt
764	sbias_rlx(:,:,:) = t_rlx_mem * sbias_rlx(:,:,:) + &
765	& t_rlx_upd * strdmp(:,:,:) * rdt
766	zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_rlx + &
767	& fbcoef(:,:) * fb_t_rlx_max
768	zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_rlx
769	DO jk = 1, jpkm1
770	tbias(:,:,jk) = tbias(:,:,jk) + tbias_rlx(:,:,jk) * zcof1(:,:)
771	sbias(:,:,jk) = sbias(:,:,jk) + sbias_rlx(:,:,jk) * zcof1(:,:)
772	tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_rlx(:,:,jk) * zcof2(:,:)
773	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_rlx(:,:,jk) * zcof2(:,:)
774	ENDDO
775	IF ( kt == nn_bias_itwrt ) THEN
776	tbias_rlx_out(:,:,:) = tbias_rlx(:,:,:)
777	sbias_rlx_out(:,:,:) = sbias_rlx(:,:,:)
778	ENDIF
779	ENDIF
780	#endif
781	! ofline bias
782	IF ( kt == nit000 ) THEN
783	IF(lwp) WRITE(numout,*) ' tra_bias: ln_bias_ofl = ',ln_bias_ofl
784	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~'
785	ENDIF
786	IF ( ln_bias_ofl ) THEN
787	IF(lwp) WRITE(numout,*) 'reading offline bias'
788	CALL fld_read( kt, 1, sf_tbias_ofl )
789	CALL fld_read( kt, 1, sf_sbias_ofl )
790
791	zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_ofl + &
792	& fbcoef(:,:) * fb_t_ofl_max
793	zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_ofl
794	DO jk = 1, jpkm1
795	tbias(:,:,jk) = tbias(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
796	sbias(:,:,jk) = sbias(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
797	tbias_p(:,:,jk) = tbias_p(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
798	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
799	ENDDO
800	ENDIF
801
802
803	!apply bias on tracers if needed
804	IF ( ln_bias_ts_app ) THEN
805
806	! Add the bias directely to T/s
807	tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + tmask(:,:,:) * tbias(:,:,:) / rdt
808	tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + tmask(:,:,:) * sbias(:,:,:) / rdt
809
810	! lateral boundary conditions (is this needed?)
811	CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1.0_wp )
812	CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1.0_wp )
813
814	ENDIF
815
816	IF ( kt == nn_bias_itwrt ) THEN
817	tbias_p_out(:,:,:) = tbias_p(:,:,:)
818	sbias_p_out(:,:,:) = sbias_p(:,:,:)
819	ENDIF
820
821	END SUBROUTINE tra_bias
822
823	SUBROUTINE dyn_bias( kt )
824	!!----------------------------------------------------------------------
825	!! * ROUTINE dyn_bias *
826	!!
827	!! ** Purpose : Computes rhd_pc, gru/v_pc bias corrected
828	!! for hydrostatic pressure gradient
829	!! depending on time step (semi-implicit/centered)
830	!! If partial steps computes bottom pressure gradient.
831	!! These correction terms will affect only the dynamical
832	!! component (pressure gradient calculation), but not
833	!! the isopycnal calculation for the lateral mixing.
834	!!
835	!! ** Method : At this stage of the computation, ta and sa are the
836	!! after temperature and salinity. If semi-implicit, these
837	!! are used to compute rho and bottom pressure gradient.
838	!! If centered, tb,sb are used instead.
839	!! If bias key is activated, the temperature,salinity are
840	!! bias corrected in the calculation of the density fields
841	!! used in the pressure gradient calculation.
842	!!
843	!!
844	!! ** Action : - rhd_pc ready. rhop will be overwriten later
845	!! - if ln_zps, bottom density gradients gru/v_pc ready.
846	!!----------------------------------------------------------------------
847	!!
848	!! * Arguments
849	INTEGER, INTENT(IN) :: kt ! ocean time-step index
850	!! * Local variables
851	REAL(wp) :: tsw(jpi,jpj,jpk,jpts)
852	!!
853	!!----------------------------------------------------------------------
854	!
855	! gtu,gsu,gtv,gsv rhop will be overwritten later in step.
856	!
857	IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg
858	tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:)
859	tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:)
860	ELSE
861	tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:)
862	tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:)
863	ENDIF
864
865	CALL eos( tsw, rhd_pc, rhop )
866
867	! is this needed?
868	CALL lbc_lnk( rhd_pc, 'T', 1.0_wp )
869
870	! Partial steps: now horizontal gradient of t,s,rd
871	! at the bottom ocean level
872	IF( ln_zps ) THEN
873	CALL zps_hde( kt, jpts, tsw, gtsu, gtsv, &
874	& rhd_pc, gru_pc , grv_pc )
875	ENDIF
876
877	END SUBROUTINE dyn_bias
878
879	SUBROUTINE bias_opn( kt )
880	!!---------------------------------------------------------------------
881	!! * ROUTINE bias_opn *
882	!!
883	!! ** Purpose : open bias restart file following the same logic as the
884	!! standard restarts.
885	!!----------------------------------------------------------------------
886	!! * Arguments
887	INTEGER, INTENT(IN) :: kt ! ocean time-step
888	!! * Local variables
889	CHARACTER(LEN=20) :: clbkt ! ocean time-step deine as a character
890	CHARACTER(LEN=50) :: clbias_tot ! total bias restart file name
891	!!----------------------------------------------------------------------
892	!
893	IF( lrst_oce .AND. .NOT.lrst_bias ) THEN ! create bias file
894	IF( nitend > 999999999 ) THEN ; WRITE(clbkt, * ) nitend
895	ELSE ; WRITE(clbkt, '(i8.8)') nitend
896	ENDIF
897	clbias_tot = TRIM(cexper)//"_"//TRIM(ADJUSTL(clbkt))//"_"//TRIM(cn_bias_tot)
898	IF(lwp) THEN
899	WRITE(numout,*)
900	SELECT CASE ( jprstlib )
901	CASE ( jprstdimg ) ; WRITE(numout,*) ' open tot bias restart binary file: '//clbias_tot
902	CASE DEFAULT ; WRITE(numout,*) ' open tot bias restart NetCDF file: '//clbias_tot
903	END SELECT
904	ENDIF
905	CALL iom_open( clbias_tot, numbias_tot , ldwrt = .TRUE., kiolib = jprstlib )
906	lrst_bias=.TRUE.
907
908	ENDIF
909	!
910	END SUBROUTINE bias_opn
911
912	SUBROUTINE bias_wrt( kt )
913	!!---------------------------------------------------------------------
914	!! * ROUTINE bias_wrt *
915	!!
916	!! ** Purpose : Write bias restart fields in the format corresponding to jprstlib
917	!!
918	!! ** Method : Write in numbias_tot when kt == nitend in output
919	!! file, save fields which are necessary for restart.
920	!!
921	!! Changed the timestep for writing to nitend rather than nitrst as this causes a
922	!! problem if the nstock list is used to determine the restart output times and
923	!! means that the bias is not output at all. M. Martin. 08/2011.
924	!! Need to check with Magdalena that this is ok for ECMWF.
925	!!----------------------------------------------------------------------
926	!! * Arguments
927	INTEGER, INTENT(IN) :: kt ! ocean time-step
928	!!----------------------------------------------------------------------
929	! ! the begining of the run [s]
930
931	IF ( ln_bias_rlx ) THEN
932	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_rlx' , tbias_rlx_out )
933	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_rlx' , sbias_rlx_out )
934	ENDIF
935
936	IF ( ln_bias_asm ) THEN
937	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm' , tbias_asm_out )
938	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm' , sbias_asm_out )
939	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_p' , tbias_p_out )
940	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p' , sbias_p_out )
941	ENDIF
942
943	IF( kt == nitend ) THEN
944	CALL iom_close( numbias_tot ) ! close the restart file (only at last time step)
945	lrst_bias = .FALSE.
946	ENDIF
947	!
948	END SUBROUTINE bias_wrt
949
950	END MODULE bias

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

Download in other formats: