Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8446)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8447)
@@ -1363,3 +1363,4 @@
    nn_bias_itwrt  = 15
    ln_itdecay     = .FALSE.
-/
+   ln_incpc       = .FALSE.
+/
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90	(revision 8446)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90	(revision 8447)
@@ -76,5 +76,5 @@
    USE oce, ONLY: &
       & tsb, tsn, tsa, &
-      & rhop,          &
+      & rhop,  &
       & gtsu, gtsv
    USE dynhpg, ONLY:   &
@@ -134,6 +134,8 @@
       & 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
+      & rn_maxlat_bias,       &  !: Max lat for latitudinal ramp
+      & rn_minlat_bias,       &  !: Min lat for latitudinal ramp
+      & zwgt,                 &  !: weight for IPC
+      & ztscale                  !: decay rate for IPC
 
    LOGICAL,  PRIVATE :: lalloc
@@ -148,5 +150,11 @@
       & 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
+      & sbias_p_out,   &   !: Output salinity bias field for P correction
+      & tbias_i_out,   &   !: Output temperature bias field for incremental P correction
+      & sbias_i_out,   &   !: Output salinity bias field for incremental P correction
+      & tbias_asm_stscale, &   !: Short time scale temperature bias field
+      & sbias_asm_stscale, &   !: Short time scale salinity bias field
+      & tbias_asm_stscale_out, &   !: Short time scale temperature bias output field
+      & sbias_asm_stscale_out   !: Short time scale salinity bias output field
 
    INTEGER, PRIVATE :: nn_lat_ramp     ! choice of latitude dependent ramp
@@ -154,7 +162,10 @@
 				       ! 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  
+   LOGICAL, PRIVATE :: ln_itdecay      ! evolve bias correction at every time step.  
+   LOGICAL, PRIVATE :: ln_incpc        ! incremental pressure correction plus pressure correction
+   LOGICAL, PRIVATE :: ln_incpc_only        ! incremental pressure correction only
+
+   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef
+   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef_stscale
 
    INTEGER, PRIVATE ::  &
@@ -200,11 +211,13 @@
       !! * Local declarations
 
-      CHARACTER(len=100) ::  cn_dir       ! dir for location ofline bias
+      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,  &  !     "      "
+      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.
+         &                   lenscl_bias, &  ! lengthscale pressure bias decay between minlat and maxlat.
+         &                   minlat_bias, &  ! used in ipc
+         &                   maxlat_bias     ! used in ipc
       
       NAMELIST/nambias/ ln_bias, ln_bias_asm, ln_bias_rlx, ln_bias_ofl,   &
@@ -215,5 +228,5 @@
          & 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
+         & nn_bias_itwrt, ln_itdecay, ln_incpc,ln_incpc_only, ztscale, zwgt
  
 
@@ -274,4 +287,8 @@
          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,*) '     incremental press. correction  ln_incpc       = ',ln_incpc
+	 WRITE(numout,*) '     incremental press. correction only  ln_incpc_only       = ',ln_incpc_only
+	 WRITE(numout,*) '     incremental press. correction  zwgt           = ',zwgt
+	 WRITE(numout,*) '     incremental press. correction  ztscale        = ',ztscale
          WRITE(numout,*) '  Parameters for parition of bias term '
          WRITE(numout,*) '                                    fb_t_rlx       = ',fb_t_rlx
@@ -296,4 +313,11 @@
             &   CALL ctl_stop (' lk_dtatem, lk_dtasal and lk_tradmp need to be true with ln_bias_rlx' )
 
+         IF ( (.NOT. ln_itdecay) .AND. ln_incpc) &   
+            &   CALL ctl_stop (' if you set ln_incpc to .true. then you need to set ln_itdecay to .true. as well' )
+
+         IF ( (.NOT. ln_incpc) .AND. ln_incpc_only) &   
+            &   CALL ctl_stop (' if you set ln_incpc_only to .true. then you need to set ln_incpc to .true. as well' )
+         
+         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 
       ENDIF
       IF( .NOT. ln_bias ) RETURN
@@ -305,9 +329,11 @@
             &      tbias_p(jpi,jpj,jpk), &
             &      sbias_p(jpi,jpj,jpk), &
