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

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

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

Download in other formats: