Changeset 11990


Ignore:
Timestamp:
2019-11-27T16:43:05+01:00 (8 weeks ago)
Author:
dford
Message:

Get the nitrogen balancing working with 3D chlorophyll increments.

Location:
branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO
Files:
7 edited
1 copied

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90

    r9862 r11990  
    5959   USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA 
    6060      & asm_phyto2d_bal_medusa 
     61   USE asmphyto3dbal_medusa, ONLY: & ! phyto3d balancing for MEDUSA 
     62      & asm_phyto3d_bal_medusa 
    6163   USE asmpco2bal, ONLY:    & ! pCO2 balancing for MEDUSA 
    6264      & asm_pco2_bal 
     
    7173      & ploss_avg,          & 
    7274      & phyt_avg,           & 
     75      & pgrow_avg_3d,       & 
     76      & ploss_avg_3d,       & 
     77      & phyt_avg_3d,        & 
    7378      & mld_max 
    7479#elif defined key_hadocc 
     
    172177   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: ploss_avg_bkg   ! Background phyto loss 
    173178   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: phyt_avg_bkg    ! Background phyto conc 
     179   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pgrow_avg_3d_bkg   ! Background phyto growth 
     180   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: ploss_avg_3d_bkg   ! Background phyto loss 
     181   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: phyt_avg_3d_bkg    ! Background phyto conc 
    174182   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: mld_max_bkg     ! Background max MLD 
    175183   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg      ! Background tracer state 
     
    213221         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. & 
    214222         & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. & 
    215          & ( .NOT. ln_slphynoninc ) ) THEN 
     223         & ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. & 
     224         & ( .NOT. ln_pchltotinc  ) ) THEN 
    216225         CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', & 
    217             &           ' if not assimilating ocean colour,',                   & 
     226            &           ' if not assimilating phytoplankton,',                   & 
    218227            &           ' so ln_phytobal will be set to .false.') 
    219228         ln_phytobal = .FALSE. 
     
    525534         ALLOCATE( ploss_avg_bkg(jpi,jpj)        ) 
    526535         ALLOCATE( phyt_avg_bkg(jpi,jpj)         ) 
     536         ALLOCATE( pgrow_avg_3d_bkg(jpi,jpj,jpk) ) 
     537         ALLOCATE( ploss_avg_3d_bkg(jpi,jpj,jpk) ) 
     538         ALLOCATE( phyt_avg_3d_bkg(jpi,jpj,jpk)  ) 
    527539         ALLOCATE( mld_max_bkg(jpi,jpj)          ) 
    528540         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) ) 
     
    530542         ploss_avg_bkg(:,:)  = 0.0 
    531543         phyt_avg_bkg(:,:)   = 0.0 
     544         pgrow_avg_3d_bkg(:,:,:)  = 0.0 
     545         ploss_avg_3d_bkg(:,:,:)  = 0.0 
     546         phyt_avg_3d_bkg(:,:,:)   = 0.0 
    532547         mld_max_bkg(:,:)    = 0.0 
    533548         tracer_bkg(:,:,:,:) = 0.0 
     
    565580            CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg ) 
    566581            CALL iom_get( inum, jpdom_autoglo, 'phyt_avg',  phyt_avg_bkg  ) 
     582            CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg_3d', pgrow_avg_3d_bkg ) 
     583            CALL iom_get( inum, jpdom_autoglo, 'ploss_avg_3d', ploss_avg_3d_bkg ) 
     584            CALL iom_get( inum, jpdom_autoglo, 'phyt_avg_3d',  phyt_avg_3d_bkg  ) 
    567585            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   ) 
    568586            pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1) 
    569587            ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1) 
    570588            phyt_avg_bkg(:,:)  = phyt_avg_bkg(:,:)  * tmask(:,:,1) 
     589            pgrow_avg_3d_bkg(:,:,:) = pgrow_avg_3d_bkg(:,:,:) * tmask(:,:,:) 
     590            ploss_avg_3d_bkg(:,:,:) = ploss_avg_3d_bkg(:,:,:) * tmask(:,:,:) 
     591            phyt_avg_3d_bkg(:,:,:)  = phyt_avg_3d_bkg(:,:,:)  * tmask(:,:,:) 
    571592            mld_max_bkg(:,:)   = mld_max_bkg(:,:)   * tmask(:,:,1) 
    572593 
     
    727748            CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) ) 
    728749            CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) ) 
     750            IF ( ln_phytobal ) THEN 
     751               CALL iom_rstput( kt, kt, inum, 'phy3d_zmi', phyto3d_balinc(:,:,:,jpzmi) ) 
     752               CALL iom_rstput( kt, kt, inum, 'phy3d_zme', phyto3d_balinc(:,:,:,jpzme) ) 
     753               CALL iom_rstput( kt, kt, inum, 'phy3d_din', phyto3d_balinc(:,:,:,jpdin) ) 
     754               CALL iom_rstput( kt, kt, inum, 'phy3d_sil', phyto3d_balinc(:,:,:,jpsil) ) 
     755               CALL iom_rstput( kt, kt, inum, 'phy3d_fer', phyto3d_balinc(:,:,:,jpfer) ) 
     756               CALL iom_rstput( kt, kt, inum, 'phy3d_det', phyto3d_balinc(:,:,:,jpdet) ) 
     757               CALL iom_rstput( kt, kt, inum, 'phy3d_dtc', phyto3d_balinc(:,:,:,jpdtc) ) 
     758               CALL iom_rstput( kt, kt, inum, 'phy3d_dic', phyto3d_balinc(:,:,:,jpdic) ) 
     759               CALL iom_rstput( kt, kt, inum, 'phy3d_alk', phyto3d_balinc(:,:,:,jpalk) ) 
     760               CALL iom_rstput( kt, kt, inum, 'phy3d_oxy', phyto3d_balinc(:,:,:,jpoxy) ) 
     761            ENDIF 
    729762#elif defined key_hadocc 
    730763            CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) ) 
     
    810843      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg        ) 
    811844      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg         ) 
     845      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg_3d'   , pgrow_avg_3d        ) 
     846      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg_3d'   , ploss_avg_3d        ) 
     847      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg_3d'    , phyt_avg_3d         ) 
    812848      CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max          ) 
    813849      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn'  , trn(:,:,:,jpchn) ) 
     
    11631199      INTEGER                          :: ji, jj, jk   ! Loop counters 
    11641200      INTEGER                          :: it           ! Index 
     1201      REAL(wp)                         :: zincper      ! IAU interval in seconds 
    11651202      REAL(wp)                         :: zincwgt      ! IAU weight for timestep 
    11661203      REAL(wp)                         :: zfrac_chn    ! Fraction of jpchn 
     
    12151252            END DO 
    12161253         ENDIF 
    1217  
     1254          
    12181255#if defined key_medusa && defined key_foam_medusa 
    1219          ! Loop over each grid point partioning the increments based on existing ratios 
    1220          DO jk = 1, jpk 
    1221             DO jj = 1, jpj 
    1222                DO ji = 1, jpi 
    1223                   IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN 
    1224                      zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd)) 
    1225                      zfrac_chd = 1.0 - zfrac_chn 
    1226                      phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn 
    1227                      phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd 
    1228                      zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn) 
    1229                      zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd) 
    1230                      zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd) 
    1231                      phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn 
    1232                      phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd 
    1233                      phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd 
    1234                   ENDIF 
     1256         IF ( ln_phytobal ) THEN 
     1257            zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt 
     1258            CALL asm_phyto3d_bal_medusa( chl_inc,                             & 
     1259               &                         zincper,                             & 
     1260               &                         rn_maxchlinc,                        & 
     1261               &                         pgrow_avg_3d_bkg, ploss_avg_3d_bkg,  & 
     1262               &                         phyt_avg_3d_bkg,                     & 
     1263               &                         tracer_bkg, phyto3d_balinc ) 
     1264         ELSE 
     1265            ! Loop over each grid point partioning the increments based on existing ratios 
     1266            DO jk = 1, jpk 
     1267               DO jj = 1, jpj 
     1268                  DO ji = 1, jpi 
     1269                     IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN 
     1270                        zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd)) 
     1271                        zfrac_chd = 1.0 - zfrac_chn 
     1272                        phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn 
     1273                        phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd 
     1274                        zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn) 
     1275                        zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd) 
     1276                        zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd) 
     1277                        phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn 
     1278                        phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd 
     1279                        phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd 
     1280                     ENDIF 
     1281                  END DO 
    12351282               END DO 
    12361283            END DO 
    1237          END DO 
     1284         ENDIF 
    12381285#elif defined key_hadocc 
    12391286         phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:) 
     
    14201467         ! Account for phytoplankton balancing if required 
    14211468         IF ( ln_phytobal ) THEN 
    1422             dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic) 
    1423             alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk) 
     1469            IF ( ln_slchltotinc ) THEN 
     1470               dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic) 
     1471               alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk) 
     1472            ENDIF 
     1473            IF ( ln_plchltotinc ) THEN 
     1474               dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto3d_balinc(:,:,1,jpdic) 
     1475               alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto3d_balinc(:,:,1,jpalk) 
     1476            ENDIF 
    14241477         ELSE 
    14251478            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) 
     
    16661719         ! Account for phytoplankton balancing if required 
    16671720         IF ( ln_phytobal ) THEN 
    1668             dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic) 
    1669             alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk) 
    1670             din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin) 
    1671             sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil) 
     1721            IF ( ln_slchltotinc ) THEN 
     1722               dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic) 
     1723               alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk) 
     1724               din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin) 
     1725               sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil) 
     1726            ENDIF 
     1727            IF ( ln_plchltotinc ) THEN 
     1728               dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto3d_balinc(:,:,:,jpdic) 
     1729               alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto3d_balinc(:,:,:,jpalk) 
     1730               din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto3d_balinc(:,:,:,jpdin) 
     1731               sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto3d_balinc(:,:,:,jpsil) 
     1732            ENDIF 
    16721733         ELSE 
    16731734            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) 
     
    18611922            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) 
    18621923#endif 
    1863             IF ( ln_phytobal ) THEN 
     1924            IF ( ln_slchltotinc .OR. ln_schltotinc ) THEN 
    18641925               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) 
    18651926            ENDIF 
     
    18811942            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) 
    18821943#endif 
    1883             IF ( ln_phytobal ) THEN 
     1944            IF ( ln_slchltotinc .OR. ln_schltotinc ) THEN 
    18841945               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) 
    18851946            ENDIF 
     
    19031964            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) 
    19041965#endif 
    1905             IF ( ln_phytobal ) THEN 
     1966            IF ( ln_slchltotinc .OR. ln_schltotinc ) THEN 
    19061967               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) 
    19071968            ENDIF 
     
    19251986            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) 
    19261987#endif 
    1927             IF ( ln_phytobal ) THEN 
     1988            IF ( ln_slchltotinc .OR. ln_schltotinc ) THEN 
    19281989               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) 
    19291990            ENDIF 
     
    19452006            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' ) 
    19462007#endif 
    1947             IF ( ln_phytobal ) THEN 
     2008            IF ( ln_slchltotinc .OR. ln_schltotinc ) THEN 
    19482009               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it) 
    19492010            ENDIF 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto3dbal_medusa.F90

    r11958 r11990  
    1 MODULE asmphyto2dbal_medusa 
     1MODULE asmphyto3dbal_medusa 
    22   !!====================================================================== 
    33   !!                       ***  MODULE asmphyto2dbal_medusa  *** 
     
    3434   PRIVATE                    
    3535 
    36    PUBLIC asm_phyto2d_bal_medusa 
     36   PUBLIC asm_phyto3d_bal_medusa 
    3737 
    3838   ! Default values for biological assimilation parameters 
     
    6969CONTAINS 
    7070 
    71    SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      & 
    72       &                               pinc_chltot,                    & 
    73       &                               ld_chldia,                      & 
    74       &                               pinc_chldia,                    & 
    75       &                               ld_chlnon,                      & 
    76       &                               pinc_chlnon,                    & 
    77       &                               ld_phytot,                      & 
    78       &                               pinc_phytot,                    & 
    79       &                               ld_phydia,                      & 
    80       &                               pinc_phydia,                    & 
    81       &                               ld_phynon,                      & 
    82       &                               pinc_phynon,                    & 
     71   SUBROUTINE asm_phyto3d_bal_medusa( pinc_chltot,                    & 
    8372      &                               pincper,                        & 
    84       &                               p_maxchlinc, ld_phytobal, pmld, & 
     73      &                               p_maxchlinc,                   & 
    8574      &                               pgrow_avg_bkg, ploss_avg_bkg,   & 
    86       &                               phyt_avg_bkg, mld_max_bkg,      & 
    87       &                               tracer_bkg, phyto2d_balinc ) 
     75      &                               phyt_avg_bkg,                   & 
     76      &                               tracer_bkg, phyto3d_balinc ) 
    8877      !!--------------------------------------------------------------------------- 
    8978      !!                    ***  ROUTINE asm_phyto2d_bal_medusa  *** 
     
    10190      !!--------------------------------------------------------------------------- 
    10291      !! 
    103       LOGICAL,  INTENT(in   )                               :: ld_chltot      ! Assim chltot y/n 
    104       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chltot    ! chltot increments 
    105       LOGICAL,  INTENT(in   )                               :: ld_chldia      ! Assim chldia y/n 
    106       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chldia    ! chldia increments 
    107       LOGICAL,  INTENT(in   )                               :: ld_chlnon      ! Assim chlnon y/n 
    108       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chlnon    ! chlnon increments 
    109       LOGICAL,  INTENT(in   )                               :: ld_phytot      ! Assim phytot y/n 
    110       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phytot    ! phytot increments 
    111       LOGICAL,  INTENT(in   )                               :: ld_phydia      ! Assim phydia y/n 
    112       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phydia    ! phydia increments 
    113       LOGICAL,  INTENT(in   )                               :: ld_phynon      ! Assim phynon y/n 
    114       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phynon    ! phynon increments 
     92      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk)       :: pinc_chltot    ! chltot increments 
    11593      REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period 
    11694      REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment 
    117       LOGICAL,  INTENT(in   )                               :: ld_phytobal    ! Balancing y/n 
    118       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld           ! Mixed layer depth 
    119       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth 
    120       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss 
    121       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto 
    122       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD 
     95      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)       :: pgrow_avg_bkg  ! Avg phyto growth 
     96      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)       :: ploss_avg_bkg  ! Avg phyto loss 
     97      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk)       :: phyt_avg_bkg   ! Avg phyto 
    12398      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables 
    124       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments 
     99      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto3d_balinc ! Balancing increments 
    125100      !! 
    126101      INTEGER                                               :: ji, jj, jk, jn ! Loop counters 
     
    144119      REAL(wp)                                              :: zrat_pds_chd   ! Ratio of jppds:jpchd 
    145120      REAL(wp)                                              :: zrat_dtc_det   ! Ratio of jpdtc:jpdet 
    146       REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy 
     121      REAL(wp),                DIMENSION(jpi,jpj,jpk)       :: cchl_p         ! C:Chl for total phy 
     122      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p_2d      ! C:Chl for total phy 
    147123      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters 
    148124      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters 
    149125      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state 
     126      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: bstate_2d      ! Background state 
    150127      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments 
     128      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: outincs_2d     ! Balancing increments 
    151129      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics 
    152       REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics 
     130      REAL(wp),                DIMENSION(jpi,jpj,1,22)      :: diag_fulldepth ! Full-depth diagnostics 
     131      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chltot_2d 
     132      REAL(wp),                DIMENSION(jpi,jpj)           :: pgrow_avg_bkg_2d 
     133      REAL(wp),                DIMENSION(jpi,jpj)           :: ploss_avg_bkg_2d 
     134      REAL(wp),                DIMENSION(jpi,jpj)           :: phyt_avg_bkg_2d 
     135      REAL(wp),                DIMENSION(jpi,jpj,1)         :: tmask_2d 
     136      REAL(wp),                DIMENSION(jpi,jpj)           :: pmld_2d 
     137      REAL(wp),                DIMENSION(jpi,jpj)           :: mld_max_bkg_2d 
    153138      !!--------------------------------------------------------------------------- 
    154139 
    155140      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value 
    156141      IF ( p_maxchlinc > 0.0 ) THEN 
    157          IF ( ld_chltot ) THEN 
    158             DO jj = 1, jpj 
    159                DO ji = 1, jpi 
    160                   pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) ) 
    161                END DO 
    162             END DO 
    163          ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN 
    164             DO jj = 1, jpj 
    165                DO ji = 1, jpi 
    166                   pinc_chltot(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) 
    167                   pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) ) 
    168                   IF ( pinc_chltot(ji,jj) .NE. ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) ) ) THEN 
    169                      zfrac = pinc_chltot(ji,jj) / ( pinc_chldia(ji,jj) + pinc_chlnon(ji,jj) ) 
    170                      pinc_chldia(ji,jj) = pinc_chldia(ji,jj) * zfrac 
    171                      pinc_chlnon(ji,jj) = pinc_chlnon(ji,jj) * zfrac 
    172                   ENDIF 
    173                END DO 
    174             END DO 
    175          ELSE IF ( ld_chldia ) THEN 
    176             DO jj = 1, jpj 
    177                DO ji = 1, jpi 
    178                   pinc_chldia(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chldia(ji,jj), p_maxchlinc ) ) 
    179                   pinc_chltot(ji,jj) = pinc_chldia(ji,jj) 
    180                END DO 
    181             END DO 
    182          ELSE IF ( ld_chlnon ) THEN 
    183             DO jj = 1, jpj 
    184                DO ji = 1, jpi 
    185                   pinc_chlnon(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chlnon(ji,jj), p_maxchlinc ) ) 
    186                   pinc_chltot(ji,jj) = pinc_chlnon(ji,jj) 
    187                END DO 
    188             END DO 
    189          ENDIF 
    190       ENDIF 
    191  
    192       IF ( ld_phytot .OR. ld_phydia .OR. ld_phynon ) THEN 
    193          CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' ) 
    194       ENDIF 
    195        
    196       IF ( ld_phytobal ) THEN   ! Nitrogen balancing 
    197  
    198          ! Set up model parameters to be passed into Hemmings balancing routine. 
    199          ! For now these are hardwired to the standard HadOCC parameter values 
    200          ! (except C:N ratios) as this is what the scheme was developed for. 
    201          ! Obviously, HadOCC and MEDUSA are rather different models, so this 
    202          ! isn't ideal, but there's not always direct analogues between the two 
    203          ! parameter sets, so it's the easiest way to get something running. 
    204          ! In the longer term, some serious MarMOT-based development is required. 
    205          modparm(1)  = 0.1                               ! grow_sat 
    206          modparm(2)  = 2.0                               ! psmax 
    207          modparm(3)  = 0.845                             ! par 
    208          modparm(4)  = 0.02                              ! alpha 
    209          modparm(5)  = 0.05                              ! resp_rate 
    210          modparm(6)  = 0.05                              ! pmort_rate 
    211          modparm(7)  = 0.01                              ! phyto_min 
    212          modparm(8)  = 0.05                              ! z_mort_1 
    213          modparm(9)  = 1.0                               ! z_mort_2 
    214          modparm(10) = ( xthetapn  + xthetapd  ) / 2.0   ! c2n_p 
    215          modparm(11) = ( xthetazmi + xthetazme ) / 2.0   ! c2n_z 
    216          modparm(12) = xthetad                           ! c2n_d 
    217          modparm(13) = 0.01                              ! graze_threshold 
    218          modparm(14) = 2.0                               ! holling_coef 
    219          modparm(15) = 0.5                               ! graze_sat 
    220          modparm(16) = 2.0                               ! graze_max 
    221  
    222          ! Set up assimilation parameters to be passed into balancing routine 
    223          ! Not sure what assimparm(1) is meant to be, but it doesn't get used 
    224          assimparm(2)  = balnutext 
    225          assimparm(3)  = balnutmin 
    226          assimparm(4)  = r 
    227          assimparm(5)  = beta_g 
    228          assimparm(6)  = beta_l 
    229          assimparm(7)  = beta_m 
    230          assimparm(8)  = a_g 
    231          assimparm(9)  = a_l 
    232          assimparm(10) = a_m 
    233          assimparm(11) = zfracb0 
    234          assimparm(12) = zfracb1 
    235          assimparm(13) = qrfmax 
    236          assimparm(14) = qafmax 
    237          assimparm(15) = zrfmax 
    238          assimparm(16) = zafmax 
    239          assimparm(17) = prfmax 
    240          assimparm(18) = incphymin 
    241          assimparm(19) = integnstep 
    242          assimparm(20) = pthreshold 
    243  
    244          ! Set up external tracer indices array bstate 
    245          i_tracer(1) = 1   ! nutrient 
    246          i_tracer(2) = 2   ! phytoplankton 
    247          i_tracer(3) = 3   ! zooplankton 
    248          i_tracer(4) = 4   ! detritus 
    249          i_tracer(5) = 5   ! DIC 
    250          i_tracer(6) = 6   ! Alkalinity 
    251  
    252          ! Set background state 
    253          bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jpdin) 
    254          bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jpphn) + tracer_bkg(:,:,:,jpphd) 
    255          bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jpzmi) + tracer_bkg(:,:,:,jpzme) 
    256          bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jpdet) 
    257          bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jpdic) 
    258          bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jpalk) 
    259  
    260          ! Calculate carbon to chlorophyll ratio for combined phytoplankton 
    261          ! and nitrogen to biomass equivalent for PZD 
    262          ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA 
    263          cchl_p(:,:) = 0.0 
    264          DO jj = 1, jpj 
    265             DO ji = 1, jpi 
    266                IF ( ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) ) .GT. 0.0 ) THEN 
    267                   cchl_p(ji,jj) = xmassc * ( ( tracer_bkg(ji,jj,1,jpphn) * xthetapn ) +      & 
    268                      &                       ( tracer_bkg(ji,jj,1,jpphd) * xthetapd )   ) /  & 
    269                      &            ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) ) 
    270                ENDIF 
    271             END DO 
    272          END DO 
    273          n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) ) 
    274          n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) ) 
    275          n2be_d = 14.01 + ( xmassc * xthetad ) 
    276  
    277          ! Call nitrogen balancing routine 
    278          CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   & 
    279             &               n2be_p, n2be_z, n2be_d, assimparm,                      & 
    280             &               INT(pincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       & 
    281             &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot(:,:), cchl_p(:,:), & 
    282             &               nbal_active, phyt_avg_bkg(:,:),                         & 
    283             &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      & 
    284             &               subsurf_active, deepneg_active,                         & 
    285             &               deeppos_active, nutprof_active,                         & 
    286             &               bstate, outincs,                                        & 
    287             &               diag_active, diag,                                      & 
    288             &               diag_fulldepth_active, diag_fulldepth ) 
    289           
    290          ! Loop over each grid point partioning the increments 
    291          phyto2d_balinc(:,:,:,:) = 0.0 
    292142         DO jk = 1, jpk 
    293143            DO jj = 1, jpj 
    294144               DO ji = 1, jpi 
    295  
    296                   ! Phytoplankton 
    297                   IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. & 
    298                      & ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) .AND. & 
    299                      & ( pinc_chltot(ji,jj) /= 0.0 ) ) THEN 
    300                      IF ( ld_chltot ) THEN 
    301                         ! Phytoplankton nitrogen split up based on existing ratios 
    302                         zfrac_phn = tracer_bkg(ji,jj,jk,jpphn) / & 
    303                            &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd)) 
    304                         zfrac_phd = tracer_bkg(ji,jj,jk,jpphd) / & 
    305                            &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd)) 
    306                      ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN 
    307                         ! Phytoplankton nitrogen split up based on assimilation increments 
    308                         zfrac_phn = pinc_chlnon(ji,jj) / pinc_chltot(ji,jj) 
    309                         zfrac_phd = pinc_chldia(ji,jj) / pinc_chltot(ji,jj) 
    310                      ENDIF 
    311  
    312                      ! Phytoplankton silicate split up based on existing ratios 
    313                      zrat_pds_phd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpphd) 
    314                       
    315                      ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen 
    316                      ! Not using pinc_chltot directly as it's only 2D 
    317                      ! This method should give same results at surface as splitting pinc_chltot would 
    318                      zrat_chn_phn = tracer_bkg(ji,jj,jk,jpchn) / tracer_bkg(ji,jj,jk,jpphn) 
    319                      zrat_chd_phd = tracer_bkg(ji,jj,jk,jpchd) / tracer_bkg(ji,jj,jk,jpphd) 
    320                       
    321                      phyto2d_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn 
    322                      phyto2d_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd 
    323                      phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_pds_phd 
    324                      phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,jk,jpphn) * zrat_chn_phn 
    325                      phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,jk,jpphd) * zrat_chd_phd 
    326                   ENDIF 
    327  
    328                   ! Zooplankton nitrogen split up based on existing ratios 
    329                   IF ( ( tracer_bkg(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpzme) > 0.0 ) ) THEN 
    330                      zfrac_zmi = tracer_bkg(ji,jj,jk,jpzmi) / & 
    331                         &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme)) 
    332                      zfrac_zme = tracer_bkg(ji,jj,jk,jpzme) / & 
    333                         &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme)) 
    334                      phyto2d_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi 
    335                      phyto2d_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme 
    336                   ENDIF 
    337  
    338                   ! Nitrogen nutrient straight from balancing scheme 
    339                   phyto2d_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1)) 
    340  
    341                   ! Nitrogen detritus straight from balancing scheme 
    342                   phyto2d_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4)) 
    343  
    344                   ! DIC straight from balancing scheme 
    345                   phyto2d_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5)) 
    346  
    347                   ! Alkalinity straight from balancing scheme 
    348                   phyto2d_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6)) 
    349  
    350                   ! Remove diatom silicate increment from nutrient silicate to conserve mass 
    351                   IF ( ( tracer_bkg(ji,jj,jk,jpsil) - phyto2d_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN 
    352                      phyto2d_balinc(ji,jj,jk,jpsil) = phyto2d_balinc(ji,jj,jk,jppds) * (-1.0) 
    353                   ENDIF 
    354  
    355                   ! Carbon detritus based on existing ratios 
    356                   IF ( ( tracer_bkg(ji,jj,jk,jpdet) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpdtc) > 0.0 ) ) THEN 
    357                      zrat_dtc_det = tracer_bkg(ji,jj,jk,jpdtc) / tracer_bkg(ji,jj,jk,jpdet) 
    358                      phyto2d_balinc(ji,jj,jk,jpdtc) = phyto2d_balinc(ji,jj,jk,jpdet) * zrat_dtc_det 
    359                   ENDIF 
    360  
    361                   ! Do nothing with iron or oxygen for the time being 
    362                   phyto2d_balinc(ji,jj,jk,jpfer) = 0.0 
    363                   phyto2d_balinc(ji,jj,jk,jpoxy) = 0.0 
    364                    
     145                  pinc_chltot(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj,jk), p_maxchlinc ) ) 
    365146               END DO 
    366147            END DO 
    367148         END DO 
    368        
    369       ELSE   ! No nitrogen balancing 
    370        
    371          ! Initialise individual chlorophyll increments to zero 
    372          phyto2d_balinc(:,:,:,jpchn) = 0.0 
    373          phyto2d_balinc(:,:,:,jpchd) = 0.0 
    374           
    375          ! Split up total surface chlorophyll increments 
     149      ENDIF 
     150 
     151      ! Set up model parameters to be passed into Hemmings balancing routine. 
     152      ! For now these are hardwired to the standard HadOCC parameter values 
     153      ! (except C:N ratios) as this is what the scheme was developed for. 
     154      ! Obviously, HadOCC and MEDUSA are rather different models, so this 
     155      ! isn't ideal, but there's not always direct analogues between the two 
     156      ! parameter sets, so it's the easiest way to get something running. 
     157      ! In the longer term, some serious MarMOT-based development is required. 
     158      modparm(1)  = 0.1                               ! grow_sat 
     159      modparm(2)  = 2.0                               ! psmax 
     160      modparm(3)  = 0.845                             ! par 
     161      modparm(4)  = 0.02                              ! alpha 
     162      modparm(5)  = 0.05                              ! resp_rate 
     163      modparm(6)  = 0.05                              ! pmort_rate 
     164      modparm(7)  = 0.01                              ! phyto_min 
     165      modparm(8)  = 0.05                              ! z_mort_1 
     166      modparm(9)  = 1.0                               ! z_mort_2 
     167      modparm(10) = ( xthetapn  + xthetapd  ) / 2.0   ! c2n_p 
     168      modparm(11) = ( xthetazmi + xthetazme ) / 2.0   ! c2n_z 
     169      modparm(12) = xthetad                           ! c2n_d 
     170      modparm(13) = 0.01                              ! graze_threshold 
     171      modparm(14) = 2.0                               ! holling_coef 
     172      modparm(15) = 0.5                               ! graze_sat 
     173      modparm(16) = 2.0                               ! graze_max 
     174 
     175      ! Set up assimilation parameters to be passed into balancing routine 
     176      ! Not sure what assimparm(1) is meant to be, but it doesn't get used 
     177      assimparm(2)  = balnutext 
     178      assimparm(3)  = balnutmin 
     179      assimparm(4)  = r 
     180      assimparm(5)  = beta_g 
     181      assimparm(6)  = beta_l 
     182      assimparm(7)  = beta_m 
     183      assimparm(8)  = a_g 
     184      assimparm(9)  = a_l 
     185      assimparm(10) = a_m 
     186      assimparm(11) = zfracb0 
     187      assimparm(12) = zfracb1 
     188      assimparm(13) = qrfmax 
     189      assimparm(14) = qafmax 
     190      assimparm(15) = zrfmax 
     191      assimparm(16) = zafmax 
     192      assimparm(17) = prfmax 
     193      assimparm(18) = incphymin 
     194      assimparm(19) = integnstep 
     195      assimparm(20) = pthreshold 
     196 
     197      ! Set up external tracer indices array bstate 
     198      i_tracer(1) = 1   ! nutrient 
     199      i_tracer(2) = 2   ! phytoplankton 
     200      i_tracer(3) = 3   ! zooplankton 
     201      i_tracer(4) = 4   ! detritus 
     202      i_tracer(5) = 5   ! DIC 
     203      i_tracer(6) = 6   ! Alkalinity 
     204 
     205      ! Set background state 
     206      bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jpdin) 
     207      bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jpphn) + tracer_bkg(:,:,:,jpphd) 
     208      bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jpzmi) + tracer_bkg(:,:,:,jpzme) 
     209      bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jpdet) 
     210      bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jpdic) 
     211      bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jpalk) 
     212 
     213      ! Calculate carbon to chlorophyll ratio for combined phytoplankton 
     214      ! and nitrogen to biomass equivalent for PZD 
     215      ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA 
     216      cchl_p(:,:,:) = 0.0 
     217      DO jk = 1, jpk 
    376218         DO jj = 1, jpj 
    377219            DO ji = 1, jpi 
    378                IF ( ( tracer_bkg(ji,jj,1,jpchn) > 0.0 ) .AND. & 
    379                   & ( tracer_bkg(ji,jj,1,jpchd) > 0.0 ) ) THEN 
    380                   IF ( ld_chltot ) THEN 
    381                      ! Chlorophyll split up based on existing ratios 
    382                      zfrac_chn = tracer_bkg(ji,jj,1,jpchn) / & 
    383                         &        ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) ) 
    384                      zfrac_chd = tracer_bkg(ji,jj,1,jpchd) / & 
    385                         &        ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd) ) 
    386                      phyto2d_balinc(ji,jj,1,jpchn) = pinc_chltot(ji,jj) * zfrac_chn 
    387                      phyto2d_balinc(ji,jj,1,jpchd) = pinc_chltot(ji,jj) * zfrac_chd 
    388                   ENDIF 
    389                   IF( ld_chldia ) THEN 
    390                      phyto2d_balinc(ji,jj,1,jpchd) = pinc_chldia(ji,jj) 
    391                   ENDIF 
    392                   IF( ld_chlnon ) THEN 
    393                      phyto2d_balinc(ji,jj,1,jpchn) = pinc_chlnon(ji,jj) 
    394                   ENDIF 
    395                    
    396                   ! Maintain stoichiometric ratios of nitrogen and silicate 
    397                   IF ( ld_chltot .OR. ld_chlnon ) THEN 
    398                      zrat_phn_chn = tracer_bkg(ji,jj,1,jpphn) / tracer_bkg(ji,jj,1,jpchn) 
    399                      phyto2d_balinc(ji,jj,1,jpphn) = phyto2d_balinc(ji,jj,1,jpchn) * zrat_phn_chn 
    400                   ENDIF 
    401                   IF ( ld_chltot .OR. ld_chldia ) THEN 
    402                      zrat_phd_chd = tracer_bkg(ji,jj,1,jpphd) / tracer_bkg(ji,jj,1,jpchd) 
    403                      phyto2d_balinc(ji,jj,1,jpphd) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_phd_chd 
    404                      zrat_pds_chd = tracer_bkg(ji,jj,1,jppds) / tracer_bkg(ji,jj,1,jpchd) 
    405                      phyto2d_balinc(ji,jj,1,jppds) = phyto2d_balinc(ji,jj,1,jpchd) * zrat_pds_chd 
    406                   ENDIF 
     220               IF ( ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd ) ) .GT. 0.0 ) THEN 
     221                  cchl_p(ji,jj,jk) = xmassc * ( ( tracer_bkg(ji,jj,jk,jpphn) * xthetapn ) +      & 
     222                     &                       ( tracer_bkg(ji,jj,jk,jpphd) * xthetapd )   ) /  & 
     223                     &            ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd ) ) 
    407224               ENDIF 
    408225            END DO 
    409226         END DO 
     227      END DO 
     228      n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) ) 
     229      n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) ) 
     230      n2be_d = 14.01 + ( xmassc * xthetad ) 
     231 
     232      pmld_2d(:,:)        = 0.5 
     233      mld_max_bkg_2d(:,:) = 0.5 
     234 
     235      DO jk = 1, jpk 
     236 
     237         tmask_2d(:,:,1)       = tmask(:,:,jk) 
     238         pinc_chltot_2d(:,:)   = pinc_chltot(:,:,jk) 
     239         cchl_p_2d(:,:)        = cchl_p(:,:,jk) 
     240         phyt_avg_bkg_2d(:,:)  = phyt_avg_bkg(:,:,jk) 
     241         pgrow_avg_bkg_2d(:,:) = pgrow_avg_bkg(:,:,jk) 
     242         ploss_avg_bkg_2d(:,:) = ploss_avg_bkg(:,:,jk) 
     243         bstate_2d(:,:,1,:)    = bstate(:,:,jk,:) 
     244         outincs_2d(:,:,:,:)   = 0.0 
    410245          
    411          ! Propagate through mixed layer 
     246         ! Call nitrogen balancing routine 
     247         CALL bio_analysis( jpi, jpj, 1, gdepw_n(:,:,2), i_tracer, modparm,   & 
     248            &               n2be_p, n2be_z, n2be_d, assimparm,                      & 
     249            &               INT(pincper), 1, INT(SUM(tmask_2d,3)), tmask_2d(:,:,:),       & 
     250            &               pmld_2d(:,:), mld_max_bkg_2d(:,:), pinc_chltot_2d(:,:), cchl_p_2d(:,:), & 
     251            &               nbal_active, phyt_avg_bkg_2d(:,:),                         & 
     252            &               gl_active, pgrow_avg_bkg_2d(:,:), ploss_avg_bkg_2d(:,:),      & 
     253            &               subsurf_active, deepneg_active,                         & 
     254            &               deeppos_active, nutprof_active,                         & 
     255            &               bstate_2d, outincs_2d,                                        & 
     256            &               diag_active, diag,                                      & 
     257            &               diag_fulldepth_active, diag_fulldepth ) 
     258 
     259         outincs(:,:,jk,:) = outincs_2d(:,:,1,:) 
     260 
     261      END DO 
     262 
     263      ! Loop over each grid point partioning the increments 
     264      phyto3d_balinc(:,:,:,:) = 0.0 
     265      DO jk = 1, jpk 
    412266         DO jj = 1, jpj 
    413267            DO ji = 1, jpi 
    414                ! 
    415                jkmax = jpk-1 
    416                DO jk = jpk-1, 1, -1 
    417                   IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. & 
    418                      & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN 
    419                      pmld(ji,jj) = gdepw_n(ji,jj,jk+1) 
    420                      jkmax = jk 
    421                   ENDIF 
    422                END DO 
    423                ! 
    424                DO jk = 2, jkmax 
    425                   phyto2d_balinc(ji,jj,jk,jpchn) = phyto2d_balinc(ji,jj,1,jpchn) 
    426                   phyto2d_balinc(ji,jj,jk,jpchd) = phyto2d_balinc(ji,jj,1,jpchd) 
    427                   phyto2d_balinc(ji,jj,jk,jpphn) = phyto2d_balinc(ji,jj,1,jpphn) 
    428                   phyto2d_balinc(ji,jj,jk,jpphd) = phyto2d_balinc(ji,jj,1,jpphd) 
    429                   phyto2d_balinc(ji,jj,jk,jppds) = phyto2d_balinc(ji,jj,1,jppds) 
    430                END DO 
    431                ! 
     268 
     269               ! Phytoplankton 
     270               IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. & 
     271                  & ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) .AND. & 
     272                  & ( pinc_chltot(ji,jj,jk) /= 0.0 ) ) THEN 
     273                  ! Phytoplankton nitrogen split up based on existing ratios 
     274                  zfrac_phn = tracer_bkg(ji,jj,jk,jpphn) / & 
     275                     &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd)) 
     276                  zfrac_phd = tracer_bkg(ji,jj,jk,jpphd) / & 
     277                     &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd)) 
     278 
     279                  ! Phytoplankton silicate split up based on existing ratios 
     280                  zrat_pds_phd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpphd) 
     281 
     282                  ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen 
     283                  ! Not using pinc_chltot directly as it's only 2D 
     284                  ! This method should give same results at surface as splitting pinc_chltot would 
     285                  zrat_chn_phn = tracer_bkg(ji,jj,jk,jpchn) / tracer_bkg(ji,jj,jk,jpphn) 
     286                  zrat_chd_phd = tracer_bkg(ji,jj,jk,jpchd) / tracer_bkg(ji,jj,jk,jpphd) 
     287 
     288                  phyto3d_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn 
     289                  phyto3d_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd 
     290                  phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpphd) * zrat_pds_phd 
     291                  phyto3d_balinc(ji,jj,jk,jpchn) = phyto3d_balinc(ji,jj,jk,jpphn) * zrat_chn_phn 
     292                  phyto3d_balinc(ji,jj,jk,jpchd) = phyto3d_balinc(ji,jj,jk,jpphd) * zrat_chd_phd 
     293               ENDIF 
     294 
     295               ! Zooplankton nitrogen split up based on existing ratios 
     296               IF ( ( tracer_bkg(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpzme) > 0.0 ) ) THEN 
     297                  zfrac_zmi = tracer_bkg(ji,jj,jk,jpzmi) / & 
     298                     &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme)) 
     299                  zfrac_zme = tracer_bkg(ji,jj,jk,jpzme) / & 
     300                     &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme)) 
     301                  phyto3d_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi 
     302                  phyto3d_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme 
     303               ENDIF 
     304 
     305               ! Nitrogen nutrient straight from balancing scheme 
     306               phyto3d_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1)) 
     307 
     308               ! Nitrogen detritus straight from balancing scheme 
     309               phyto3d_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4)) 
     310 
     311               ! DIC straight from balancing scheme 
     312               phyto3d_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5)) 
     313 
     314               ! Alkalinity straight from balancing scheme 
     315               phyto3d_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6)) 
     316 
     317               ! Remove diatom silicate increment from nutrient silicate to conserve mass 
     318               IF ( ( tracer_bkg(ji,jj,jk,jpsil) - phyto3d_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN 
     319                  phyto3d_balinc(ji,jj,jk,jpsil) = phyto3d_balinc(ji,jj,jk,jppds) * (-1.0) 
     320               ENDIF 
     321 
     322               ! Carbon detritus based on existing ratios 
     323               IF ( ( tracer_bkg(ji,jj,jk,jpdet) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpdtc) > 0.0 ) ) THEN 
     324                  zrat_dtc_det = tracer_bkg(ji,jj,jk,jpdtc) / tracer_bkg(ji,jj,jk,jpdet) 
     325                  phyto3d_balinc(ji,jj,jk,jpdtc) = phyto3d_balinc(ji,jj,jk,jpdet) * zrat_dtc_det 
     326               ENDIF 
     327 
     328               ! Do nothing with iron or oxygen for the time being 
     329               phyto3d_balinc(ji,jj,jk,jpfer) = 0.0 
     330               phyto3d_balinc(ji,jj,jk,jpoxy) = 0.0 
     331 
    432332            END DO 
    433333         END DO 
    434  
    435          ! Set other balancing increments to zero 
    436          phyto2d_balinc(:,:,:,jpzmi) = 0.0 
    437          phyto2d_balinc(:,:,:,jpzme) = 0.0 
    438          phyto2d_balinc(:,:,:,jpdin) = 0.0 
    439          phyto2d_balinc(:,:,:,jpsil) = 0.0 
    440          phyto2d_balinc(:,:,:,jpfer) = 0.0 
    441          phyto2d_balinc(:,:,:,jpdet) = 0.0 
    442          phyto2d_balinc(:,:,:,jpdtc) = 0.0 
    443          phyto2d_balinc(:,:,:,jpdic) = 0.0 
    444          phyto2d_balinc(:,:,:,jpalk) = 0.0 
    445          phyto2d_balinc(:,:,:,jpoxy) = 0.0 
    446  
    447       ENDIF 
     334      END DO 
    448335       
    449336      ! If performing extra tidal mixing in the Indonesian Throughflow, 
     
    453340         DO jn = 1, jptra 
    454341            DO jk = 1, jpk 
    455                phyto2d_balinc(:,:,jk,jn) = phyto2d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) ) 
     342               phyto3d_balinc(:,:,jk,jn) = phyto3d_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) ) 
    456343            END DO 
    457344         END DO 
    458345      ENDIF 
    459346 
    460    END SUBROUTINE asm_phyto2d_bal_medusa 
     347   END SUBROUTINE asm_phyto3d_bal_medusa 
    461348 
    462349#else 
     
    465352   !!---------------------------------------------------------------------- 
    466353CONTAINS 
    467    SUBROUTINE asm_phyto2d_bal_medusa( ld_chltot,                      & 
    468       &                              pinc_chltot,                    & 
    469       &                              ld_chldia,                      & 
    470       &                              pinc_chldia,                    & 
    471       &                              ld_chlnon,                      & 
    472       &                              pinc_chlnon,                    & 
    473       &                              ld_phytot,                      & 
    474       &                              pinc_phytot,                    & 
    475       &                              ld_phydia,                      & 
    476       &                              pinc_phydia,                    & 
    477       &                              ld_phynon,                      & 
    478       &                              pinc_phynon,                    & 
     354   SUBROUTINE asm_phyto3d_bal_medusa(pinc_chltot,                    & 
    479355      &                              pincper,                        & 
    480       &                              p_maxchlinc, ld_phytobal, pmld, & 
     356      &                              p_maxchlinc,                   & 
    481357      &                              pgrow_avg_bkg, ploss_avg_bkg,   & 
    482       &                              phyt_avg_bkg, mld_max_bkg,      & 
    483       &                              tracer_bkg, phyto2d_balinc ) 
    484       LOGICAL :: ld_chltot 
    485       REAL    :: pinc_chltot(:,:) 
    486       LOGICAL :: ld_chldia 
    487       REAL    :: pinc_chldia(:,:) 
    488       LOGICAL :: ld_chlnon 
    489       REAL    :: pinc_chlnon(:,:) 
    490       LOGICAL :: ld_phytot 
    491       REAL    :: pinc_phytot(:,:) 
    492       LOGICAL :: ld_phydia 
    493       REAL    :: pinc_phydia(:,:) 
    494       LOGICAL :: ld_phynon 
    495       REAL    :: pinc_phynon(:,:) 
     358      &                              phyt_avg_bkg,                   & 
     359      &                              tracer_bkg, phyto3d_balinc ) 
     360      REAL    :: pinc_chltot(:,:,:) 
    496361      REAL    :: pincper 
    497362      REAL    :: p_maxchlinc 
    498       LOGICAL :: ld_phytobal 
    499       REAL    :: pmld(:,:) 
    500       REAL    :: pgrow_avg_bkg(:,:) 
    501       REAL    :: ploss_avg_bkg(:,:) 
    502       REAL    :: phyt_avg_bkg(:,:) 
    503       REAL    :: mld_max_bkg(:,:) 
     363      REAL    :: pgrow_avg_bkg(:,:,:) 
     364      REAL    :: ploss_avg_bkg(:,:,:) 
     365      REAL    :: phyt_avg_bkg(:,:,:) 
    504366      REAL    :: tracer_bkg(:,:,:,:) 
    505367      REAL    :: phyto2d_balinc(:,:,:,:) 
    506       WRITE(*,*) 'asm_phyto2d_bal_medusa: You should not have seen this print! error?' 
    507    END SUBROUTINE asm_phyto2d_bal_medusa 
     368      WRITE(*,*) 'asm_phyto3d_bal_medusa: You should not have seen this print! error?' 
     369   END SUBROUTINE asm_phyto3d_bal_medusa 
    508370#endif 
    509371 
    510372   !!====================================================================== 
    511 END MODULE asmphyto2dbal_medusa 
     373END MODULE asmphyto3dbal_medusa 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_fin.F90

    r9292 r11990  
    5151                                   mld_max, pgrow_avg,                  & 
    5252                                   ploss_avg, phyt_avg,                 & 
     53                                   ploss_avg_3d, phyt_avg_3d,           & 
     54                                   pgrow_avg_3d,                        & 
    5355# endif 
    5456                                   za_sed_c, za_sed_ca, za_sed_fe,      & 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/bio_medusa_init.F90

    r9292 r11990  
    3636      USE par_oce,           ONLY: jpi, jpj, jpk 
    3737# if defined key_foam_medusa 
    38       USE sms_medusa,        ONLY: jdms, pgrow_avg, ploss_avg, phyt_avg, mld_max 
     38      USE sms_medusa,        ONLY: jdms, pgrow_avg, ploss_avg, phyt_avg, mld_max, & 
     39         &                         pgrow_avg_3d, ploss_avg_3d, phyt_avg_3d 
    3940# else 
    4041      USE sms_medusa,        ONLY: jdms 
     
    201202      ploss_avg(:,:) = 0.0 
    202203      phyt_avg(:,:)  = 0.0 
     204      pgrow_avg_3d(:,:,:) = 0.0 
     205      ploss_avg_3d(:,:,:) = 0.0 
     206      phyt_avg_3d(:,:,:)  = 0.0 
    203207      IF( kt == nittrc000 ) THEN 
    204208         mld_max(:,:) = 0.0 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/plankton.F90

    r9292 r11990  
    5252# if defined key_foam_medusa 
    5353                                   pgrow_avg, ploss_avg, phyt_avg,         & 
     54                                   pgrow_avg_3d, ploss_avg_3d, phyt_avg_3d, & 
    5455# endif 
    5556                                   xkphd, xkphn, xkzme, xkzmi,             & 
     
    242243                                  ((zphn(ji,jj) + zphd(ji,jj)) *       & 
    243244                                   fse3t(ji,jj,jk) * fq0) 
     245               !! 
     246               pgrow_avg_3d(ji,jj,jk) = (fprn(ji,jj) * zphn(ji,jj)) +  & 
     247                                        (fprd(ji,jj) * zphd(ji,jj)) 
     248               ploss_avg_3d(ji,jj,jk) = fgmepd(ji,jj) + fdpd(ji,jj) +  & 
     249                                        fdpd2(ji,jj)                +  & 
     250                                        fgmepn(ji,jj) + fdpn(ji,jj) +  & 
     251                                        fdpn2(ji,jj)  + fgmipn(ji,jj) 
     252               phyt_avg_3d(ji,jj,jk)  = zphn(ji,jj) + zphd(ji,jj) 
    244253            ENDIF 
    245254         ENDDO 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90

    r9292 r11990  
    458458   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ploss_avg  !: Mixed layer average phytoplankton loss 
    459459   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: phyt_avg   !: Mixed layer average phytoplankton 
     460   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: pgrow_avg_3d  !: Mixed layer average phytoplankton growth 
     461   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ploss_avg_3d  !: Mixed layer average phytoplankton loss 
     462   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: phyt_avg_3d   !: Mixed layer average phytoplankton 
    460463   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_max    !: Maximum mixed layer depth 
    461464!! 
     
    537540# if defined key_foam_medusa 
    538541      ALLOCATE( pgrow_avg(jpi,jpj)   , ploss_avg(jpi,jpj)   ,       & 
     542         &      pgrow_avg_3d(jpi,jpj,jpk) , ploss_avg_3d(jpi,jpj,jpk) , & 
     543         &      phyt_avg_3d(jpi,jpj,jpk)  ,                         & 
    539544         &      phyt_avg(jpi,jpj)    , mld_max(jpi,jpj)     ,    STAT=ierr(9) ) 
    540545# endif 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r9292 r11990  
    105105      USE sms_medusa,                 ONLY: hist_pco2, xobs_xco2a,          & 
    106106                                            pgrow_avg, ploss_avg,           & 
     107                                            pgrow_avg_3d, ploss_avg_3d,     & 
     108                                            phyt_avg_3d,                    & 
    107109                                            phyt_avg, mld_max 
    108110# else 
  • branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_3dnitbal/NEMOGCM/NEMO/TOP_SRC/trcrst.F90

    r9292 r11990  
    370370         mld_max(:,:)   = 0.0 
    371371      ENDIF 
     372      IF( iom_varid( numrtr, 'pgrow_avg_3d', ldstop = .FALSE. ) > 0 ) THEN 
     373         IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg present - reading in ...' 
     374         CALL iom_get( numrtr, jpdom_autoglo, 'pgrow_avg_3d',  pgrow_avg_3d(:,:,:)  ) 
     375         CALL iom_get( numrtr, jpdom_autoglo, 'ploss_avg_3d',  ploss_avg_3d(:,:,:)  ) 
     376         CALL iom_get( numrtr, jpdom_autoglo, 'phyt_avg_3d',   phyt_avg_3d(:,:,:)   ) 
     377      ELSE 
     378         IF(lwp) WRITE(numout,*) ' MEDUSA pgrow_avg_3d absent - setting to zero ...' 
     379         pgrow_avg_3d(:,:,:) = 0.0 
     380         ploss_avg_3d(:,:,:) = 0.0 
     381         phyt_avg_3d(:,:,:)  = 0.0 
     382      ENDIF 
    372383# endif 
    373384 
     
    555566      CALL iom_rstput( kt, nitrst, numrtw, 'phyt_avg',  phyt_avg(:,:)  ) 
    556567      CALL iom_rstput( kt, nitrst, numrtw, 'mld_max',   mld_max(:,:)   ) 
     568      CALL iom_rstput( kt, nitrst, numrtw, 'pgrow_avg_3d', pgrow_avg_3d(:,:,:) ) 
     569      CALL iom_rstput( kt, nitrst, numrtw, 'ploss_avg_3d', ploss_avg_3d(:,:,:) ) 
     570      CALL iom_rstput( kt, nitrst, numrtw, 'phyt_avg_3d',  phyt_avg_3d(:,:,:)  ) 
    557571# endif 
    558572!! 
Note: See TracChangeset for help on using the changeset viewer.