+            &      tbias_i(jpi,jpj,jpk), &
+            &      sbias_i(jpi,jpj,jpk), &
             &      rhd_pc(jpi,jpj,jpk) , &
             &      gru_pc(jpi,jpj)     , &
             &      grv_pc(jpi,jpj)       )
 
-         ALLOCATE( fbcoef(jpi,jpj)       )
+         ALLOCATE( fbcoef(jpi,jpj), fbcoef_stscale(jpi,jpj) )
          
          IF( ln_bias_asm ) ALLOCATE(  tbias_asm(jpi,jpj,jpk),     &
@@ -322,4 +348,12 @@
                                       tbias_rlx_out(jpi,jpj,jpk), &  
                                       sbias_rlx_out(jpi,jpj,jpk)  )
+
+         IF( ln_incpc )    ALLOCATE(  tbias_asm_stscale(jpi,jpj,jpk), &
+            &                         sbias_asm_stscale(jpi,jpj,jpk), &
+                                      tbias_asm_stscale_out(jpi,jpj,jpk), &
+                                      sbias_asm_stscale_out(jpi,jpj,jpk), &
+                                      tbias_i_out(jpi,jpj,jpk),   &  
+                                      sbias_i_out(jpi,jpj,jpk)   )
+
          lalloc = .TRUE.
 
@@ -360,14 +394,29 @@
          tbias_p(:,:,:) = 0.0_wp
          sbias_p(:,:,:) = 0.0_wp
+         tbias_i(:,:,:) = 0.0_wp
+         sbias_i(:,:,:) = 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
+            tbias_asm(:,:,:) = 0.0_wp
+            sbias_asm(:,:,:) = 0.0_wp
+            tbias_asm_out(:,:,:) = 0.0_wp
+            sbias_asm_out(:,:,:) = 0.0_wp
+         ENDIF
+
+         IF ( ln_incpc ) THEN   !incr pressure correction
+            tbias_asm_stscale(:,:,:) = 0.0_wp
+            sbias_asm_stscale(:,:,:) = 0.0_wp
+            tbias_asm_stscale_out(:,:,:) = 0.0_wp
+            sbias_asm_stscale_out(:,:,:) = 0.0_wp
+         ENDIF
+
+
          numbias_tot    = 0
          ! Get bias from file and prevent fail if the file does not exist
@@ -404,4 +453,17 @@
                ENDIF
             ENDIF
+
+            IF ( ln_incpc .and. .not.ln_bsyncro ) THEN
+               IF(lwp) WRITE(numout,*) 'Reading short time scale bias correction fields for bias asm from file ',cn_bias_tot
+               IF( iom_varid( numbias_tot, 'tbias_asm_stscale' ) > 0 ) THEN
+                  ! Get the T and S bias data
+                  CALL iom_get( numbias_tot, jpdom_autoglo, 'tbias_asm_stscale', tbias_asm_stscale )
+                  CALL iom_get( numbias_tot, jpdom_autoglo, 'sbias_asm_stscale', sbias_asm_stscale )
+               ELSE
+                  CALL ctl_stop( 'Short time scale bias assim variables not found in ',cn_bias_tot )
+               ENDIF
+            ENDIF  
+
+
 
             ! Close the file
@@ -484,5 +546,19 @@
          fbcoef(:,:) = 0.0_wp
          fctamp      = 0.0_wp
-      ENDIF
+         fbcoef_stscale(:,:) = 0.0_wp
+      ENDIF
+
+
+      IF ( ln_incpc) THEN
+         minlat_bias = 3.0_wp
+         maxlat_bias = 8.0_wp  
+         WHERE ( abs( gphit(:,:) ) <= minlat_bias )
+            fbcoef_stscale(:,:)=0._wp
+         ELSEWHERE ( abs( gphit(:,:) ) >= maxlat_bias ) 
+            fbcoef_stscale(:,:)=1._wp
+         ELSEWHERE
+            fbcoef_stscale(:,:)=1._wp - ((maxlat_bias - abs( gphit(:,:)))/(maxlat_bias-minlat_bias))
+         ENDWHERE   
+      ENDIF 
 
 
@@ -517,4 +593,6 @@
       REAL(wp)            ::   fb_t_asm_max, fb_t_rlx_max, fb_t_ofl_max
       REAL(wp)            ::   ztfrac, ztsday
+      REAL(wp)            ::   zfrac, zfrac1 ! temporal weights for inst pcbias (names could be changed)
+      REAL(wp)            ::   zdecay        ! used in inst pcorr
       REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2
 
@@ -528,4 +606,7 @@
       tbias_p(:,:,:) = 0.0_wp
       sbias_p(:,:,:) = 0.0_wp
+      tbias_i(:,:,:) = 0.0_wp
+      sbias_i(:,:,:) = 0.0_wp
+
       IF ( ln_bias_asm ) THEN
       
@@ -535,16 +616,18 @@
             &                              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)
+	    ztsday  = rday / REAL(rdt,wp)
+
+            zdecay = (1-ztscale)**(1/REAL(ztsday,wp)) ! used in ipc
+            zfrac1 = zdecay**REAL(kt,wp) ! used in ipc
+            IF ( zfrac1 <= 0.0 ) zfrac1 = 0.0_wp
 
             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
-
+                  ztfrac = REAL(kt,wp) / REAL(nitiaufin,wp)  ! nudging factor during the IAU
 	       
                   IF (lwp) THEN
@@ -552,31 +635,63 @@
                      WRITE(numout,*) '~~~~~~~~~~~~'
                      WRITE(numout,* ) "proportion of  increment applied in pcbias ",ztfrac
-                     WRITE(numout,* ) "proportion of  historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
+                     WRITE(numout,* ) "proportion of  historical pcbias applied ",t_asm_mem**(REAL(kt,wp)/ztsday)
                   ENDIF
-	        
-                  DO jk = 1, jpkm1
+
+                  DO jk = 1, jpkm1  
                      tbias(:,:,jk) = tbias(:,:,jk) +                            &
-                     &                ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +                    &
+                     &                ( t_asm_mem**(REAL(kt,wp)/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)  +                     &
+                     &               ( t_asm_mem**(REAL(kt,wp)/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)  +                     &
+                     &               ( t_asm_mem**(REAL(kt,wp)/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)  +                     &
+                     &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk)  +                     &
                      &               ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk) ) * zcof2(:,:)
                   ENDDO
 
+                  IF (ln_incpc) THEN
+
+                     IF (lwp) THEN
+                       WRITE(numout,*) 'tra_bias : bias weights'
+                       WRITE(numout,*) '~~~~~~~~~~~~'
+                       WRITE(numout,* ) "IPC - proportion of  increment applied in pcbias ",ztfrac
+                       WRITE(numout,* ) "IPC - proportion of  historical pcbias applied ",zfrac1
+                     ENDIF
+
+                     DO jk = 1, jpkm1
+
+                        tbias_i(:,:,jk) =  ( t_bkginc(:,:,jk) * tmask(:,:,jk)* zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) )         &
+                        &                + ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
+                        sbias_i(:,:,jk) =  ( s_bkginc(:,:,jk) * tmask(:,:,jk)* zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) )         &
+                        &                + ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
+
+                     ENDDO
+
+                  ENDIF
+
                   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)  +       &
+                           tbias_asm_out(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/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)  +       &
+                           sbias_asm_out(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk)  +       &
                            &                     ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
                          END DO
+
+                        IF ( ln_incpc) THEN
+                           DO jk = 1, jpkm1
+                              tbias_asm_stscale_out(:,:,jk) = ( t_bkginc(:,:,jk) * tmask(:,:,jk) *  zwgt * ztfrac ) + ( tbias_asm_stscale(:,:,jk) * zfrac1 )
+                              sbias_asm_stscale_out(:,:,jk) = ( s_bkginc(:,:,jk) * tmask(:,:,jk) *  zwgt * ztfrac ) + ( sbias_asm_stscale(:,:,jk) * zfrac1 )
+                           ENDDO
+                        ENDIF
+
                      ENDIF
+
                   ENDIF
 
@@ -584,15 +699,23 @@
                   IF ( kt == nitiaufin ) THEN
                      DO jk = 1, jpkm1
-                        tbias_asm(:,:,jk) =  t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk)  +       &
+                        tbias_asm(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/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)  +       &
+                        sbias_asm(:,:,jk) =  t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk)  +       &
                         &                     ztfrac * t_asm_upd * s_bkginc(:,:,jk) * tmask(:,:,jk)
                      END DO
+
+                     IF (ln_incpc) THEN
+                        DO jk = 1, jpkm1
+                           tbias_asm_stscale(:,:,jk) = ( t_bkginc(:,:,jk) * tmask(:,:,jk) * zwgt * ztfrac ) + ( tbias_asm_stscale(:,:,jk) * zfrac1 )
+                           sbias_asm_stscale(:,:,jk) = ( s_bkginc(:,:,jk) * tmask(:,:,jk) * zwgt * ztfrac ) + ( sbias_asm_stscale(:,:,jk) * zfrac1 )
+                        ENDDO
+                     ENDIF
+     
                   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
-
+                  ztfrac=t_asm_mem**(REAL(kt-nitiaufin,wp)/REAL(nitiaufin,wp))  ! decay from end of IAU
+                  
                   DO jk = 1, jpkm1
                      tbias(:,:,jk) = tbias(:,:,jk) +                            &
@@ -606,33 +729,58 @@
                   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
-
+                 IF (ln_incpc) THEN
+
+                   zfrac  = zdecay**REAL(kt-nitiaufin,wp) 
+                   IF ( zfrac <= 0.0 ) zfrac = 0.0_wp
+   
+                   DO jk = 1, jpkm1
+                      tbias_i(:,:,jk) =  tbias_asm_stscale(:,:,jk) * zfrac * (1.0 - fbcoef_stscale(:,:))
+                      sbias_i(:,:,jk) =  sbias_asm_stscale(:,:,jk) * zfrac * (1.0 - fbcoef_stscale(:,:))
+                   ENDDO
+
+                   IF (lwp) THEN
+                      WRITE(numout,*) 'tra_bias : bias weights'
+                      WRITE(numout,*) '~~~~~~~~~~~~'
+                      WRITE(numout,* ) "IPC - proportion of increments + historical pcbias applied ",zfrac
+                   ENDIF
+
+                 ENDIF
+
+                 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
+
+                       IF ( ln_incpc ) THEN
+                          IF (lwp) WRITE(numout,*) 'after end of IAU - IPC - saving s/tbias_asm_stscale'
+                             DO jk = 1, jpkm1
+                                tbias_asm_stscale_out(:,:,jk) =  tbias_asm_stscale(:,:,jk) * zfrac
+                                sbias_asm_stscale_out(:,:,jk) =  sbias_asm_stscale(:,:,jk) * zfrac
+                             ENDDO
+                          ENDIF
+                       ENDIF
+                    ENDIF
+                 ENDIF
 
             ELSE ! no assimilation increments, simply decay pressure correction (e.g for forecasts)
+                 ! this occurs at obsoper step as well ( ln_asmiau is false ) 
 
                DO jk = 1, jpkm1
                   tbias(:,:,jk) = tbias(:,:,jk) +                                                         &
-                  &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof1(:,:)
+                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk) ) * zcof1(:,:)
                   sbias(:,:,jk) = sbias(:,:,jk) +                                                         &
-                  &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof1(:,:)
+                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk) ) * zcof1(:,:)
                   tbias_p(:,:,jk) = tbias_p(:,:,jk) +                                                     &
-                  &               ( t_asm_mem**(real(kt)/ztsday) * tbias_asm(:,:,jk) ) * zcof2(:,:)
+                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * tbias_asm(:,:,jk) ) * zcof2(:,:)
                   sbias_p(:,:,jk) = sbias_p(:,:,jk) +                                                     &  
-                  &               ( t_asm_mem**(real(kt)/ztsday) * sbias_asm(:,:,jk) ) * zcof2(:,:)
+                  &               ( t_asm_mem**(REAL(kt,wp)/ztsday) * sbias_asm(:,:,jk) ) * zcof2(:,:)
                ENDDO
 
@@ -640,7 +788,16 @@
                   WRITE(numout,*) 'tra_bias : bias weights'
                   WRITE(numout,*) '~~~~~~~~~~~~'
-                  WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(real(kt)/ztsday)
+                  WRITE(numout,* ) "proportion of historical pcbias applied ",t_asm_mem**(REAL(kt,wp)/ztsday)
                ENDIF
 
+               IF (ln_incpc) THEN
+                   IF (lwp) WRITE(numout,*) 'obsoper or forecast mode - IPC - computing tbias_i and sbias_i'
+                   DO jk = 1, jpkm1
+                      tbias_i(:,:,jk) = tbias_i(:,:,jk) + ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
+                      sbias_i(:,:,jk) = sbias_i(:,:,jk) + ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
+                   ENDDO
+                   tbias_i(:,:,:) =tbias_i(:,:,:)*tmask(:,:,:)
+                   sbias_i(:,:,:) =sbias_i(:,:,:)*tmask(:,:,:)
+               ENDIF
             ENDIF
  
@@ -731,6 +888,10 @@
 
       IF ( kt == nn_bias_itwrt ) THEN
-         tbias_p_out(:,:,:)   = tbias_p(:,:,:)
-         sbias_p_out(:,:,:)   = sbias_p(:,:,:)
+            tbias_p_out(:,:,:)   = tbias_p(:,:,:)
+            sbias_p_out(:,:,:)   = sbias_p(:,:,:)
+            IF (ln_incpc) THEN
+               tbias_i_out(:,:,:)   = tbias_i(:,:,:)
+               sbias_i_out(:,:,:)   = sbias_i(:,:,:)            
+            ENDIF
       ENDIF
 
@@ -774,16 +935,36 @@
          tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:)
          tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:)
+         IF ( ln_incpc ) THEN
+            tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:)
+            tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:)
+            IF ( ln_incpc_only ) THEN
+               IF(lwp) WRITE(numout,*) 'ln_incpc_only =', ln_incpc_only, 'tsw updated with IPC only and ln_dynhpg_imp = ',ln_dynhpg_imp
+               tsw(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - tbias_i(:,:,:)
+               tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - sbias_i(:,:,:)
+            ENDIF
+         ENDIF
       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(:,:,:))
+         IF ( ln_incpc ) THEN
+            tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:) - tbias_i(:,:,:)
+            tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:) - sbias_i(:,:,:)
+            IF ( ln_incpc_only ) THEN
+               IF(lwp) WRITE(numout,*) 'ln_incpc_only =', ln_incpc_only, 'tsw updated with IPC only and ln_dynhpg_imp = ',ln_dynhpg_imp
+               tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_i(:,:,:)
+               tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_i(:,:,:)
+            ENDIF
+         ENDIF
+      ENDIF
+
+      IF(lwp) WRITE(numout,*) 'dyn_bias( kt ) calculating rhd_pc, kt =', kt
+      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,  &
@@ -856,4 +1037,11 @@
          CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_p'   , sbias_p_out )         
       ENDIF
+
+      IF ( ln_incpc ) THEN
+         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_asm_stscale' , tbias_asm_stscale_out )   
+         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_asm_stscale' , sbias_asm_stscale_out ) 
+         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'tbias_i'   , tbias_i_out )
+         CALL iom_rstput( nn_bias_itwrt, nn_bias_itwrt, numbias_tot, 'sbias_i'   , sbias_i_out )   
+      ENDIF
       
       IF( kt == nitend ) THEN
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/biaspar.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/biaspar.F90	(revision 8446)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/biaspar.F90	(revision 8447)
@@ -27,6 +27,8 @@
       & tbias, &       !: Temperature bias field for T correction
       & tbias_p, &     !:  "           "    "    "   P correction
+      & tbias_i, &     !:  "           "    "    "   incremental P correction
       & sbias, &       !: Salinity bias field    for S correction
       & sbias_p, &     !:  "           "    "    "   P correction
+      & sbias_i, &     !:  "           "    "    "   incremental P correction
       & rhd_pc         !: Press corrtd density from online to use in dyn_hpg
 
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 8446)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 8447)
@@ -113,5 +113,5 @@
 
       ENDIF
-      !
+
       SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation
       CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate
@@ -134,4 +134,7 @@
       !
       IF ( ln_bias .AND. ln_bias_pc_app )  THEN
+         IF(lwp) THEN 
+         WRITE(numout,*) " ! restore original density"
+         ENDIF
          rhd(:,:,:) = z_rhd_st(:,:,:)     ! restore original density
          gru(:,:)   = z_gru_st(:,:)
Index: /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 8446)
+++ /branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 8447)
@@ -457,4 +457,6 @@
       END SELECT
       !
+      CALL lbc_lnk( prd, 'T', 1.0_wp )
+      !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk )
       !
