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

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

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

Download in other formats: