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

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

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

Download in other formats: