Index: /branches/UKMO/dev_r5518_pcbias_ipc/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90
===================================================================
--- /branches/UKMO/dev_r5518_pcbias_ipc/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90	(revision 7080)
+++ /branches/UKMO/dev_r5518_pcbias_ipc/NEMOGCM/NEMO/OPA_SRC/ASM/bias.F90	(revision 7081)
@@ -133,6 +133,6 @@
       & 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
 
    LOGICAL,  PRIVATE :: lalloc
@@ -153,7 +153,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
+   LOGICAL, PRIVATE :: lat_ramp_incpc  ! ramp for incremental pressure correction
+
+   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef
+   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: fbcoef_stscale
 
    INTEGER, PRIVATE ::  &
@@ -214,5 +217,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, lat_ramp_incpc
  
 
@@ -273,4 +276,6 @@
          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,*) '     lat ramp for inc.press.corr    lat_ramp_incpc = ',lat_ramp_incpc
          WRITE(numout,*) '  Parameters for parition of bias term '
          WRITE(numout,*) '                                    fb_t_rlx       = ',fb_t_rlx
@@ -304,9 +309,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),     &
@@ -321,4 +328,8 @@
                                       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))
+
          lalloc = .TRUE.
 
@@ -369,4 +380,11 @@
            sbias_asm(:,:,:) = 0.0_wp
          ENDIF
+
+         IF ( ln_incpc ) THEN   !incr pressure correction
+           tbias_asm_stscale(:,:,:) = 0.0_wp
+           sbias_asm_stscale(:,:,:) = 0.0_wp
+         ENDIF
+
+
          numbias_tot    = 0
          ! Get bias from file and prevent fail if the file does not exist
@@ -402,5 +420,20 @@
                   CALL ctl_stop( 'Bias assim variables not found in ',cn_bias_tot )
                ENDIF
+
+
+               IF ( ln_incpc ) 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  
+
             ENDIF
+
+
 
             ! Close the file
@@ -483,5 +516,20 @@
          fbcoef(:,:) = 0.0_wp
          fctamp      = 0.0_wp
-      ENDIF
+         fbcoef_stscale(:,:) = 0.0_wp
+      ENDIF
+
+
+      IF ( lat_ramp_incpc) THEN
+       ! not sure if this should be a special case of nn_lat_ramp == 2
+         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 
 
 
@@ -516,4 +564,8 @@
       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)            ::   ztscale       ! reduce the inst pcbias by this amount per 24 hours
+      REAL(wp)            ::   zwgt          ! Weight for the inst pcorr term
+      REAL(wp)            ::   zdecay        ! used in inst pcorr
       REAL(wp), DIMENSION(jpi,jpj) :: zcof1, zcof2
 
@@ -527,4 +579,10 @@
       tbias_p(:,:,:) = 0.0_wp
       sbias_p(:,:,:) = 0.0_wp
+
+      ztscale = 0.1_wp
+      zwgt    = 1.0_wp
+      tbias_i(:,:,:) = 0.0_wp
+      sbias_i(:,:,:) = 0.0_wp
+
       IF ( ln_bias_asm ) THEN
       
@@ -539,4 +597,8 @@
 	    ztsday  = rday / real(rdt)
 
+            zdecay = (1-ztscale)**(1/real(ztsday)) ! used in ipc
+            zfrac1 = max(0.0_wp, zdecay**real(kt)) ! used in ipc
+
+
             IF( ln_asmiau .and. ln_trainc ) THEN  ! nudge in increments and decay historical contribution
 
@@ -545,5 +607,5 @@
 
                   ztfrac = real(kt) / real(nitiaufin)  ! nudging factor during the IAU
-
+                  
 	       
                   IF (lwp) THEN
@@ -569,4 +631,20 @@
                   ENDDO
 
+
+                  IF (ln_incpc) THEN
+
+                     DO jk = 1, jpkm1
+
+                        tbias_i(:,:,jk) = tbias_i(:,:,jk) +                         &             ! not sure if it should sum to previous              
+                        &               ( t_bkginc(:,:,jk) * zwgt * ztfrac * (1.0 - fbcoef_stscale(:,:)) )         &
+                        &                - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
+                        sbias_i(:,:,jk) = sbias_i(:,:,jk) +                         &                            
+                        &               ( s_bkginc(:,:,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
@@ -593,5 +671,6 @@
 
                   ztfrac=t_asm_mem**(real(kt-nitiaufin)/real(nitiaufin))  ! decay from end of IAU
-
+                  
+                  
                   DO jk = 1, jpkm1
                      tbias(:,:,jk) = tbias(:,:,jk) +                            &
@@ -603,5 +682,21 @@
                      sbias_p(:,:,jk) = sbias_p(:,:,jk) +                        &  
                      &               ( ztfrac * sbias_asm(:,:,jk) ) * zcof2(:,:)
+
                   ENDDO
+
+                 IF (ln_incpc) THEN
+
+                   zfrac  = max(0.0_wp, zdecay**real(kt-nitiaufin) )
+                   
+                   DO jk = 1, jpkm1
+                      tbias_i(:,:,jk) = tbias_i(:,:,jk) +                         &                ! not sure if it should sum to previous          
+                      &               ( t_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) )         &
+                      &                - ( tbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
+                      sbias_i(:,:,jk) = sbias_i(:,:,jk) +                         &                            ! not sure if it should sum to previous
+                      &               ( s_bkginc(:,:,jk) * zwgt * zfrac * (1.0 - fbcoef_stscale(:,:)) )         &
+                      &                - ( sbias_asm_stscale(:,:,jk) * zfrac1 * (1.0 - fbcoef_stscale(:,:)) )
+                    ENDDO
+
+                 ENDIF
 
                   IF (lwp) THEN
@@ -773,7 +868,19 @@
          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_i(:,:,:) - tbias_p(:,:,:)
+            tsw(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) -sbias_i(:,:,:) - sbias_p(:,:,:)
+         ENDIF
+
       ELSE
          tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_p(:,:,:)
          tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_p(:,:,:)
+
+         IF (ln_incpc) THEN
+            tsw(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) - tbias_i(:,:,:) - tbias_p(:,:,:)
+            tsw(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) - sbias_i(:,:,:) - sbias_p(:,:,:)
+         ENDIF
+
       ENDIF
 
