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

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

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

Download in other formats: