Changeset 10622


Ignore:
Timestamp:
2019-02-01T17:27:20+01:00 (14 months ago)
Author:
dford
Message:

Implement biogeochemistry assimilation for FABM-ERSEM.

Location:
branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM
Files:
8 edited
1 copied

Legend:

Unmodified
Added
Removed
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/CONFIG/SHARED/namelist_ref

    r10574 r10622  
    12681268    ln_slchldiainc = .false. !  Logical switch for applying slchldia increments 
    12691269    ln_slchlnoninc = .false. !  Logical switch for applying slchlnon increments 
     1270    ln_slchlnaninc = .false. !  Logical switch for applying slchlnan increments 
     1271    ln_slchlpicinc = .false. !  Logical switch for applying slchlpic increments 
     1272    ln_slchldininc = .false. !  Logical switch for applying slchldin increments 
    12701273    ln_schltotinc  = .false. !  Logical switch for applying schltot increments 
    12711274    ln_slphytotinc = .false. !  Logical switch for applying slphytot increments 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90

    r10574 r10622  
    1111   !! 'key_medusa'          : MEDUSA model 
    1212   !! 'key_roam'            : MEDUSA extras for carbonate cycle 
     13   !! 'key_fabm'            : ERSEM model coupled via FABM 
    1314   !! 'key_karaml'          : Kara mixed layer depth 
    1415   !!--------------------------------------------------------------------------- 
     
    3839      & ln_kara,            & 
    3940#endif    
    40       & hmld,               &  
     41      & hmld,               & 
     42      & hmld_tref,          & 
    4143      & hmlp,               & 
    4244      & hmlpt  
     
    4547      & c_asmbal,           & 
    4648      & nitbkg_r,           & 
     49      & nitavgbkg_r,        & 
    4750      & nitdin_r,           & 
    4851      & nitiaustr_r,        & 
    4952      & nitiaufin_r 
     53   USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA 
     54      & asm_phyto2d_bal_medusa 
     55   USE asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC 
     56      & asm_phyto2d_bal_hadocc 
     57   USE asmphyto2dbal_ersem,  ONLY: & ! phyto2d balancing for ERSEM 
     58      & asm_phyto2d_bal_ersem 
     59   USE asmpco2bal, ONLY:    & ! pCO2 balancing 
     60      & asm_pco2_bal 
    5061#if defined key_top 
    5162   USE trc, ONLY:           & ! passive tracer variables 
    5263      & trn,                & 
    53       & trb 
     64      & trb,                & 
     65      & nittrc000 
    5466   USE par_trc, ONLY:       & ! passive tracer parameters 
    5567      & jptra 
    5668#endif 
    5769#if defined key_medusa 
    58    USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA 
    59       & asm_phyto2d_bal_medusa 
    60    USE asmpco2bal, ONLY:    & ! pCO2 balancing for MEDUSA 
    61       & asm_pco2_bal 
    6270   USE sms_medusa             ! MEDUSA parameters 
    6371   USE par_medusa             ! MEDUSA parameters 
     
    7280      & mld_max 
    7381#elif defined key_hadocc 
    74    USE asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC 
    75       & asm_phyto2d_bal_hadocc 
    76    USE asmpco2bal, ONLY:    & ! pCO2 balancing for HadOCC 
    77       & asm_pco2_bal 
    7882   USE par_hadocc             ! HadOCC parameters 
    7983   USE had_bgc_const          ! HadOCC parameters 
     
    8690   USE had_bgc_const, ONLY: & ! HadOCC C:Chl ratio 
    8791      & cchl_p 
     92#elif defined key_fabm 
     93   USE par_fabm               ! FABM-ERSEM parameters 
     94   USE trc, ONLY:           & ! FABM-ERSEM diagnostic variables 
     95      & pgrow_avg,          & 
     96      & ploss_avg,          & 
     97      & phyt_avg,           & 
     98      & mld_max 
    8899#endif 
    89100 
     
    97108   PUBLIC  asm_bgc_init_bkg       ! called by asm_inc_init in asminc.F90 
    98109   PUBLIC  asm_bgc_bal_wri        ! called by nemo_gcm in nemogcm.F90 
     110   PUBLIC  asm_bgc_bkg_alloc      ! called by asm_bkg_wri in asmbkg.F90 
     111   PUBLIC  asm_bgc_bkg_tavg       ! called by asm_bkg_wri in asmbkg.F90 
    99112   PUBLIC  asm_bgc_bkg_wri        ! called by asm_bkg_wri in asmbkg.F90 
    100113   PUBLIC  phyto2d_asm_inc        ! called by bgc_asm_inc in asminc.F90 
     
    109122   LOGICAL, PUBLIC :: ln_slchldiainc = .FALSE. !: No surface diatom     log10(chlorophyll) increment 
    110123   LOGICAL, PUBLIC :: ln_slchlnoninc = .FALSE. !: No surface non-diatom log10(chlorophyll) increment 
     124   LOGICAL, PUBLIC :: ln_slchlnaninc = .FALSE. !: No surface nanophyto  log10(chlorophyll) increment 
     125   LOGICAL, PUBLIC :: ln_slchlpicinc = .FALSE. !: No surface picophyto  log10(chlorophyll) increment 
     126   LOGICAL, PUBLIC :: ln_slchldininc = .FALSE. !: No surface dinoflag   log10(chlorophyll) increment 
    111127   LOGICAL, PUBLIC :: ln_schltotinc  = .FALSE. !: No surface total      chlorophyll        increment 
    112128   LOGICAL, PUBLIC :: ln_slphytotinc = .FALSE. !: No surface total      log10(phyto C)     increment 
     
    146162   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchldia_bkginc ! slchldia inc 
    147163   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchlnon_bkginc ! slchlnon inc 
     164   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchlnan_bkginc ! slchlnan inc 
     165   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchlpic_bkginc ! slchlpic inc 
     166   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchldin_bkginc ! slchldin inc 
    148167   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: schltot_bkginc  ! schltot inc 
    149168   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slphytot_bkginc ! slphytot inc 
     
    172191   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: chl_bkg         ! Background chl 
    173192   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: cchl_p_bkg      ! Background c:chl 
     193#if defined key_fabm 
     194   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: totalk_bkg      ! Background total alkalinity 
     195#endif 
     196 
     197   REAL(wp), DIMENSION(:,:),     SAVE, ALLOCATABLE :: pgrow_avg_tavg  ! Time average of pgrow_avg 
     198   REAL(wp), DIMENSION(:,:),     SAVE, ALLOCATABLE :: ploss_avg_tavg  ! Time average of ploss_avg 
     199   REAL(wp), DIMENSION(:,:),     SAVE, ALLOCATABLE :: phyt_avg_tavg   ! Time average of phyt_avg 
     200   REAL(wp), DIMENSION(:,:,:,:), SAVE, ALLOCATABLE :: trn_tavg        ! Time average of trn 
     201#if defined key_hadocc 
     202   REAL(wp), DIMENSION(:,:,:),   SAVE, ALLOCATABLE :: HADOCC_CHL_tavg ! Time average of HADOCC_CHL 
     203   REAL(wp), DIMENSION(:,:,:),   SAVE, ALLOCATABLE :: cchl_p_tavg     ! Time average of cchl_p 
     204#elif defined key_fabm 
     205   REAL(wp), DIMENSION(:,:,:),   SAVE, ALLOCATABLE :: totalk_tavg     ! Time average of O3_TA 
     206#endif 
    174207 
    175208#  include "domzgr_substitute.h90" 
     
    190223      !!------------------------------------------------------------------------ 
    191224 
    192 #if ! defined key_top || ( ! defined key_hadocc && ! defined key_medusa ) 
     225#if ! defined key_top || ( ! defined key_hadocc && ! defined key_medusa && ! defined key_fabm ) 
    193226      CALL ctl_stop( ' Attempting to assimilate biogeochemical observations', & 
    194227         &           ' but no compatible biogeochemical model is available' ) 
     
    208241#endif 
    209242 
     243#if defined key_fabm 
     244      IF ( ln_pphinc ) THEN 
     245         CALL ctl_stop( ' Cannot currently assimilate pH into FABM-ERSEM' ) 
     246      ENDIF 
     247#endif 
     248 
    210249      IF ( ( ln_phytobal ).AND.                                       & 
    211250         & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. & 
    212          & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. & 
    213          & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. & 
    214          & ( .NOT. ln_slphynoninc ) ) THEN 
     251         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. & 
     252         & ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. & 
     253         & ( .NOT. ln_schltotinc  ).AND.( .NOT. ln_slphytotinc ).AND. & 
     254         & ( .NOT. ln_slphydiainc ).AND.( .NOT. ln_slphynoninc ) ) THEN 
    215255         CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', & 
    216256            &           ' if not assimilating ocean colour,',                   & 
     
    221261      IF ( ( ln_balwri ).AND.                                         & 
    222262         & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. & 
    223          & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. & 
     263         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. & 
     264         & ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. & 
     265         & ( .NOT. ln_schltotinc  ).AND.                              & 
    224266         & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. & 
    225267         & ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. & 
     
    244286 
    245287      IF ( ( ln_slchltotinc .OR. ln_schltotinc  ) .AND. & 
    246          & ( ln_slchldiainc .OR. ln_slchlnoninc ) ) THEN 
     288         & ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slchlnaninc .OR. & 
     289         &   ln_slchlpicinc .OR. ln_slchldininc ) ) THEN 
    247290         CALL ctl_stop( ' Can only assimilate total or PFT surface chlorophyll, not both' ) 
    248291      ENDIF 
     
    250293      IF ( ln_phytobal .AND.                                      & 
    251294         & ( ( ln_slchlnoninc .AND.( .NOT. ln_slchldiainc ) ).OR. & 
     295#if defined key_fabm 
     296         &   (  ( ln_slchldiainc .OR. ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc ) .AND. & 
     297         &      ( ( .NOT. ln_slchldiainc ) .OR. ( .NOT. ln_slchlnaninc ) .OR. & 
     298         &        ( .NOT. ln_slchlpicinc ) .OR. ( .NOT. ln_slchldininc ) ) ) .OR. & 
     299#endif 
    252300         &   ( ln_slchldiainc .AND.( .NOT. ln_slchlnoninc ) ) ) ) THEN 
    253301         CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', & 
     
    306354      ENDIF 
    307355       
     356      IF ( ln_slchlnaninc ) THEN 
     357         ALLOCATE( slchlnan_bkginc(jpi,jpj) ) 
     358         CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnan', slchlnan_bkginc ) 
     359      ENDIF 
     360       
     361      IF ( ln_slchlpicinc ) THEN 
     362         ALLOCATE( slchlpic_bkginc(jpi,jpj) ) 
     363         CALL asm_bgc_read_incs_2d( knum, 'bckinslchlpic', slchlpic_bkginc ) 
     364      ENDIF 
     365       
     366      IF ( ln_slchldininc ) THEN 
     367         ALLOCATE( slchldin_bkginc(jpi,jpj) ) 
     368         CALL asm_bgc_read_incs_2d( knum, 'bckinslchldin', slchldin_bkginc ) 
     369      ENDIF 
     370       
    308371      IF ( ln_schltotinc ) THEN 
    309372         ALLOCATE( schltot_bkginc(jpi,jpj) ) 
     
    379442       
    380443      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & 
     444         & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. & 
    381445         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & 
    382446         & ln_slphynoninc ) THEN 
     
    516580      !!------------------------------------------------------------------------ 
    517581 
    518 #if defined key_top && ( defined key_hadocc || defined key_medusa ) 
     582#if defined key_top && ( defined key_hadocc || defined key_medusa || defined key_fabm ) 
    519583      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & 
     584         & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. & 
    520585         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & 
    521586         & ln_slphynoninc .OR. ln_plchltotinc .OR. ln_pchltotinc ) THEN 
     
    557622         CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) ) 
    558623         CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) ) 
     624#elif defined key_fabm 
     625         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl1', tracer_bkg(:,:,:,jp_fabm_chl1) ) 
     626         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl2', tracer_bkg(:,:,:,jp_fabm_chl2) ) 
     627         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl3', tracer_bkg(:,:,:,jp_fabm_chl3) ) 
     628         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl4', tracer_bkg(:,:,:,jp_fabm_chl4) ) 
     629         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1c',  tracer_bkg(:,:,:,jp_fabm_p1c)  ) 
     630         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1n',  tracer_bkg(:,:,:,jp_fabm_p1n)  ) 
     631         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1p',  tracer_bkg(:,:,:,jp_fabm_p1p)  ) 
     632         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1s',  tracer_bkg(:,:,:,jp_fabm_p1s)  ) 
     633         CALL iom_get( inum, jpdom_autoglo, 'ersem_p2c',  tracer_bkg(:,:,:,jp_fabm_p2c)  ) 
     634         CALL iom_get( inum, jpdom_autoglo, 'ersem_p2n',  tracer_bkg(:,:,:,jp_fabm_p2n)  ) 
     635         CALL iom_get( inum, jpdom_autoglo, 'ersem_p2p',  tracer_bkg(:,:,:,jp_fabm_p2p)  ) 
     636         CALL iom_get( inum, jpdom_autoglo, 'ersem_p3c',  tracer_bkg(:,:,:,jp_fabm_p3c)  ) 
     637         CALL iom_get( inum, jpdom_autoglo, 'ersem_p3n',  tracer_bkg(:,:,:,jp_fabm_p3n)  ) 
     638         CALL iom_get( inum, jpdom_autoglo, 'ersem_p3p',  tracer_bkg(:,:,:,jp_fabm_p3p)  ) 
     639         CALL iom_get( inum, jpdom_autoglo, 'ersem_p4c',  tracer_bkg(:,:,:,jp_fabm_p4c)  ) 
     640         CALL iom_get( inum, jpdom_autoglo, 'ersem_p4n',  tracer_bkg(:,:,:,jp_fabm_p4n)  ) 
     641         CALL iom_get( inum, jpdom_autoglo, 'ersem_p4p',  tracer_bkg(:,:,:,jp_fabm_p4p)  ) 
    559642#endif 
    560643          
     
    588671            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) 
    589672            CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) ) 
     673#elif defined key_fabm 
     674            CALL iom_get( inum, jpdom_autoglo, 'ersem_z4c',  tracer_bkg(:,:,:,jp_fabm_z4c)  ) 
     675            CALL iom_get( inum, jpdom_autoglo, 'ersem_z5c',  tracer_bkg(:,:,:,jp_fabm_z5c)  ) 
     676            CALL iom_get( inum, jpdom_autoglo, 'ersem_z5n',  tracer_bkg(:,:,:,jp_fabm_z5n)  ) 
     677            CALL iom_get( inum, jpdom_autoglo, 'ersem_z5p',  tracer_bkg(:,:,:,jp_fabm_z5p)  ) 
     678            CALL iom_get( inum, jpdom_autoglo, 'ersem_z6c',  tracer_bkg(:,:,:,jp_fabm_z6c)  ) 
     679            CALL iom_get( inum, jpdom_autoglo, 'ersem_z6n',  tracer_bkg(:,:,:,jp_fabm_z6n)  ) 
     680            CALL iom_get( inum, jpdom_autoglo, 'ersem_z6p',  tracer_bkg(:,:,:,jp_fabm_z6p)  ) 
     681            CALL iom_get( inum, jpdom_autoglo, 'ersem_n1p',  tracer_bkg(:,:,:,jp_fabm_n1p)  ) 
     682            CALL iom_get( inum, jpdom_autoglo, 'ersem_n3n',  tracer_bkg(:,:,:,jp_fabm_n3n)  ) 
     683            CALL iom_get( inum, jpdom_autoglo, 'ersem_n4n',  tracer_bkg(:,:,:,jp_fabm_n4n)  ) 
     684            CALL iom_get( inum, jpdom_autoglo, 'ersem_n5s',  tracer_bkg(:,:,:,jp_fabm_n5s)  ) 
     685            CALL iom_get( inum, jpdom_autoglo, 'ersem_o2o',  tracer_bkg(:,:,:,jp_fabm_o2o)  ) 
     686            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c',  tracer_bkg(:,:,:,jp_fabm_o3c)  ) 
     687            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) ) 
     688            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:)              ) 
     689            totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:) 
    590690#endif 
    591691         ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN 
     
    596696            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) 
    597697            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) 
     698#elif defined key_fabm 
     699            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c',  tracer_bkg(:,:,:,jp_fabm_o3c)  ) 
     700            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) ) 
     701            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:)              ) 
     702            totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:) 
    598703#endif 
    599704            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   ) 
     
    622727         CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) 
    623728         CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) 
     729#elif defined key_fabm 
     730         CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c',  tracer_bkg(:,:,:,jp_fabm_o3c)  ) 
     731         CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_o3ba) ) 
     732         CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:)              ) 
     733         totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:) 
    624734#endif 
    625735         CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   ) 
     
    687797 
    688798         IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & 
     799            & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. & 
    689800            & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & 
    690801            & ln_slphynoninc ) THEN 
     
    716827               CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) ) 
    717828            ENDIF 
     829#elif defined key_fabm 
     830            CALL iom_rstput( kt, kt, inum, 'phy2d_chl1', phyto2d_balinc(:,:,:,jp_fabm_chl1) ) 
     831            CALL iom_rstput( kt, kt, inum, 'phy2d_chl2', phyto2d_balinc(:,:,:,jp_fabm_chl2) ) 
     832            CALL iom_rstput( kt, kt, inum, 'phy2d_chl3', phyto2d_balinc(:,:,:,jp_fabm_chl3) ) 
     833            CALL iom_rstput( kt, kt, inum, 'phy2d_chl4', phyto2d_balinc(:,:,:,jp_fabm_chl4) ) 
     834            CALL iom_rstput( kt, kt, inum, 'phy2d_p1c',  phyto2d_balinc(:,:,:,jp_fabm_p1c)  ) 
     835            CALL iom_rstput( kt, kt, inum, 'phy2d_p1n',  phyto2d_balinc(:,:,:,jp_fabm_p1n)  ) 
     836            CALL iom_rstput( kt, kt, inum, 'phy2d_p1p',  phyto2d_balinc(:,:,:,jp_fabm_p1p)  ) 
     837            CALL iom_rstput( kt, kt, inum, 'phy2d_p1s',  phyto2d_balinc(:,:,:,jp_fabm_p1s)  ) 
     838            CALL iom_rstput( kt, kt, inum, 'phy2d_p2c',  phyto2d_balinc(:,:,:,jp_fabm_p2c)  ) 
     839            CALL iom_rstput( kt, kt, inum, 'phy2d_p2n',  phyto2d_balinc(:,:,:,jp_fabm_p2n)  ) 
     840            CALL iom_rstput( kt, kt, inum, 'phy2d_p2p',  phyto2d_balinc(:,:,:,jp_fabm_p2p)  ) 
     841            CALL iom_rstput( kt, kt, inum, 'phy2d_p3c',  phyto2d_balinc(:,:,:,jp_fabm_p3c)  ) 
     842            CALL iom_rstput( kt, kt, inum, 'phy2d_p3n',  phyto2d_balinc(:,:,:,jp_fabm_p3n)  ) 
     843            CALL iom_rstput( kt, kt, inum, 'phy2d_p3p',  phyto2d_balinc(:,:,:,jp_fabm_p3p)  ) 
     844            CALL iom_rstput( kt, kt, inum, 'phy2d_p4c',  phyto2d_balinc(:,:,:,jp_fabm_p4c)  ) 
     845            CALL iom_rstput( kt, kt, inum, 'phy2d_p4n',  phyto2d_balinc(:,:,:,jp_fabm_p4n)  ) 
     846            CALL iom_rstput( kt, kt, inum, 'phy2d_p4p',  phyto2d_balinc(:,:,:,jp_fabm_p4p)  ) 
     847            IF ( ln_phytobal ) THEN 
     848               CALL iom_rstput( kt, kt, inum, 'phy2d_z4c',   phyto2d_balinc(:,:,:,jp_fabm_z4c)  ) 
     849               CALL iom_rstput( kt, kt, inum, 'phy2d_z5c',   phyto2d_balinc(:,:,:,jp_fabm_z5c)  ) 
     850               CALL iom_rstput( kt, kt, inum, 'phy2d_z5n',   phyto2d_balinc(:,:,:,jp_fabm_z5n)  ) 
     851               CALL iom_rstput( kt, kt, inum, 'phy2d_z5p',   phyto2d_balinc(:,:,:,jp_fabm_z5p)  ) 
     852               CALL iom_rstput( kt, kt, inum, 'phy2d_z6c',   phyto2d_balinc(:,:,:,jp_fabm_z6c)  ) 
     853               CALL iom_rstput( kt, kt, inum, 'phy2d_z6n',   phyto2d_balinc(:,:,:,jp_fabm_z6n)  ) 
     854               CALL iom_rstput( kt, kt, inum, 'phy2d_z6p',   phyto2d_balinc(:,:,:,jp_fabm_z6p)  ) 
     855               CALL iom_rstput( kt, kt, inum, 'phy2d_n1p',   phyto2d_balinc(:,:,:,jp_fabm_n1p)  ) 
     856               CALL iom_rstput( kt, kt, inum, 'phy2d_n3n',   phyto2d_balinc(:,:,:,jp_fabm_n3n)  ) 
     857               CALL iom_rstput( kt, kt, inum, 'phy2d_n4n',   phyto2d_balinc(:,:,:,jp_fabm_n4n)  ) 
     858               CALL iom_rstput( kt, kt, inum, 'phy2d_n5s',   phyto2d_balinc(:,:,:,jp_fabm_n5s)  ) 
     859               CALL iom_rstput( kt, kt, inum, 'phy2d_o2o',   phyto2d_balinc(:,:,:,jp_fabm_o2o)  ) 
     860               CALL iom_rstput( kt, kt, inum, 'phy2d_o3c',   phyto2d_balinc(:,:,:,jp_fabm_o3c)  ) 
     861               CALL iom_rstput( kt, kt, inum, 'phy2d_o3ba',  phyto2d_balinc(:,:,:,jp_fabm_o3ba) ) 
     862            ENDIF 
    718863#endif 
    719864         ENDIF 
     
    728873#elif defined key_hadocc 
    729874            CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) ) 
     875#elif defined key_fabm 
     876            CALL iom_rstput( kt, kt, inum, 'phy3d_chl1', phyto3d_balinc(:,:,:,jp_fabm_chl1) ) 
     877            CALL iom_rstput( kt, kt, inum, 'phy3d_chl2', phyto3d_balinc(:,:,:,jp_fabm_chl2) ) 
     878            CALL iom_rstput( kt, kt, inum, 'phy3d_chl3', phyto3d_balinc(:,:,:,jp_fabm_chl3) ) 
     879            CALL iom_rstput( kt, kt, inum, 'phy3d_chl4', phyto3d_balinc(:,:,:,jp_fabm_chl4) ) 
     880            CALL iom_rstput( kt, kt, inum, 'phy3d_p1c',  phyto3d_balinc(:,:,:,jp_fabm_p1c)  ) 
     881            CALL iom_rstput( kt, kt, inum, 'phy3d_p1n',  phyto3d_balinc(:,:,:,jp_fabm_p1n)  ) 
     882            CALL iom_rstput( kt, kt, inum, 'phy3d_p1p',  phyto3d_balinc(:,:,:,jp_fabm_p1p)  ) 
     883            CALL iom_rstput( kt, kt, inum, 'phy3d_p1s',  phyto3d_balinc(:,:,:,jp_fabm_p1s)  ) 
     884            CALL iom_rstput( kt, kt, inum, 'phy3d_p2c',  phyto3d_balinc(:,:,:,jp_fabm_p2c)  ) 
     885            CALL iom_rstput( kt, kt, inum, 'phy3d_p2n',  phyto3d_balinc(:,:,:,jp_fabm_p2n)  ) 
     886            CALL iom_rstput( kt, kt, inum, 'phy3d_p2p',  phyto3d_balinc(:,:,:,jp_fabm_p2p)  ) 
     887            CALL iom_rstput( kt, kt, inum, 'phy3d_p3c',  phyto3d_balinc(:,:,:,jp_fabm_p3c)  ) 
     888            CALL iom_rstput( kt, kt, inum, 'phy3d_p3n',  phyto3d_balinc(:,:,:,jp_fabm_p3n)  ) 
     889            CALL iom_rstput( kt, kt, inum, 'phy3d_p3p',  phyto3d_balinc(:,:,:,jp_fabm_p3p)  ) 
     890            CALL iom_rstput( kt, kt, inum, 'phy3d_p4c',  phyto3d_balinc(:,:,:,jp_fabm_p4c)  ) 
     891            CALL iom_rstput( kt, kt, inum, 'phy3d_p4n',  phyto3d_balinc(:,:,:,jp_fabm_p4n)  ) 
     892            CALL iom_rstput( kt, kt, inum, 'phy3d_p4p',  phyto3d_balinc(:,:,:,jp_fabm_p4p)  ) 
    730893#endif 
    731894         ENDIF 
     
    738901            CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) ) 
    739902            CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) ) 
     903#elif defined key_fabm 
     904            CALL iom_rstput( kt, kt, inum, 'pco2_o3c',  pco2_balinc(:,:,:,jp_fabm_o3c)  ) 
     905            CALL iom_rstput( kt, kt, inum, 'pco2_o3ba', pco2_balinc(:,:,:,jp_fabm_o3ba) ) 
    740906#endif 
    741907         ELSE IF ( ln_sfco2inc ) THEN 
     
    746912            CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) ) 
    747913            CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) ) 
     914#elif defined key_fabm 
     915            CALL iom_rstput( kt, kt, inum, 'fco2_o3c',  pco2_balinc(:,:,:,jp_fabm_o3c)  ) 
     916            CALL iom_rstput( kt, kt, inum, 'fco2_o3ba', pco2_balinc(:,:,:,jp_fabm_o3ba) ) 
    748917#endif 
    749918         ENDIF 
     
    756925            CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) ) 
    757926            CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) ) 
     927#elif defined key_fabm 
     928            CALL iom_rstput( kt, kt, inum, 'ph_o3c',  ph_balinc(:,:,:,jp_fabm_o3c)  ) 
     929            CALL iom_rstput( kt, kt, inum, 'ph_o3ba', ph_balinc(:,:,:,jp_fabm_o3ba) ) 
    758930#endif 
    759931         ENDIF 
     
    774946   !!=========================================================================== 
    775947 
    776    SUBROUTINE asm_bgc_bkg_wri( kt, knum ) 
     948   SUBROUTINE asm_bgc_bkg_alloc 
     949      !!------------------------------------------------------------------------ 
     950      !!                    ***  ROUTINE asm_bgc_bkg_alloc  *** 
     951      !! 
     952      !! ** Purpose :   allocate time-average arrays for background 
     953      !! 
     954      !! ** Method  :   allocate time-average arrays for background 
     955      !! 
     956      !! ** Action  :   allocate time-average arrays for background 
     957      !! 
     958      !! References :   asm_bkg_wri 
     959      !!------------------------------------------------------------------------ 
     960      !! 
     961      INTEGER :: ierror 
     962      !! 
     963      !!------------------------------------------------------------------------ 
     964       
     965      ALLOCATE( pgrow_avg_tavg(jpi,jpj), STAT=ierror ) 
     966      IF( ierror > 0 ) THEN 
     967         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate pgrow_avg_tavg' ) 
     968      ENDIF 
     969      pgrow_avg_tavg(:,:) = 0.0 
     970       
     971      ALLOCATE( ploss_avg_tavg(jpi,jpj), STAT=ierror ) 
     972      IF( ierror > 0 ) THEN 
     973         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate ploss_avg_tavg' ) 
     974      ENDIF 
     975      ploss_avg_tavg(:,:) = 0.0 
     976       
     977      ALLOCATE( phyt_avg_tavg(jpi,jpj), STAT=ierror ) 
     978      IF( ierror > 0 ) THEN 
     979         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate phyt_avg_tavg' ) 
     980      ENDIF 
     981      phyt_avg_tavg(:,:) = 0.0 
     982       
     983      ALLOCATE( trn_tavg(jpi,jpj,jpk,jptra), STAT=ierror ) 
     984      IF( ierror > 0 ) THEN 
     985         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate trn_tavg' ) 
     986      ENDIF 
     987      trn_tavg(:,:,:,:) = 0.0 
     988 
     989#if defined key_hadocc 
     990      ALLOCATE( HADOCC_CHL_tavg(jpi,jpj,jpk), STAT=ierror ) 
     991      IF( ierror > 0 ) THEN 
     992         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate HADOCC_CHL_tavg' ) 
     993      ENDIF 
     994      HADOCC_CHL_tavg(:,:,:) = 0.0 
     995 
     996      ALLOCATE( cchl_p_tavg(jpi,jpj,jpk), STAT=ierror ) 
     997      IF( ierror > 0 ) THEN 
     998         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate cchl_p_tavg' ) 
     999      ENDIF 
     1000      cchl_p_tavg(:,:,:) = 0.0 
     1001 
     1002#elif defined key_fabm 
     1003      ALLOCATE( totalk_tavg(jpi,jpj,jpk), STAT=ierror ) 
     1004      IF( ierror > 0 ) THEN 
     1005         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate totalk_tavg' ) 
     1006      ENDIF 
     1007      totalk_tavg(:,:,:) = 0.0 
     1008#endif 
     1009                                  
     1010   END SUBROUTINE asm_bgc_bkg_alloc 
     1011 
     1012   !!=========================================================================== 
     1013   !!=========================================================================== 
     1014   !!=========================================================================== 
     1015 
     1016   SUBROUTINE asm_bgc_bkg_tavg(kt, pnumtimes_tavg) 
    7771017      !!------------------------------------------------------------------------ 
    7781018      !!                    ***  ROUTINE asm_bgc_bkg_wri  *** 
     
    7871027      !!------------------------------------------------------------------------ 
    7881028      !! 
    789       INTEGER, INTENT(in   ) :: kt     ! Current time-step 
    790       INTEGER, INTENT(in   ) :: knum   ! i/o unit of increments file 
    791       !! 
    792       !!------------------------------------------------------------------------ 
    793  
     1029      INTEGER,  INTENT(in   ) :: kt               ! Current time-step 
     1030      !! 
     1031      REAL(wp), INTENT(in   ) :: pnumtimes_tavg   ! No of times to average over 
     1032      !! 
     1033      !!------------------------------------------------------------------------ 
     1034       
     1035      IF (kt == nittrc000) THEN 
     1036         pgrow_avg_tavg(:,:) = 0.0 
     1037         ploss_avg_tavg(:,:) = 0.0 
     1038         phyt_avg_tavg(:,:)  = 0.0 
     1039         trn_tavg(:,:,:,:)   = 0.0 
    7941040#if defined key_hadocc 
    795       CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg             ) 
    796       CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg             ) 
    797       CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg              ) 
    798       CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max               ) 
    799       CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_nut'  , trn(:,:,:,jp_had_nut) ) 
    800       CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_phy'  , trn(:,:,:,jp_had_phy) ) 
    801       CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_zoo'  , trn(:,:,:,jp_had_zoo) ) 
    802       CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_pdn'  , trn(:,:,:,jp_had_pdn) ) 
    803       CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_dic'  , trn(:,:,:,jp_had_dic) ) 
    804       CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_alk'  , trn(:,:,:,jp_had_alk) ) 
    805       CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_chl'  , HADOCC_CHL(:,:,:)     ) 
    806       CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:)         ) 
     1041         HADOCC_CHL_tavg(:,:,:) = 0.0 
     1042         cchl_p_tavg(:,:,:)     = 0.0 
     1043#elif defined key_fabm 
     1044         totalk_tavg(:,:,:)  = 0.0 
     1045#endif 
     1046      ENDIF 
     1047       
     1048      pgrow_avg_tavg(:,:) = pgrow_avg_tavg(:,:) + pgrow_avg(:,:) / pnumtimes_tavg 
     1049      ploss_avg_tavg(:,:) = ploss_avg_tavg(:,:) + ploss_avg(:,:) / pnumtimes_tavg 
     1050      phyt_avg_tavg(:,:)  = phyt_avg_tavg(:,:)  + phyt_avg(:,:)  / pnumtimes_tavg 
     1051      trn_tavg(:,:,:,:)   = trn_tavg(:,:,:,:)   + trn(:,:,:,:)   / pnumtimes_tavg 
     1052#if defined key_hadocc 
     1053      HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL_tavg(:,:,:) + HADOCC_CHL(:,:,:) / pnumtimes_tavg 
     1054      cchl_p_tavg(:,:,:)     = cchl_p_tavg(:,:,:)     + cchl_p(:,:,:)     / pnumtimes_tavg 
     1055#elif defined key_fabm 
     1056      totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) + & 
     1057         &                  fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) / pnumtimes_tavg 
     1058      totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:) 
     1059#endif 
     1060                                  
     1061   END SUBROUTINE asm_bgc_bkg_tavg 
     1062 
     1063   !!=========================================================================== 
     1064   !!=========================================================================== 
     1065   !!=========================================================================== 
     1066 
     1067   SUBROUTINE asm_bgc_bkg_wri( kt, knum, ld_avgbkg ) 
     1068      !!------------------------------------------------------------------------ 
     1069      !!                    ***  ROUTINE asm_bgc_bkg_wri  *** 
     1070      !! 
     1071      !! ** Purpose :   write out bgc background 
     1072      !! 
     1073      !! ** Method  :   write out bgc background 
     1074      !! 
     1075      !! ** Action  :   write out bgc background 
     1076      !! 
     1077      !! References :   asm_bkg_wri 
     1078      !!------------------------------------------------------------------------ 
     1079      !! 
     1080      INTEGER, INTENT(in   ) :: kt          ! Current time-step 
     1081      INTEGER, INTENT(in   ) :: knum        ! i/o unit of increments file 
     1082      LOGICAL, INTENT(in   ) :: ld_avgbkg   ! Averaged background? 
     1083      !! 
     1084      INTEGER                :: nitbgcbkg_r ! Period referenced to nit000 
     1085      !! 
     1086      !!------------------------------------------------------------------------ 
     1087       
     1088      IF (ld_avgbkg) THEN 
     1089         nitbgcbkg_r = nitavgbkg_r 
     1090      ELSE 
     1091         nitbgcbkg_r = nitbkg_r 
     1092         pgrow_avg_tavg(:,:) = pgrow_avg(:,:) 
     1093         ploss_avg_tavg(:,:) = ploss_avg(:,:) 
     1094         phyt_avg_tavg(:,:)  = phyt_avg(:,:) 
     1095         trn_tavg(:,:,:,:)   = trn(:,:,:,:) 
     1096#if defined key_hadocc 
     1097         HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL(:,:,:) 
     1098         cchl_p_tavg(:,:,:)     = cchl_p(:,:,:) 
     1099#elif defined key_fabm 
     1100         totalk_tavg(:,:,:)  = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) 
     1101         totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:) 
     1102#endif 
     1103      ENDIF 
     1104 
     1105#if defined key_hadocc 
     1106      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg'   , pgrow_avg_tavg             ) 
     1107      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg'   , ploss_avg_tavg             ) 
     1108      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg'    , phyt_avg_tavg              ) 
     1109      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max'     , mld_max                    ) 
     1110      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_nut'  , trn_tavg(:,:,:,jp_had_nut) ) 
     1111      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_phy'  , trn_tavg(:,:,:,jp_had_phy) ) 
     1112      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_zoo'  , trn_tavg(:,:,:,jp_had_zoo) ) 
     1113      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_pdn'  , trn_tavg(:,:,:,jp_had_pdn) ) 
     1114      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_dic'  , trn_tavg(:,:,:,jp_had_dic) ) 
     1115      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_alk'  , trn_tavg(:,:,:,jp_had_alk) ) 
     1116      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_chl'  , HADOCC_CHL(:,:,:)     ) 
     1117      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:)         ) 
    8071118#elif defined key_medusa 
    808       CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg        ) 
    809       CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg        ) 
    810       CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg         ) 
    811       CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max          ) 
    812       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn'  , trn(:,:,:,jpchn) ) 
    813       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chd'  , trn(:,:,:,jpchd) ) 
    814       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phn'  , trn(:,:,:,jpphn) ) 
    815       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phd'  , trn(:,:,:,jpphd) ) 
    816       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_pds'  , trn(:,:,:,jppds) ) 
    817       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zmi'  , trn(:,:,:,jpzmi) ) 
    818       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zme'  , trn(:,:,:,jpzme) ) 
    819       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_din'  , trn(:,:,:,jpdin) ) 
    820       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_sil'  , trn(:,:,:,jpsil) ) 
    821       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_fer'  , trn(:,:,:,jpfer) ) 
    822       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_det'  , trn(:,:,:,jpdet) ) 
    823       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dtc'  , trn(:,:,:,jpdtc) ) 
    824       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dic'  , trn(:,:,:,jpdic) ) 
    825       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_alk'  , trn(:,:,:,jpalk) ) 
    826       CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_oxy'  , trn(:,:,:,jpoxy) ) 
     1119      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg'   , pgrow_avg_tavg        ) 
     1120      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg'   , ploss_avg_tavg        ) 
     1121      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg'    , phyt_avg_tavg         ) 
     1122      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max'     , mld_max               ) 
     1123      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chn'  , trn_tavg(:,:,:,jpchn) ) 
     1124      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chd'  , trn_tavg(:,:,:,jpchd) ) 
     1125      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phn'  , trn_tavg(:,:,:,jpphn) ) 
     1126      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phd'  , trn_tavg(:,:,:,jpphd) ) 
     1127      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_pds'  , trn_tavg(:,:,:,jppds) ) 
     1128      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zmi'  , trn_tavg(:,:,:,jpzmi) ) 
     1129      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zme'  , trn_tavg(:,:,:,jpzme) ) 
     1130      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_din'  , trn_tavg(:,:,:,jpdin) ) 
     1131      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_sil'  , trn_tavg(:,:,:,jpsil) ) 
     1132      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_fer'  , trn_tavg(:,:,:,jpfer) ) 
     1133      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_det'  , trn_tavg(:,:,:,jpdet) ) 
     1134      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dtc'  , trn_tavg(:,:,:,jpdtc) ) 
     1135      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dic'  , trn_tavg(:,:,:,jpdic) ) 
     1136      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_alk'  , trn_tavg(:,:,:,jpalk) ) 
     1137      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_oxy'  , trn_tavg(:,:,:,jpoxy) ) 
     1138#elif defined key_fabm 
     1139      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg'   , pgrow_avg_tavg               ) 
     1140      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg'   , ploss_avg_tavg               ) 
     1141      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg'    , phyt_avg_tavg                ) 
     1142      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max'     , mld_max                      ) 
     1143      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl1'  , trn_tavg(:,:,:,jp_fabm_chl1) ) 
     1144      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl2'  , trn_tavg(:,:,:,jp_fabm_chl2) ) 
     1145      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl3'  , trn_tavg(:,:,:,jp_fabm_chl3) ) 
     1146      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl4'  , trn_tavg(:,:,:,jp_fabm_chl4) ) 
     1147      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1c'   , trn_tavg(:,:,:,jp_fabm_p1c)  ) 
     1148      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1n'   , trn_tavg(:,:,:,jp_fabm_p1n)  ) 
     1149      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1p'   , trn_tavg(:,:,:,jp_fabm_p1p)  ) 
     1150      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1s'   , trn_tavg(:,:,:,jp_fabm_p1s)  ) 
     1151      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2c'   , trn_tavg(:,:,:,jp_fabm_p2c)  ) 
     1152      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2n'   , trn_tavg(:,:,:,jp_fabm_p2n)  ) 
     1153      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2p'   , trn_tavg(:,:,:,jp_fabm_p2p)  ) 
     1154      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3c'   , trn_tavg(:,:,:,jp_fabm_p3c)  ) 
     1155      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3n'   , trn_tavg(:,:,:,jp_fabm_p3n)  ) 
     1156      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3p'   , trn_tavg(:,:,:,jp_fabm_p3p)  ) 
     1157      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4c'   , trn_tavg(:,:,:,jp_fabm_p4c)  ) 
     1158      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4n'   , trn_tavg(:,:,:,jp_fabm_p4n)  ) 
     1159      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4p'   , trn_tavg(:,:,:,jp_fabm_p4p)  ) 
     1160      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z4c'   , trn_tavg(:,:,:,jp_fabm_z4c)  ) 
     1161      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5c'   , trn_tavg(:,:,:,jp_fabm_z5c)  ) 
     1162      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5n'   , trn_tavg(:,:,:,jp_fabm_z5n)  ) 
     1163      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5p'   , trn_tavg(:,:,:,jp_fabm_z5p)  ) 
     1164      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6c'   , trn_tavg(:,:,:,jp_fabm_z6c)  ) 
     1165      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6n'   , trn_tavg(:,:,:,jp_fabm_z6n)  ) 
     1166      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6p'   , trn_tavg(:,:,:,jp_fabm_z6p)  ) 
     1167      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n1p'   , trn_tavg(:,:,:,jp_fabm_n1p)  ) 
     1168      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n3n'   , trn_tavg(:,:,:,jp_fabm_n3n)  ) 
     1169      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n4n'   , trn_tavg(:,:,:,jp_fabm_n4n)  ) 
     1170      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n5s'   , trn_tavg(:,:,:,jp_fabm_n5s)  ) 
     1171      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o2o'   , trn_tavg(:,:,:,jp_fabm_o2o)  ) 
     1172      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3c'   , trn_tavg(:,:,:,jp_fabm_o3c)  ) 
     1173      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ba'  , trn_tavg(:,:,:,jp_fabm_o3ba) ) 
     1174      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ta'  , totalk_tavg                  ) 
    8271175#endif 
    8281176 
     
    9111259      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phynon ! Local phynon incs 
    9121260      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg 
     1261#elif defined key_fabm 
     1262      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs 
     1263      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg 
     1264      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldin ! Local chldin incs 
     1265      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldin ! Local chldin bkg 
     1266      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnan ! Local chlnan incs 
     1267      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnan ! Local chlnan bkg 
     1268      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlpic ! Local chlpic incs 
     1269      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlpic ! Local chlpic bkg 
    9131270#endif 
    9141271      !!------------------------------------------------------------------------ 
     
    9261283#elif defined key_hadocc 
    9271284            zbkg_chltot(:,:) = chl_bkg(:,:,1) 
     1285#elif defined key_fabm 
     1286            zbkg_chltot(:,:) = tracer_bkg(:,:,1,jp_fabm_chl1) + & 
     1287               &               tracer_bkg(:,:,1,jp_fabm_chl2) + & 
     1288               &               tracer_bkg(:,:,1,jp_fabm_chl3) + & 
     1289               &               tracer_bkg(:,:,1,jp_fabm_chl4) 
    9281290#endif 
    9291291            CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot ) 
     
    9341296         ENDIF 
    9351297 
    936 #if defined key_medusa 
     1298#if defined key_medusa || defined key_fabm 
    9371299         ! Diatom chlorophyll 
    9381300         IF ( ln_slchldiainc ) THEN 
     1301#if defined key_medusa 
    9391302            zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd) 
     1303#elif defined key_fabm 
     1304            zbkg_chldia(:,:) = tracer_bkg(:,:,1,jp_fabm_chl1) 
     1305#endif 
    9401306            CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia ) 
    9411307         ELSE 
     
    9541320#endif 
    9551321 
     1322#if defined key_fabm 
     1323         ! Nanophytoplankton chlorophyll 
     1324         IF ( ln_slchlnaninc ) THEN 
     1325            zbkg_chlnan(:,:) = tracer_bkg(:,:,1,jp_fabm_chl2) 
     1326            CALL asm_bgc_unlog_2d( zbkg_chlnan, slchlnan_bkginc, zinc_chlnan ) 
     1327         ELSE 
     1328            zinc_chlnan(:,:) = 0.0 
     1329         ENDIF 
     1330 
     1331         ! Picophytoplankton chlorophyll 
     1332         IF ( ln_slchlpicinc ) THEN 
     1333            zbkg_chlpic(:,:) = tracer_bkg(:,:,1,jp_fabm_chl3) 
     1334            CALL asm_bgc_unlog_2d( zbkg_chlpic, slchlpic_bkginc, zinc_chlpic ) 
     1335         ELSE 
     1336            zinc_chlpic(:,:) = 0.0 
     1337         ENDIF 
     1338 
     1339         ! Dinoflagellate chlorophyll 
     1340         IF ( ln_slchldininc ) THEN 
     1341            zbkg_chldin(:,:) = tracer_bkg(:,:,1,jp_fabm_chl4) 
     1342            CALL asm_bgc_unlog_2d( zbkg_chldin, slchldin_bkginc, zinc_chldin ) 
     1343         ELSE 
     1344            zinc_chldin(:,:) = 0.0 
     1345         ENDIF 
     1346#endif 
     1347 
    9561348         ! Total phytoplankton carbon 
    9571349         IF ( ln_slphytotinc ) THEN 
     
    9881380         ! Select mixed layer 
    9891381         IF ( ll_asmdin ) THEN 
    990 #if defined key_top && ( defined key_hadocc || defined key_medusa ) 
     1382#if defined key_top && ( defined key_hadocc || defined key_medusa || defined key_fabm ) 
    9911383            CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', & 
    9921384               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' ) 
     
    10141406#endif 
    10151407            CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points] 
    1016                !zmld(:,:) = hmld_tref(:,:) 
    1017                CALL ctl_stop( ' hmld_tref mixed layer requested for phyto2d assimilation,', & 
    1018                   &           ' but is not available in this version' ) 
     1408               zmld(:,:) = hmld_tref(:,:) 
    10191409            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points] 
    10201410               zmld(:,:) = hmlpt(:,:) 
     
    10531443            &                         cchl_p_bkg(:,:,1),                   & 
    10541444            &                         tracer_bkg, phyto2d_balinc ) 
     1445#elif defined key_fabm 
     1446         CALL asm_phyto2d_bal_ersem(  (ln_slchltotinc .OR. ln_schltotinc), & 
     1447            &                         zinc_chltot,                         & 
     1448            &                         ln_slchldiainc,                      & 
     1449            &                         zinc_chldia,                         & 
     1450            &                         ln_slchlnaninc,                      & 
     1451            &                         zinc_chlnan,                         & 
     1452            &                         ln_slchlpicinc,                      & 
     1453            &                         zinc_chlpic,                         & 
     1454            &                         ln_slchldininc,                      & 
     1455            &                         zinc_chldin,                         & 
     1456            &                         zincper,                             & 
     1457            &                         rn_maxchlinc, ln_phytobal, zmld,     & 
     1458            &                         pgrow_avg_bkg, ploss_avg_bkg,        & 
     1459            &                         phyt_avg_bkg, mld_max_bkg,           & 
     1460            &                         totalk_bkg,                          & 
     1461            &                         tracer_bkg, phyto2d_balinc ) 
    10551462#else 
    10561463         CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', & 
     
    10981505                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt 
    10991506            END WHERE 
     1507#elif defined key_fabm 
     1508            WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & 
     1509               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp ) 
     1510               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & 
     1511                  &                           phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt 
     1512               trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + & 
     1513                  &                           phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt 
     1514            END WHERE 
    11001515#endif 
    11011516 
     
    11321547               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 
    11331548            END WHERE 
     1549#elif defined key_fabm 
     1550            WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & 
     1551               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp ) 
     1552               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & 
     1553                  &                           phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) 
     1554               trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) 
     1555            END WHERE 
    11341556#endif 
    11351557  
     
    11631585      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights 
    11641586      ! 
    1165       INTEGER                          :: ji, jj, jk   ! Loop counters 
    1166       INTEGER                          :: it           ! Index 
    1167       REAL(wp)                         :: zincwgt      ! IAU weight for timestep 
    1168       REAL(wp)                         :: zfrac_chn    ! Fraction of jpchn 
    1169       REAL(wp)                         :: zfrac_chd    ! Fraction of jpchd 
    1170       REAL(wp)                         :: zrat_phn_chn ! jpphn:jpchn ratio 
    1171       REAL(wp)                         :: zrat_phd_chd ! jpphd:jpchd ratio 
    1172       REAL(wp)                         :: zrat_pds_chd ! jppds:jpchd ratio 
    1173       REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc      ! Chlorophyll increments 
    1174       REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl      ! Chlorophyll background 
     1587      INTEGER                          :: ji, jj, jk    ! Loop counters 
     1588      INTEGER                          :: it            ! Index 
     1589      REAL(wp)                         :: zincwgt       ! IAU weight for timestep 
     1590      REAL(wp)                         :: zfrac_chn     ! Fraction of jpchn 
     1591      REAL(wp)                         :: zfrac_chd     ! Fraction of jpchd 
     1592      REAL(wp)                         :: zfrac_chl1    ! Fraction of jp_fabm_chl1 
     1593      REAL(wp)                         :: zfrac_chl2    ! Fraction of jp_fabm_chl2 
     1594      REAL(wp)                         :: zfrac_chl3    ! Fraction of jp_fabm_chl3 
     1595      REAL(wp)                         :: zfrac_chl4    ! Fraction of jp_fabm_chl4 
     1596      REAL(wp)                         :: zrat_phn_chn  ! jpphn:jpchn ratio 
     1597      REAL(wp)                         :: zrat_phd_chd  ! jpphd:jpchd ratio 
     1598      REAL(wp)                         :: zrat_pds_chd  ! jppds:jpchd ratio 
     1599      REAL(wp)                         :: zrat_p1c_chl1 ! jp_fabm_p1c:jp_fabm_chl1 ratio 
     1600      REAL(wp)                         :: zrat_p1n_chl1 ! jp_fabm_p1n:jp_fabm_chl1 ratio 
     1601      REAL(wp)                         :: zrat_p1p_chl1 ! jp_fabm_p1p:jp_fabm_chl1 ratio 
     1602      REAL(wp)                         :: zrat_p1s_chl1 ! jp_fabm_p1s:jp_fabm_chl1 ratio 
     1603      REAL(wp)                         :: zrat_p2c_chl2 ! jp_fabm_p2c:jp_fabm_chl2 ratio 
     1604      REAL(wp)                         :: zrat_p2n_chl2 ! jp_fabm_p2n:jp_fabm_chl2 ratio 
     1605      REAL(wp)                         :: zrat_p2p_chl2 ! jp_fabm_p2p:jp_fabm_chl2 ratio 
     1606      REAL(wp)                         :: zrat_p3c_chl3 ! jp_fabm_p3c:jp_fabm_chl3 ratio 
     1607      REAL(wp)                         :: zrat_p3n_chl3 ! jp_fabm_p3n:jp_fabm_chl3 ratio 
     1608      REAL(wp)                         :: zrat_p3p_chl3 ! jp_fabm_p3p:jp_fabm_chl3 ratio 
     1609      REAL(wp)                         :: zrat_p4c_chl4 ! jp_fabm_p4c:jp_fabm_chl4 ratio 
     1610      REAL(wp)                         :: zrat_p4n_chl4 ! jp_fabm_p4n:jp_fabm_chl4 ratio 
     1611      REAL(wp)                         :: zrat_p4p_chl4 ! jp_fabm_p4p:jp_fabm_chl4 ratio 
     1612      REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc       ! Chlorophyll increments 
     1613      REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl       ! Chlorophyll background 
    11751614      !!------------------------------------------------------------------------ 
    11761615 
     
    11891628#elif defined key_hadocc 
    11901629            bkg_chl(:,:,:) = chl_bkg(:,:,:) 
     1630#elif defined key_fabm 
     1631            bkg_chl(:,:,:) = tracer_bkg(:,:,:,jp_fabm_chl1) + & 
     1632               &             tracer_bkg(:,:,:,jp_fabm_chl2) + & 
     1633               &             tracer_bkg(:,:,:,jp_fabm_chl3) + & 
     1634               &             tracer_bkg(:,:,:,jp_fabm_chl4) 
    11911635#endif 
    11921636            DO jk = 1, jpk 
     
    12401684#elif defined key_hadocc 
    12411685         phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:) 
     1686#elif defined key_fabm 
     1687         ! Loop over each grid point partioning the increments based on existing ratios 
     1688         DO jk = 1, jpk 
     1689            DO jj = 1, jpj 
     1690               DO ji = 1, jpi 
     1691                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_chl1) > 0.0 ) .AND. & 
     1692                     & ( tracer_bkg(ji,jj,jk,jp_fabm_chl2) > 0.0 ) .AND. & 
     1693                     & ( tracer_bkg(ji,jj,jk,jp_fabm_chl3) > 0.0 ) .AND. & 
     1694                     & ( tracer_bkg(ji,jj,jk,jp_fabm_chl4) > 0.0 ) ) THEN 
     1695                     zfrac_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_chl1) / bkg_chl(ji,jj,jk) 
     1696                     zfrac_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_chl2) / bkg_chl(ji,jj,jk) 
     1697                     zfrac_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_chl3) / bkg_chl(ji,jj,jk) 
     1698                     zfrac_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_chl4) / bkg_chl(ji,jj,jk) 
     1699                     phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) = chl_inc(ji,jj,jk) * zfrac_chl1 
     1700                     phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) = chl_inc(ji,jj,jk) * zfrac_chl2 
     1701                     phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) = chl_inc(ji,jj,jk) * zfrac_chl3 
     1702                     phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) = chl_inc(ji,jj,jk) * zfrac_chl4 
     1703                     zrat_p1c_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1c) / tracer_bkg(ji,jj,jk,jp_fabm_chl1) 
     1704                     zrat_p1n_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1n) / tracer_bkg(ji,jj,jk,jp_fabm_chl1) 
     1705                     zrat_p1p_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1p) / tracer_bkg(ji,jj,jk,jp_fabm_chl1) 
     1706                     zrat_p1s_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_p1s) / tracer_bkg(ji,jj,jk,jp_fabm_chl1) 
     1707                     zrat_p2c_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2c) / tracer_bkg(ji,jj,jk,jp_fabm_chl2) 
     1708                     zrat_p2n_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2n) / tracer_bkg(ji,jj,jk,jp_fabm_chl2) 
     1709                     zrat_p2p_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_p2p) / tracer_bkg(ji,jj,jk,jp_fabm_chl2) 
     1710                     zrat_p3c_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3c) / tracer_bkg(ji,jj,jk,jp_fabm_chl3) 
     1711                     zrat_p3n_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3n) / tracer_bkg(ji,jj,jk,jp_fabm_chl3) 
     1712                     zrat_p3p_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_p3p) / tracer_bkg(ji,jj,jk,jp_fabm_chl3) 
     1713                     zrat_p4c_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4c) / tracer_bkg(ji,jj,jk,jp_fabm_chl4) 
     1714                     zrat_p4n_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4n) / tracer_bkg(ji,jj,jk,jp_fabm_chl4) 
     1715                     zrat_p4p_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_p4p) / tracer_bkg(ji,jj,jk,jp_fabm_chl4) 
     1716                     phyto3d_balinc(ji,jj,jk,jp_fabm_p1c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1c_chl1 
     1717                     phyto3d_balinc(ji,jj,jk,jp_fabm_p1n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1n_chl1 
     1718                     phyto3d_balinc(ji,jj,jk,jp_fabm_p1p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1p_chl1 
     1719                     phyto3d_balinc(ji,jj,jk,jp_fabm_p1s) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl1) * zrat_p1s_chl1 
     1720                     phyto3d_balinc(ji,jj,jk,jp_fabm_p2c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2c_chl2 
     1721                     phyto3d_balinc(ji,jj,jk,jp_fabm_p2n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2n_chl2 
     1722                     phyto3d_balinc(ji,jj,jk,jp_fabm_p2p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl2) * zrat_p2p_chl2 
     1723                     phyto3d_balinc(ji,jj,jk,jp_fabm_p3c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3c_chl3 
     1724                     phyto3d_balinc(ji,jj,jk,jp_fabm_p3n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3n_chl3 
     1725                     phyto3d_balinc(ji,jj,jk,jp_fabm_p3p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl3) * zrat_p3p_chl3 
     1726                     phyto3d_balinc(ji,jj,jk,jp_fabm_p4c) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4c_chl4 
     1727                     phyto3d_balinc(ji,jj,jk,jp_fabm_p4n) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4n_chl4 
     1728                     phyto3d_balinc(ji,jj,jk,jp_fabm_p4p) = phyto3d_balinc(ji,jj,jk,jp_fabm_chl4) * zrat_p4p_chl4 
     1729                  ENDIF 
     1730               END DO 
     1731            END DO 
     1732         END DO 
    12421733#else 
    12431734         CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', & 
     
    12851776                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt 
    12861777            END WHERE 
     1778#elif defined key_fabm 
     1779            WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & 
     1780               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp ) 
     1781               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & 
     1782                  &                           phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt 
     1783               trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + & 
     1784                  &                           phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt 
     1785            END WHERE 
    12871786#endif 
    12881787 
     
    13181817                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1) 
    13191818               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 
     1819            END WHERE 
     1820#elif defined key_fabm 
     1821            WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & 
     1822               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp ) 
     1823               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & 
     1824                  &                           phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) 
     1825               trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) 
    13201826            END WHERE 
    13211827#endif 
     
    14471953            &               pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) ) 
    14481954 
     1955#elif defined key_fabm 
     1956         ! Account for phytoplankton balancing if required 
     1957         IF ( ln_phytobal ) THEN 
     1958            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_o3c) + phyto2d_balinc(:,:,1,jp_fabm_o3c) 
     1959            alk_bkg_temp(:,:) = totalk_bkg(:,:,1)             + phyto2d_balinc(:,:,1,jp_fabm_o3ba) 
     1960         ELSE 
     1961            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_o3c) 
     1962            alk_bkg_temp(:,:) = totalk_bkg(:,:,1) 
     1963         ENDIF 
     1964 
     1965         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & 
     1966            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        & 
     1967            &               pco2_balinc(:,:,1,jp_fabm_o3c), pco2_balinc(:,:,1,jp_fabm_o3ba) ) 
     1968 
    14491969#else 
    14501970         CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', & 
     
    14541974         ! Select mixed layer 
    14551975         IF ( ll_asmdin ) THEN 
    1456 #if defined key_hadocc || defined key_medusa 
     1976#if defined key_hadocc || defined key_medusa || defined key_fabm 
    14571977            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', & 
    14581978               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' ) 
     
    14802000#endif 
    14812001            CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points] 
    1482                !zmld(:,:) = hmld_tref(:,:) 
    1483                CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', & 
    1484                   &           ' but is not available in this version' ) 
     2002               zmld(:,:) = hmld_tref(:,:) 
    14852003            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points] 
    14862004               zmld(:,:) = hmlpt(:,:) 
     
    15552073                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt 
    15562074            END WHERE 
     2075#elif defined key_fabm 
     2076            WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & 
     2077               &   trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp ) 
     2078               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & 
     2079                  &                           pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt 
     2080               trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + & 
     2081                  &                           pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt 
     2082            END WHERE 
    15572083#endif 
    15582084 
     
    15892115                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1) 
    15902116               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 
     2117            END WHERE 
     2118#elif defined key_fabm 
     2119            WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. & 
     2120               &   trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp ) 
     2121               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + & 
     2122                  &                           pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) 
     2123               trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) 
    15912124            END WHERE 
    15922125#endif 
     
    18542387 
    18552388         IF ( ln_pno3inc ) THEN 
    1856 #if defined key_hadocc || defined key_medusa 
     2389#if defined key_hadocc || defined key_medusa || defined key_fabm 
    18572390#if defined key_hadocc 
    18582391            it = jp_had_nut 
    18592392#elif defined key_medusa 
    18602393            it = jpdin 
     2394#elif defined key_fabm 
     2395            it = jp_fabm_n3n 
    18612396#endif 
    18622397            IF ( ln_phytobal ) THEN 
     
    18782413 
    18792414         IF ( ln_psi4inc ) THEN 
     2415#if defined key_medusa || defined key_fabm 
    18802416#if defined key_medusa 
    18812417            it = jpsil 
     2418#elif defined key_fabm 
     2419            it = jp_fabm_n5s 
     2420#endif 
    18822421            IF ( ln_phytobal ) THEN 
    18832422               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) 
     
    18982437 
    18992438         IF ( ln_pdicinc ) THEN 
    1900 #if defined key_hadocc || defined key_medusa 
     2439#if defined key_hadocc || defined key_medusa || defined key_fabm 
    19012440#if defined key_hadocc 
    19022441            it = jp_had_dic 
    19032442#elif defined key_medusa 
    19042443            it = jpdic 
     2444#elif defined key_fabm 
     2445            it = jp_fabm_o3c 
    19052446#endif 
    19062447            IF ( ln_phytobal ) THEN 
     
    19222463 
    19232464         IF ( ln_palkinc ) THEN 
    1924 #if defined key_hadocc || defined key_medusa 
     2465#if defined key_hadocc || defined key_medusa || defined key_fabm 
    19252466#if defined key_hadocc 
    19262467            it = jp_had_alk 
    19272468#elif defined key_medusa 
    19282469            it = jpalk 
     2470#elif defined key_fabm 
     2471            it = jp_fabm_o3ba 
    19292472#endif 
    19302473            IF ( ln_phytobal ) THEN 
     
    19462489 
    19472490         IF ( ln_po2inc ) THEN 
     2491#if defined key_medusa || defined key_fabm 
    19482492#if defined key_medusa 
    19492493            it = jpoxy 
     2494#elif defined key_fabm 
     2495            it = jp_fabm_o2o 
     2496#endif 
    19502497            IF ( ln_phytobal ) THEN 
    19512498               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) 
     
    20042551                  trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt 
    20052552               END WHERE 
     2553#elif defined key_fabm 
     2554               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. & 
     2555                  &   trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp ) 
     2556                  trn(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt 
     2557                  trb(:,:,:,jp_fabm_n3n) = trb(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt 
     2558               END WHERE 
    20062559#else 
    20072560               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) 
     
    20152568                  trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt 
    20162569                  trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt 
     2570               END WHERE 
     2571#elif defined key_fabm 
     2572               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. & 
     2573                  &   trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp ) 
     2574                  trn(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt 
     2575                  trb(:,:,:,jp_fabm_n5s) = trb(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt 
    20172576               END WHERE 
    20182577#else 
     
    20342593                  trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt 
    20352594               END WHERE 
     2595#elif defined key_fabm 
     2596               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. & 
     2597                  &   trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp ) 
     2598                  trn(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt 
     2599                  trb(:,:,:,jp_fabm_o3c) = trb(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt 
     2600               END WHERE 
    20362601#else 
    20372602               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) 
     
    20522617                  trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt 
    20532618               END WHERE 
     2619#elif defined key_fabm 
     2620               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. & 
     2621                  &   trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp ) 
     2622                  trn(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt 
     2623                  trb(:,:,:,jp_fabm_o3ba) = trb(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt 
     2624               END WHERE 
    20542625#else 
    20552626               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) 
     
    20632634                  trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt 
    20642635                  trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt 
     2636               END WHERE 
     2637#elif defined key_fabm 
     2638               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. & 
     2639                  &   trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp ) 
     2640                  trn(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt 
     2641                  trb(:,:,:,jp_fabm_o2o) = trb(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt 
    20652642               END WHERE 
    20662643#else 
     
    20912668            ! Initialize the now fields with the background + increment 
    20922669            ! Background currently is what the model is initialised with 
    2093 #if defined key_hadocc 
    2094             CALL ctl_warn( ' Doing direct initialisation of HadOCC with 3D BGC assimilation', & 
     2670            CALL ctl_warn( ' Doing direct initialisation with 3D BGC assimilation', & 
    20952671               &           ' Background state is taken from model rather than background file' ) 
    2096 #elif defined key_medusa 
    2097             CALL ctl_warn( ' Doing direct initialisation of MEDUSA with 3D BGC assimilation', & 
    2098                &           ' Background state is taken from model rather than background file' ) 
    2099 #endif 
    21002672 
    21012673            IF ( ln_pno3inc ) THEN 
     
    21122684                  trb(:,:,:,jpdin) = trn(:,:,:,jpdin) 
    21132685               END WHERE 
     2686#elif defined key_fabm 
     2687               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. & 
     2688                  &   trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) > 0.0_wp ) 
     2689                  trn(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) + pno3_bkginc(:,:,:) 
     2690                  trb(:,:,:,jp_fabm_n3n) = trn(:,:,:,jp_fabm_n3n) 
     2691               END WHERE 
    21142692#else 
    21152693               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) 
     
    21232701                  trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) 
    21242702                  trb(:,:,:,jpsil) = trn(:,:,:,jpsil) 
     2703               END WHERE 
     2704#elif defined key_fabm 
     2705               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. & 
     2706                  &   trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) > 0.0_wp ) 
     2707                  trn(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) + psi4_bkginc(:,:,:) 
     2708                  trb(:,:,:,jp_fabm_n5s) = trn(:,:,:,jp_fabm_n5s) 
    21252709               END WHERE 
    21262710#else 
     
    21422726                  trb(:,:,:,jpdic) = trn(:,:,:,jpdic) 
    21432727               END WHERE 
     2728#elif defined key_fabm 
     2729               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. & 
     2730                  &   trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) > 0.0_wp ) 
     2731                  trn(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) + pdic_bkginc(:,:,:) 
     2732                  trb(:,:,:,jp_fabm_o3c) = trn(:,:,:,jp_fabm_o3c) 
     2733               END WHERE 
    21442734#else 
    21452735               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) 
     
    21602750                  trb(:,:,:,jpalk) = trn(:,:,:,jpalk) 
    21612751               END WHERE 
     2752#elif defined key_fabm 
     2753               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. & 
     2754                  &   trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) > 0.0_wp ) 
     2755                  trn(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) + palk_bkginc(:,:,:) 
     2756                  trb(:,:,:,jp_fabm_o3ba) = trn(:,:,:,jp_fabm_o3ba) 
     2757               END WHERE 
    21622758#else 
    21632759               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) 
     
    21712767                  trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) 
    21722768                  trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy) 
     2769               END WHERE 
     2770#elif defined key_fabm 
     2771               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. & 
     2772                  &   trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) > 0.0_wp ) 
     2773                  trn(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) + po2_bkginc(:,:,:) 
     2774                  trb(:,:,:,jp_fabm_o2o) = trn(:,:,:,jp_fabm_o2o) 
    21732775               END WHERE 
    21742776#else 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90

    r10574 r10622  
    5252   USE asminc, ONLY: ln_avgbkg 
    5353#if defined key_top 
    54    USE asmbgc, ONLY: asm_bgc_bkg_wri 
     54   USE asmbgc, ONLY: asm_bgc_bkg_alloc, & 
     55      &              asm_bgc_bkg_tavg,  & 
     56      &              asm_bgc_bkg_wri 
    5557#endif 
    5658   IMPLICIT NONE 
     
    140142          
    141143         numtimes_tavg = REAL ( nitavgbkg_r -  nn_it000 + 1 ) 
    142       ENDIF    
     144      ENDIF 
     145       
     146#if defined key_top 
     147      ! Allocate BGC average arrays whatever, to save code repetition later 
     148      IF ( kt == ( nn_it000 - 1) ) THEN 
     149         CALL asm_bgc_bkg_alloc 
     150      ENDIF 
     151#endif 
    143152 
    144153      ! If creating an averaged assim bkg, sum the contribution every timestep 
     
    157166#if defined key_zdftke 
    158167         en_tavg(:,:,:)       = en_tavg(:,:,:) + en(:,:,:) / numtimes_tavg 
     168#endif 
     169#if defined key_top 
     170         CALL asm_bgc_bkg_tavg( kt, numtimes_tavg ) 
    159171#endif 
    160172      ENDIF 
     
    226238             
    227239#if defined key_top 
    228             CALL asm_bgc_bkg_wri( kt, inum ) 
     240            CALL asm_bgc_bkg_wri( kt, inum, ln_avgbkg ) 
    229241#endif 
    230242            CALL iom_close( inum ) 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r10574 r10622  
    167167         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
    168168         &                 ln_phytobal, ln_slchltotinc, ln_slchldiainc,    & 
     169         &                 ln_slchlnaninc, ln_slchlpicinc, ln_slchldininc, & 
    169170         &                 ln_slchlnoninc, ln_schltotinc, ln_slphytotinc,  & 
    170171         &                 ln_slphydiainc, ln_slphynoninc, ln_spco2inc,    & 
     
    228229         WRITE(numout,*) '      Logical switch for applying slchldia increments     ln_slchldiainc = ', ln_slchldiainc 
    229230         WRITE(numout,*) '      Logical switch for applying slchlnon increments     ln_slchlnoninc = ', ln_slchlnoninc 
     231         WRITE(numout,*) '      Logical switch for applying slchlnan increments     ln_slchlnaninc = ', ln_slchlnaninc 
     232         WRITE(numout,*) '      Logical switch for applying slchlpic increments     ln_slchlpicinc = ', ln_slchlpicinc 
     233         WRITE(numout,*) '      Logical switch for applying slchldin increments     ln_slchldininc = ', ln_slchldininc 
    230234         WRITE(numout,*) '      Logical switch for applying schltot increments       ln_schltotinc = ', ln_schltotinc 
    231235         WRITE(numout,*) '      Logical switch for applying slphytot increments     ln_slphytotinc = ', ln_slphytotinc 
     
    295299      ENDIF 
    296300      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & 
     301         & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. & 
    297302         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & 
    298303         & ln_slphynoninc .OR. ln_spco2inc    .OR. ln_sfco2inc    .OR. & 
     
    13531358      ! Ocean colour variables first 
    13541359      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & 
     1360         & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. & 
    13551361         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & 
    13561362         & ln_slphynoninc ) THEN 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_ersem.F90

    r10574 r10622  
    1 MODULE asmphyto2dbal_medusa 
     1MODULE asmphyto2dbal_ersem 
    22   !!====================================================================== 
    3    !!                       ***  MODULE asmphyto2dbal_medusa  *** 
    4    !! Calculate increments to MEDUSA based on surface phyto2d increments 
     3   !!                       ***  MODULE asmphyto2dbal_ersem  *** 
     4   !! Calculate increments to ERSEM based on surface phyto2d increments 
    55   !! 
    66   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al. 
     
    1010   !! 
    1111   !!====================================================================== 
    12    !! History :  3.6  ! 2017-08 (D. Ford)  Adapted from asmphyto2dbal_hadocc 
     12   !! History :  3.6  ! 2019-01 (D. Ford)  Adapted from asmphyto2dbal_medusa 
    1313   !!---------------------------------------------------------------------- 
    14 #if defined key_asminc && defined key_medusa 
     14#if defined key_asminc && defined key_fabm 
    1515   !!---------------------------------------------------------------------- 
    1616   !! 'key_asminc'          : assimilation increment interface 
    17    !! 'key_medusa'          : MEDUSA model 
     17   !! 'key_fabm'            : FABM-ERSEM model 
    1818   !!---------------------------------------------------------------------- 
    19    !! asm_phyto2d_bal_medusa : routine to calculate increments to MEDUSA 
     19   !! asm_phyto2d_bal_ersem : routine to calculate increments to ERSEM 
    2020   !!---------------------------------------------------------------------- 
    2121   USE par_kind,      ONLY: wp             ! kind parameters 
    2222   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes 
    2323   USE dom_oce,       ONLY: gdepw_n        ! domain information 
    24    USE zdftmx,        ONLY: ln_tmx_itf, &  ! Indonesian Throughflow 
    25       &                     mask_itf       ! tidal mixing mask 
    2624   USE iom                                 ! i/o 
    27    USE sms_medusa                          ! MEDUSA parameters 
    28    USE par_medusa                          ! MEDUSA parameters 
     25   USE par_fabm                            ! FABM-ERSEM parameters 
    2926   USE par_trc,       ONLY: jptra          ! Tracer parameters 
    3027   USE bioanalysis                         ! Nitrogen balancing 
     
    3330   PRIVATE                    
    3431 
    35    PUBLIC asm_phyto2d_bal_medusa 
     32   PUBLIC asm_phyto2d_bal_ersem 
    3633 
    3734   ! Default values for biological assimilation parameters 
     
    6865CONTAINS 
    6966 
    70    SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      & 
     67   SUBROUTINE asm_phyto2d_bal_ersem( ld_chltot,                      & 
    7168      &                               pinc_chltot,                    & 
    7269      &                               ld_chldia,                      & 
    7370      &                               pinc_chldia,                    & 
    74       &                               ld_chlnon,                      & 
    75       &                               pinc_chlnon,                    & 
    76       &                               ld_phytot,                      & 
    77       &                               pinc_phytot,                    & 
    78       &                               ld_phydia,                      & 
    79       &                               pinc_phydia,                    & 
    80       &                               ld_phynon,                      & 
    81       &                               pinc_phynon,                    & 
     71      &                               ld_chlnan,                      & 
     72      &                               pinc_chlnan,                    & 
     73      &                               ld_chlpic,                      & 
     74      &                               pinc_chlpic,                    & 
     75      &                               ld_chldin,                      & 
     76      &                               pinc_chldin,                    & 
    8277      &                               pincper,                        & 
    8378      &                               p_maxchlinc, ld_phytobal, pmld, & 
    8479      &                               pgrow_avg_bkg, ploss_avg_bkg,   & 
    8580      &                               phyt_avg_bkg, mld_max_bkg,      & 
     81      &                               totalk_bkg,                     & 
    8682      &                               tracer_bkg, phyto2d_balinc ) 
    8783      !!--------------------------------------------------------------------------- 
    88       !!                    ***  ROUTINE asm_phyto2d_bal_medusa  *** 
     84      !!                    ***  ROUTINE asm_phyto2d_bal_ersem  *** 
    8985      !! 
    90       !! ** Purpose :   calculate increments to MEDUSA from 2d phytoplankton increments 
     86      !! ** Purpose :   calculate increments to ERSEM from 2d phytoplankton increments 
    9187      !! 
    92       !! ** Method  :   average up MEDUSA to look like HadOCC 
     88      !! ** Method  :   EITHER (ld_phytobal == .TRUE.): 
     89      !!                average up ERSEM to look like HadOCC 
    9390      !!                call nitrogen balancing scheme 
    9491      !!                separate back out to MEDUSA 
     92      !!                OR (ld_phytobal == .FALSE.): 
     93      !!                calculate increments to maintain background stoichiometry 
    9594      !! 
    9695      !! ** Action  :   populate phyto2d_balinc 
     
    9897      !! References :   Hemmings et al., 2008, J. Mar. Res. 
    9998      !!                Ford et al., 2012, Ocean Sci. 
     99      !!                Skakala et al., 2018, JGR 
    100100      !!--------------------------------------------------------------------------- 
    101101      !! 
     
    104104      LOGICAL,  INTENT(in   )                               :: ld_chldia      ! Assim chldia y/n 
    105105      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chldia    ! chldia increments 
    106       LOGICAL,  INTENT(in   )                               :: ld_chlnon      ! Assim chlnon y/n 
    107       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chlnon    ! chlnon increments 
    108       LOGICAL,  INTENT(in   )                               :: ld_phytot      ! Assim phytot y/n 
    109       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phytot    ! phytot increments 
    110       LOGICAL,  INTENT(in   )                               :: ld_phydia      ! Assim phydia y/n 
    111       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phydia    ! phydia increments 
    112       LOGICAL,  INTENT(in   )                               :: ld_phynon      ! Assim phynon y/n 
    113       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phynon    ! phynon increments 
     106      LOGICAL,  INTENT(in   )                               :: ld_chlnan      ! Assim chlnan y/n 
     107      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chlnan    ! chlnan increments 
     108      LOGICAL,  INTENT(in   )                               :: ld_chlpic      ! Assim chlpic y/n 
     109      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chlpic    ! chlpic increments 
     110      LOGICAL,  INTENT(in   )                               :: ld_chldin      ! Assim chldin y/n 
     111      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chldin    ! chldin increments 
    114112      REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period 
    115113      REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment 
     
    120118      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto 
    121119      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD 
     120      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)       :: totalk_bkg     ! Total alkalinity 
    122121      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables 
    123122      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments 
     
    126125      INTEGER                                               :: jkmax          ! Loop index 
    127126      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices 
     127      REAL(wp)                                              :: zmassc         ! Carbon molar mass 
     128      REAL(wp)                                              :: zmassn         ! Nitrogen molar mass 
     129      REAL(wp)                                              :: z4qnc          ! Z4/qnc (mesozoo N:C) 
    128130      REAL(wp)                                              :: n2be_p         ! N:biomass for total phy 
    129131      REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo 
    130132      REAL(wp)                                              :: n2be_d         ! N:biomass for detritus 
    131       REAL(wp)                                              :: zfrac          ! Fraction 
    132       REAL(wp)                                              :: zfrac_chn      ! Fraction of jpchn 
    133       REAL(wp)                                              :: zfrac_chd      ! Fraction of jpchd 
    134       REAL(wp)                                              :: zfrac_phn      ! Fraction of jpphn 
    135       REAL(wp)                                              :: zfrac_phd      ! Fraction of jpphd 
    136       REAL(wp)                                              :: zfrac_zmi      ! Fraction of jpzmi 
    137       REAL(wp)                                              :: zfrac_zme      ! Fraction of jpzme 
    138       REAL(wp)                                              :: zrat_pds_phd   ! Ratio of jppds:jpphd 
    139       REAL(wp)                                              :: zrat_chd_phd   ! Ratio of jpchd:jpphd 
    140       REAL(wp)                                              :: zrat_chn_phn   ! Ratio of jpchn:jpphn 
    141       REAL(wp)                                              :: zrat_phn_chn   ! Ratio of jpphn:jpchn 
    142       REAL(wp)                                              :: zrat_phd_chd   ! Ratio of jpphd:jpchd 
    143       REAL(wp)                                              :: zrat_pds_chd   ! Ratio of jppds:jpchd 
    144       REAL(wp)                                              :: zrat_dtc_det   ! Ratio of jpdtc:jpdet 
     133      REAL(wp)                                              :: zfrac          ! Fractions 
     134      REAL(wp)                                              :: zfrac_chl1     !  
     135      REAL(wp)                                              :: zfrac_chl2     !  
     136      REAL(wp)                                              :: zfrac_chl3     !  
     137      REAL(wp)                                              :: zfrac_chl4     !  
     138      REAL(wp)                                              :: zfrac_p1n      !  
     139      REAL(wp)                                              :: zfrac_p2n      !  
     140      REAL(wp)                                              :: zfrac_p3n      !  
     141      REAL(wp)                                              :: zfrac_p4n      !  
     142      REAL(wp)                                              :: zfrac_z4n      !  
     143      REAL(wp)                                              :: zfrac_z5n      !  
     144      REAL(wp)                                              :: zfrac_z6n      !  
     145      REAL(wp)                                              :: zfrac_n3n      !  
     146      REAL(wp)                                              :: zfrac_n4n      !  
     147      REAL(wp)                                              :: zfrac_r4n      !  
     148      REAL(wp)                                              :: zfrac_r6n      !  
     149      REAL(wp)                                              :: zfrac_r8n      !  
     150      REAL(wp)                                              :: zrat_chl1_p1n  ! Ratios 
     151      REAL(wp)                                              :: zrat_p1c_p1n   !  
     152      REAL(wp)                                              :: zrat_p1p_p1n   !  
     153      REAL(wp)                                              :: zrat_p1s_p1n   !  
     154      REAL(wp)                                              :: zrat_chl2_p2n  !  
     155      REAL(wp)                                              :: zrat_p2c_p2n   !  
     156      REAL(wp)                                              :: zrat_p2p_p2n   !  
     157      REAL(wp)                                              :: zrat_chl3_p3n  !  
     158      REAL(wp)                                              :: zrat_p3c_p3n   !  
     159      REAL(wp)                                              :: zrat_p3p_p3n   !  
     160      REAL(wp)                                              :: zrat_chl4_p4n  !  
     161      REAL(wp)                                              :: zrat_p4c_p4n   !  
     162      REAL(wp)                                              :: zrat_p4p_p4n   !  
     163      REAL(wp)                                              :: zrat_z4c_z4n   !  
     164      REAL(wp)                                              :: zrat_z5c_z5n   !  
     165      REAL(wp)                                              :: zrat_z5p_z5n   !  
     166      REAL(wp)                                              :: zrat_z6c_z6n   !  
     167      REAL(wp)                                              :: zrat_z6p_z6n   !  
     168      REAL(wp)                                              :: zrat_r4c_r4n   !  
     169      REAL(wp)                                              :: zrat_r4p_r4n   !  
     170      REAL(wp)                                              :: zrat_r6c_r6n   !  
     171      REAL(wp)                                              :: zrat_r6p_r6n   !  
     172      REAL(wp)                                              :: zrat_r6s_r6n   !  
     173      REAL(wp)                                              :: zrat_r8c_r8n   !  
     174      REAL(wp)                                              :: zrat_r8p_r8n   !  
     175      REAL(wp)                                              :: zrat_r8s_r8n   !  
     176      REAL(wp)                                              :: zrat_p1c_chl1  !  
     177      REAL(wp)                                              :: zrat_p1n_chl1  !  
     178      REAL(wp)                                              :: zrat_p1p_chl1  !  
     179      REAL(wp)                                              :: zrat_p1s_chl1  !  
     180      REAL(wp)                                              :: zrat_p2c_chl2  !  
     181      REAL(wp)                                              :: zrat_p2n_chl2  !  
     182      REAL(wp)                                              :: zrat_p2p_chl2  !  
     183      REAL(wp)                                              :: zrat_p3c_chl3  !  
     184      REAL(wp)                                              :: zrat_p3n_chl3  !  
     185      REAL(wp)                                              :: zrat_p3p_chl3  !  
     186      REAL(wp)                                              :: zrat_p4c_chl4  !  
     187      REAL(wp)                                              :: zrat_p4n_chl4  !  
     188      REAL(wp)                                              :: zrat_p4p_chl4  !  
    145189      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy 
    146190      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters 
     
    150194      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics 
    151195      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics 
     196      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chltot_temp 
    152197      !!--------------------------------------------------------------------------- 
     198       
     199      ! Set parameters 
     200      zmassc = 12.01 
     201      zmassn = 14.01 
     202      z4qnc  = 0.0126 
     203      !z4qnc  = model%state_variables(jp_fabm_z4c)%parameters%qnc%value 
     204      !z4qnc  = get_property_by_name(model%state_variables(jp_fabm_z4c)%parameters, 'qnc') 
     205      IF (lwp) WRITE(numout,*) 'z4qnc = ', z4qnc 
    153206 
    154207      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value 
     
    160213               END DO 
    161214            END DO 
    162          ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN 
     215         ELSE 
    163216            DO jj = 1, jpj 
    164217               DO ji = 1, jpi 
    165                   pinc_chltot(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) 
    166                   pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) ) 
    167                   IF ( pinc_chltot(ji,jj) .NE. ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) ) ) THEN 
    168                      zfrac = pinc_chltot(ji,jj) / ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) ) 
    169                      pinc_chldia(ji,jj) = pinc_chldia(ji,jj) * zfrac 
    170                      pinc_chlnon(ji,jj) = pinc_chlnon(ji,jj) * zfrac 
    171                   ENDIF 
    172                END DO 
    173             END DO 
    174          ELSE IF ( ld_chldia ) THEN 
    175             DO jj = 1, jpj 
    176                DO ji = 1, jpi 
    177                   pinc_chldia(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chldia(ji,jj), p_maxchlinc ) ) 
    178                   pinc_chltot(ji,jj) = pinc_chldia(ji,jj) 
    179                END DO 
    180             END DO 
    181          ELSE IF ( ld_chlnon ) THEN 
    182             DO jj = 1, jpj 
    183                DO ji = 1, jpi 
    184                   pinc_chlnon(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chlnon(ji,jj), p_maxchlinc ) ) 
    185                   pinc_chltot(ji,jj) = pinc_chlnon(ji,jj) 
     218                  IF ( ld_chldia .AND. ld_chlnan .AND. ld_chlpic .AND. ld_chldin ) THEN 
     219                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + & 
     220                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj) 
     221                  ELSE IF ( ld_chldia .AND. ld_chlnan .AND. ld_chlpic ) THEN 
     222                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + & 
     223                        &                      pinc_chlpic(ji,jj) 
     224                  ELSE IF ( ld_chldia .AND. ld_chlnan .AND. ld_chldin ) THEN 
     225                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + & 
     226                        &                      pinc_chldin(ji,jj) 
     227                  ELSE IF ( ld_chldia .AND. ld_chlpic .AND. ld_chldin ) THEN 
     228                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + & 
     229                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj) 
     230                  ELSE IF ( ld_chlnan .AND. ld_chlpic .AND. ld_chldin ) THEN 
     231                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + & 
     232                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj) 
     233                  ELSE IF ( ld_chldia .AND. ld_chlnan ) THEN 
     234                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) 
     235                  ELSE IF ( ld_chldia .AND. ld_chlpic ) THEN 
     236                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlpic(ji,jj) 
     237                  ELSE IF ( ld_chldia .AND. ld_chldin ) THEN 
     238                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chldin(ji,jj) 
     239                  ELSE IF ( ld_chlnan .AND. ld_chlpic ) THEN 
     240                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + pinc_chlpic(ji,jj) 
     241                  ELSE IF ( ld_chlnan .AND. ld_chldin ) THEN 
     242                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + pinc_chldin(ji,jj) 
     243                  ELSE IF ( ld_chlpic .AND. ld_chldin ) THEN 
     244                     pinc_chltot_temp(ji,jj) = pinc_chlpic(ji,jj) + pinc_chldin(ji,jj) 
     245                  ELSE IF ( ld_chldia ) THEN 
     246                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) 
     247                  ELSE IF ( ld_chlnan ) THEN 
     248                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) 
     249                  ELSE IF ( ld_chlpic ) THEN 
     250                     pinc_chltot_temp(ji,jj) = pinc_chlpic(ji,jj) 
     251                  ELSE IF ( ld_chldin ) THEN 
     252                     pinc_chltot_temp(ji,jj) = pinc_chldin(ji,jj) 
     253                  ENDIF 
     254                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot_temp(ji,jj), p_maxchlinc ) ) 
     255                  IF ( pinc_chltot(ji,jj) .NE. pinc_chltot_temp(ji,jj) ) THEN 
     256                     zfrac = pinc_chltot(ji,jj) / pinc_chltot_temp(ji,jj) 
     257                     IF ( ld_chldia ) THEN 
     258                        pinc_chldia(ji,jj) = pinc_chldia(ji,jj) * zfrac 
     259                     ENDIF 
     260                     IF ( ld_chlnan ) THEN 
     261                        pinc_chlnan(ji,jj) = pinc_chlnan(ji,jj) * zfrac 
     262                     ENDIF 
     263                     IF ( ld_chlpic ) THEN 
     264                        pinc_chlpic(ji,jj) = pinc_chlpic(ji,jj) * zfrac 
     265                     ENDIF 
     266                     IF ( ld_chldin ) THEN 
     267                        pinc_chldin(ji,jj) = pinc_chldin(ji,jj) * zfrac 
     268                     ENDIF 
     269                  ENDIF 
    186270               END DO 
    187271            END DO 
    188272         ENDIF 
    189273      ENDIF 
    190  
    191       IF ( ld_phytot .OR. ld_phydia .OR. ld_phynon ) THEN 
    192          CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' ) 
    193       ENDIF 
     274       
     275      ! Initialise balancing increments 
     276      phyto2d_balinc(:,:,:,:) = 0.0 
    194277       
    195278      IF ( ld_phytobal ) THEN   ! Nitrogen balancing 
     
    197280         ! Set up model parameters to be passed into Hemmings balancing routine. 
    198281         ! For now these are hardwired to the standard HadOCC parameter values 
    199          ! (except C:N ratios) as this is what the scheme was developed for. 
    200          ! Obviously, HadOCC and MEDUSA are rather different models, so this 
     282         ! as this is what the scheme was developed for. 
     283         ! Obviously, HadOCC and ERSEM are rather different models, so this 
    201284         ! isn't ideal, but there's not always direct analogues between the two 
    202285         ! parameter sets, so it's the easiest way to get something running. 
     
    211294         modparm(8)  = 0.05                              ! z_mort_1 
    212295         modparm(9)  = 1.0                               ! z_mort_2 
    213          modparm(10) = ( xthetapn  + xthetapd  ) / 2.0   ! c2n_p 
    214          modparm(11) = ( xthetazmi + xthetazme ) / 2.0   ! c2n_z 
    215          modparm(12) = xthetad                           ! c2n_d 
     296         modparm(10) = 6.625                             ! c2n_p 
     297         modparm(11) = 5.625                             ! c2n_z 
     298         modparm(12) = 7.5                               ! c2n_d 
    216299         modparm(13) = 0.01                              ! graze_threshold 
    217300         modparm(14) = 2.0                               ! holling_coef 
     
    250333 
    251334         ! Set background state 
    252          bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jpdin) 
    253          bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jpphn) + tracer_bkg(:,:,:,jpphd) 
    254          bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jpzmi) + tracer_bkg(:,:,:,jpzme) 
    255          bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jpdet) 
    256          bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jpdic) 
    257          bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jpalk) 
     335         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jp_fabm_n3n) + & 
     336            &                        tracer_bkg(:,:,:,jp_fabm_n4n) 
     337         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jp_fabm_p1n) + & 
     338            &                        tracer_bkg(:,:,:,jp_fabm_p2n) + & 
     339            &                        tracer_bkg(:,:,:,jp_fabm_p3n) + & 
     340            &                        tracer_bkg(:,:,:,jp_fabm_p4n) 
     341         bstate(:,:,:,i_tracer(3)) = (tracer_bkg(:,:,:,jp_fabm_z4c) * z4qnc) + & 
     342            &                        tracer_bkg(:,:,:,jp_fabm_z5n) + & 
     343            &                        tracer_bkg(:,:,:,jp_fabm_z6n) 
     344         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jp_fabm_r4n) + & 
     345            &                        tracer_bkg(:,:,:,jp_fabm_r6n) + & 
     346            &                        tracer_bkg(:,:,:,jp_fabm_r8n) 
     347         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jp_fabm_o3c) 
     348         bstate(:,:,:,i_tracer(6)) = totalk_bkg(:,:,:) 
    258349 
    259350         ! Calculate carbon to chlorophyll ratio for combined phytoplankton 
    260          ! and nitrogen to biomass equivalent for PZD 
    261          ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA 
     351         ! and nitrogen to biomass equivalent for PZD (hardwire as per HadOCC) 
    262352         cchl_p(:,:) = 0.0 
    263353         DO jj = 1, jpj 
    264354            DO ji = 1, jpi 
    265                IF ( ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) ) .GT. 0.0 ) THEN 
    266                   cchl_p(ji,jj) = xmassc * ( ( tracer_bkg(ji,jj,1,jpphn) * xthetapn ) +      & 
    267                      &                       ( tracer_bkg(ji,jj,1,jpphd) * xthetapd )   ) /  & 
    268                      &            ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) ) 
     355               IF ( ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + tracer_bkg(ji,jj,1,jp_fabm_chl2) + & 
     356                  &   tracer_bkg(ji,jj,1,jp_fabm_chl3) + tracer_bkg(ji,jj,1,jp_fabm_chl4) ) .GT. 0.0 ) THEN 
     357                  cchl_p(ji,jj) = zmassc * ( tracer_bkg(ji,jj,1,jp_fabm_p1c) +      & 
     358                     &                       tracer_bkg(ji,jj,1,jp_fabm_p2c) +      & 
     359                     &                       tracer_bkg(ji,jj,1,jp_fabm_p3c) +      & 
     360                     &                       tracer_bkg(ji,jj,1,jp_fabm_p4c)   ) /  & 
     361                     &            ( tracer_bkg(ji,jj,1,jp_fabm_chl1) +             & 
     362                     &              tracer_bkg(ji,jj,1,jp_fabm_chl2) +             & 
     363                     &              tracer_bkg(ji,jj,1,jp_fabm_chl3) +             & 
     364                     &              tracer_bkg(ji,jj,1,jp_fabm_chl4)   ) 
    269365               ENDIF 
    270366            END DO 
    271367         END DO 
    272          n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) ) 
    273          n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) ) 
    274          n2be_d = 14.01 + ( xmassc * xthetad ) 
     368         n2be_p = zmassn + ( zmassc * modparm(10) ) 
     369         n2be_z = zmassn + ( zmassc * modparm(11) ) 
     370         n2be_d = zmassn + ( zmassc * modparm(12) ) 
    275371 
    276372         ! Call nitrogen balancing routine 
     
    288384          
    289385         ! Loop over each grid point partioning the increments 
    290          phyto2d_balinc(:,:,:,:) = 0.0 
    291386         DO jk = 1, jpk 
    292387            DO jj = 1, jpj 
     
    294389 
    295390                  ! Phytoplankton 
    296                   IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. & 
    297                      & ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) .AND. & 
     391                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) > 0.0 ) .AND. & 
     392                     & ( tracer_bkg(ji,jj,jk,jp_fabm_p2n) > 0.0 ) .AND. & 
     393                     & ( tracer_bkg(ji,jj,jk,jp_fabm_p3n) > 0.0 ) .AND. & 
     394                     & ( tracer_bkg(ji,jj,jk,jp_fabm_p4n) > 0.0 ) .AND. & 
    298395                     & ( pinc_chltot(ji,jj) /= 0.0 ) ) THEN 
    299396                     IF ( ld_chltot ) THEN 
    300397                        ! Phytoplankton nitrogen split up based on existing ratios 
    301                         zfrac_phn = tracer_bkg(ji,jj,jk,jpphn) / & 
    302                            &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd)) 
    303                         zfrac_phd = tracer_bkg(ji,jj,jk,jpphd) / & 
    304                            &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd)) 
    305                      ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN 
     398                        zfrac_p1n = tracer_bkg(ji,jj,jk,jp_fabm_p1n) / & 
     399                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + & 
     400                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + & 
     401                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + & 
     402                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) ) 
     403                        zfrac_p2n = tracer_bkg(ji,jj,jk,jp_fabm_p2n) / & 
     404                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + & 
     405                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + & 
     406                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + & 
     407                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) ) 
     408                        zfrac_p3n = tracer_bkg(ji,jj,jk,jp_fabm_p3n) / & 
     409                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + & 
     410                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + & 
     411                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + & 
     412                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) ) 
     413                        zfrac_p4n = tracer_bkg(ji,jj,jk,jp_fabm_p4n) / & 
     414                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_p1n) + & 
     415                           &          tracer_bkg(ji,jj,jk,jp_fabm_p2n) + & 
     416                           &          tracer_bkg(ji,jj,jk,jp_fabm_p3n) + & 
     417                           &          tracer_bkg(ji,jj,jk,jp_fabm_p4n) ) 
     418                     ELSE 
    306419                        ! Phytoplankton nitrogen split up based on assimilation increments 
    307                         zfrac_phn = pinc_chlnon(ji,jj) / pinc_chltot(ji,jj) 
    308                         zfrac_phd = pinc_chldia(ji,jj) / pinc_chltot(ji,jj) 
     420                        zfrac_p1n = pinc_chldia(ji,jj) / pinc_chltot(ji,jj) 
     421                        zfrac_p2n = pinc_chlnan(ji,jj) / pinc_chltot(ji,jj) 
     422                        zfrac_p3n = pinc_chlpic(ji,jj) / pinc_chltot(ji,jj) 
     423                        zfrac_p4n = pinc_chldin(ji,jj) / pinc_chltot(ji,jj) 
    309424                     ENDIF 
    310  
    311                      ! Phytoplankton silicate split up based on existing ratios 
    312                      zrat_pds_phd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpphd) 
    313425                      
    314                      ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen 
    315                      ! Not using pinc_chltot directly as it's only 2D 
    316                      ! This method should give same results at surface as splitting pinc_chltot would 
    317                      zrat_chn_phn = tracer_bkg(ji,jj,jk,jpchn) / tracer_bkg(ji,jj,jk,jpphn) 
    318                      zrat_chd_phd = tracer_bkg(ji,jj,jk,jpchd) / tracer_bkg(ji,jj,jk,jpphd) 
     426                     ! Other phytoplankton variables split up based on existing ratios with nitrogen 
     427                     zrat_chl1_p1n = tracer_bkg(ji,jj,jk,jp_fabm_chl1) / tracer_bkg(ji,jj,jk,jp_fabm_p1n) 
     428                     zrat_p1c_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_p1c)  / tracer_bkg(ji,jj,jk,jp_fabm_p1n) 
     429                     zrat_p1p_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_p1p)  / tracer_bkg(ji,jj,jk,jp_fabm_p1n) 
     430                     zrat_p1s_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_p1s)  / tracer_bkg(ji,jj,jk,jp_fabm_p1n) 
     431                     zrat_chl2_p2n = tracer_bkg(ji,jj,jk,jp_fabm_chl2) / tracer_bkg(ji,jj,jk,jp_fabm_p2n) 
     432                     zrat_p2c_p2n  = tracer_bkg(ji,jj,jk,jp_fabm_p2c)  / tracer_bkg(ji,jj,jk,jp_fabm_p2n) 
     433                     zrat_p2p_p2n  = tracer_bkg(ji,jj,jk,jp_fabm_p2p)  / tracer_bkg(ji,jj,jk,jp_fabm_p2n) 
     434                     zrat_chl3_p3n = tracer_bkg(ji,jj,jk,jp_fabm_chl3) / tracer_bkg(ji,jj,jk,jp_fabm_p3n) 
     435                     zrat_p3c_p3n  = tracer_bkg(ji,jj,jk,jp_fabm_p3c)  / tracer_bkg(ji,jj,jk,jp_fabm_p3n) 
     436                     zrat_p3p_p3n  = tracer_bkg(ji,jj,jk,jp_fabm_p3p)  / tracer_bkg(ji,jj,jk,jp_fabm_p3n) 
     437                     zrat_chl4_p4n = tracer_bkg(ji,jj,jk,jp_fabm_chl4) / tracer_bkg(ji,jj,jk,jp_fabm_p4n) 
     438                     zrat_p4c_p4n  = tracer_bkg(ji,jj,jk,jp_fabm_p4c)  / tracer_bkg(ji,jj,jk,jp_fabm_p4n) 
     439                     zrat_p4p_p4n  = tracer_bkg(ji,jj,jk,jp_fabm_p4p)  / tracer_bkg(ji,jj,jk,jp_fabm_p4n) 
    319440                      
    320                      phyto2d_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn 
    321                      phyto2d_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd 
    322                      phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_pds_phd 
    323                      phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,jk,jpphn) * zrat_chn_phn 
    324                      phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_chd_phd 
     441                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p1n 
     442                     phyto2d_balinc(ji,jj,jk,jp_fabm_p2n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p2n 
     443                     phyto2d_balinc(ji,jj,jk,jp_fabm_p3n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p3n 
     444                     phyto2d_balinc(ji,jj,jk,jp_fabm_p4n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p4n 
     445                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl1) = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_chl1_p1n 
     446                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_p1c_p1n 
     447                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_p1p_p1n 
     448                     phyto2d_balinc(ji,jj,jk,jp_fabm_p1s)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p1n) * zrat_p1s_p1n 
     449                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl2) = phyto2d_balinc(ji,jj,jk,jp_fabm_p2n) * zrat_chl2_p2n 
     450                     phyto2d_balinc(ji,jj,jk,jp_fabm_p2c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p2n) * zrat_p2c_p2n 
     451                     phyto2d_balinc(ji,jj,jk,jp_fabm_p2p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p2n) * zrat_p2p_p2n 
     452                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl3) = phyto2d_balinc(ji,jj,jk,jp_fabm_p3n) * zrat_chl3_p3n 
     453                     phyto2d_balinc(ji,jj,jk,jp_fabm_p3c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p3n) * zrat_p3c_p3n 
     454                     phyto2d_balinc(ji,jj,jk,jp_fabm_p3p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p3n) * zrat_p3p_p3n 
     455                     phyto2d_balinc(ji,jj,jk,jp_fabm_chl4) = phyto2d_balinc(ji,jj,jk,jp_fabm_p4n) * zrat_chl4_p4n 
     456                     phyto2d_balinc(ji,jj,jk,jp_fabm_p4c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p4n) * zrat_p4c_p4n 
     457                     phyto2d_balinc(ji,jj,jk,jp_fabm_p4p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_p4n) * zrat_p4p_p4n 
    325458                  ENDIF 
    326459 
    327460                  ! Zooplankton nitrogen split up based on existing ratios 
    328                   IF ( ( tracer_bkg(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpzme) > 0.0 ) ) THEN 
    329                      zfrac_zmi = tracer_bkg(ji,jj,jk,jpzmi) / & 
    330                         &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme)) 
    331                      zfrac_zme = tracer_bkg(ji,jj,jk,jpzme) / & 
    332                         &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme)) 
    333                      phyto2d_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi 
    334                      phyto2d_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme 
    335                   ENDIF 
    336  
    337                   ! Nitrogen nutrient straight from balancing scheme 
    338                   phyto2d_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1)) 
    339  
    340                   ! Nitrogen detritus straight from balancing scheme 
    341                   phyto2d_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4)) 
     461                  ! Update carbon and phosphorus according to existing ratios 
     462                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) > 0.0 ) .AND. & 
     463                     & ( tracer_bkg(ji,jj,jk,jp_fabm_z5n) > 0.0 ) .AND. & 
     464                     & ( tracer_bkg(ji,jj,jk,jp_fabm_z6n) > 0.0 ) ) THEN 
     465                     zfrac_z4n = ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) / & 
     466                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) + & 
     467                        &          tracer_bkg(ji,jj,jk,jp_fabm_z5n) + & 
     468                        &          tracer_bkg(ji,jj,jk,jp_fabm_z6n) ) 
     469                     zfrac_z5n = tracer_bkg(ji,jj,jk,jp_fabm_z5n) / & 
     470                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) + & 
     471                        &          tracer_bkg(ji,jj,jk,jp_fabm_z5n) + & 
     472                        &          tracer_bkg(ji,jj,jk,jp_fabm_z6n) ) 
     473                     zfrac_z6n = tracer_bkg(ji,jj,jk,jp_fabm_z6n) / & 
     474                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_z4c) * z4qnc ) + & 
     475                        &          tracer_bkg(ji,jj,jk,jp_fabm_z5n) + & 
     476                        &          tracer_bkg(ji,jj,jk,jp_fabm_z6n) ) 
     477                     zrat_z4c_z4n = 1.0 / z4qnc 
     478                     zrat_z5c_z5n = tracer_bkg(ji,jj,jk,jp_fabm_z5c) / tracer_bkg(ji,jj,jk,jp_fabm_z5n) 
     479                     zrat_z5p_z5n = tracer_bkg(ji,jj,jk,jp_fabm_z5p) / tracer_bkg(ji,jj,jk,jp_fabm_z5n) 
     480                     zrat_z6c_z6n = tracer_bkg(ji,jj,jk,jp_fabm_z6c) / tracer_bkg(ji,jj,jk,jp_fabm_z6n) 
     481                     zrat_z6p_z6n = tracer_bkg(ji,jj,jk,jp_fabm_z6p) / tracer_bkg(ji,jj,jk,jp_fabm_z6n) 
     482                     phyto2d_balinc(ji,jj,jk,jp_fabm_z5n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z5n 
     483                     phyto2d_balinc(ji,jj,jk,jp_fabm_z6n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z6n 
     484                     phyto2d_balinc(ji,jj,jk,jp_fabm_z4c) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z4n * zrat_z4c_z4n 
     485                     phyto2d_balinc(ji,jj,jk,jp_fabm_z5c) = phyto2d_balinc(ji,jj,jk,jp_fabm_z5n) * zrat_z5c_z5n 
     486                     phyto2d_balinc(ji,jj,jk,jp_fabm_z6c) = phyto2d_balinc(ji,jj,jk,jp_fabm_z6n) * zrat_z6c_z6n 
     487                     phyto2d_balinc(ji,jj,jk,jp_fabm_z5p) = phyto2d_balinc(ji,jj,jk,jp_fabm_z5n) * zrat_z5p_z5n 
     488                     phyto2d_balinc(ji,jj,jk,jp_fabm_z6p) = phyto2d_balinc(ji,jj,jk,jp_fabm_z6n) * zrat_z6p_z6n 
     489                  ENDIF 
     490 
     491                  ! Nitrogen nutrient split between nitrate and ammonium based on existing ratios 
     492                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_n3n) > 0.0 ) .AND. & 
     493                     & ( tracer_bkg(ji,jj,jk,jp_fabm_n4n) > 0.0 ) ) THEN 
     494                     zfrac_n3n = tracer_bkg(ji,jj,jk,jp_fabm_n3n) / & 
     495                        &        (tracer_bkg(ji,jj,jk,jp_fabm_n3n) + tracer_bkg(ji,jj,jk,jp_fabm_n4n)) 
     496                     zfrac_n4n = tracer_bkg(ji,jj,jk,jp_fabm_n4n) / & 
     497                        &        (tracer_bkg(ji,jj,jk,jp_fabm_n3n) + tracer_bkg(ji,jj,jk,jp_fabm_n4n)) 
     498                     phyto2d_balinc(ji,jj,jk,jp_fabm_n3n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n3n 
     499                     phyto2d_balinc(ji,jj,jk,jp_fabm_n4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n4n 
     500                  ENDIF 
     501 
     502                  ! Detritus nitrogen split up based on existing ratios 
     503                  ! Update carbon and phosphorus according to existing ratios 
     504                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_r4n) > 0.0 ) .AND. & 
     505                     & ( tracer_bkg(ji,jj,jk,jp_fabm_r6n) > 0.0 ) .AND. & 
     506                     & ( tracer_bkg(ji,jj,jk,jp_fabm_r8n) > 0.0 ) ) THEN 
     507                     zfrac_r4n = tracer_bkg(ji,jj,jk,jp_fabm_r4n) / & 
     508                        &        (tracer_bkg(ji,jj,jk,jp_fabm_r4n) + & 
     509                        &         tracer_bkg(ji,jj,jk,jp_fabm_r6n) + & 
     510                        &         tracer_bkg(ji,jj,jk,jp_fabm_r8n)) 
     511                     zfrac_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6n) / & 
     512                        &        (tracer_bkg(ji,jj,jk,jp_fabm_r4n) + & 
     513                        &         tracer_bkg(ji,jj,jk,jp_fabm_r6n) + & 
     514                        &         tracer_bkg(ji,jj,jk,jp_fabm_r8n)) 
     515                     zfrac_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8n) / & 
     516                        &        (tracer_bkg(ji,jj,jk,jp_fabm_r4n) + & 
     517                        &         tracer_bkg(ji,jj,jk,jp_fabm_r6n) + & 
     518                        &         tracer_bkg(ji,jj,jk,jp_fabm_r8n)) 
     519                     zrat_r4c_r4n = tracer_bkg(ji,jj,jk,jp_fabm_r4c) / tracer_bkg(ji,jj,jk,jp_fabm_r4n) 
     520                     zrat_r4p_r4n = tracer_bkg(ji,jj,jk,jp_fabm_r4p) / tracer_bkg(ji,jj,jk,jp_fabm_r4n) 
     521                     zrat_r6c_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6c) / tracer_bkg(ji,jj,jk,jp_fabm_r6n) 
     522                     zrat_r6p_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6p) / tracer_bkg(ji,jj,jk,jp_fabm_r6n) 
     523                     zrat_r6s_r6n = tracer_bkg(ji,jj,jk,jp_fabm_r6s) / tracer_bkg(ji,jj,jk,jp_fabm_r6n) 
     524                     zrat_r8c_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8c) / tracer_bkg(ji,jj,jk,jp_fabm_r8n) 
     525                     zrat_r8p_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8p) / tracer_bkg(ji,jj,jk,jp_fabm_r8n) 
     526                     zrat_r8s_r8n = tracer_bkg(ji,jj,jk,jp_fabm_r8s) / tracer_bkg(ji,jj,jk,jp_fabm_r8n) 
     527                     phyto2d_balinc(ji,jj,jk,jp_fabm_r4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r4n 
     528                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r6n 
     529                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r8n 
     530                     phyto2d_balinc(ji,jj,jk,jp_fabm_r4c) = phyto2d_balinc(ji,jj,jk,jp_fabm_r4n) * zrat_r4c_r4n 
     531                     phyto2d_balinc(ji,jj,jk,jp_fabm_r4p) = phyto2d_balinc(ji,jj,jk,jp_fabm_r4n) * zrat_r4p_r4n 
     532                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6c) = phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) * zrat_r6c_r6n 
     533                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6p) = phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) * zrat_r6p_r6n 
     534                     phyto2d_balinc(ji,jj,jk,jp_fabm_r6s) = phyto2d_balinc(ji,jj,jk,jp_fabm_r6n) * zrat_r6s_r6n 
     535                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8c) = phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) * zrat_r8c_r8n 
     536                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8p) = phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) * zrat_r8p_r8n 
     537                     phyto2d_balinc(ji,jj,jk,jp_fabm_r8s) = phyto2d_balinc(ji,jj,jk,jp_fabm_r8n) * zrat_r8s_r8n 
     538                  ENDIF 
    342539 
    343540                  ! DIC straight from balancing scheme 
    344                   phyto2d_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5)) 
     541                  phyto2d_balinc(ji,jj,jk,jp_fabm_o3c) = outincs(ji,jj,jk,i_tracer(5)) 
    345542 
    346543                  ! Alkalinity straight from balancing scheme 
    347                   phyto2d_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6)) 
    348  
    349                   ! Remove diatom silicate increment from nutrient silicate to conserve mass 
    350                   IF ( ( tracer_bkg(ji,jj,jk,jpsil) - phyto2d_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN 
    351                      phyto2d_balinc(ji,jj,jk,jpsil) = phyto2d_balinc(ji,jj,jk,jppds) * (-1.0) 
    352                   ENDIF 
    353  
    354                   ! Carbon detritus based on existing ratios 
    355                   IF ( ( tracer_bkg(ji,jj,jk,jpdet) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpdtc) > 0.0 ) ) THEN 
    356                      zrat_dtc_det = tracer_bkg(ji,jj,jk,jpdtc) / tracer_bkg(ji,jj,jk,jpdet) 
    357                      phyto2d_balinc(ji,jj,jk,jpdtc) = phyto2d_balinc(ji,jj,jk,jpdet) * zrat_dtc_det 
    358                   ENDIF 
    359  
    360                   ! Do nothing with iron or oxygen for the time being 
    361                   phyto2d_balinc(ji,jj,jk,jpfer) = 0.0 
    362                   phyto2d_balinc(ji,jj,jk,jpoxy) = 0.0 
     544                  phyto2d_balinc(ji,jj,jk,jp_fabm_o3ba) = outincs(ji,jj,jk,i_tracer(6)) 
     545 
     546                  ! Remove P/R silicon increments from silicate to conserve mass 
     547                  zfrac = phyto2d_balinc(ji,jj,jk,jp_fabm_p1s) + & 
     548                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r6s) + & 
     549                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r8s) 
     550                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_n5s) - zfrac ) > 0.0 ) THEN 
     551                     phyto2d_balinc(ji,jj,jk,jp_fabm_n5s) = zfrac * (-1.0) 
     552                  ENDIF 
     553 
     554                  ! Remove P/Z/R phosphorus increments from phosphate to conserve mass 
     555                  zfrac = phyto2d_balinc(ji,jj,jk,jp_fabm_p1p) + & 
     556                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_p2p) + & 
     557                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_p3p) + & 
     558                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_p4p) + & 
     559                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_z5p) + & 
     560                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_z6p) + & 
     561                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r4p) + & 
     562                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r6p) + & 
     563                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_r8p) 
     564                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_n1p) - zfrac ) > 0.0 ) THEN 
     565                     phyto2d_balinc(ji,jj,jk,jp_fabm_n1p) = zfrac * (-1.0) 
     566                  ENDIF 
    363567                   
    364568               END DO 
     
    366570         END DO 
    367571       
    368       ELSE   ! No nitrogen balancing 
    369        
    370          ! Initialise individual chlorophyll increments to zero 
    371          phyto2d_balinc(:,:,:,jpchn) = 0.0 
    372          phyto2d_balinc(:,:,:,jpchd) = 0.0 
     572      ELSE   ! No nitrogen balancing - just update phytoplankton 
    373573          
    374574         ! Split up total surface chlorophyll increments 
    375575         DO jj = 1, jpj 
    376576            DO ji = 1, jpi 
    377                IF ( ( tracer_bkg(ji,jj,1,jpchn) > 0.0 ) .AND. & 
    378                   & ( tracer_bkg(ji,jj,1,jpchd) > 0.0 ) ) THEN 
     577               IF ( ( tracer_bkg(ji,jj,1,jp_fabm_chl1) > 0.0 ) .AND. & 
     578                  & ( tracer_bkg(ji,jj,1,jp_fabm_chl2) > 0.0 ) .AND. & 
     579                  & ( tracer_bkg(ji,jj,1,jp_fabm_chl3) > 0.0 ) .AND. & 
     580                  & ( tracer_bkg(ji,jj,1,jp_fabm_chl4) > 0.0 ) ) THEN 
    379581                  IF ( ld_chltot ) THEN 
    380582                     ! Chlorophyll split up based on existing ratios 
    381                      zfrac_chn = tracer_bkg(ji,jj,1,jpchn) / & 
    382                         &        ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) ) 
    383                      zfrac_chd = tracer_bkg(ji,jj,1,jpchd) / & 
    384                         &        ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) ) 
    385                      phyto2d_balinc(ji,jj,1,jpchn) = pinc_chltot(ji,jj) * zfrac_chn 
    386                      phyto2d_balinc(ji,jj,1,jpchd) = pinc_chltot(ji,jj) * zfrac_chd 
     583                     zfrac_chl1 = tracer_bkg(ji,jj,1,jp_fabm_chl1) /  & 
     584                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + & 
     585                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + & 
     586                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + & 
     587                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) ) 
     588                     zfrac_chl2 = tracer_bkg(ji,jj,1,jp_fabm_chl2) /  & 
     589                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + & 
     590                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + & 
     591                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + & 
     592                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) ) 
     593                     zfrac_chl3 = tracer_bkg(ji,jj,1,jp_fabm_chl3) /  & 
     594                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + & 
     595                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + & 
     596                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + & 
     597                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) ) 
     598                     zfrac_chl4 = tracer_bkg(ji,jj,1,jp_fabm_chl4) /  & 
     599                        &        ( tracer_bkg(ji,jj,1,jp_fabm_chl1) + & 
     600                        &          tracer_bkg(ji,jj,1,jp_fabm_chl2) + & 
     601                        &          tracer_bkg(ji,jj,1,jp_fabm_chl3) + & 
     602                        &          tracer_bkg(ji,jj,1,jp_fabm_chl4) ) 
     603                     phyto2d_balinc(ji,jj,1,jp_fabm_chl1) = pinc_chltot(ji,jj) * zfrac_chl1 
     604                     phyto2d_balinc(ji,jj,1,jp_fabm_chl2) = pinc_chltot(ji,jj) * zfrac_chl2 
     605                     phyto2d_balinc(ji,jj,1,jp_fabm_chl3) = pinc_chltot(ji,jj) * zfrac_chl3 
     606                     phyto2d_balinc(ji,jj,1,jp_fabm_chl4) = pinc_chltot(ji,jj) * zfrac_chl4 
    387607                  ENDIF 
    388608                  IF( ld_chldia ) THEN 
    389                      phyto2d_balinc(ji,jj,1,jpchd) = pinc_chldia(ji,jj) 
    390                   ENDIF 
    391                   IF( ld_chlnon ) THEN 
    392                      phyto2d_balinc(ji,jj,1,jpchn) = pinc_chlnon(ji,jj) 
     609                     phyto2d_balinc(ji,jj,1,jp_fabm_chl1) = pinc_chldia(ji,jj) 
     610                  ENDIF 
     611                  IF( ld_chlnan ) THEN 
     612                     phyto2d_balinc(ji,jj,1,jp_fabm_chl2) = pinc_chlnan(ji,jj) 
     613                  ENDIF 
     614                  IF( ld_chlpic ) THEN 
     615                     phyto2d_balinc(ji,jj,1,jp_fabm_chl3) = pinc_chlpic(ji,jj) 
     616                  ENDIF 
     617                  IF( ld_chldin ) THEN 
     618                     phyto2d_balinc(ji,jj,1,jp_fabm_chl4) = pinc_chldin(ji,jj) 
    393619                  ENDIF 
    394620                   
    395                   ! Maintain stoichiometric ratios of nitrogen and silicate 
    396                   IF ( ld_chltot .OR. ld_chlnon ) THEN 
    397                      zrat_phn_chn = tracer_bkg(ji,jj,1,jpphn) / tracer_bkg(ji,jj,1,jpchn) 
    398                      phyto2d_balinc(ji,jj,1,jpphn) = phyto2d_balinc(ji,jj,1,jpchn) * zrat_phn_chn 
    399                   ENDIF 
     621                  ! Maintain stoichiometric ratios of carbon, nitrogen, phosphorus and silicon 
    400622                  IF ( ld_chltot .OR. ld_chldia ) THEN 
    401                      zrat_phd_chd = tracer_bkg(ji,jj,1,jpphd) / tracer_bkg(ji,jj,1,jpchd) 
    402                      phyto2d_balinc(ji,jj,1,jpphd) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_phd_chd 
    403                      zrat_pds_chd = tracer_bkg(ji,jj,1,jppds) / tracer_bkg(ji,jj,1,jpchd) 
    404                      phyto2d_balinc(ji,jj,1,jppds) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_pds_chd 
     623                     zrat_p1c_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1c) / tracer_bkg(ji,jj,1,jp_fabm_chl1) 
     624                     zrat_p1n_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1n) / tracer_bkg(ji,jj,1,jp_fabm_chl1) 
     625                     zrat_p1p_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1p) / tracer_bkg(ji,jj,1,jp_fabm_chl1) 
     626                     zrat_p1s_chl1 = tracer_bkg(ji,jj,1,jp_fabm_p1s) / tracer_bkg(ji,jj,1,jp_fabm_chl1) 
     627                     phyto2d_balinc(ji,jj,1,jp_fabm_p1c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1c_chl1 
     628                     phyto2d_balinc(ji,jj,1,jp_fabm_p1n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1n_chl1 
     629                     phyto2d_balinc(ji,jj,1,jp_fabm_p1p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1p_chl1 
     630                     phyto2d_balinc(ji,jj,1,jp_fabm_p1s) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) * zrat_p1s_chl1 
     631                  ENDIF 
     632                  IF ( ld_chltot .OR. ld_chlnan ) THEN 
     633                     zrat_p2c_chl2 = tracer_bkg(ji,jj,1,jp_fabm_p2c) / tracer_bkg(ji,jj,1,jp_fabm_chl2) 
     634                     zrat_p2n_chl2 = tracer_bkg(ji,jj,1,jp_fabm_p2n) / tracer_bkg(ji,jj,1,jp_fabm_chl2) 
     635                     zrat_p2p_chl2 = tracer_bkg(ji,jj,1,jp_fabm_p2p) / tracer_bkg(ji,jj,1,jp_fabm_chl2) 
     636                     phyto2d_balinc(ji,jj,1,jp_fabm_p2c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2) * zrat_p2c_chl2 
     637                     phyto2d_balinc(ji,jj,1,jp_fabm_p2n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2) * zrat_p2n_chl2 
     638                     phyto2d_balinc(ji,jj,1,jp_fabm_p2p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2) * zrat_p2p_chl2 
     639                  ENDIF 
     640                  IF ( ld_chltot .OR. ld_chlpic ) THEN 
     641                     zrat_p3c_chl3 = tracer_bkg(ji,jj,1,jp_fabm_p3c) / tracer_bkg(ji,jj,1,jp_fabm_chl3) 
     642                     zrat_p3n_chl3 = tracer_bkg(ji,jj,1,jp_fabm_p3n) / tracer_bkg(ji,jj,1,jp_fabm_chl3) 
     643                     zrat_p3p_chl3 = tracer_bkg(ji,jj,1,jp_fabm_p3p) / tracer_bkg(ji,jj,1,jp_fabm_chl3) 
     644                     phyto2d_balinc(ji,jj,1,jp_fabm_p3c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3) * zrat_p3c_chl3 
     645                     phyto2d_balinc(ji,jj,1,jp_fabm_p3n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3) * zrat_p3n_chl3 
     646                     phyto2d_balinc(ji,jj,1,jp_fabm_p3p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3) * zrat_p3p_chl3 
     647                  ENDIF 
     648                  IF ( ld_chltot .OR. ld_chldin ) THEN 
     649                     zrat_p4c_chl4 = tracer_bkg(ji,jj,1,jp_fabm_p4c) / tracer_bkg(ji,jj,1,jp_fabm_chl4) 
     650                     zrat_p4n_chl4 = tracer_bkg(ji,jj,1,jp_fabm_p4n) / tracer_bkg(ji,jj,1,jp_fabm_chl4) 
     651                     zrat_p4p_chl4 = tracer_bkg(ji,jj,1,jp_fabm_p4p) / tracer_bkg(ji,jj,1,jp_fabm_chl4) 
     652                     phyto2d_balinc(ji,jj,1,jp_fabm_p4c) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4) * zrat_p4c_chl4 
     653                     phyto2d_balinc(ji,jj,1,jp_fabm_p4n) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4) * zrat_p4n_chl4 
     654                     phyto2d_balinc(ji,jj,1,jp_fabm_p4p) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4) * zrat_p4p_chl4 
    405655                  ENDIF 
    406656               ENDIF 
     
    422672               ! 
    423673               DO jk = 2, jkmax 
    424                   phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,1,jpchn) 
    425                   phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,1,jpchd) 
    426                   phyto2d_balinc(ji,jj,jk,jpphn) = phyto2d_balinc(ji,jj,1,jpphn) 
    427                   phyto2d_balinc(ji,jj,jk,jpphd) = phyto2d_balinc(ji,jj,1,jpphd) 
    428                   phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,1,jppds) 
     674                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl1) = phyto2d_balinc(ji,jj,1,jp_fabm_chl1) 
     675                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1c) 
     676                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1n) 
     677                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1p) 
     678                  phyto2d_balinc(ji,jj,jk,jp_fabm_p1s)  = phyto2d_balinc(ji,jj,1,jp_fabm_p1s) 
     679                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl2) = phyto2d_balinc(ji,jj,1,jp_fabm_chl2) 
     680                  phyto2d_balinc(ji,jj,jk,jp_fabm_p2c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p2c) 
     681                  phyto2d_balinc(ji,jj,jk,jp_fabm_p2n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p2n) 
     682                  phyto2d_balinc(ji,jj,jk,jp_fabm_p2p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p2p) 
     683                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl3) = phyto2d_balinc(ji,jj,1,jp_fabm_chl3) 
     684                  phyto2d_balinc(ji,jj,jk,jp_fabm_p3c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p3c) 
     685                  phyto2d_balinc(ji,jj,jk,jp_fabm_p3n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p3n) 
     686                  phyto2d_balinc(ji,jj,jk,jp_fabm_p3p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p3p) 
     687                  phyto2d_balinc(ji,jj,jk,jp_fabm_chl4) = phyto2d_balinc(ji,jj,1,jp_fabm_chl4) 
     688                  phyto2d_balinc(ji,jj,jk,jp_fabm_p4c)  = phyto2d_balinc(ji,jj,1,jp_fabm_p4c) 
     689                  phyto2d_balinc(ji,jj,jk,jp_fabm_p4n)  = phyto2d_balinc(ji,jj,1,jp_fabm_p4n) 
     690                  phyto2d_balinc(ji,jj,jk,jp_fabm_p4p)  = phyto2d_balinc(ji,jj,1,jp_fabm_p4p) 
    429691               END DO 
    430692               ! 
     
    432694         END DO 
    433695 
    434          ! Set other balancing increments to zero 
    435          phyto2d_balinc(:,:,:,jpzmi) = 0.0 
    436          phyto2d_balinc(:,:,:,jpzme) = 0.0 
    437          phyto2d_balinc(:,:,:,jpdin) = 0.0 
    438          phyto2d_balinc(:,:,:,jpsil) = 0.0 
    439          phyto2d_balinc(:,:,:,jpfer) = 0.0 
    440          phyto2d_balinc(:,:,:,jpdet) = 0.0 
    441          phyto2d_balinc(:,:,:,jpdtc) = 0.0 
    442          phyto2d_balinc(:,:,:,jpdic) = 0.0 
    443          phyto2d_balinc(:,:,:,jpalk) = 0.0 
    444          phyto2d_balinc(:,:,:,jpoxy) = 0.0 
    445  
    446696      ENDIF 
    447        
    448       ! If performing extra tidal mixing in the Indonesian Throughflow, 
    449       ! increments have been found to make the carbon cycle unstable 
    450       ! Therefore, mask these out 
    451       IF ( ln_tmx_itf ) THEN 
    452          DO jn = 1, jptra 
    453             DO jk = 1, jpk 
    454                phyto2d_balinc(:,:,jk,jn) = phyto2d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) ) 
    455             END DO 
    456          END DO 
    457       ENDIF 
    458  
    459    END SUBROUTINE asm_phyto2d_bal_medusa 
     697 
     698   END SUBROUTINE asm_phyto2d_bal_ersem 
    460699 
    461700#else 
     
    464703   !!---------------------------------------------------------------------- 
    465704CONTAINS 
    466    SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      & 
    467       &                              pinc_chltot,                    & 
    468       &                              ld_chldia,                      & 
    469       &                              pinc_chldia,                    & 
    470       &                              ld_chlnon,                      & 
    471       &                              pinc_chlnon,                    & 
    472       &                              ld_phytot,                      & 
    473       &                              pinc_phytot,                    & 
    474       &                              ld_phydia,                      & 
    475       &                              pinc_phydia,                    & 
    476       &                              ld_phynon,                      & 
    477       &                              pinc_phynon,                    & 
    478       &                              pincper,                        & 
    479       &                              p_maxchlinc, ld_phytobal, pmld, & 
    480       &                              pgrow_avg_bkg, ploss_avg_bkg,   & 
    481       &                              phyt_avg_bkg, mld_max_bkg,      & 
    482       &                              tracer_bkg, phyto2d_balinc ) 
     705   SUBROUTINE asm_phyto2d_bal_ersem(  ld_chltot,                      & 
     706      &                               pinc_chltot,                    & 
     707      &                               ld_chldia,                      & 
     708      &                               pinc_chldia,                    & 
     709      &                               ld_chlnan,                      & 
     710      &                               pinc_chlnan,                    & 
     711      &                               ld_chlpic,                      & 
     712      &                               pinc_chlpic,                    & 
     713      &                               ld_chldin,                      & 
     714      &                               pinc_chldin,                    & 
     715      &                               pincper,                        & 
     716      &                               p_maxchlinc, ld_phytobal, pmld, & 
     717      &                               pgrow_avg_bkg, ploss_avg_bkg,   & 
     718      &                               phyt_avg_bkg, mld_max_bkg,      & 
     719      &                               totalk_bkg,                     & 
     720      &                               tracer_bkg, phyto2d_balinc ) 
    483721      LOGICAL :: ld_chltot 
    484722      REAL    :: pinc_chltot(:,:) 
    485723      LOGICAL :: ld_chldia 
    486724      REAL    :: pinc_chldia(:,:) 
    487       LOGICAL :: ld_chlnon 
    488       REAL    :: pinc_chlnon(:,:) 
    489       LOGICAL :: ld_phytot 
    490       REAL    :: pinc_phytot(:,:) 
    491       LOGICAL :: ld_phydia 
    492       REAL    :: pinc_phydia(:,:) 
    493       LOGICAL :: ld_phynon 
    494       REAL    :: pinc_phynon(:,:) 
     725      LOGICAL :: ld_chlnan 
     726      REAL    :: pinc_chlnan(:,:) 
     727      LOGICAL :: ld_chlpic 
     728      REAL    :: pinc_chlpic(:,:) 
     729      LOGICAL :: ld_chldin 
     730      REAL    :: pinc_chldin(:,:) 
    495731      REAL    :: pincper 
    496732      REAL    :: p_maxchlinc 
     
    501737      REAL    :: phyt_avg_bkg(:,:) 
    502738      REAL    :: mld_max_bkg(:,:) 
     739      REAL    :: totalk_bkg(:,:,:) 
    503740      REAL    :: tracer_bkg(:,:,:,:) 
    504741      REAL    :: phyto2d_balinc(:,:,:,:) 
    505       WRITE(*,*) 'asm_phyto2d_bal_medusa: You should not have seen this print! error?' 
    506    END SUBROUTINE asm_phyto2d_bal_medusa 
     742      WRITE(*,*) 'asm_phyto2d_bal_ersem: You should not have seen this print! error?' 
     743   END SUBROUTINE asm_phyto2d_bal_ersem 
    507744#endif 
    508745 
    509746   !!====================================================================== 
    510 END MODULE asmphyto2dbal_medusa 
     747END MODULE asmphyto2dbal_ersem 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/par_fabm.F90

    r10390 r10622  
    1212 
    1313   ! Variables needed for OBS/ASM 
    14    INTEGER, PUBLIC :: jp_fabm_chl1, jp_fabm_chl2, & 
    15                       jp_fabm_chl3, jp_fabm_chl4, & 
    16                       jp_fabm_p1c,  jp_fabm_p1n,  & 
    17                       jp_fabm_p1p,  jp_fabm_p1s,  & 
    18                       jp_fabm_p2c,  jp_fabm_p2n,  & 
    19                       jp_fabm_p2p,  jp_fabm_p3c,  & 
    20                       jp_fabm_p3n,  jp_fabm_p3p,  & 
    21                       jp_fabm_p4c,  jp_fabm_p4n,  & 
    22                       jp_fabm_p4p,  jp_fabm_z4c,  & 
    23                       jp_fabm_z5c,  jp_fabm_z5n,  & 
    24                       jp_fabm_z5p,  jp_fabm_z6c,  & 
    25                       jp_fabm_z6n,  jp_fabm_z6p,  & 
    26                       jp_fabm_n1p,  jp_fabm_n3n,  & 
    27                       jp_fabm_n4n,  jp_fabm_n5s,  & 
    28                       jp_fabm_o2o,  jp_fabm_o3c,  & 
    29                       jp_fabm_o3a,  jp_fabm_o3ph, & 
    30                       jp_fabm_o3pc 
     14   INTEGER, PUBLIC :: jp_fabm_chl1,  jp_fabm_chl2, & 
     15                      jp_fabm_chl3,  jp_fabm_chl4, & 
     16                      jp_fabm_p1c,   jp_fabm_p1n,  & 
     17                      jp_fabm_p1p,   jp_fabm_p1s,  & 
     18                      jp_fabm_p2c,   jp_fabm_p2n,  & 
     19                      jp_fabm_p2p,   jp_fabm_p3c,  & 
     20                      jp_fabm_p3n,   jp_fabm_p3p,  & 
     21                      jp_fabm_p4c,   jp_fabm_p4n,  & 
     22                      jp_fabm_p4p,   jp_fabm_z4c,  & 
     23                      jp_fabm_z5c,   jp_fabm_z5n,  & 
     24                      jp_fabm_z5p,   jp_fabm_z6c,  & 
     25                      jp_fabm_z6n,   jp_fabm_z6p,  & 
     26                      jp_fabm_n1p,   jp_fabm_n3n,  & 
     27                      jp_fabm_n4n,   jp_fabm_n5s,  & 
     28                      jp_fabm_o2o,   jp_fabm_o3c,  & 
     29                      jp_fabm_o3ta,  jp_fabm_o3ba, & 
     30                      jp_fabm_o3pc,  jp_fabm_o3ph, & 
     31                      jp_fabm_r4n,   jp_fabm_r4c,  & 
     32                      jp_fabm_r4p,   jp_fabm_r6n,  & 
     33                      jp_fabm_r6c,   jp_fabm_r6p,  & 
     34                      jp_fabm_r6s,   jp_fabm_r8n,  & 
     35                      jp_fabm_r8c,   jp_fabm_r8p,  & 
     36                      jp_fabm_r8s,                 & 
     37                      jp_fabm_pgrow, jp_fabm_ploss 
    3138 
    3239#if defined key_fabm 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/trcini_fabm.F90

    r10390 r10622  
    105105      jp_fabm_o2o  = fabm_state_index( 'O2_o' ) 
    106106      jp_fabm_o3c  = fabm_state_index( 'O3_c' ) 
    107       jp_fabm_o3a  = fabm_state_index( 'O3_bioalk' ) 
     107      jp_fabm_o3ba = fabm_state_index( 'O3_bioalk' ) 
     108      jp_fabm_r4n  = fabm_state_index( 'R4_n' ) 
     109      jp_fabm_r4c  = fabm_state_index( 'R4_c' ) 
     110      jp_fabm_r4p  = fabm_state_index( 'R4_p' ) 
     111      jp_fabm_r6n  = fabm_state_index( 'R6_n' ) 
     112      jp_fabm_r6c  = fabm_state_index( 'R6_c' ) 
     113      jp_fabm_r6p  = fabm_state_index( 'R6_p' ) 
     114      jp_fabm_r6s  = fabm_state_index( 'R6_s' ) 
     115      jp_fabm_r8n  = fabm_state_index( 'R8_n' ) 
     116      jp_fabm_r8c  = fabm_state_index( 'R8_c' ) 
     117      jp_fabm_r8p  = fabm_state_index( 'R8_p' ) 
     118      jp_fabm_r8s  = fabm_state_index( 'R8_s' ) 
    108119 
    109120      ! Get indexes for select diagnostic variables 
    110       jp_fabm_o3ph = fabm_diag_index( 'O3_pH' ) 
    111       jp_fabm_o3pc = fabm_diag_index( 'O3_pCO2' ) 
     121      jp_fabm_o3ta  = fabm_diag_index( 'O3_TA' ) 
     122      jp_fabm_o3ph  = fabm_diag_index( 'O3_pH' ) 
     123      jp_fabm_o3pc  = fabm_diag_index( 'O3_pCO2' ) 
     124      jp_fabm_pgrow = fabm_diag_index( 'p_grow_sum_result' ) 
     125      jp_fabm_ploss = fabm_diag_index( 'p_loss_sum_result' ) 
     126       
     127      MLD_MAX(:,:)   = 0.0 
     128      PGROW_AVG(:,:) = 0.0 
     129      PLOSS_AVG(:,:) = 0.0 
     130      PHYT_AVG(:,:)  = 0.0 
    112131 
    113132      IF (lwp) THEN 
     
    445464      END DO 
    446465      IF (fabm_state_index == -1) THEN 
    447          CALL ctl_stop( 'Could not find '//TRIM(state_name)//' state variable' ) 
     466         CALL ctl_warn( 'Could not find '//TRIM(state_name)//' state variable' ) 
    448467      ELSE 
    449468         IF (lwp) WRITE(numout,*) 'Index for '//TRIM(state_name)//' is: ', fabm_state_index 
     
    477496      END DO 
    478497      IF (fabm_diag_index == -1) THEN 
    479          CALL ctl_stop( 'Could not find '//TRIM(diag_name)//' diagnostic' ) 
     498         CALL ctl_warn( 'Could not find '//TRIM(diag_name)//' diagnostic' ) 
    480499      ELSE 
    481500         IF (lwp) WRITE(numout,*) 'Index for '//TRIM(diag_name)//' is: ', fabm_diag_index 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/FABM/trcsms_fabm.F90

    r10156 r10622  
    3333   USE inputs_fabm 
    3434   USE vertical_movement_fabm 
     35   USE zdfmxl 
     36   USE asmbgc, ONLY: mld_choice_bgc 
     37   USE lbclnk 
    3538 
    3639   !USE fldread         !  time interpolation 
     
    113116 
    114117      CALL st2d_fabm_nxt( kt ) 
     118       
     119      CALL asmdiags_fabm( kt ) 
    115120 
    116121      IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrfabm ) 
     
    130135 
    131136   END SUBROUTINE trc_sms_fabm 
     137    
     138   SUBROUTINE asmdiags_fabm( kt ) 
     139      INTEGER, INTENT(IN) :: kt 
     140      INTEGER :: ji,jj,jk,jkmax 
     141      REAL(wp), DIMENSION(jpi,jpj,jpk) :: pgrow_3d, ploss_3d, zmld 
     142       
     143      IF (kt == nittrc000) THEN 
     144         MLD_MAX(:,:) = 0.0 
     145      ENDIF 
     146      PGROW_AVG(:,:) = 0.0 
     147      PLOSS_AVG(:,:) = 0.0 
     148      PHYT_AVG(:,:)  = 0.0 
     149         
     150      pgrow_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_pgrow) 
     151      ploss_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_ploss) 
     152       
     153      SELECT CASE( mld_choice_bgc ) 
     154      CASE ( 1 )                   ! Turbocline/mixing depth [W points] 
     155         zmld(:,:) = hmld(:,:) 
     156      CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points] 
     157         zmld(:,:) = hmlp(:,:) 
     158      CASE ( 3 )                   ! Kara MLD [Interpolated] 
     159#if defined key_karaml 
     160         IF ( ln_kara ) THEN 
     161            zmld(:,:) = hmld_kara(:,:) 
     162         ELSE 
     163            CALL ctl_stop( ' Kara mixed layer requested for BGC assimilation,', & 
     164               &           ' but ln_kara=.false.' ) 
     165         ENDIF 
     166#else 
     167         CALL ctl_stop( ' Kara mixed layer requested for BGC assimilation,', & 
     168            &           ' but is not defined' ) 
     169#endif 
     170      CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points] 
     171         zmld(:,:) = hmld_tref(:,:) 
     172      CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points] 
     173         zmld(:,:) = hmlpt(:,:) 
     174      END SELECT 
     175    
     176      DO jj = 2, jpjm1 
     177         DO ji = 2, jpim1 
     178            ! 
     179            jkmax = jpk-1 
     180            DO jk = jpk-1, 1, -1 
     181               IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. & 
     182                  & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN 
     183                  zmld(ji,jj) = gdepw_n(ji,jj,jk+1) 
     184                  jkmax = jk 
     185               ENDIF 
     186            END DO 
     187            ! 
     188            DO jk = 1, jkmax 
     189               PHYT_AVG(ji,jj) = PHYT_AVG(ji,jj) + & 
     190                  &              trn(ji,jj,jk,jp_fabm_p1n) + & 
     191                  &              trn(ji,jj,jk,jp_fabm_p2n) + & 
     192                  &              trn(ji,jj,jk,jp_fabm_p3n) + & 
     193                  &              trn(ji,jj,jk,jp_fabm_p4n) 
     194               IF ( pgrow_3d(ji,jj,jk) .GT. 0.0 ) THEN 
     195                  PGROW_AVG(ji,jj) = PGROW_AVG(ji,jj) + & 
     196                     &               pgrow_3d(ji,jj,jk) 
     197               ENDIF 
     198               IF ( ploss_3d(ji,jj,jk) .GT. 0.0 ) THEN 
     199                  PLOSS_AVG(ji,jj) = PLOSS_AVG(ji,jj) + & 
     200                     &               ploss_3d(ji,jj,jk) 
     201               ENDIF 
     202            END DO 
     203            
     204            PHYT_AVG(ji,jj)  = PHYT_AVG(ji,jj)  / REAL(jkmax) 
     205            PGROW_AVG(ji,jj) = PGROW_AVG(ji,jj) / REAL(jkmax) 
     206            PLOSS_AVG(ji,jj) = PLOSS_AVG(ji,jj) / REAL(jkmax) 
     207    
     208            IF ( zmld(ji,jj) .GT. MLD_MAX(ji,jj) ) THEN 
     209               MLD_MAX(ji,jj) = zmld(ji,jj) 
     210            ENDIF 
     211            ! 
     212         END DO 
     213      END DO 
     214       
     215      PHYT_AVG(:,:)  = PHYT_AVG(:,:)  * tmask(:,:,1) 
     216      PGROW_AVG(:,:) = PGROW_AVG(:,:) * tmask(:,:,1) 
     217      PLOSS_AVG(:,:) = PLOSS_AVG(:,:) * tmask(:,:,1) 
     218      MLD_MAX(:,:)   = MLD_MAX(:,:)   * tmask(:,:,1) 
     219       
     220   END SUBROUTINE asmdiags_fabm 
    132221 
    133222   SUBROUTINE compute_fabm() 
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/TOP_SRC/trc.F90

    r10162 r10622  
    225225#endif 
    226226 
     227#if defined key_fabm 
     228   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  PGROW_AVG  !: Phytoplankton growth for use in ASM code 
     229   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  PLOSS_AVG  !: Phytoplankton loss   for use in ASM code 
     230   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  PHYT_AVG   !: Phytoplankton        for use in ASM code 
     231   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  MLD_MAX    !: Maximum MLD          for use in ASM code 
     232#endif 
     233 
    227234   !!---------------------------------------------------------------------- 
    228235   !! NEMO/TOP 3.3.1 , NEMO Consortium (2010) 
     
    253260! FABM <<<+++ 
    254261         &      ln_trc_sbc(jptra)     , ln_trc_cbc(jptra)     , ln_trc_obc(jptra)     ,       & 
     262         &      PGROW_AVG(jpi,jpj)    , PLOSS_AVG(jpi,jpj)    , PHYT_AVG(jpi,jpj)     ,       & 
     263         &      MLD_MAX(jpi,jpj)      ,                                                       & 
    255264#endif 
    256265#if defined key_bdy 
Note: See TracChangeset for help on using the changeset viewer.