MODULE bias !! Is used by OPA and STEP !!====================================================================== !! *** Module bias *** !! Code to estimate and apply bias correction. !! The bias is in T/S and Pressure. It follows the generalized !! bias algorithm presented in Balmaseda et al 2007. !! !! It can be read from a file offline, estimated online from relaxation !! terms or from assimilation increments (this latter estimtd in inner loop) !! !! The parameters for the time evolution of the bias can be different !! from the relaxation terms and for the assim increments. Only the !! parameter for the relaxtion terms are considered here. !! !! The offline bias term can contain the seasonal cycle. !! !! The time evolution of the bias for relaxtion is estimated as followed !! bias_rlx(t+1)=t_rlx_mem*bias_rlx(t)+t_rlx_upd*(t/s)trdmp. !! !! The total bias in T/S is partion between the correction to T/S only !! (fb_t) and the correction applied to the pressure gradient (fb_p). !! We impose that (fb_p=1.-fb_t). These factors can be different for the !! different terms(fb_t_rxl,fb_t_asm,fb_t_ofl) !! !! (t/s)bias = fb_t_ofl * (t/s)bias_ofl + !! fb_t_rlx * (t/s)bias_rlx + !! fb_t_asm * (t/s)bias_asm !! (t/s)bias_p =fb_p_ofl * (t/s)bias_ofl+ !! fb_p_rlx * (t/s)bias_rlx_p + !! fb_p_asm * (t/s)bias_asm_p !! (t/s)bias is applied directely to correct T and S !! (t/s)bias_p is used to compute rhd_pc and gru/v_pc !! !! Note: the above is an adhoc /simple way of making the partition !! between bias in pressure and bias in T/S. It would be better !! if the partition is done at the time of estimating the time !! evolution of the bias. That would mean doubling the number of !! 3D arrays. !! !! New addtion: (not in Balmaseda et al 2007): !! A more physical alternative to the partition of the bias can be !! done in terms of the inertial frequency: when the time scales of !! adjustment are shorter that >1/f (Equator), the bias correction should !! be in the the pressure term. Otherwise, it can act directly in T/S. !! NOTE that the bias correction in the pressure term here (following !! (Bell et al 2007) is just the "opposite" as the semi-prognostic method !! in Greatbatch et al 2004. !! The use of this partition is controlled by ln_inertial=.true. !! !! !! 2009-03 (M.A. Balmaseda ECMWF) !!====================================================================== !!---------------------------------------------------------------------- !! bias_init : Read in the bias namelist and the bias arrays !! tra_bias : Apply the bias fields on T/S directly !! dyn_bias : Compute density correction for dynamic hpg !! bias_opn : open bias files for restart capabilities !! bias_wrt : write bias fies " " " !!---------------------------------------------------------------------- !! * Modules used USE par_kind, ONLY: & & wp USE par_oce, ONLY: & & jpi, & & jpj, & & jpk USE dom_oce, ONLY: & & rdt, & & ln_zps, & & gphit USE phycst, ONLY: & & rday, & & rad USE oce, ONLY: & & tsb, tsn, tsa, & & rhop, & & gtsu, gtsv USE dynhpg, ONLY: & & ln_dynhpg_imp USE tradmp USE dtatsd, ONLY: & & ln_tsd_tradmp USE in_out_manager, ONLY: & & lwp, & & numnam_ref, & & numnam_cfg, & & numond, & & numout, & & lrst_oce, & & nit000 USE iom USE eosbn2 USE zpshde ! partial step: hor. derivative (zps_hde routine) USE biaspar USE fldread ! read input fields USE lbclnk ! lateral boundary conditions (or mpp link) USE asmpar USE asminc USE lib_mpp, ONLY: & & ctl_stop, & & ctl_nam IMPLICIT NONE !! * Routine accessibility PRIVATE PUBLIC bias_init, & !: Read in the bias namelist and the bias arrays & tra_bias, & !: Estimate/apply bias on traces & dyn_bias, & !: " density correction for pressure gradient. & bias_opn, & & bias_wrt !! * Substitutions needed to have the variable fsdept_n # include "domzgr_substitute.h90" !! * Shared variables !! * Private module variables REAL(wp), PRIVATE :: & & bias_time_unit_asm, & !: bias_asm units in s ( per day = 86400 s) & bias_time_unit_rlx, & !: bias_rlx units in s ( 1 second) & bias_time_unit_ofl, & !: bias_ofl units in s ( 1 second) & t_rlx_mem, & !: time param for mem in bias_rlx model & t_rlx_upd, & !: time param for update in bias_rlx model !: (pct of incr for computation of bias) & t_asm_mem, & !: time param for mem in bias_asm model & t_asm_upd, & !: time param for update in bias_asm model !: (pct of incr for computation of bias) & fb_t_rlx, & !: parition of bias in T for rlx bias term & fb_t_asm, & !: parition of bias in T for asm bias term & fb_t_ofl, & !: parition of bias in T for ofl bias term & fb_p_rlx, & !: parition of bias in P for rlx bias term & fb_p_asm, & !: parition of bias in P for asm bias term & fb_p_ofl, & !: parition of bias in P for ofl bias term & fctamp, & !: amplification factor for T if inertial & rn_maxlat_bias, & !: Max lat for latitudinal ramp & rn_minlat_bias !: Min lat for latitudinal ramp LOGICAL, PRIVATE :: lalloc REAL(wp), PRIVATE, DIMENSION(:,:,:), ALLOCATABLE :: & & tbias_asm, & !: Temperature bias field & sbias_asm, & !: Salinity bias field & tbias_rlx, & !: Temperature bias field & sbias_rlx, & !: Salinity bias field & tbias_asm_out, & !: Output temperature bias field & sbias_asm_out, & !: Output salinity bias field & tbias_rlx_out, & !: Output temperature bias field & sbias_rlx_out, & !: Output salinity bias field & tbias_p_out, & !: Output temperature bias field for P correction & sbias_p_out !: Output salinity bias field for P correction INTEGER, PRIVATE :: nn_lat_ramp ! choice of latitude dependent ramp ! for the pressure correction. ! 1:inertial ramp, 2:linear ramp, else:no ramp LOGICAL, PRIVATE :: ln_bsyncro ! syncronous or assincrous bias correction LOGICAL, PRIVATE :: ln_itdecay ! evolve bias correction at every time step. REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef INTEGER, PRIVATE :: & & numbias_asm, & ! netcdf id of bias file from assim & numbias_tot, & ! netcdf id of bias file with total bias & nn_bias_itwrt ! time step for outputting bias pressure corr CHARACTER(LEN=128), PRIVATE :: & & cn_bias_asm, & ! name of bias file from assim & cn_bias_tot ! name of bias with total/rlx bias ! Structure of input T and S bias offline (file informations, fields read) TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) :: sf_tbias_ofl TYPE(FLD), PRIVATE, ALLOCATABLE, DIMENSION(:) :: sf_sbias_ofl TYPE(FLD_N), PRIVATE ::& ! information about the fields to be read & sn_tbias_ofl, sn_sbias_ofl CONTAINS SUBROUTINE bias_init !!---------------------------------------------------------------------- !! *** ROUTINE bias_init *** !! !! ** Purpose : Read in the bias namelist and read in the bias fields. !! !! ** Method : Read in the bias namelist and read in the bias fields. !! !! ** Action : !! !! History : !! ! 08-05 (D. Lea) Initial version !! ! 08-10 (M. Martin) Tidy !! ! 09-03 (M. Balmaseda). Generalize to estimate the bias !! from relax and offline bias term. !! Introduce parameters to control the !! model for the bias !! (variables and time evolution) !!---------------------------------------------------------------------- IMPLICIT NONE !! * Local declarations CHARACTER(len=100) :: cn_dir ! dir for location ofline bias INTEGER :: ierror INTEGER :: ios ! Local integer output status for namelist read REAL(wp) :: eft_rlx, & ! efolding time (bias memory) in days & eft_asm, & ! " " & log2, & & lenscl_bias !lengthscale of the pressure bias decay between minlat and maxlat. NAMELIST/nambias/ ln_bias, ln_bias_asm, ln_bias_rlx, ln_bias_ofl, & & ln_bias_ts_app, ln_bias_pc_app, & & fb_t_asm, fb_t_rlx, fb_t_ofl, fb_p_asm, fb_p_rlx, fb_p_ofl, & & eft_rlx, t_rlx_upd, eft_asm, t_asm_upd, nn_lat_ramp, & & bias_time_unit_asm, bias_time_unit_rlx, bias_time_unit_ofl, & & cn_bias_tot, cn_bias_asm, cn_dir, sn_tbias_ofl, sn_sbias_ofl, & & ln_bsyncro, fctamp, rn_maxlat_bias, rn_minlat_bias, & & nn_bias_itwrt, ln_itdecay !----------------------------------------------------------------------- ! Read Namelist : bias interface !----------------------------------------------------------------------- REWIND( numnam_ref ) ! Namelist nambias in reference namelist : Bias pressure correction READ ( numnam_ref, nambias, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in reference namelist', lwp ) ! Set additional default values (note that most values are set in the reference namelist) IF ( ln_asmiau ) nn_bias_itwrt = nitiaufin ! ... default values (NB: frequency positive => hours, negative => months) ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! sn_tbias_ofl = FLD_N( 'tbias_ofl' , -1. , 'tbias' , .TRUE. , .FALSE. , 'yearly', '', '', '' ) sn_sbias_ofl = FLD_N( 'sbias_ofl' , -1. , 'sbias' , .TRUE. , .FALSE. , 'yearly', '', '', '' ) REWIND( numnam_cfg ) ! Namelist nambias in configuration namelist : Bias pressure correction READ ( numnam_cfg, nambias, IOSTAT = ios, ERR = 902 ) 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambias in configuration namelist', lwp ) IF(lwm) WRITE ( numond, nambias ) IF ( ( .NOT. ln_bias_asm ) .AND. ( .NOT. ln_bias_ofl ) .AND. ( .NOT. ln_bias_rlx ) ) THEN ln_bias_ts_app = .FALSE. ln_bias_pc_app = .FALSE. ln_bias = .FALSE. ENDIF ! set up decay scales log2 = LOG( 2.0_wp ) t_rlx_mem = EXP( - log2 * rdt / ( eft_rlx * rday ) ) t_asm_mem = EXP( - log2 * bias_time_unit_asm/ ( eft_asm * rday ) ) ! Control print IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'bias_init : ' WRITE(numout,*) '~~~~~~~~~~~ ' WRITE(numout,*) ' Namelist nambias : ' WRITE(numout,*) ' Bias switches/options/variables ' WRITE(numout,*) ' bias main switch ln_bias = ',ln_bias WRITE(numout,*) ' bias from assim ln_bias_asm = ',ln_bias_asm WRITE(numout,*) ' bias from relax ln_bias_rlx = ',ln_bias_rlx WRITE(numout,*) ' bias from offln ln_bias_ofl = ',ln_bias_ofl WRITE(numout,*) ' bias T and S apply ln_bias_ts_app = ',ln_bias_ts_app WRITE(numout,*) ' bias pressure correctn apply ln_bias_pc_app = ',ln_bias_pc_app WRITE(numout,*) ' bias pressure correctn apply ln_bias_pc_app = ',ln_bias_pc_app WRITE(numout,*) ' lat ramp for bias correction nn_lat_ramp = ',nn_lat_ramp WRITE(numout,*) ' time step for writing bias fld nn_bias_itwrt = ',nn_bias_itwrt WRITE(numout,*) ' evolve pcbias at each timestep ln_itdecay = ',ln_itdecay WRITE(numout,*) ' Parameters for parition of bias term ' WRITE(numout,*) ' fb_t_rlx = ',fb_t_rlx WRITE(numout,*) ' fb_t_asm = ',fb_t_asm WRITE(numout,*) ' fb_t_ofl = ',fb_t_ofl WRITE(numout,*) ' fb_p_rlx = ',fb_p_rlx WRITE(numout,*) ' fb_p_asm = ',fb_p_asm WRITE(numout,*) ' fb_p_ofl = ',fb_p_ofl WRITE(numout,*) ' Parameters for time evolution of bias ' WRITE(numout,*) ' Rlx efolding time (mem) eft_rlx,t_rlx_mem = ', eft_rlx, t_rlx_mem, 1. - log2 * rdt / (eft_rlx * rday) WRITE(numout,*) ' uptdate factor t_rlx_upd = ',t_rlx_upd WRITE(numout,*) ' Asm efolding time (mem) eft_asm,t_asm_mem = ', eft_asm, t_asm_mem, 1. - log2 * rdt / (eft_asm * rday) WRITE(numout,*) ' uptdate factor t_asm_upd = ',t_asm_upd WRITE(numout,*) ' Filenames and input structures' WRITE(numout,*) ' bias_tot filename cn_bias_to = ',cn_bias_tot WRITE(numout,*) ' bias_asm filename cn_bias_asm = ',cn_bias_asm WRITE(numout,*) ' bias_asm time unit (secs) bias_time_unit_asm = ',bias_time_unit_asm WRITE(numout,*) ' structure Tem bias ofl sn_tbias_ofl = ',sn_tbias_ofl WRITE(numout,*) ' structure Sal bias ofl sn_sbias_ofl = ',sn_sbias_ofl IF ( ( (.NOT. ln_tsd_tradmp) .OR. (.NOT. ln_tradmp) ) .AND. ln_bias_rlx ) & & CALL ctl_stop (' lk_dtatem, lk_dtasal and lk_tradmp need to be true with ln_bias_rlx' ) ENDIF IF( .NOT. ln_bias ) RETURN IF( .NOT. lalloc ) THEN ALLOCATE( tbias(jpi,jpj,jpk) , & & sbias(jpi,jpj,jpk) , & & tbias_p(jpi,jpj,jpk), & & sbias_p(jpi,jpj,jpk), & & rhd_pc(jpi,jpj,jpk) , & & gru_pc(jpi,jpj) , & & grv_pc(jpi,jpj) ) ALLOCATE( fbcoef(jpi,jpj) ) IF( ln_bias_asm ) ALLOCATE( tbias_asm(jpi,jpj,jpk), & & sbias_asm(jpi,jpj,jpk), & tbias_asm_out(jpi,jpj,jpk), & sbias_asm_out(jpi,jpj,jpk), & tbias_p_out(jpi,jpj,jpk), & sbias_p_out(jpi,jpj,jpk) ) IF( ln_bias_rlx ) ALLOCATE( tbias_rlx(jpi,jpj,jpk), & & sbias_rlx(jpi,jpj,jpk), & tbias_rlx_out(jpi,jpj,jpk), & sbias_rlx_out(jpi,jpj,jpk) ) lalloc = .TRUE. ENDIF IF( ln_bias_ofl ) THEN ! set sf_tbias_ofl and sf_sbias_ofl strctrs ! ! tbias ! ALLOCATE( sf_tbias_ofl(1), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'bias_init: unable to allocate sf_tbias_ofl structure' ) ; RETURN ENDIF ALLOCATE( sf_tbias_ofl(1)%fnow(jpi,jpj,jpk) ) ALLOCATE( sf_tbias_ofl(1)%fdta(jpi,jpj,jpk,2) ) ! fill structure with values and control print CALL fld_fill( sf_tbias_ofl, (/ sn_tbias_ofl /), cn_dir, 'bias_init', 'Offline T bias term ', 'nam_tbias_ofl' ) ! ! salinity bias ! ALLOCATE( sf_sbias_ofl(1), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'bias_init: unable to allocate sf_sbias_ofl structure' ) ; RETURN ENDIF ALLOCATE( sf_sbias_ofl(1)%fnow(jpi,jpj,jpk) ) ALLOCATE( sf_sbias_ofl(1)%fdta(jpi,jpj,jpk,2) ) ! fill structure with values and control print CALL fld_fill( sf_sbias_ofl, (/ sn_sbias_ofl /), cn_dir, 'bias_init', 'Offline S bias term ', 'nam_sbias_ofl' ) ENDIF ! Read total bias IF ( ln_bias ) THEN tbias(:,:,:) = 0.0_wp sbias(:,:,:) = 0.0_wp tbias_p(:,:,:) = 0.0_wp sbias_p(:,:,:) = 0.0_wp gru_pc(:,:) = 0.0_wp grv_pc(:,:) = 0.0_wp IF ( ln_bias_rlx ) THEN tbias_rlx(:,:,:) = 0.0_wp sbias_rlx(:,:,:) = 0.0_wp ENDIF IF ( ln_bias_asm ) THEN !now rlx and asm bias in same file tbias_asm(:,:,:) = 0.0_wp sbias_asm(:,:,:) = 0.0_wp ENDIF numbias_tot = 0 ! Get bias from file and prevent fail if the file does not exist IF(lwp) WRITE(numout,*) 'Opening ',TRIM( cn_bias_tot ) CALL iom_open( cn_bias_tot, numbias_tot, ldstop=.FALSE. ) IF ( numbias_tot > 0 ) THEN ! Could check validity time of bias fields here... ! Get the T and S bias data IF(lwp) WRITE(numout,*) 'Reading bias fields from tot...' !Search for bias from relaxation term if needed. Use same file IF ( ln_bias_rlx ) THEN IF(lwp) WRITE(numout,*) 'Reading bias fields for bias rlx from file ',cn_bias_tot IF( iom_varid( numbias_tot, 'tbias_rlx' ) > 0 ) THEN ! Get the T and S bias data CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_rlx', tbias_rlx ) CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_rlx', sbias_rlx ) ELSE CALL ctl_stop( 'Bias relaxation variables not found in ',cn_bias_tot ) ENDIF ENDIF !Search for bias from assim term if needed. Use same file IF ( ln_bias_asm .and. .not.ln_bsyncro ) THEN IF(lwp) WRITE(numout,*) 'Reading a-syncro bias fields for bias asm from file ',cn_bias_tot IF( iom_varid( numbias_tot, 'tbias_asm' ) > 0 ) THEN ! Get the T and S bias data CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm', tbias_asm ) CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm', sbias_asm ) ELSE CALL ctl_stop( 'Bias assim variables not found in ',cn_bias_tot ) ENDIF ENDIF ! Close the file CALL iom_close(numbias_tot) ELSE IF(lwp) WRITE(numout,*) 'No bias file found so T and S bias fields are set to zero' ENDIF ENDIF ! for the time being, the bias_asm is read in the same file as ! bias_rlx ! Implications: bias_asm is estimated/evolved in time in the second outer ! loop only, when the assimilation increments are ready. ! bias_asm is kept cte during the first outer loop. ! => Assyncronous bias correction. ! Alternative: Syncronous bias correction: ! bias_asm estimated/evolved in the first outer loop ! with the asm incrments of the previous cycle. ! bias_asm kept cte during the second outer loop. ! Implication: bias_asm should be estimated really in the ! inner loop. IF ( ln_bsyncro ) THEN ! Read bias from assimilation from a separate file IF ( ln_bias_asm ) THEN tbias_asm(:,:,:) = 0.0_wp sbias_asm(:,:,:) = 0.0_wp numbias_asm = 0 ! Get bias from file and prevent fail if the file does not exist IF(lwp) WRITE(numout,*) 'Opening file for syncro assim bias ',TRIM( cn_bias_asm ) CALL iom_open( cn_bias_asm, numbias_asm, ldstop=.FALSE. ) IF ( numbias_asm > 0 ) THEN ! Could check validity time of bias fields here... ! Get the T and S bias data IF(lwp) WRITE(numout,*) 'Reading syncro bias fields from asm from file ',cn_bias_asm CALL iom_get( numbias_asm, jpdom_autoglo, 'tbias_asm', tbias_asm ) CALL iom_get( numbias_asm, jpdom_autoglo, 'sbias_asm', sbias_asm ) ! this is only applicable if tbias_asm were to be calculated in the inner loop tbias_asm(:,:,:) = tbias_asm(:,:,:) * rdt / bias_time_unit_asm sbias_asm(:,:,:) = sbias_asm(:,:,:) * rdt / bias_time_unit_asm ! Close the file CALL iom_close(numbias_asm) ELSE IF(lwp) WRITE(numout,*) 'No bias file found from asm so T and S bias fields are set to zero' ENDIF ENDIF ENDIF !latitudinal dependence of partition coeficients. Adhoc IF ( nn_lat_ramp == 1 ) THEN ! Use the inertial ramp. lenscl_bias = ( rn_maxlat_bias - rn_minlat_bias )*2._wp WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias ) fbcoef(:,:) = 0._wp ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias ) fbcoef(:,:) = 1._wp ELSEWHERE fbcoef(:,:) = 1._wp - exp( -( abs( gphit(:,:) ) - rn_minlat_bias ) & * ( abs( gphit(:,:) ) - rn_minlat_bias ) / lenscl_bias ) ENDWHERE ELSEIF ( nn_lat_ramp == 2 ) THEN ! Use a linear ramp consist with the geostrophic velocity balance ramp in NEMOVAR WHERE ( abs( gphit(:,:) ) <= rn_minlat_bias ) fbcoef(:,:) = 0._wp ELSEWHERE ( abs( gphit(:,:) ) >= rn_maxlat_bias ) fbcoef(:,:) = 1._wp ELSEWHERE fbcoef(:,:) = 1._wp - ((rn_maxlat_bias - abs( gphit(:,:)))/(rn_maxlat_bias - rn_minlat_bias)) ENDWHERE ELSE fbcoef(:,:) = 0.0_wp fctamp = 0.0_wp ENDIF END SUBROUTINE bias_init SUBROUTINE tra_bias ( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_bias *** !! !! ** Purpose : Update bias field every time step !! !! ** Method : add contributions to bias from 3 terms !! !! ** Action : Bias from assimilation (read in bias_init) !! Bias from relaxation term is estimated according to !! the prescribed time evolution of the bias !! Bias from ofl term is read from external file !! The difference contributions are added and the partition !! into direct bias in T/S and pressure perfomed. !! !! History : 09-03 (M. Balmaseda) !!---------------------------------------------------------------------- !! called every timestep after dta_sst if ln_bias true. IMPLICIT NONE !! * Arguments INTEGER, INTENT(IN) :: kt ! ocean time-step index !! * Local variables INTEGER :: ji,jj,jk, it ! local loop index REAL(wp) :: tsclf ! time saling factor REAL(wp) :: fb_t_asm_max, fb_t_rlx_max, fb_t_ofl_max REAL(wp) :: ztfrac, ztsday REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2 IF ( .NOT. ln_bias ) RETURN fb_t_rlx_max = MIN(fb_t_rlx*fctamp,1.0_wp) fb_t_asm_max = MIN(fb_t_asm*fctamp,1.0_wp) fb_t_ofl_max = MIN(fb_t_ofl*fctamp,1.0_wp) tbias(:,:,:) = 0.0_wp sbias(:,:,:) = 0.0_wp tbias_p(:,:,:) = 0.0_wp sbias_p(:,:,:) = 0.0_wp IF ( ln_bias_asm ) THEN tsclf = 1 IF ( .NOT.ln_bsyncro ) tsclf = rdt / bias_time_unit_asm zcof1(:,:) = tsclf * ( ( 1.0_wp - fbcoef(:,:) ) * fb_t_asm + & & fbcoef(:,:) * fb_t_asm_max ) zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_asm IF ( ln_itdecay ) THEN !decay the pressure correction at each time step ztsday = rday / real(rdt) IF( ln_asmiau .and. ln_trainc ) THEN ! nudge in increments and decay historical contribution IF ( kt <= nitiaufin ) THEN ! During IAU calculate the fraction of increments to apply at each time step ztfrac = real(kt) / real(nitiaufin) ! nudging factor during the IAU IF (lwp) THEN WRITE(numout,*) 'tra_bias : bias weights' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,* ) "proportion of increment applied in pcbias ",ztfrac WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday) ENDIF DO jk = 1, jpkm1 tbias(:,:,jk) = tbias(:,:,jk) + & & ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) + & & ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:) sbias(:,:,jk) = sbias(:,:,jk) + & & ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) + & & ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof1(:,:) tbias_p(:,:,jk) = tbias_p(:,:,jk) + & & ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) + & & ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:) sbias_p(:,:,jk) = sbias_p(:,:,jk) + & & ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) + & & ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:) ENDDO IF ( .not.ln_bsyncro ) THEN IF ( kt == nn_bias_itwrt ) THEN DO jk = 1, jpkm1 tbias_asm_out(:,:,jk) = t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) + & & ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) sbias_asm_out(:,:,jk) = t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) + & & ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) END DO ENDIF ENDIF ! update the historical component with the increments at the end of the IAU IF ( kt == nitiaufin ) THEN DO jk = 1, jpkm1 tbias_asm(:,:,jk) = t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) + & & ztfrac * t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) sbias_asm(:,:,jk) = t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) + & & ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) END DO ENDIF ELSE ! decay pressure correction from combined historical component and increments after IAU ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin)) ! decay from end of IAU DO jk = 1, jpkm1 tbias(:,:,jk) = tbias(:,:,jk) + & & ( ztfrac * tbias_asm(:,:,jk) ) * zcof1(:,:) sbias(:,:,jk) = sbias(:,:,jk) + & & ( ztfrac * sbias_asm(:,:,jk) ) * zcof1(:,:) tbias_p(:,:,jk) = tbias_p(:,:,jk) + & & ( ztfrac * tbias_asm(:,:,jk) ) * zcof2(:,:) sbias_p(:,:,jk) = sbias_p(:,:,jk) + & & ( ztfrac * sbias_asm(:,:,jk) ) * zcof2(:,:) ENDDO IF (lwp) THEN WRITE(numout,*) 'tra_bias : bias weights' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,* ) "proportion of increments + historical pcbias applied ",ztfrac ENDIF IF ( ln_trainc .and. .not.ln_bsyncro ) THEN IF ( kt == nn_bias_itwrt ) THEN DO jk = 1, jpkm1 tbias_asm_out(:,:,jk) = ztfrac * tbias_asm(:,:,jk) sbias_asm_out(:,:,jk) = ztfrac * sbias_asm(:,:,jk) END DO ENDIF ENDIF ENDIF ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts) DO jk = 1, jpkm1 tbias(:,:,jk) = tbias(:,:,jk) + & & ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof1(:,:) sbias(:,:,jk) = sbias(:,:,jk) + & & ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof1(:,:) tbias_p(:,:,jk) = tbias_p(:,:,jk) + & & ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof2(:,:) sbias_p(:,:,jk) = sbias_p(:,:,jk) + & & ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof2(:,:) ENDDO IF (lwp) THEN WRITE(numout,*) 'tra_bias : bias weights' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday) ENDIF ENDIF ELSE ! maintain a constant pressure correction DO jk = 1, jpkm1 tbias(:,:,jk) = tbias(:,:,jk) + tbias_asm(:,:,jk) * zcof1(:,:) sbias(:,:,jk) = sbias(:,:,jk) + sbias_asm(:,:,jk) * zcof1(:,:) tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_asm(:,:,jk) * zcof2(:,:) sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_asm(:,:,jk) * zcof2(:,:) END DO IF( ln_asmiau .and. ln_trainc .and. .not.ln_bsyncro ) THEN ! if last outer loop (ln_asmiau=true and ln_trainc=true). t/sbias_asm ! is updated, only once (end of run) taking into account units. IF ( kt == nn_bias_itwrt ) THEN IF(lwp) WRITE(numout,*)' estimating asm bias at timestep: ',kt DO jk = 1, jpkm1 tbias_asm_out(:,:,jk) = t_asm_mem * tbias_asm(:,:,jk) + & & t_asm_upd * t_bkginc(:,:,jk) * tmask(:,:,jk) sbias_asm_out(:,:,jk) = t_asm_mem * sbias_asm(:,:,jk) + & & t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) END DO ENDIF ENDIF ENDIF ENDIF #if defined key_tradmp ! Time evolution of bias from relaxation IF ( ln_bias_rlx ) THEN tbias_rlx(:,:,:) = t_rlx_mem * tbias_rlx(:,:,:) + & & t_rlx_upd * ttrdmp(:,:,:) * rdt sbias_rlx(:,:,:) = t_rlx_mem * sbias_rlx(:,:,:) + & & t_rlx_upd * strdmp(:,:,:) * rdt zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_rlx + & & fbcoef(:,:) * fb_t_rlx_max zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_rlx DO jk = 1, jpkm1 tbias(:,:,jk) = tbias(:,:,jk) + tbias_rlx(:,:,jk) * zcof1(:,:) sbias(:,:,jk) = sbias(:,:,jk) + sbias_rlx(:,:,jk) * zcof1(:,:) tbias_p(:,:,jk) = tbias_p(:,:,jk) + tbias_rlx(:,:,jk) * zcof2(:,:) sbias_p(:,:,jk) = sbias_p(:,:,jk) + sbias_rlx(:,:,jk) * zcof2(:,:) ENDDO IF ( kt == nn_bias_itwrt ) THEN tbias_rlx_out(:,:,:) = tbias_rlx(:,:,:) sbias_rlx_out(:,:,:) = sbias_rlx(:,:,:) ENDIF ENDIF #endif ! ofline bias IF ( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) ' tra_bias: ln_bias_ofl = ',ln_bias_ofl IF(lwp) WRITE(numout,*) ' ~~~~~~~~~' ENDIF IF ( ln_bias_ofl ) THEN IF(lwp) WRITE(numout,*) 'reading offline bias' CALL fld_read( kt, 1, sf_tbias_ofl ) CALL fld_read( kt, 1, sf_sbias_ofl ) zcof1(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_t_ofl + & & fbcoef(:,:) * fb_t_ofl_max zcof2(:,:) = ( 1.0_wp - fbcoef(:,:) ) * fb_p_ofl DO jk = 1, jpkm1 tbias(:,:,jk) = tbias(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:) sbias(:,:,jk) = sbias(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof1(:,:) tbias_p(:,:,jk) = tbias_p(:,:,jk) + sf_tbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:) sbias_p(:,:,jk) = sbias_p(:,:,jk) + sf_sbias_ofl(1)%fnow(:,:,jk) * zcof2(:,:) ENDDO ENDIF !apply bias on tracers if needed IF ( ln_bias_ts_app ) THEN ! Add the bias directely to T/s tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + tmask(:,:,:) * tbias(:,:,:) / rdt tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + tmask(:,:,:) * sbias(:,:,:) / rdt ! lateral boundary conditions (is this needed?) CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1.0_wp ) CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1.0_wp ) ENDIF IF ( kt == nn_bias_itwrt ) THEN tbias_p_out(:,:,:) = tbias_p(:,:,:) sbias_p_out(:,:,:) = sbias_p(:,:,:) ENDIF END SUBROUTINE tra_bias SUBROUTINE dyn_bias( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_bias *** !! !! ** Purpose : Computes rhd_pc, gru/v_pc bias corrected !! for hydrostatic pressure gradient !! depending on time step (semi-implicit/centered) !! If partial steps computes bottom pressure gradient. !! These correction terms will affect only the dynamical !! component (pressure gradient calculation), but not !! the isopycnal calculation for the lateral mixing. !! !! ** Method : At this stage of the computation, ta and sa are the !! after temperature and salinity. If semi-implicit, these !! are used to compute rho and bottom pressure gradient. !! If centered, tb,sb are used instead. !! If bias key is activated, the temperature,salinity are !! bias corrected in the calculation of the density fields !! used in the pressure gradient calculation. !! !! !! ** Action : - rhd_pc ready. rhop will be overwriten later !! - if ln_zps, bottom density gradients gru/v_pc ready. !!---------------------------------------------------------------------- !! !! * Arguments INTEGER, INTENT(IN) :: kt ! ocean time-step index !! * Local variables REAL(wp) :: tsw(jpi,jpj,jpk,jpts) !! !!---------------------------------------------------------------------- ! ! gtu,gsu,gtv,gsv rhop will be overwritten later in step. ! IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:) tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:) ELSE tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:) tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:) ENDIF CALL eos( tsw, rhd_pc, rhop , fsdept_n(:,:,:)) ! is this needed? CALL lbc_lnk( rhd_pc, 'T', 1.0_wp ) ! Partial steps: now horizontal gradient of t,s,rd ! at the bottom ocean level IF( ln_zps ) THEN CALL zps_hde( kt, jpts, tsw, gtsu, gtsv, & & rhd_pc, gru_pc , grv_pc ) ENDIF END SUBROUTINE dyn_bias SUBROUTINE bias_opn( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE bias_opn *** !! !! ** Purpose : open bias restart file following the same logic as the !! standard restarts. !!---------------------------------------------------------------------- !! * Arguments INTEGER, INTENT(IN) :: kt ! ocean time-step !! * Local variables CHARACTER(LEN=20) :: clbkt ! ocean time-step deine as a character CHARACTER(LEN=50) :: clbias_tot ! total bias restart file name !!---------------------------------------------------------------------- ! IF( lrst_oce .AND. .NOT.lrst_bias ) THEN ! create bias file IF( nitend > 999999999 ) THEN ; WRITE(clbkt, * ) nitend ELSE ; WRITE(clbkt, '(i8.8)') nitend ENDIF clbias_tot = TRIM(cexper)//"_"//TRIM(ADJUSTL(clbkt))//"_"//TRIM(cn_bias_tot) IF(lwp) THEN WRITE(numout,*) SELECT CASE ( jprstlib ) CASE ( jprstdimg ) ; WRITE(numout,*) ' open tot bias restart binary file: '//clbias_tot CASE DEFAULT ; WRITE(numout,*) ' open tot bias restart NetCDF file: '//clbias_tot END SELECT ENDIF CALL iom_open( clbias_tot, numbias_tot , ldwrt = .TRUE., kiolib = jprstlib ) lrst_bias=.TRUE. ENDIF ! END SUBROUTINE bias_opn SUBROUTINE bias_wrt( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE bias_wrt *** !! !! ** Purpose : Write bias restart fields in the format corresponding to jprstlib !! !! ** Method : Write in numbias_tot when kt == nitend in output !! file, save fields which are necessary for restart. !! !! Changed the timestep for writing to nitend rather than nitrst as this causes a !! problem if the nstock list is used to determine the restart output times and !! means that the bias is not output at all. M. Martin. 08/2011. !! Need to check with Magdalena that this is ok for ECMWF. !!---------------------------------------------------------------------- !! * Arguments INTEGER, INTENT(IN) :: kt ! ocean time-step !!---------------------------------------------------------------------- ! ! the begining of the run [s] IF ( ln_bias_rlx ) THEN CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_rlx' , tbias_rlx_out ) CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_rlx' , sbias_rlx_out ) ENDIF IF ( ln_bias_asm ) THEN CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm' , tbias_asm_out ) CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm' , sbias_asm_out ) CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_p' , tbias_p_out ) CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p' , sbias_p_out ) ENDIF IF( kt == nitend ) THEN CALL iom_close( numbias_tot ) ! close the restart file (only at last time step) lrst_bias = .FALSE. ENDIF ! END SUBROUTINE bias_wrt END MODULE bias