Index: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_medusa.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_hadocc.F90	(revision 8428)
+++ branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_medusa.F90	(revision 8436)
@@ -1,6 +1,6 @@
-MODULE asmlogchlbal_hadocc
+MODULE asmlogchlbal_medusa
    !!======================================================================
-   !!                       ***  MODULE asmlogchlbal_hadocc  ***
-   !! Calculate increments to HadOCC based on surface logchl increments
+   !!                       ***  MODULE asmlogchlbal_medusa  ***
+   !! Calculate increments to MEDUSA based on surface logchl increments
    !!
    !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
@@ -10,12 +10,13 @@
    !!
    !!======================================================================
-   !! History :  3.6  ! 2017-08  (D. Ford)     Adapted from bioanal.F90
-   !!----------------------------------------------------------------------
-#if defined key_asminc && defined key_hadocc
+   !! History :  3.6  ! 2017-08 (D. Ford)  Adapted from asmlogchlbal_hadocc
+   !!----------------------------------------------------------------------
+#if defined key_asminc && defined key_medusa && defined key_foam_medusa
    !!----------------------------------------------------------------------
    !! 'key_asminc'          : assimilation increment interface
-   !! 'key_hadocc'          : HadOCC model
-   !!----------------------------------------------------------------------
-   !! asm_logchl_bal_hadocc : routine to calculate increments to HadOCC
+   !! 'key_medusa'          : MEDUSA model
+   !! 'key_foam_medusa'     : MEDUSA extras for FOAM OBS and ASM
+   !!----------------------------------------------------------------------
+   !! asm_logchl_bal_medusa : routine to calculate increments to MEDUSA
    !!----------------------------------------------------------------------
    USE par_kind,      ONLY: wp             ! kind parameters
@@ -24,13 +25,7 @@
    USE zdfmxl                              ! mixed layer depth
    USE iom                                 ! i/o
-   USE trc,           ONLY: trn, trb,    & ! HadOCC variables
-      &                     HADOCC_CHL,  &
-      &                     pgrow_avg,   &
-      &                     ploss_avg,   &
-      &                     phyt_avg,    &
-      &                     mld_max
-   USE par_hadocc                          ! HadOCC parameters
-   USE had_bgc_stnd,  ONLY: kmt            ! HadOCC parameters
-   USE had_bgc_const                       ! HadOCC parameters
+   USE trc,           ONLY: trn, trb       ! MEDUSA variables
+   USE sms_medusa                          ! MEDUSA parameters
+   USE par_medusa                          ! MEDUSA parameters
    USE par_trc,       ONLY: jptra          ! Tracer parameters
    USE bioanalysis                         ! Nitrogen balancing
@@ -39,5 +34,5 @@
    PRIVATE                   
 
-   PUBLIC asm_logchl_bal_hadocc
+   PUBLIC asm_logchl_bal_medusa
 
    ! Default values for biological assimilation parameters
@@ -74,13 +69,18 @@
 CONTAINS
 
-   SUBROUTINE asm_logchl_bal_hadocc( logchl_bkginc, aincper, mld_choice_bgc, &
-      &                              k_maxchlinc, ld_logchlbal, logchl_balinc )
+   SUBROUTINE asm_logchl_bal_medusa( logchl_bkginc, aincper, mld_choice_bgc, &
+      &                              k_maxchlinc, ld_logchlbal,              &
+      &                              pgrow_avg_bkg, ploss_avg_bkg,           &
+      &                              phyt_avg_bkg, mld_max_bkg,              &
+      &                              logchl_balinc )
       !!---------------------------------------------------------------------------
-      !!                    ***  ROUTINE asm_logchl_bal_hadocc  ***
-      !!
-      !! ** Purpose :   calculate increments to HadOCC from logchl increments
+      !!                    ***  ROUTINE asm_logchl_bal_medusa  ***
+      !!
+      !! ** Purpose :   calculate increments to MEDUSA from logchl increments
       !!
       !! ** Method  :   convert logchl increments to chl increments
+      !!                average up MEDUSA to look like HadOCC
       !!                call nitrogen balancing scheme
+      !!                separate back out to MEDUSA
       !!
       !! ** Action  :   populate logchl_balinc
@@ -95,10 +95,29 @@
       REAL(wp), INTENT(in   )                               :: k_maxchlinc    ! Max chl increment
       LOGICAL,  INTENT(in   )                               :: ld_logchlbal   ! Balancing y/n
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc  ! Balancing increments
       !!
-      INTEGER                                               :: ji, jj, jk     ! Loop counters
+      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
       INTEGER                                               :: jkmax          ! Loop index
       INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
+      REAL(wp)                                              :: n2be_p         ! N:biomass for total phy
+      REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo
+      REAL(wp)                                              :: n2be_d         ! N:biomass for detritus
+      REAL(wp)                                              :: zfrac_chn      ! Fraction of jpchn
+      REAL(wp)                                              :: zfrac_chd      ! Fraction of jpchd
+      REAL(wp)                                              :: zfrac_phn      ! Fraction of jpphn
+      REAL(wp)                                              :: zfrac_phd      ! Fraction of jpphd
+      REAL(wp)                                              :: zfrac_zmi      ! Fraction of jpzmi
+      REAL(wp)                                              :: zfrac_zme      ! Fraction of jpzme
+      REAL(wp)                                              :: zrat_pds_phd   ! Ratio of jppds:jpphd
+      REAL(wp)                                              :: zrat_chd_phd   ! Ratio of jpchd:jpphd
+      REAL(wp)                                              :: zrat_chn_phn   ! Ratio of jpchn:jpphn
+      REAL(wp)                                              :: zrat_dtc_det   ! Ratio of jpdtc:jpdet
       REAL(wp),                DIMENSION(jpi,jpj)           :: chl_inc        ! Chlorophyll increments
+      REAL(wp),                DIMENSION(jpi,jpj)           :: medusa_chl     ! MEDUSA total chlorophyll
+      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy
       REAL(wp),                DIMENSION(jpi,jpj)           :: zmld           ! Mixed layer depth
       REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
@@ -117,8 +136,9 @@
       ! 3) Subtract background from analysis to get chl incs
       ! If k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
+      medusa_chl(:,:) = trb(:,:,1,jpchn) + trb(:,:,1,jpchd)
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF ( HADOCC_CHL(ji,jj,1) > 0.0 ) THEN
-               chl_inc(ji,jj) = 10**( LOG10( HADOCC_CHL(ji,jj,1) ) + logchl_bkginc(ji,jj) ) - HADOCC_CHL(ji,jj,1)
+            IF ( medusa_chl(ji,jj) > 0.0 ) THEN
+               chl_inc(ji,jj) = 10**( LOG10( medusa_chl(ji,jj) ) + logchl_bkginc(ji,jj) ) - medusa_chl(ji,jj)
                IF ( k_maxchlinc > 0.0 ) THEN
                   chl_inc(ji,jj) = MAX( -1.0 * k_maxchlinc, MIN( chl_inc(ji,jj), k_maxchlinc ) )
@@ -158,21 +178,27 @@
       IF ( ld_logchlbal ) THEN   ! Nitrogen balancing
 
-         ! Set up model parameters to be passed into Hemmings balancing routine
-         modparm(1)  = grow_sat
-         modparm(2)  = psmax
-         modparm(3)  = par
-         modparm(4)  = alpha
-         modparm(5)  = resp_rate
-         modparm(6)  = pmort_rate
-         modparm(7)  = phyto_min
-         modparm(8)  = z_mort_1
-         modparm(9)  = z_mort_2
-         modparm(10) = c2n_p
-         modparm(11) = c2n_z
-         modparm(12) = c2n_d
-         modparm(13) = graze_threshold
-         modparm(14) = holling_coef
-         modparm(15) = graze_sat
-         modparm(16) = graze_max
+         ! Set up model parameters to be passed into Hemmings balancing routine.
+         ! For now these are hardwired to the standard HadOCC parameter values
+         ! (except C:N ratios) as this is what the scheme was developed for.
+         ! Obviously, HadOCC and MEDUSA are rather different models, so this
+         ! isn't ideal, but there's not always direct analogues between the two
+         ! parameter sets, so it's the easiest way to get something running.
+         ! In the longer term, some serious MarMOT-based development is required.
+         modparm(1)  = 0.1                               ! grow_sat
+         modparm(2)  = 2.0                               ! psmax
+         modparm(3)  = 0.845                             ! par
+         modparm(4)  = 0.02                              ! alpha
+         modparm(5)  = 0.05                              ! resp_rate
+         modparm(6)  = 0.05                              ! pmort_rate
+         modparm(7)  = 0.01                              ! phyto_min
+         modparm(8)  = 0.05                              ! z_mort_1
+         modparm(9)  = 1.0                               ! z_mort_2
+         modparm(10) = ( xthetapn  + xthetapd  ) / 2.0   ! c2n_p
+         modparm(11) = ( xthetazmi + xthetazme ) / 2.0   ! c2n_z
+         modparm(12) = xthetad                           ! c2n_d
+         modparm(13) = 0.01                              ! graze_threshold
+         modparm(14) = 2.0                               ! holling_coef
+         modparm(15) = 0.5                               ! graze_sat
+         modparm(16) = 2.0                               ! graze_max
 
          ! Set up assimilation parameters to be passed into balancing routine
