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

Last change on this file since 7172 was 7172, checked in by isabella, 4 years ago

adding updates to variables tbias_asm_stscale, sbias_asm_stscale

File size: 46.0 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
673                     IF (ln_incpc) THEN
674                        DO jk = 1, jpk
675                           tbias_asm_stscale(:,:,jk) = ( t_bkginc(:,:,jk) * zwgt * ztfrac ) + ( tbias_asm_stscale(:,:,jk) * zfrac1 )
676                           sbias_asm_stscale(:,:,jk) = ( s_bkginc(:,:,jk) * zwgt * ztfrac ) + ( sbias_asm_stscale(:,:,jk) * zfrac1 )
677                        ENDDO
678                     ENDIF
679                  ENDIF
680               
681               ELSE ! decay pressure correction from combined historical component and increments after IAU
682
683                  ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin))  ! decay from end of IAU
684                 
685                 
686                  DO jk = 1, jpkm1
687                     tbias(:,:,jk) = tbias(:,:,jk) +                            &
688                     &                ( ztfrac * tbias_asm(:,:,jk) ) * zcof1(:,:)
689                     sbias(:,:,jk) = sbias(:,:,jk) +                            &
690                     &               (  ztfrac * sbias_asm(:,:,jk) ) * zcof1(:,:)
691                     tbias_p(:,:,jk) = tbias_p(:,:,jk) +                        &
692                     &               (  ztfrac * tbias_asm(:,:,jk) ) * zcof2(:,:)
693                     sbias_p(:,:,jk) = sbias_p(:,:,jk) +                        & 
694                     &               ( ztfrac * sbias_asm(:,:,jk) ) * zcof2(:,:)
695
696                  ENDDO
697
698                 IF (ln_incpc) THEN
699
700                   zfrac  = max(0.0_wp, zdecay**real(kt-nitiaufin) )
701                   
702                   DO jk = 1, jpkm1
703                      tbias_i(:,:,jk) = ( t_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) )         &
704                      &                - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
705                      sbias_i(:,:,jk) =  ( s_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) )         &
706                      &                - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
707                   ENDDO
708
709                 ENDIF
710
711                  IF (lwp) THEN
712                     WRITE(numout,*) 'tra_bias : bias weights'
713                     WRITE(numout,*) '~~~~~~~~~~~~'
714                     WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac
715                  ENDIF
716
717                  IF ( ln_trainc .and. .not.ln_bsyncro ) THEN
718                     IF ( kt == nn_bias_itwrt ) THEN
719                        DO jk = 1, jpkm1
720                           tbias_asm_out(:,:,jk) =  ztfrac * tbias_asm(:,:,jk) 
721                           sbias_asm_out(:,:,jk) =  ztfrac * sbias_asm(:,:,jk) 
722                         END DO
723                     ENDIF
724                  ENDIF
725
726               ENDIF
727
728
729            ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts)
730
731               DO jk = 1, jpkm1
732                  tbias(:,:,jk) = tbias(:,:,jk) +                                                         &
733                  &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof1(:,:)
734                  sbias(:,:,jk) = sbias(:,:,jk) +                                                         &
735                  &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof1(:,:)
736                  tbias_p(:,:,jk) = tbias_p(:,:,jk) +                                                     &
737                  &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof2(:,:)
738                  sbias_p(:,:,jk) = sbias_p(:,:,jk) +                                                     & 
739                  &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof2(:,:)
740               ENDDO
741
742               IF (lwp) THEN
743                  WRITE(numout,*) 'tra_bias : bias weights'
744                  WRITE(numout,*) '~~~~~~~~~~~~'
745                  WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
746               ENDIF
747
748            ENDIF
749 
750         ELSE ! maintain a constant pressure correction 
751
752            DO jk = 1, jpkm1
753               tbias(:,:,jk) = tbias(:,:,jk) + tbias_asm(:,:,jk) * zcof1(:,:)
754               sbias(:,:,jk) = sbias(:,:,jk) + sbias_asm(:,:,jk) * zcof1(:,:)
755               tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_asm(:,:,jk) * zcof2(:,:)
756               sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_asm(:,:,jk) * zcof2(:,:)
757            END DO     
758
759            IF( ln_asmiau .and. ln_trainc .and. .not.ln_bsyncro ) THEN   
760            ! if last outer loop (ln_asmiau=true and ln_trainc=true). t/sbias_asm
761            ! is updated, only once (end of run) taking into account units.
762               IF ( kt == nn_bias_itwrt ) THEN
763                 IF(lwp) WRITE(numout,*)' estimating asm bias at timestep: ',kt
764                 DO jk = 1, jpkm1
765                   tbias_asm_out(:,:,jk) = t_asm_mem * tbias_asm(:,:,jk)  +             &
766                   &                      t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk)
767                   sbias_asm_out(:,:,jk) = t_asm_mem * sbias_asm(:,:,jk) +             &
768                   &                      t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)               
769                 END DO
770               ENDIF
771             ENDIF
772 
773         ENDIF
774
775      ENDIF
776
777
778#if   defined key_tradmp 
779      ! Time evolution of bias from relaxation
780      IF ( ln_bias_rlx ) THEN
781         tbias_rlx(:,:,:) = t_rlx_mem * tbias_rlx(:,:,:) + &
782            &               t_rlx_upd * ttrdmp(:,:,:) * rdt
783         sbias_rlx(:,:,:) = t_rlx_mem * sbias_rlx(:,:,:) + &
784            &               t_rlx_upd * strdmp(:,:,:) * rdt
785         zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_rlx +      &
786            &                    fbcoef(:,:)   * fb_t_rlx_max
787         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_rlx
788         DO jk = 1, jpkm1
789            tbias(:,:,jk)   = tbias(:,:,jk) + tbias_rlx(:,:,jk) * zcof1(:,:)
790            sbias(:,:,jk)   = sbias(:,:,jk) + sbias_rlx(:,:,jk) * zcof1(:,:)
791            tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_rlx(:,:,jk) * zcof2(:,:)
792            sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_rlx(:,:,jk) * zcof2(:,:)
793         ENDDO
794         IF ( kt == nn_bias_itwrt ) THEN
795            tbias_rlx_out(:,:,:) = tbias_rlx(:,:,:)
796            sbias_rlx_out(:,:,:) = sbias_rlx(:,:,:)
797         ENDIF
798      ENDIF
799#endif
800      ! ofline bias
801      IF ( kt == nit000 ) THEN
802         IF(lwp) WRITE(numout,*) ' tra_bias: ln_bias_ofl = ',ln_bias_ofl
803         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~'
804      ENDIF
805      IF ( ln_bias_ofl ) THEN
806         IF(lwp) WRITE(numout,*) 'reading offline bias'
807         CALL fld_read( kt, 1, sf_tbias_ofl ) 
808         CALL fld_read( kt, 1, sf_sbias_ofl ) 
809
810         zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_ofl +           &
811            &                    fbcoef(:,:)   * fb_t_ofl_max
812         zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_ofl
813         DO jk = 1, jpkm1
814            tbias(:,:,jk)   = tbias(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
815            sbias(:,:,jk)   = sbias(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:)
816            tbias_p(:,:,jk) = tbias_p(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
817            sbias_p(:,:,jk) = sbias_p(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:)
818         ENDDO
819      ENDIF
820
821
822      !apply bias on tracers if needed     
823      IF ( ln_bias_ts_app ) THEN
824         
825         ! Add the bias directely to T/s
826         tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + tmask(:,:,:) * tbias(:,:,:) / rdt
827         tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + tmask(:,:,:) * sbias(:,:,:) / rdt
828
829         ! lateral boundary conditions (is this needed?)
830         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1.0_wp )
831         CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1.0_wp )   
832
833      ENDIF
834
835      IF ( kt == nn_bias_itwrt ) THEN
836         tbias_p_out(:,:,:)   = tbias_p(:,:,:)
837         sbias_p_out(:,:,:)   = sbias_p(:,:,:)
838         tbias_i_out(:,:,:)   = tbias_i(:,:,:)
839         sbias_i_out(:,:,:)   = sbias_i(:,:,:)
840      ENDIF
841
842   END SUBROUTINE tra_bias
843
844   SUBROUTINE dyn_bias( kt )
845      !!----------------------------------------------------------------------
846      !!                   ***  ROUTINE dyn_bias  ***
847      !!
848      !! ** Purpose :   Computes rhd_pc, gru/v_pc bias corrected
849      !!                for hydrostatic pressure gradient
850      !!                depending on time step (semi-implicit/centered)
851      !!                If partial steps computes bottom pressure gradient.
852      !!                These correction terms will affect only the dynamical
853      !!                component (pressure gradient calculation), but not
854      !!                the isopycnal calculation for the lateral mixing.
855      !!
856      !! ** Method  :   At this stage of the computation, ta and sa are the
857      !!                after temperature and salinity. If semi-implicit, these
858      !!                are used to compute rho and bottom pressure gradient.
859      !!                If centered, tb,sb are used instead.
860      !!                If bias key is activated, the temperature,salinity are
861      !!                bias corrected in the calculation of the density fields
862      !!                used in the pressure gradient calculation.
863      !!
864      !!
865      !! ** Action  : - rhd_pc ready. rhop will be overwriten later
866      !!              - if ln_zps, bottom density gradients gru/v_pc ready.
867      !!----------------------------------------------------------------------
868      !!
869      !! * Arguments
870      INTEGER, INTENT(IN) ::   kt    ! ocean time-step index
871      !! * Local variables
872      REAL(wp) :: tsw(jpi,jpj,jpk,jpts)
873      !!
874      !!----------------------------------------------------------------------
875      !
876      ! gtu,gsu,gtv,gsv rhop will be overwritten later in step.
877      !
878      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg
879         tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:)
880         tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:)
881         IF ( ln_incpc ) THEN
882            tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:)
883            tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:)
884         ENDIF
885      ELSE
886         tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:)
887         tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:)
888         IF ( ln_incpc ) THEN
889            tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:)
890            tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:)
891         ENDIF
892      ENDIF
893
894      CALL eos( tsw, rhd_pc, rhop )
895     
896      ! is this needed?
897      CALL lbc_lnk( rhd_pc, 'T', 1.0_wp )
898
899      ! Partial steps: now horizontal gradient of t,s,rd
900      ! at the bottom ocean level
901      IF( ln_zps    )  THEN
902         CALL zps_hde( kt, jpts, tsw, gtsu, gtsv,  &
903            &          rhd_pc, gru_pc , grv_pc  )
904      ENDIF
905
906   END SUBROUTINE dyn_bias
907   
908   SUBROUTINE bias_opn( kt )
909      !!---------------------------------------------------------------------
910      !!                   ***  ROUTINE bias_opn  ***
911      !!                     
912      !! ** Purpose :  open bias restart file following the same logic as the
913      !!               standard restarts.
914      !!----------------------------------------------------------------------
915      !! * Arguments
916      INTEGER, INTENT(IN) ::   kt     ! ocean time-step
917      !! * Local variables
918      CHARACTER(LEN=20)   ::   clbkt    ! ocean time-step deine as a character
919      CHARACTER(LEN=50)   ::   clbias_tot   ! total bias restart file name
920      !!----------------------------------------------------------------------
921      !
922      IF( lrst_oce .AND. .NOT.lrst_bias ) THEN       ! create bias file
923         IF( nitend > 999999999 ) THEN   ;   WRITE(clbkt, *       ) nitend
924         ELSE                            ;   WRITE(clbkt, '(i8.8)') nitend
925         ENDIF
926         clbias_tot = TRIM(cexper)//"_"//TRIM(ADJUSTL(clbkt))//"_"//TRIM(cn_bias_tot)
927         IF(lwp) THEN
928            WRITE(numout,*)
929            SELECT CASE ( jprstlib )
930            CASE ( jprstdimg )   ;   WRITE(numout,*) '             open tot bias restart binary file: '//clbias_tot
931            CASE DEFAULT         ;   WRITE(numout,*) '             open tot bias restart NetCDF file: '//clbias_tot
932            END SELECT
933         ENDIF
934         CALL iom_open( clbias_tot, numbias_tot , ldwrt = .TRUE., kiolib = jprstlib )
935         lrst_bias=.TRUE.
936
937      ENDIF
938      !
939   END SUBROUTINE bias_opn
940
941   SUBROUTINE bias_wrt( kt )
942      !!---------------------------------------------------------------------
943      !!                   ***  ROUTINE bias_wrt  ***
944      !!                     
945      !! ** Purpose :   Write bias restart fields in the format corresponding to jprstlib
946      !!
947      !! ** Method  :   Write in numbias_tot when kt == nitend in output
948      !!                file, save fields which are necessary for restart.
949      !!
950      !! Changed the timestep for writing to nitend rather than nitrst as this causes a
951      !! problem if the nstock list is used to determine the restart output times and
952      !! means that the bias is not output at all. M. Martin. 08/2011.
953      !! Need to check with Magdalena that this is ok for ECMWF.
954      !!----------------------------------------------------------------------
955      !! * Arguments
956      INTEGER, INTENT(IN) ::   kt   ! ocean time-step
957      !!----------------------------------------------------------------------
958      !                                                                     ! the begining of the run [s]
959
960      IF ( ln_bias_rlx ) THEN
961         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_rlx' , tbias_rlx_out )   
962         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_rlx' , sbias_rlx_out )   
963      ENDIF
964     
965      IF ( ln_bias_asm ) THEN
966         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm' , tbias_asm_out )   
967         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm' , sbias_asm_out )   
968         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_p'   , tbias_p_out )
969         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p'   , sbias_p_out )
970         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_i'   , tbias_i_out )
971         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_i'   , sbias_i_out )         
972      ENDIF
973
974      IF ( ln_incpc ) THEN
975         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm_stscale' , tbias_asm_stscale )   
976         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm_stscale' , sbias_asm_stscale )   
977      ENDIF
978     
979      IF( kt == nitend ) THEN
980         CALL iom_close( numbias_tot )     ! close the restart file (only at last time step)
981         lrst_bias = .FALSE.
982      ENDIF
983      !
984   END SUBROUTINE bias_wrt
985
986END MODULE bias
Note: See TracBrowser for help on using the repository browser.