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

Last change on this file since 7081 was 7081, checked in by isabella, 7 years ago
Initial changes for ipc
File size: 45.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
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) = tbias_i(:,:,jk) + & ! not sure if it should sum to previous
639	& ( t_bkginc(:,:,jk) * zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) ) &
640	& - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
641	sbias_i(:,:,jk) = sbias_i(:,:,jk) + &
642	& ( s_bkginc(:,:,jk) * zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) ) &
643	& - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
644
645	ENDDO
646
647	ENDIF
648
649	IF ( .not.ln_bsyncro ) THEN
650	IF ( kt == nn_bias_itwrt ) THEN
651	DO jk = 1, jpkm1
652	tbias_asm_out(:,:,jk) = t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
653	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
654	sbias_asm_out(:,:,jk) = t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
655	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
656	END DO
657	ENDIF
658	ENDIF
659
660	! update the historical component with the increments at the end of the IAU
661	IF ( kt == nitiaufin ) THEN
662	DO jk = 1, jpkm1
663	tbias_asm(:,:,jk) = t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) + &
664	& ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
665	sbias_asm(:,:,jk) = t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) + &
666	& ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
667	END DO
668	ENDIF
669
670	ELSE ! decay pressure correction from combined historical component and increments after IAU
671
672	ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin)) ! decay from end of IAU
673
674
675	DO jk = 1, jpkm1
676	tbias(:,:,jk) = tbias(:,:,jk) + &
677	& ( ztfrac * tbias_asm(:,:,jk) ) * zcof1(:,:)
678	sbias(:,:,jk) = sbias(:,:,jk) + &
679	& ( ztfrac * sbias_asm(:,:,jk) ) * zcof1(:,:)
680	tbias_p(:,:,jk) = tbias_p(:,:,jk) + &
681	& ( ztfrac * tbias_asm(:,:,jk) ) * zcof2(:,:)
682	sbias_p(:,:,jk) = sbias_p(:,:,jk) + &
683	& ( ztfrac * sbias_asm(:,:,jk) ) * zcof2(:,:)
684
685	ENDDO
686
687	IF (ln_incpc) THEN
688
689	zfrac = max(0.0_wp, zdecay**real(kt-nitiaufin) )
690
691	DO jk = 1, jpkm1
692	tbias_i(:,:,jk) = tbias_i(:,:,jk) + & ! not sure if it should sum to previous
693	& ( t_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) ) &
694	& - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
695	sbias_i(:,:,jk) = sbias_i(:,:,jk) + & ! not sure if it should sum to previous
696	& ( s_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) ) &
697	& - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
698	ENDDO
699
700	ENDIF
701
702	IF (lwp) THEN
703	WRITE(numout,*) 'tra_bias : bias weights'
704	WRITE(numout,*) '~~~~~~~~~~~~'
705	WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac
706	ENDIF
707
708	IF ( ln_trainc .and. .not.ln_bsyncro ) THEN
709	IF ( kt == nn_bias_itwrt ) THEN
710	DO jk = 1, jpkm1
711	tbias_asm_out(:,:,jk) = ztfrac * tbias_asm(:,:,jk)
712	sbias_asm_out(:,:,jk) = ztfrac * sbias_asm(:,:,jk)
713	END DO
714	ENDIF
715	ENDIF
716
717	ENDIF
718
719
720	ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts)
721
722	DO jk = 1, jpkm1
723	tbias(:,:,jk) = tbias(:,:,jk) + &
724	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) ) * zcof1(:,:)
725	sbias(:,:,jk) = sbias(:,:,jk) + &
726	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) ) * zcof1(:,:)
727	tbias_p(:,:,jk) = tbias_p(:,:,jk) + &
728	& ( t_asm_mem*(real(kt)/ztsday) tbias_asm(:,:,jk) ) * zcof2(:,:)
729	sbias_p(:,:,jk) = sbias_p(:,:,jk) + &
730	& ( t_asm_mem*(real(kt)/ztsday) sbias_asm(:,:,jk) ) * zcof2(:,:)
731	ENDDO
732
733	IF (lwp) THEN
734	WRITE(numout,*) 'tra_bias : bias weights'
735	WRITE(numout,*) '~~~~~~~~~~~~'
736	WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
737	ENDIF
738
739	ENDIF
740
741	ELSE ! maintain a constant pressure correction
742
743	DO jk = 1, jpkm1
744	tbias(:,:,jk) = tbias(:,:,jk) + tbias_asm(:,:,jk) * zcof1(:,:)
745	sbias(:,:,jk) = sbias(:,:,jk) + sbias_asm(:,:,jk) * zcof1(:,:)
746	tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_asm(:,:,jk) * zcof2(:,:)
747	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_asm(:,:,jk) * zcof2(:,:)
748	END DO
749
750	IF( ln_asmiau .and. ln_trainc .and. .not.ln_bsyncro ) THEN
751	! if last outer loop (ln_asmiau=true and ln_trainc=true). t/sbias_asm
752	! is updated, only once (end of run) taking into account units.
753	IF ( kt == nn_bias_itwrt ) THEN
754	IF(lwp) WRITE(numout,*)' estimating asm bias at timestep: ',kt
755	DO jk = 1, jpkm1
756	tbias_asm_out(:,:,jk) = t_asm_mem * tbias_asm(:,:,jk) + &
757	& t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
758	sbias_asm_out(:,:,jk) = t_asm_mem * sbias_asm(:,:,jk) + &
759	& t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
760	END DO
761	ENDIF
762	ENDIF
763
764	ENDIF
765
766	ENDIF
767
768
769	#if defined key_tradmp
770	! Time evolution of bias from relaxation
771	IF ( ln_bias_rlx ) THEN
772	tbias_rlx(:,:,:) = t_rlx_mem * tbias_rlx(:,:,:) + &
773	& t_rlx_upd * ttrdmp(:,:,:) * rdt
774	sbias_rlx(:,:,:) = t_rlx_mem * sbias_rlx(:,:,:) + &
775	& t_rlx_upd * strdmp(:,:,:) * rdt
776	zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_rlx + &
777	& fbcoef(:,:) * fb_t_rlx_max
778	zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_rlx
779	DO jk = 1, jpkm1
780	tbias(:,:,jk) = tbias(:,:,jk) + tbias_rlx(:,:,jk) * zcof1(:,:)
781	sbias(:,:,jk) = sbias(:,:,jk) + sbias_rlx(:,:,jk) * zcof1(:,:)
782	tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_rlx(:,:,jk) * zcof2(:,:)
783	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_rlx(:,:,jk) * zcof2(:,:)
784	ENDDO
785	IF ( kt == nn_bias_itwrt ) THEN
786	tbias_rlx_out(:,:,:) = tbias_rlx(:,:,:)
787	sbias_rlx_out(:,:,:) = sbias_rlx(:,:,:)
788	ENDIF
789	ENDIF
790	#endif
791	! ofline bias
792	IF ( kt == nit000 ) THEN
793	IF(lwp) WRITE(numout,*) ' tra_bias: ln_bias_ofl = ',ln_bias_ofl
794	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~'
795	ENDIF
796	IF ( ln_bias_ofl ) THEN
797	IF(lwp) WRITE(numout,*) 'reading offline bias'
798	CALL fld_read( kt, 1, sf_tbias_ofl )
799	CALL fld_read( kt, 1, sf_sbias_ofl )
800
801	zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_ofl + &
802	& fbcoef(:,:) * fb_t_ofl_max
803	zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_ofl
804	DO jk = 1, jpkm1
805	tbias(:,:,jk) = tbias(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
806	sbias(:,:,jk) = sbias(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
807	tbias_p(:,:,jk) = tbias_p(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
808	sbias_p(:,:,jk) = sbias_p(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
809	ENDDO
810	ENDIF
811
812
813	!apply bias on tracers if needed
814	IF ( ln_bias_ts_app ) THEN
815
816	! Add the bias directely to T/s
817	tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + tmask(:,:,:) * tbias(:,:,:) / rdt
818	tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + tmask(:,:,:) * sbias(:,:,:) / rdt
819
820	! lateral boundary conditions (is this needed?)
821	CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1.0_wp )
822	CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1.0_wp )
823
824	ENDIF
825
826	IF ( kt == nn_bias_itwrt ) THEN
827	tbias_p_out(:,:,:) = tbias_p(:,:,:)
828	sbias_p_out(:,:,:) = sbias_p(:,:,:)
829	ENDIF
830
831	END SUBROUTINE tra_bias
832
833	SUBROUTINE dyn_bias( kt )
834	!!----------------------------------------------------------------------
835	!! * ROUTINE dyn_bias *
836	!!
837	!! ** Purpose : Computes rhd_pc, gru/v_pc bias corrected
838	!! for hydrostatic pressure gradient
839	!! depending on time step (semi-implicit/centered)
840	!! If partial steps computes bottom pressure gradient.
841	!! These correction terms will affect only the dynamical
842	!! component (pressure gradient calculation), but not
843	!! the isopycnal calculation for the lateral mixing.
844	!!
845	!! ** Method : At this stage of the computation, ta and sa are the
846	!! after temperature and salinity. If semi-implicit, these
847	!! are used to compute rho and bottom pressure gradient.
848	!! If centered, tb,sb are used instead.
849	!! If bias key is activated, the temperature,salinity are
850	!! bias corrected in the calculation of the density fields
851	!! used in the pressure gradient calculation.
852	!!
853	!!
854	!! ** Action : - rhd_pc ready. rhop will be overwriten later
855	!! - if ln_zps, bottom density gradients gru/v_pc ready.
856	!!----------------------------------------------------------------------
857	!!
858	!! * Arguments
859	INTEGER, INTENT(IN) :: kt ! ocean time-step index
860	!! * Local variables
861	REAL(wp) :: tsw(jpi,jpj,jpk,jpts)
862	!!
863	!!----------------------------------------------------------------------
864	!
865	! gtu,gsu,gtv,gsv rhop will be overwritten later in step.
866	!
867	IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg
868	tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:)
869	tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:)
870
871	IF (ln_incpc) THEN
872	tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) -tbias_i(:,:,:) - tbias_p(:,:,:)
873	tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) -sbias_i(:,:,:) - sbias_p(:,:,:)
874	ENDIF
875
876	ELSE
877	tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:)
878	tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:)
879
880	IF (ln_incpc) THEN
881	tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_i(:,:,:) - tbias_p(:,:,:)
882	tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_i(:,:,:) - sbias_p(:,:,:)
883	ENDIF
884
885	ENDIF
886
887	CALL eos( tsw, rhd_pc, rhop )
888
889	! is this needed?
890	CALL lbc_lnk( rhd_pc, 'T', 1.0_wp )
891
892	! Partial steps: now horizontal gradient of t,s,rd
893	! at the bottom ocean level
894	IF( ln_zps ) THEN
895	CALL zps_hde( kt, jpts, tsw, gtsu, gtsv, &
896	& rhd_pc, gru_pc , grv_pc )
897	ENDIF
898
899	END SUBROUTINE dyn_bias
900
901	SUBROUTINE bias_opn( kt )
902	!!---------------------------------------------------------------------
903	!! * ROUTINE bias_opn *
904	!!
905	!! ** Purpose : open bias restart file following the same logic as the
906	!! standard restarts.
907	!!----------------------------------------------------------------------
908	!! * Arguments
909	INTEGER, INTENT(IN) :: kt ! ocean time-step
910	!! * Local variables
911	CHARACTER(LEN=20) :: clbkt ! ocean time-step deine as a character
912	CHARACTER(LEN=50) :: clbias_tot ! total bias restart file name
913	!!----------------------------------------------------------------------
914	!
915	IF( lrst_oce .AND. .NOT.lrst_bias ) THEN ! create bias file
916	IF( nitend > 999999999 ) THEN ; WRITE(clbkt, * ) nitend
917	ELSE ; WRITE(clbkt, '(i8.8)') nitend
918	ENDIF
919	clbias_tot = TRIM(cexper)//"_"//TRIM(ADJUSTL(clbkt))//"_"//TRIM(cn_bias_tot)
920	IF(lwp) THEN
921	WRITE(numout,*)
922	SELECT CASE ( jprstlib )
923	CASE ( jprstdimg ) ; WRITE(numout,*) ' open tot bias restart binary file: '//clbias_tot
924	CASE DEFAULT ; WRITE(numout,*) ' open tot bias restart NetCDF file: '//clbias_tot
925	END SELECT
926	ENDIF
927	CALL iom_open( clbias_tot, numbias_tot , ldwrt = .TRUE., kiolib = jprstlib )
928	lrst_bias=.TRUE.
929
930	ENDIF
931	!
932	END SUBROUTINE bias_opn
933
934	SUBROUTINE bias_wrt( kt )
935	!!---------------------------------------------------------------------
936	!! * ROUTINE bias_wrt *
937	!!
938	!! ** Purpose : Write bias restart fields in the format corresponding to jprstlib
939	!!
940	!! ** Method : Write in numbias_tot when kt == nitend in output
941	!! file, save fields which are necessary for restart.
942	!!
943	!! Changed the timestep for writing to nitend rather than nitrst as this causes a
944	!! problem if the nstock list is used to determine the restart output times and
945	!! means that the bias is not output at all. M. Martin. 08/2011.
946	!! Need to check with Magdalena that this is ok for ECMWF.
947	!!----------------------------------------------------------------------
948	!! * Arguments
949	INTEGER, INTENT(IN) :: kt ! ocean time-step
950	!!----------------------------------------------------------------------
951	! ! the begining of the run [s]
952
953	IF ( ln_bias_rlx ) THEN
954	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_rlx' , tbias_rlx_out )
955	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_rlx' , sbias_rlx_out )
956	ENDIF
957
958	IF ( ln_bias_asm ) THEN
959	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm' , tbias_asm_out )
960	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm' , sbias_asm_out )
961	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_p' , tbias_p_out )
962	CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p' , sbias_p_out )
963	ENDIF
964
965	IF( kt == nitend ) THEN
966	CALL iom_close( numbias_tot ) ! close the restart file (only at last time step)
967	lrst_bias = .FALSE.
968	ENDIF
969	!
970	END SUBROUTINE bias_wrt
971
972	END MODULE bias

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

Download in other formats: