Index: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 8406)
+++ branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 8428)
@@ -22,4 +22,5 @@
    !!   ssh_asm_inc    : Apply the SSH increment
    !!   seaice_asm_inc : Apply the seaice increment
+   !!   logchl_asm_inc : Apply the logchl increment
    !!----------------------------------------------------------------------
    USE wrk_nemo         ! Memory Allocation
@@ -44,4 +45,29 @@
 #endif
    USE sbc_oce          ! Surface boundary condition variables.
+   USE zdfmxl, ONLY :  &  
+!   &  hmld_tref,       &   
+#if defined key_karaml
+   &  hmld_kara,       &
+   &  ln_kara,         &
+#endif   
+   &  hmld,            & 
+   &  hmlp,            &
+   &  hmlpt 
+#if defined key_top
+   USE trc, ONLY: &
+      & trn,      &
+      & trb
+   USE par_trc, ONLY: &
+      & jptra
+#endif
+#if defined key_medusa && defined key_foam_medusa
+   USE asmlogchlbal_medusa, ONLY: &
+      & asm_logchl_bal_medusa
+   USE par_medusa
+#elif defined key_hadocc
+   USE asmlogchlbal_hadocc, ONLY: &
+      & asm_logchl_bal_hadocc
+   USE par_hadocc
+#endif
 
    IMPLICIT NONE
@@ -54,4 +80,5 @@
    PUBLIC   ssh_asm_inc    !: Apply the SSH increment
    PUBLIC   seaice_asm_inc !: Apply the seaice increment
+   PUBLIC   logchl_asm_inc !: Apply the logchl increment
 
 #if defined key_asminc
@@ -61,4 +88,5 @@
 #endif
    LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields
+   LOGICAL, PUBLIC :: ln_balwri = .FALSE.      !: No output of the assimilation balancing increments
    LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment
    LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization
@@ -67,4 +95,6 @@
    LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment
    LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment
+   LOGICAL, PUBLIC :: ln_logchlinc = .FALSE.   !: No log10(chlorophyll) increment
+   LOGICAL, PUBLIC :: ln_logchlbal = .FALSE.   !: Don't apply multivariate log10(chlorophyll) balancing
    LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check
    LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
@@ -91,4 +121,19 @@
    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), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: logchl_bkginc !: Increment to background logchl
+#if defined key_top
+   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: logchl_balinc  !: Increment to BGC variables from logchl assim
+#endif
+   REAL(wp) :: rn_maxchlinc = -999.0  !: maximum absolute non-log chlorophyll increment from logchl assimilation
+                                      !: <= 0 implies no maximum applied (switch turned off)
+                                      !:  > 0 implies maximum absolute chl increment capped at this value
+
+   INTEGER :: mld_choice_bgc    = 5   !: choice of mld criteria to use for biogeochemistry assimilation
+                                      !: 1) hmld      - Turbocline/mixing depth                           [W points]
+                                      !: 2) hmlp      - Density criterion (0.01 kg/m^3 change from 10m)   [W points]
+                                      !: 3) hmld_kara - Kara MLD                                          [Interpolated]
+                                      !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points]
+                                      !: 5) hmlpt     - Density criterion (0.01 kg/m^3 change from 10m)   [T points]
 
    !! * Substitutions
@@ -133,10 +178,12 @@
       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_logchlinc, ln_logchlbal,                     &
          &                 ln_asmdin, ln_asmiau,                           &
          &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
          &                 ln_salfix, salfixmin, nn_divdmp,                &
-         &                 ln_seaiceinc, ln_temnofreeze
+         &                 ln_seaiceinc, ln_temnofreeze,                   &
+         &                 mld_choice_bgc, rn_maxchlinc
       !!----------------------------------------------------------------------
 
@@ -163,4 +210,5 @@
          WRITE(numout,*) '   Namelist namasm : 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
@@ -168,4 +216,6 @@
          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 logchl increments         ln_logchlinc = ', ln_logchlinc
+         WRITE(numout,*) '      Logical switch for logchl multivariate balancing      ln_logchlbal = ', ln_logchlbal
          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
@@ -176,4 +226,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
 
@@ -223,6 +275,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_logchlinc ) )) &
+         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' and ln_logchlinc is set to .true.', &
          &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
          &                ' Inconsistent options')
@@ -233,6 +287,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_logchlinc ) )  &
+         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
+         &                ' and ln_logchlinc are set to .false. :', &
          &                ' The assimilation increments are not applied')
 
@@ -254,4 +309,16 @@
          &                ' Background time step for Direct Initialization is outside', &
          &                ' the cycle interval')
+
+      IF ( ( ln_balwri ).AND.( .NOT. ( ln_logchlinc ) ) ) THEN
+         CALL ctl_warn( ' Balancing increments are only calculated for logchl', &
+            &           ' Not assimilating logchl, so ln_balwri will be set to .false.')
+         ln_balwri = .FALSE.
+      ENDIF
+
+      IF ( ( ln_logchlbal ).AND.( .NOT. ( ln_logchlinc ) ) ) THEN
+         CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', &
+            &           ' Not assimilating logchl, so ln_logchlbal will be set to .false.')
+         ln_logchlbal = .FALSE.
+      ENDIF
 
       IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
@@ -350,5 +417,14 @@
       ssh_iau(:,:)    = 0.0
 #endif
-      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN
+      IF ( ln_logchlinc ) THEN
+         ALLOCATE( logchl_bkginc(jpi,jpj) )
+         logchl_bkginc(:,:) = 0.0
+#if defined key_top
+         ALLOCATE( logchl_balinc(jpi,jpj,jpk,jptra) )
+         logchl_balinc(:,:,:,:) = 0.0
+#endif
+      ENDIF
+      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) &
+         &  .OR.( ln_logchlinc ) ) THEN
 
          !--------------------------------------------------------------------
@@ -422,4 +498,13 @@
             ! to allow for differences in masks
             WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
+         ENDIF
+
+         IF ( ln_logchlinc ) THEN
+            CALL iom_get( inum, jpdom_autoglo, 'bckinlogchl', logchl_bkginc(:,:), 1 )
+            ! Apply the masks
+            logchl_bkginc(:,:) = logchl_bkginc(:,:) * tmask(:,:,1)
+            ! Set missing increments to 0.0 rather than 1e+20
+            ! to allow for differences in masks
+            WHERE( ABS( logchl_bkginc(:,:) ) > 1.0e+10 ) logchl_bkginc(:,:) = 0.0
          ENDIF
 
@@ -1155,4 +1240,119 @@
 
    END SUBROUTINE seaice_asm_inc
+
+   SUBROUTINE logchl_asm_inc( kt )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE logchl_asm_inc  ***
+      !!          
+      !! ** Purpose : Apply the chlorophyll assimilation increments.
+      !!
+      !! ** Method  : Calculate increments to state variables using nitrogen
+      !!              balancing.
+      !!              Direct initialization or Incremental Analysis Updating.
+      !!
+      !! ** Action  : 
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(IN) :: kt   ! Current time step
+      !
+      INTEGER  :: jk              ! Loop counter
+      INTEGER  :: it              ! Index
+      REAL(wp) :: zincwgt         ! IAU weight for current time step
+      REAL(wp) :: zincper         ! IAU interval in seconds
+      !!----------------------------------------------------------------------
+
+      IF ( kt <= nit000 ) THEN
+
+         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
+
+#if defined key_medusa && defined key_foam_medusa
+         !CALL asm_logchl_bal_medusa()
+         CALL ctl_stop( 'Attempting to assimilate logchl into MEDUSA, ', &
+            &           'but not fully implemented yet' )
+#elif defined key_hadocc
+         CALL asm_logchl_bal_hadocc( logchl_bkginc, zincper, mld_choice_bgc, &
+            &                        rn_maxchlinc, ln_logchlbal, logchl_balinc )
+#else
+         CALL ctl_stop( 'Attempting to assimilate logchl, ', &
+            &           'but not defined a biogeochemical model' )
+#endif
+
+      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
+
+            IF(lwp) THEN
+               WRITE(numout,*) 
+               WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', &
+                  &  kt,' with IAU weight = ', wgtiau(it)
+               WRITE(numout,*) '~~~~~~~~~~~~'
+            ENDIF
+
+            ! Update the biogeochemical variables
+            ! Add directly to trn and trb, rather than to tra, as not a tendency
+#if defined key_medusa && defined key_foam_medusa
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
+                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
+            END DO
+#elif defined key_hadocc
+            DO jk = 1, jpkm1
+               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
+                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
+                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
+            END DO
+#endif
+
+            ! Do not deallocate arrays - needed by asm_bal_wri
+            ! which is called at end of model run
+         ENDIF
+
+      ELSEIF ( ln_asmdin ) THEN 
+
+         !--------------------------------------------------------------------
+         ! Direct Initialization
+         !--------------------------------------------------------------------
+         
+         IF ( kt == nitdin_r ) THEN
+
+            neuler = 0                    ! Force Euler forward step
+
+#if defined key_medusa && defined key_foam_medusa
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
+               &                         logchl_balinc(:,:,:,jp_msa0:jp_msa1)
+            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
+#elif defined key_hadocc
+            ! Initialize the now fields with the background + increment
+            ! Background currently is what the model is initialised with
+            CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', &
+               &           ' Background state is taken from model rather than background file' )
+            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
+               &                         logchl_balinc(:,:,:,jp_had0:jp_had1)
+            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
+#endif
+ 
+            ! Do not deallocate arrays - needed by asm_bal_wri
+            ! which is called at end of model run
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE logchl_asm_inc
    
    !!======================================================================
