Index: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 11939)
+++ branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 11944)
@@ -21,5 +21,6 @@
    !!   dyn_asm_inc    : Apply the dynamic (u and v) increments
    !!   ssh_asm_inc    : Apply the SSH increment
-   !!   seaice_asm_inc : Apply the sea ice increment
+   !!   seaice_asm_inc : Apply the sea ice concentration increment
+   !!   sit_asm_inc    : Apply the sea ice thickness increment
    !!----------------------------------------------------------------------
    USE wrk_nemo         ! Memory Allocation
@@ -49,5 +50,7 @@
    PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments
    PUBLIC   ssh_asm_inc    !: Apply the SSH increment
-   PUBLIC   seaice_asm_inc !: Apply the seaice increment
+   PUBLIC   seaice_asm_inc !: Apply the seaice concentration increment
+   PUBLIC   sit_asm_inc    !: Apply the seaice thickness increment
+   PUBLIC   bgc_asm_inc    !: Apply the biogeochemistry increments
 
 #if defined key_asminc
@@ -63,4 +66,6 @@
    LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment
    LOGICAL, PUBLIC :: ln_seaiceinc = .FALSE.   !: No sea ice concentration increment
+   LOGICAL, PUBLIC :: ln_sitinc = .FALSE.      !: No sea ice thickness increment
+   LOGICAL, PUBLIC :: lk_bgcinc = .FALSE.      !: No biogeochemistry increments
    LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check
    LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
@@ -85,4 +90,5 @@
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   sit_bkginc            ! Increment to the background sea ice thickness
 
    !! * Substitutions
@@ -127,9 +133,17 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace
       !!
-      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
+      NAMELIST/nam_asminc/ ln_bkgwri, ln_balwri,                           &
          &                 ln_trainc, ln_dyninc, ln_sshinc,                &
+         &                 ln_phytobal, ln_slchltotinc, ln_slchldiainc,    &
+         &                 ln_slchlnoninc, ln_schltotinc, ln_slphytotinc,  &
+         &                 ln_slphydiainc, ln_slphynoninc, ln_spco2inc,    &
+         &                 ln_sfco2inc, ln_plchltotinc, ln_pchltotinc,     &
+         &                 ln_pno3inc, ln_psi4inc, ln_pdicinc, ln_palkinc, &
+         &                 ln_pphinc, ln_po2inc,                           &
          &                 ln_asmdin, ln_asmiau,                           &
          &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
-         &                 ln_salfix, salfixmin, nn_divdmp
+         &                 ln_salfix, salfixmin, nn_divdmp,                &
+         &                 ln_seaiceinc, ln_sitinc, ln_temnofreeze,        &
+         &                 mld_choice_bgc, rn_maxchlinc
       !!----------------------------------------------------------------------
 
@@ -138,4 +152,5 @@
       !-----------------------------------------------------------------------
       ln_seaiceinc = .FALSE.
+      ln_sitinc = .FALSE.
       ln_temnofreeze = .FALSE.
 
@@ -154,11 +169,31 @@
          WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters'
+         WRITE(numout,*) '   Namelist nam_asminc : set assimilation increment parameters'
          WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
+         WRITE(numout,*) '      Logical switch for writing out balancing increments      ln_balwri = ', ln_balwri
          WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
          WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
          WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc
          WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin
-         WRITE(numout,*) '      Logical switch for applying sea ice increments           ln_seaiceinc = ', ln_seaiceinc
+         WRITE(numout,*) '      Logical switch for applying SIC increments               ln_seaiceinc = ', ln_seaiceinc
+         WRITE(numout,*) '      Logical switch for applying SIT increments               ln_sitinc = ', ln_sitinc
+         WRITE(numout,*) '      Logical switch for phytoplankton balancing             ln_phytobal = ', ln_phytobal
+         WRITE(numout,*) '      Logical switch for applying slchltot increments     ln_slchltotinc = ', ln_slchltotinc
+         WRITE(numout,*) '      Logical switch for applying slchldia increments     ln_slchldiainc = ', ln_slchldiainc
+         WRITE(numout,*) '      Logical switch for applying slchlnon increments     ln_slchlnoninc = ', ln_slchlnoninc
+         WRITE(numout,*) '      Logical switch for applying schltot increments       ln_schltotinc = ', ln_schltotinc
+         WRITE(numout,*) '      Logical switch for applying slphytot increments     ln_slphytotinc = ', ln_slphytotinc
+         WRITE(numout,*) '      Logical switch for applying slphydia increments     ln_slphydiainc = ', ln_slphydiainc
+         WRITE(numout,*) '      Logical switch for applying slphynon increments     ln_slphynoninc = ', ln_slphynoninc
+         WRITE(numout,*) '      Logical switch for applying spco2 increments           ln_spco2inc = ', ln_spco2inc
+         WRITE(numout,*) '      Logical switch for applying sfco2 increments           ln_sfco2inc = ', ln_sfco2inc
+         WRITE(numout,*) '      Logical switch for applying plchltot increments     ln_plchltotinc = ', ln_plchltotinc
+         WRITE(numout,*) '      Logical switch for applying pchltot increments       ln_pchltotinc = ', ln_pchltotinc
+         WRITE(numout,*) '      Logical switch for applying pno3 increments             ln_pno3inc = ', ln_pno3inc
+         WRITE(numout,*) '      Logical switch for applying psi4 increments             ln_psi4inc = ', ln_psi4inc
+         WRITE(numout,*) '      Logical switch for applying pdic increments             ln_pdicinc = ', ln_pdicinc
+         WRITE(numout,*) '      Logical switch for applying palk increments             ln_palkinc = ', ln_palkinc
+         WRITE(numout,*) '      Logical switch for applying pph increments               ln_pphinc = ', ln_pphinc
+         WRITE(numout,*) '      Logical switch for applying po2 increments               ln_po2inc = ', ln_po2inc
          WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
          WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
@@ -169,4 +204,6 @@
          WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
          WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
+         WRITE(numout,*) '      Choice of MLD for BGC assimilation                  mld_choice_bgc = ', mld_choice_bgc
+         WRITE(numout,*) '      Maximum absolute chlorophyll increment (<=0 = off)    rn_maxchlinc = ', rn_maxchlinc
       ENDIF
 
@@ -216,7 +253,8 @@
 
       IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
-           .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc) )) &
-         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc,', &
-         &                ' and ln_seaiceinc is set to .true.', &
+         & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. &
+         &        ( ln_sitinc ).OR.( lk_bgcinc ) )) &
+         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' ln_sitinc and ln_(bgc-variable)inc is set to .true.', &
          &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
          &                ' Inconsistent options')
@@ -227,7 +265,7 @@
 
       IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
-         &                     )  &
-         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc', &
-         &                ' and ln_seaiceinc are set to .false. :', &
+         & .AND.( .NOT. ln_sitinc ).AND.( .NOT. lk_bgcinc ) )  &
+         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' ln_sitinc and ln_(bgc-variable)inc are set to .false. :', &
          &                ' The assimilation increments are not applied')
 
@@ -333,4 +371,5 @@
       ALLOCATE( ssh_bkginc(jpi,jpj)   )
       ALLOCATE( seaice_bkginc(jpi,jpj))
+      ALLOCATE( sit_bkginc(jpi,jpj)   )
 #if defined key_asminc
       ALLOCATE( ssh_iau(jpi,jpj)      )
@@ -342,8 +381,10 @@
       ssh_bkginc(:,:) = 0.0
       seaice_bkginc(:,:) = 0.0
+      sit_bkginc(:,:) = 0.0
 #if defined key_asminc
       ssh_iau(:,:)    = 0.0
 #endif
-      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN
+      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) &
+         &  .OR.( ln_sitinc ).OR.( lk_bgcinc ) ) THEN
 
          !--------------------------------------------------------------------
@@ -408,4 +449,13 @@
             ! to allow for differences in masks
             WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
+         ENDIF
+
+         IF ( ln_sitinc ) THEN
+            CALL iom_get( inum, jpdom_autoglo, 'bckinsit', sit_bkginc, 1 )
+            ! Apply the masks
+            sit_bkginc(:,:) = sit_bkginc(:,:) * tmask(:,:,1)
+            ! Set missing increments to 0.0 rather than 1e+20
+            ! to allow for differences in masks
+            WHERE( ABS( sit_bkginc(:,:) ) > 1.0e+10 ) sit_bkginc(:,:) = 0.0
          ENDIF
 
@@ -768,4 +818,5 @@
       ! Perhaps the following call should be in step
       IF   ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
+      IF   ( ln_sitinc  )      CALL sit_asm_inc ( kt )      ! apply sea ice thickness increment
       !
    END SUBROUTINE tra_asm_inc
@@ -925,9 +976,231 @@
    END SUBROUTINE ssh_asm_inc
 
+   SUBROUTINE sit_asm_inc( kt, kindic )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE sit_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the sea ice thickness assimilation increment.
+      !!
+      !! ** Method  : Direct initialization or Incremental Analysis Updating.
+      !!
+      !! ** Action  : 
+      !!
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      !
+      INTEGER, INTENT(in)           ::   kt   ! Current time step
+      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
+      !
+      INTEGER  ::   it
+      REAL(wp) ::   zincwgt   ! IAU weight for current time step
+! #if defined key_lim2
+!       REAL(wp), DIMENSION(jpi,jpj) ::   zofrld, zohicif, zseaicendg, zhicifinc  ! LIM
+!       REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres
+!       !!THICKNESS INCS NOT SET UP FOR LIM
+! #endif
+      !!----------------------------------------------------------------------
+
+      IF ( ln_asmiau ) THEN
+
+         !--------------------------------------------------------------------
+         ! Incremental Analysis Updating
+         !--------------------------------------------------------------------
+
+         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
+
+            it = kt - nit000 + 1
+            zincwgt = wgtiau(it)      ! IAU weight for the current time step 
+            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
+            ! EF: Actually CICE is expecting a tendency so is divided by rdt below
+
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               WRITE(numout,*) 'sit_asm_inc : sea ice thick IAU at time step = ', &
+                  &  kt,' with IAU weight = ', wgtiau(it)
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Sea-ice : LIM-3 case (to add)
+
+! #if defined key_lim2
+!             ! Sea-ice : LIM-2 case (to add if needed)
+!             zofrld (:,:) = frld(:,:)
+!             zohicif(:,:) = hicif(:,:)
+!             !
+!             frld  = MIN( MAX( frld (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
+!             pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
+!             fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction
+!             !
+!             zseaicendg(:,:) = zofrld(:,:) - frld(:,:)   ! find out actual sea ice nudge applied
+!             !
+!             ! Nudge sea ice depth to bring it up to a required minimum depth
+!             WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
+!                zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt    
+!             ELSEWHERE
+!                zhicifinc(:,:) = 0.0_wp
+!             END WHERE
+!             !
+!             ! nudge ice depth
+!             hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
+!             phicif(:,:) = phicif(:,:) + zhicifinc(:,:)       
+!             !
+!             ! seaice salinity balancing (to add)
+! #endif
+
+#if defined key_cice && defined key_asminc
+            ! Sea-ice thickness : CICE case. Pass ice thickness increment tendency into CICE
+            ndsit_da(:,:) = sit_bkginc(:,:) * zincwgt / rdt
+#endif
+
+            IF ( kt == nitiaufin_r ) THEN
+               DEALLOCATE( sit_bkginc )
+            ENDIF
+
+         ELSE
+
+#if defined key_cice && defined key_asminc
+            ! Sea-ice thickness : CICE case. Zero ice increment tendency into CICE
+            ndsit_da(:,:) = 0.0_wp
+#endif
+
+         ENDIF
+
+      ELSEIF ( ln_asmdin ) THEN
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+            ! Sea-ice : LIM-3 case (to add)
+
+! #if defined key_lim2
+!             ! Sea-ice : LIM-2 case (add if needed)
+!             zofrld(:,:)=frld(:,:)
+!             zohicif(:,:)=hicif(:,:)
+!             ! 
+!             ! Initialize the now fields the background + increment
+!             frld (:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
+!             pfrld(:,:) = frld(:,:) 
+!             fr_i (:,:) = 1.0_wp - frld(:,:)                ! adjust ice fraction
+!             zseaicendg(:,:) = zofrld(:,:) - frld(:,:)      ! find out actual sea ice nudge applied
+!             !
+!             ! Nudge sea ice depth to bring it up to a required minimum depth
+!             WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
+!                zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt    
+!             ELSEWHERE
+!                zhicifinc(:,:) = 0.0_wp
+!             END WHERE
+!             !
+!             ! nudge ice depth
+!             hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
+!             phicif(:,:) = phicif(:,:)       
+!             !
+!             ! seaice salinity balancing (to add)
+! #endif
+ 
+#if defined key_cice && defined key_asminc
+            ! Sea-ice thickness : CICE case. Pass ice thickness increment tendency into CICE
+           ndsit_da(:,:) = sit_bkginc(:,:) / rdt
+#endif
+           IF ( .NOT. PRESENT(kindic) ) THEN
+              DEALLOCATE( sit_bkginc )
+           END IF
+
+         ELSE
+
+#if defined key_cice && defined key_asminc
+            ! Sea-ice thicnkness : CICE case. Zero ice thickness increment tendency into CICE 
+            ndsit_da(:,:) = 0.0_wp
+#endif
+         
+         ENDIF
+
+!#if defined defined key_lim2 || defined key_cice
+!
+!            IF (ln_seaicebal ) THEN       
+!             !! balancing salinity increments
+!             !! simple case from limflx.F90 (doesn't include a mass flux)
+!             !! assumption is that as ice concentration is reduced or increased
+!             !! the snow and ice depths remain constant
+!             !! note that snow is being created where ice concentration is being increased
+!             !! - could be more sophisticated and
+!             !! not do this (but would need to alter h_snow)
+!
+!             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store
+!
+!             DO jj = 1, jpj
+!               DO ji = 1, jpi 
+!           ! calculate change in ice and snow mass per unit area
+!           ! positive values imply adding salt to the ocean (results from ice formation)
+!           ! fwf : ice formation and melting
+!
+!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt
+!
+!           ! change salinity down to mixed layer depth
+!                 mld=hmld_kara(ji,jj)
+!
+!           ! prevent small mld
+!           ! less than 10m can cause salinity instability 
+!                 IF (mld < 10) mld=10
+!
+!           ! set to bottom of a level 
+!                 DO jk = jpk-1, 2, -1
+!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN 
+!                     mld=gdepw(ji,jj,jk+1)
+!                     jkmax=jk
+!                   ENDIF
+!                 ENDDO
+!
+!            ! avoid applying salinity balancing in shallow water or on land
+!            ! 
+!
+!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
+!
+!                 dsal_ocn=0.0_wp
+!                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing
+!
+!                 if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
+!                              dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
+!
+!           ! put increments in for levels in the mixed layer
+!           ! but prevent salinity below a threshold value 
+!
+!                   DO jk = 1, jkmax              
+!
+!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 
+!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
+!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
+!                     ENDIF
+!
+!                   ENDDO
+!
+!      !            !  salt exchanges at the ice/ocean interface
+!      !            zpmess         = zfons / rdt_ice    ! rdt_ice is ice timestep
+!      !
+!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean
+!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt  
+!      !!               
+!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d) 
+!      !!                                                     ! E-P (kg m-2 s-2)
+!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2)
+!               ENDDO !ji
+!             ENDDO !jj!
+!
+!            ENDIF !ln_seaicebal
+!
+!#endif
+
+
+   END SUBROUTINE sit_asm_inc
+
    SUBROUTINE seaice_asm_inc( kt, kindic )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE seaice_asm_inc  ***
       !!          
-      !! ** Purpose : Apply the sea ice assimilation increment.
+      !! ** Purpose : Apply the sea ice concentration assimilation increment.
       !!
       !! ** Method  : Direct initialization or Incremental Analysis Updating.
