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

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

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

Download in other formats: