Index: /branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_logchlbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90
===================================================================
--- /branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_logchlbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90	(revision 8437)
+++ /branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_logchlbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90	(revision 8438)
@@ -22,4 +22,5 @@
    USE par_trc,     ONLY: jptra          ! Tracer parameters
 
+
    IMPLICIT NONE
    PRIVATE                   
@@ -27,49 +28,39 @@
    PUBLIC asm_logchl_bal_ersem
 
+   REAL(wp), PARAMETER :: tiny = 1e-20
+
 CONTAINS
-
-   SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, &
-      &                             k_maxchlinc, logchl_bkginc, logchl_balinc )
-      !!---------------------------------------------------------------------------
-      !!                    ***  ROUTINE asm_logchl_bal_ersem  ***
+    SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, &
+      &                             k_maxchlinc, logchl_bkginc, logchl_balinc )                     
+
+                                                                        
+                                         
+ !! INTEGER, INTENT(in   ) :: npfts
+ !!REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: inc_log_chlTot
+
+      LOGICAL,  INTENT(in   )                           :: ld_logchlpftinc
+      INTEGER,  INTENT(in   )                           :: npfts
+      INTEGER,  INTENT(in   )                           :: mld_choice_bgc
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,npfts) :: logchl_bkginc
+      REAL(wp), INTENT(in   )                               :: k_maxchlinc
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra)   :: logchl_balinc
+
       !!
-      !! ** Purpose :   calculate increments to ERSEM from logchl increments
-      !!
-      !! ** Method  :   convert logchl increments to chl increments
-      !!                split between the ERSEM PFTs
-      !!                spread through the mixed layer
-      !!                [forthcoming: calculate increments to nutrients and zooplankton]
-      !!
-      !! ** Action  :   populate logchl_balinc
-      !!
-      !! References :   forthcoming...
-      !!---------------------------------------------------------------------------
-      !!
-      LOGICAL,  INTENT(in   )                               :: ld_logchlpftinc
-      INTEGER,  INTENT(in   )                               :: npfts
-      INTEGER,  INTENT(in   )                               :: mld_choice_bgc
-      REAL(wp), INTENT(in   )                               :: k_maxchlinc
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,npfts)     :: logchl_bkginc
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc
-      !!
-      INTEGER                      :: ji, jj, jk
-      INTEGER                      :: jkmax
-      REAL(wp)                     :: chl_tot, chl_inc
-      REAL(wp), DIMENSION(jpi,jpj) :: zmld
-      !!---------------------------------------------------------------------------
-
-      ! Split surface logchl incs into surface Chl1-4 incs
-      !
-      ! In order to transform logchl incs to chl incs, need to account for the background,
-      ! cannot simply do 10^logchl_bkginc. Need to:
-      ! 1) Add logchl inc to log10(background) to get log10(analysis)
-      ! 2) Take 10^log10(analysis) to get analysis
-      ! 3) Subtract background from analysis to get chl incs
-      ! If k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
-      !
-      ! Only apply increments if all of Chl1-4 background values are > 0
-      ! In theory, they always will be, and if any are not that's a sign
-      ! that something's going wrong which the assimilation might make worse
-      !
+ 
+ INTEGER                      :: ji, jj, jk, jp
+ INTEGER                      :: jkmax
+ REAL(wp)                     :: chl_tot, chl_inc
+ REAL(wp), DIMENSION(jpi,jpj) :: zmld
+ REAL(wp), DIMENSION(jpi,jpj,npfts,jptra) :: cc_chl_field  !! CROSS-COV BETWEEN SPECIFIC SIZE CLASS AND GIVEN BIOGEOCHEM VAR
+ REAL(wp), DIMENSION(jptra) :: threshold
+ 
+
+
+         ! STEP 1: TRANSFORM ASSIMILATED LOG-FIELD INCREMENTS INTO INCREMENTS AND PRODUCE THEM FOR THE SIZE-CLASSES IF NEEDED, ALSO SAVE INCREMENTS OF LOGS IN SEPARATE VARIABLES (THIS IS IN CASE ONLY TOTAL CHLOROPHYLL WAS ASSIMILATED) 
+
+  IF (lwp) WRITE(numout,*) 'Code has gotten to point started'
+  CALL flush(numout)
+
+
       IF ( ld_logchlpftinc ) THEN
          !
@@ -151,7 +142,85 @@
       ENDIF
 
-      ! Propagate surface Chl1-4 incs through mixed layer
-      ! First, choose mixed layer definition
-      !
+      
+      
+   
+
+
+
+  IF (lwp) WRITE(numout,*) 'Code has gotten to point log_inc transformed to inc'
+  CALL flush(numout)
+   ! STEP 2: READ CROSS_COVARIANCES
+
+  CALL read_crosscovariances(cc_chl_field, npfts)
+
+ 
+ IF (lwp) WRITE(numout,*) 'Code has gotten to point cc read'
+ CALL flush(numout)
+  !! STEP 3: CALCULATE INCREMENTS OF OTHER BIOGEOCHEMICAL VARIABLES 
+
+
+  threshold(jp_fabm_m1+jp_fabm_n3n) = 100
+  threshold(jp_fabm_m1+jp_fabm_n4n) = 10
+  threshold(jp_fabm_m1+jp_fabm_n1p) = 15
+  threshold(jp_fabm_m1+jp_fabm_n5s) = 80
+  threshold(jp_fabm_m1+jp_fabm_z4c) = 20
+  threshold(jp_fabm_m1+jp_fabm_z5c) = 20
+  threshold(jp_fabm_m1+jp_fabm_p1c) = 25
+  threshold(jp_fabm_m1+jp_fabm_p2c) = 25
+  threshold(jp_fabm_m1+jp_fabm_p3c) = 25
+  threshold(jp_fabm_m1+jp_fabm_p4c) = 25
+  threshold(jp_fabm_m1+jp_fabm_p1n) = 10
+  threshold(jp_fabm_m1+jp_fabm_p2n) = 10
+  threshold(jp_fabm_m1+jp_fabm_p3n) = 10
+  threshold(jp_fabm_m1+jp_fabm_p4n) = 10
+  threshold(jp_fabm_m1+jp_fabm_p1p) = 10
+  threshold(jp_fabm_m1+jp_fabm_p2p) = 10
+  threshold(jp_fabm_m1+jp_fabm_p3p) = 10
+  threshold(jp_fabm_m1+jp_fabm_p4p) = 10
+
+ 
+
+ 
+  DO jj = 1, jpj
+         DO ji = 1, jpi
+                 DO jp=1,jptra
+
+                  
+            !CALL t;ge'e;fdgfgd
+
+                         IF ( (jp /= jp_fabm_m1+jp_fabm_chl1) .AND. &
+                            & (jp /= jp_fabm_m1+jp_fabm_chl2) .AND. &
+                            & (jp /= jp_fabm_m1+jp_fabm_chl3) .AND. &
+                            & (jp /= jp_fabm_m1+jp_fabm_chl4)) THEN 
+            
+                         CALL asm_inc_biogeo_ersem(npfts, jp, trn(ji,jj,1,jp), logchl_bkginc(ji,jj,:),  cc_chl_field(ji,jj,:,:), threshold(jp), &
+                          & logchl_balinc(ji,jj,1,jp))
+
+
+                         END IF
+
+
+         END DO
+
+      END DO
+
+    END DO
+
+
+
+ IF (lwp) WRITE(numout,*) 'Code has gotten to point surface logs calculated'
+
+
+
+
+ 
+
+
+
+
+ CALL flush(numout)
+       !! STEP 4: PROPAGATE EACH INCREMENT THROUGH THE MIXED LAYER
+
+
       SELECT CASE( mld_choice_bgc )
       CASE ( 1 )                   ! Turbocline/mixing depth [W points]
@@ -179,5 +248,6 @@
       ! Now set MLD to bottom of a level and propagate incs equally through mixed layer
       !
-      DO jj = 1, jpj
+
+       DO jj = 1, jpj
          DO ji = 1, jpi
             !
@@ -192,16 +262,150 @@
             !
             DO jk = 2, jkmax
-               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
-               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
-               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
-               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
+
+               DO jp = 1, jptra
+
+               !IF (lwp) WRITE(numout,*) logchl_balinc_ersem_chl1(ji,jj,1)
+               !CALL flush(numout)
+               logchl_balinc(ji,jj,jk,jp) = logchl_balinc(ji,jj,1,jp)
+             
+               END DO
             END DO
-            !
+
+  
          END DO
       END DO
 
-      ! Multivariate balancing forthcoming...
-
-   END SUBROUTINE asm_logchl_bal_ersem
+
+   IF (lwp) WRITE(numout,*) 'Code has gotten to point log transformed into inc and propagated'
+   CALL flush(numout)
+
+END SUBROUTINE asm_logchl_bal_ersem
+
+
+
+!CONTAINS 
+SUBROUTINE read_crosscovariances(cc_chl_field, npfts)
+
+  !!     status = nf90_open(path = filepath, cmode = nf90_nowrite, ncid = ccfile)
+
+ INTEGER, INTENT(in  ) :: npfts
+ REAL(wp), INTENT( out), DIMENSION(jpi,jpj,npfts,jptra) :: cc_chl_field  !! CROSS-COV BETWEEN SPECIFIC SIZE CLASS AND GIVEN BIOGEOCHEM VAR
+ INTEGER :: ncid
+ INTEGER :: pft
+ CHARACTER(LEN=256)  :: path_chl(npfts)
+ CHARACTER(LEN=256)  :: id_chl(npfts)
+
+ IF (npfts /= 4) THEN
+
+ IF (lwp) WRITE(numout,*) 'Wrong number of Pfts for reading crosscovs files'
+  CALL flush(numout)
+
+ ELSE
+
+ path_chl(1) = "crosscovs_chl1.nc"
+ path_chl(2) = "crosscovs_chl2.nc"
+ path_chl(3) = "crosscovs_chl3.nc"
+ path_chl(4) = "crosscovs_chl4.nc"
+
+ id_chl(1) = "cc_Chl1"
+ id_chl(2) = "cc_Chl2"
+ id_chl(3) = "cc_Chl3"
+ id_chl(4) = "cc_Chl4"
+
+
+ DO pft = 1, npfts
+
+ CALL iom_open( path_chl(pft), ncid )
+ CALL iom_get( ncid, jpdom_autoglo, 'cc_NO3', cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n3n))
+ CALL iom_get( ncid, jpdom_autoglo, 'cc_NH4', cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n4n))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_PO4", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n1p))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_SiO", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n5s))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_Z4c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_z4c))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_Z5c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_z5c))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P1c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p1c))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P2c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p2c))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P3c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p3c))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P4c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p4c))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P1n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p1n))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P2n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p2n))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P3n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p3n))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P4n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p4n))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P1p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p1p))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P2p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p2p))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P3p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p3p))
+ CALL iom_get( ncid, jpdom_autoglo, "cc_P4p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p4p))
+ CALL iom_get( ncid, jpdom_autoglo, id_chl(pft), cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_chl1))
+ CALL iom_close( ncid )
+
+ END DO
+
+ END IF
+
+
+  END SUBROUTINE read_crosscovariances
+
+
+!CONTAINS
+ SUBROUTINE asm_inc_biogeo_ersem(npfts, jp, bgf_value, chl, cc, threshold, output)
+ INTEGER, INTENT(in ) :: npfts
+ INTEGER, INTENT(in ) :: jp
+ REAL, INTENT(in  )  :: bgf_value
+ REAL, INTENT(in   ), DIMENSION(npfts) :: chl
+ REAL, INTENT(in   ), DIMENSION(npfts,jptra) :: cc  !! VARIANCE OF CHL SIZE-CLASS 
+ REAL, INTENT(in   ) :: threshold
+ REAL, INTENT(  out) :: output
+
+
+ REAL :: total, weight_1, weight_2, weight_3, weight_4
+
+
+    output = 0.0
+
+    IF ( cc(1,jp)*cc(2,jp)*cc(3,jp)*cc(4,jp)*cc(1, jp_fabm_m1+jp_fabm_chl1)*cc(2, jp_fabm_m1+jp_fabm_chl2)*cc(3, jp_fabm_m1+jp_fabm_chl3)*cc(4, jp_fabm_m1+jp_fabm_chl4) /= 0.0  ) THEN
+
+       total = ABS(cc(1,jp)/cc(1, jp_fabm_m1+jp_fabm_chl1)) + ABS(cc(2,jp)/cc(2, jp_fabm_m1+jp_fabm_chl2)) + ABS(cc(3,jp)/cc(3, jp_fabm_m1+jp_fabm_chl3)) + ABS(cc(4,jp)/cc(4, jp_fabm_m1+jp_fabm_chl4))
+       weight_1 = ABS(cc(1,jp)/cc(1, jp_fabm_m1+jp_fabm_chl1))/total
+       weight_2 = ABS(cc(2,jp)/cc(2, jp_fabm_m1+jp_fabm_chl2))/total
+       weight_3 = ABS(cc(3,jp)/cc(3, jp_fabm_m1+jp_fabm_chl3))/total
+       weight_4 = ABS(cc(4,jp)/cc(4, jp_fabm_m1+jp_fabm_chl4))/total
+
+!       PRINT *, 'weights', weight_1, weight_2, weight_3, weight_4 
+      
+       output = weight_1*(cc(1,jp)/cc(1, jp_fabm_m1+jp_fabm_chl1))*chl(1) + weight_2*(cc(2,jp)/cc(2, jp_fabm_m1+jp_fabm_chl2))*chl(2) + weight_3*(cc(3,jp)/cc(3, jp_fabm_m1+jp_fabm_chl3))*chl(3) + weight_4*(cc(4,jp)/cc(4, jp_fabm_m1+jp_fabm_chl4))*chl(4)
+
+
+       ! FIRST BREAK = DON'T INCREASE VALUE MORE THAN TENFOLD 
+       
+       output = MIN(output, LOG10(11.0))
+      
+
+       ! SECOND BREAK = IF VALUE LARGER THAN THRESHOLD DO NOT INCREASE IT
+
+       IF (bgf_value > threshold) THEN
+           output = MIN(output, 0.0)
+       ENDIF
+
+        ! THIRD BREAK = IF VALUE MORE THAN FIFTH OF THRESHOLD THAN DON'T INCREASE IT BY MORE THAN AN INCREMENT EQUAL TO TENTH OF THRESHOLD (INCREASE IT SLOWLY)
+
+       IF (bgf_value > 0.2*threshold) THEN
+           output = MIN(output, LOG10(0.1*threshold/bgf_value + 1))
+       ENDIF
+
+       output = (10**( output ) - 1)*bgf_value
+  
+       PRINT *, 'output', output
+         
+      
+    ENDIF
+
+      
+     !! calculate increments for one or all variables
+
+   END SUBROUTINE asm_inc_biogeo_ersem
+
+
+
+  
+
 
 #else
@@ -210,4 +414,5 @@
    !!----------------------------------------------------------------------
 CONTAINS
+   CONTAINS
    SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, &
       &                             k_maxchlinc, logchl_bkginc, logchl_balinc )
@@ -220,6 +425,12 @@
       WRITE(*,*) 'asm_logchl_bal_ersem: You should not have seen this print! error?'
    END SUBROUTINE asm_logchl_bal_ersem
+
+
+
+
+
+  
+
 #endif
-
    !!======================================================================
 END MODULE asmlogchlbal_ersem
