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

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

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

Download in other formats: