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, 4 years ago

correcting formulas

File size: 45.6 KB
Line 
1MODULE 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_mem*bias_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   
181CONTAINS
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)
234901   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 )
250902   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
979END MODULE bias
Note: See TracBrowser for help on using the repository browser.