source: branches/UKMO/dev_r5518_pcbias_ipc/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90 @ 7201

Last change on this file since 7201 was 7201, checked in by isabella, 5 years ago

adding output in forecast

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