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

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

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

Download in other formats: