New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
bias.F90 in branches/UKMO/dev_r5518_pcbias_ipc/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

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

Last change on this file since 7653 was 7653, checked in by isabella, 7 years ago

adding changes to reset variables following options

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