@@ -207,39 +233,126 @@
 
          ! Set background state
-         bstate(:,:,:,i_tracer(1)) = trb(:,:,:,jp_had_nut)
-         bstate(:,:,:,i_tracer(2)) = trb(:,:,:,jp_had_phy)
-         bstate(:,:,:,i_tracer(3)) = trb(:,:,:,jp_had_zoo)
-         bstate(:,:,:,i_tracer(4)) = trb(:,:,:,jp_had_pdn)
-         bstate(:,:,:,i_tracer(5)) = trb(:,:,:,jp_had_dic)
-         bstate(:,:,:,i_tracer(6)) = trb(:,:,:,jp_had_alk)
+         bstate(:,:,:,i_tracer(1)) = trb(:,:,:,jpdin)
+         bstate(:,:,:,i_tracer(2)) = trb(:,:,:,jpphn) + trb(:,:,:,jpphd)
+         bstate(:,:,:,i_tracer(3)) = trb(:,:,:,jpzmi) + trb(:,:,:,jpzme)
+         bstate(:,:,:,i_tracer(4)) = trb(:,:,:,jpdet)
+         bstate(:,:,:,i_tracer(5)) = trb(:,:,:,jpdic)
+         bstate(:,:,:,i_tracer(6)) = trb(:,:,:,jpalk)
+
+         ! Calculate carbon to chlorophyll ratio for combined phytoplankton
+         ! and nitrogen to biomass equivalent for PZD
+         ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA
+         !cchl_p(:,:) = ( trb(:,:,1,jpchn) + trb(:,:,1,jpchd ) ) / &
+         !   &          ( ( trb(:,:,1,jpphn) * xthetapn ) + ( trb(:,:,1,jpphd) * xthetapd ) )
+         cchl_p(:,:) = 0.0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( ( trb(ji,jj,1,jpchn) + trb(ji,jj,1,jpchd ) ) .GT. 0.0 ) THEN
+                  cchl_p(ji,jj) = ( ( trb(ji,jj,1,jpphn) * xthetapn ) + ( trb(ji,jj,1,jpphd) * xthetapd ) ) / &
+                     &            ( trb(ji,jj,1,jpchn) + trb(ji,jj,1,jpchd ) )
+               ENDIF
+            END DO
+         END DO
+         n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) )
+         n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) )
+         n2be_d = 14.01 + ( xmassc * xthetad )
+         
+         WRITE(numout,*) 'DAF: nproc, min/max cchl_p, min/max chl_inc = ', nproc, MINVAL(cchl_p), MAXVAL(cchl_p), MINVAL(chl_inc), MAXVAL(chl_inc)
 
          ! Call nitrogen balancing routine
-         CALL bio_analysis( jpi, jpj, jpk, ZDZ(:,:,:), i_tracer, modparm,          &
-            &               n2be_p, n2be_z, n2be_d, assimparm,                     &
-            &               INT(aincper), 1, kmt(:,:), tmask(:,:,:),               &
-            &               zmld(:,:), mld_max(:,:), chl_inc(:,:), cchl_p(:,:,1),  &
-            &               nbal_active, phyt_avg(:,:),                            &
-            &               gl_active, pgrow_avg(:,:), ploss_avg(:,:),             &
-            &               subsurf_active, deepneg_active,                        &
-            &               deeppos_active, nutprof_active,                        &
-            &               bstate, outincs,                                       &
-            &               diag_active, diag,                                     &
+         CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   &
+            &               n2be_p, n2be_z, n2be_d, assimparm,                      &
+            &               INT(aincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       &
+            &               zmld(:,:), mld_max_bkg(:,:), chl_inc(:,:), cchl_p(:,:), &
+            &               nbal_active, phyt_avg_bkg(:,:),                         &
+            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      &
+            &               subsurf_active, deepneg_active,                         &
+            &               deeppos_active, nutprof_active,                         &
+            &               bstate, outincs,                                        &
+            &               diag_active, diag,                                      &
             &               diag_fulldepth_active, diag_fulldepth )
-
-         ! Save balancing increments
-         logchl_balinc(:,:,:,jp_had_nut) = outincs(:,:,:,i_tracer(1))
-         logchl_balinc(:,:,:,jp_had_phy) = outincs(:,:,:,i_tracer(2))
-         logchl_balinc(:,:,:,jp_had_zoo) = outincs(:,:,:,i_tracer(3))
-         logchl_balinc(:,:,:,jp_had_pdn) = outincs(:,:,:,i_tracer(4))
-         logchl_balinc(:,:,:,jp_had_dic) = outincs(:,:,:,i_tracer(5))
-         logchl_balinc(:,:,:,jp_had_alk) = outincs(:,:,:,i_tracer(6))
+         
+         WRITE(numout,*) 'DAF: nproc, min/max outincs(phy)1,20 = ', nproc, MINVAL(outincs(:,:,1,i_tracer(2))), MAXVAL(outincs(:,:,1,i_tracer(2))), MINVAL(outincs(:,:,20,i_tracer(2))), MAXVAL(outincs(:,:,20,i_tracer(2)))
+
+         ! Loop over each grid point partioning the increments
+         logchl_balinc(:,:,:,:) = 0.0
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+
+                  IF ( ( trb(ji,jj,jk,jpphn) > 0.0 ) .AND. ( trb(ji,jj,jk,jpphd) > 0.0 ) ) THEN
+                     ! Phytoplankton nitrogen and silicate split up based on existing ratios
+                     zfrac_phn = trb(ji,jj,jk,jpphn) / (trb(ji,jj,jk,jpphn) + trb(ji,jj,jk,jpphd))
+                     zfrac_phd = 1.0 - zfrac_phn
+                     zrat_pds_phd = trb(ji,jj,jk,jppds) / trb(ji,jj,jk,jpphd)
+                     logchl_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn
+                     logchl_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd
+                     logchl_balinc(ji,jj,jk,jppds) = logchl_balinc(ji,jj,jk,jpphd) * zrat_pds_phd
+
+                     ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen
+                     ! Not using chl_inc directly as it's only 2D
+                     ! This method should give same results at surface as splitting chl_inc would
+                     zrat_chn_phn = trb(ji,jj,jk,jpchn) / trb(ji,jj,jk,jpphn)
+                     zrat_chd_phd = trb(ji,jj,jk,jpchd) / trb(ji,jj,jk,jpphd)
+                     logchl_balinc(ji,jj,jk,jpchn) = logchl_balinc(ji,jj,jk,jpphn) * zrat_chn_phn
+                     logchl_balinc(ji,jj,jk,jpchd) = logchl_balinc(ji,jj,jk,jpphd) * zrat_chd_phd
+                  ENDIF
+
+                  IF ( ( trb(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( trb(ji,jj,jk,jpzme) > 0.0 ) ) THEN
+                     ! Zooplankton nitrogen split up based on existing ratios
+                     zfrac_zmi = trb(ji,jj,jk,jpzmi) / (trb(ji,jj,jk,jpzmi) + trb(ji,jj,jk,jpzme))
+                     zfrac_zme = 1.0 - zfrac_zmi
+                     logchl_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi
+                     logchl_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme
+                  ENDIF
+
+                  ! Nitrogen nutrient straight from balancing scheme
+                  logchl_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1))
+
+                  ! Nitrogen detritus straight from balancing scheme
+                  logchl_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4))
+
+                  ! DIC straight from balancing scheme
+                  logchl_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5))
+
+                  ! Alkalinity straight from balancing scheme
+                  logchl_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6))
+
+                  ! Remove diatom silicate increment from nutrient silicate to conserve mass
+                  IF ( ( trb(ji,jj,jk,jpsil) - logchl_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN
+                     logchl_balinc(ji,jj,jk,jpsil) = logchl_balinc(ji,jj,jk,jppds) * (-1.0)
+                  ENDIF
+
+                  IF ( ( trb(ji,jj,jk,jpdet) > 0.0 ) .AND. ( trb(ji,jj,jk,jpdtc) > 0.0 ) ) THEN
+                     ! Carbon detritus based on existing ratios
+                     zrat_dtc_det = trb(ji,jj,jk,jpdtc) / trb(ji,jj,jk,jpdet)
+                     logchl_balinc(ji,jj,jk,jpdtc) = logchl_balinc(ji,jj,jk,jpdet) * zrat_dtc_det
+                  ENDIF
+
+                  ! Do nothing with iron or oxygen for the time being
+                  logchl_balinc(ji,jj,jk,jpfer) = 0.0
+                  logchl_balinc(ji,jj,jk,jpoxy) = 0.0
+                  
+               END DO
+            END DO
+         END DO
       
       ELSE   ! No nitrogen balancing
       
-         ! Initialise phytoplankton increment to zero
-         logchl_balinc(:,:,:,jp_had_phy) = 0.0
+         ! Initialise individual chlorophyll increments to zero
+         logchl_balinc(:,:,:,jpchn) = 0.0
+         logchl_balinc(:,:,:,jpchd) = 0.0
          
-         ! Convert surface chlorophyll increment to phytoplankton nitrogen
-         logchl_balinc(:,:,1,jp_had_phy) = ( cchl_p(:,:,1) / (mw_carbon * c2n_p) ) * chl_inc(:,:)
+         ! Split up total surface chlorophyll increments
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( medusa_chl(ji,jj) > 0.0 ) THEN
+                  zfrac_chn = trb(ji,jj,1,jpchn) / medusa_chl(ji,jj)
+                  zfrac_chd = 1.0 - zfrac_chn
+                  logchl_balinc(ji,jj,1,jpchn) = chl_inc(ji,jj) * zfrac_chn
+                  logchl_balinc(ji,jj,1,jpchd) = chl_inc(ji,jj) * zfrac_chd
+               ENDIF
+            END DO
+         END DO
          
          ! Propagate through mixed layer
@@ -257,5 +370,6 @@
                !
                DO jk = 2, jkmax
-                  logchl_balinc(ji,jj,jk,jp_had_phy) = logchl_balinc(ji,jj,1,jp_had_phy)
+                  logchl_balinc(ji,jj,jk,jpchn) = logchl_balinc(ji,jj,1,jpchn)
+                  logchl_balinc(ji,jj,jk,jpchd) = logchl_balinc(ji,jj,1,jpchd)
                END DO
                !
@@ -264,13 +378,21 @@
 
          ! Set other balancing increments to zero
-         logchl_balinc(:,:,:,jp_had_nut) = 0.0
-         logchl_balinc(:,:,:,jp_had_zoo) = 0.0
-         logchl_balinc(:,:,:,jp_had_pdn) = 0.0
-         logchl_balinc(:,:,:,jp_had_dic) = 0.0
-         logchl_balinc(:,:,:,jp_had_alk) = 0.0
-      
+         logchl_balinc(:,:,:,jpphn) = 0.0
+         logchl_balinc(:,:,:,jpphd) = 0.0
+         logchl_balinc(:,:,:,jppds) = 0.0
+         logchl_balinc(:,:,:,jpzmi) = 0.0
+         logchl_balinc(:,:,:,jpzme) = 0.0
+         logchl_balinc(:,:,:,jpdin) = 0.0
+         logchl_balinc(:,:,:,jpsil) = 0.0
+         logchl_balinc(:,:,:,jpfer) = 0.0
+         logchl_balinc(:,:,:,jpdet) = 0.0
+         logchl_balinc(:,:,:,jpdtc) = 0.0
+         logchl_balinc(:,:,:,jpdic) = 0.0
+         logchl_balinc(:,:,:,jpalk) = 0.0
+         logchl_balinc(:,:,:,jpoxy) = 0.0
+
       ENDIF
 
-   END SUBROUTINE asm_logchl_bal_hadocc
+   END SUBROUTINE asm_logchl_bal_medusa
 
 #else
@@ -279,5 +401,5 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE asm_logchl_bal_hadocc( logchl_bkginc, aincper, mld_choice_bgc, &
+   SUBROUTINE asm_logchl_bal_medusa( logchl_bkginc, aincper, mld_choice_bgc, &
       &                              k_maxchlinc, logchl_balinc )
       REAL    :: logchl_bkginc(:,:)
@@ -286,8 +408,8 @@
       REAL    :: k_maxchlinc
       REAL(   :: logchl_balinc(:,:,:,:)
-      WRITE(*,*) 'asm_logchl_bal_hadocc: You should not have seen this print! error?'
-   END SUBROUTINE asm_logchl_bal_hadocc
+      WRITE(*,*) 'asm_logchl_bal_medusa: You should not have seen this print! error?'
+   END SUBROUTINE asm_logchl_bal_medusa
 #endif
 
    !!======================================================================
-END MODULE asmlogchlbal_hadocc
+END MODULE asmlogchlbal_medusa
