New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
asmphyto2dbal_ersem.F90 in branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_ersem.F90 @ 10667

Last change on this file since 10667 was 10667, checked in by dford, 5 years ago

Reference FABM indices to jp_fabm_m1.
File size: 53.2 KB
Line 
1MODULE asmphyto2dbal_ersem
2   !!======================================================================
3   !!                       ***  MODULE asmphyto2dbal_ersem  ***
4   !! Calculate increments to ERSEM based on surface phyto2d increments
5   !!
6   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
7   !! For licensing reasons this is kept in its own internal Met Office
8   !! branch (dev/frdf/vn3.6_nitrogen_balancing) rather than in the Paris
9   !! repository, and must be merged in when building.
10   !!
11   !!======================================================================
12   !! History :  3.6  ! 2019-01 (D. Ford)  Adapted from asmphyto2dbal_medusa
13   !!----------------------------------------------------------------------
14#if defined key_asminc && defined key_fabm
15   !!----------------------------------------------------------------------
16   !! 'key_asminc'          : assimilation increment interface
17   !! 'key_fabm'            : FABM-ERSEM model
18   !!----------------------------------------------------------------------
19   !! asm_phyto2d_bal_ersem : routine to calculate increments to ERSEM
20   !!----------------------------------------------------------------------
21   USE par_kind,      ONLY: wp             ! kind parameters
22   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
23   USE dom_oce,       ONLY: gdepw_n        ! domain information
24   USE iom                                 ! i/o
25   USE par_fabm                            ! FABM-ERSEM parameters
26   USE par_trc,       ONLY: jptra          ! Tracer parameters
27   USE bioanalysis                         ! Nitrogen balancing
28
29   IMPLICIT NONE
30   PRIVATE                   
31
32   PUBLIC asm_phyto2d_bal_ersem
33
34   ! Default values for biological assimilation parameters
35   ! Should match Hemmings et al. (2008)
36   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
37   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
38   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
39   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
40   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
41   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
42   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
43   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
44   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
45   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
46   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
47   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
48   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
49   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
50   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
51   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
52   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
53   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
54   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
55   !
56   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
57   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
58   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
59   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
60   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
61   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
62   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
63   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
64
65CONTAINS
66
67   SUBROUTINE asm_phyto2d_bal_ersem(  ld_chltot,                      &
68      &                               pinc_chltot,                    &
69      &                               ld_chldia,                      &
70      &                               pinc_chldia,                    &
71      &                               ld_chlnan,                      &
72      &                               pinc_chlnan,                    &
73      &                               ld_chlpic,                      &
74      &                               pinc_chlpic,                    &
75      &                               ld_chldin,                      &
76      &                               pinc_chldin,                    &
77      &                               pincper,                        &
78      &                               p_maxchlinc, ld_phytobal, pmld, &
79      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
80      &                               phyt_avg_bkg, mld_max_bkg,      &
81      &                               totalk_bkg,                     &
82      &                               tracer_bkg, phyto2d_balinc )
83      !!---------------------------------------------------------------------------
84      !!                    ***  ROUTINE asm_phyto2d_bal_ersem  ***
85      !!
86      !! ** Purpose :   calculate increments to ERSEM from 2d phytoplankton increments
87      !!
88      !! ** Method  :   EITHER (ld_phytobal == .TRUE.):
89      !!                average up ERSEM to look like HadOCC
90      !!                call nitrogen balancing scheme
91      !!                separate back out to MEDUSA
92      !!                OR (ld_phytobal == .FALSE.):
93      !!                calculate increments to maintain background stoichiometry
94      !!
95      !! ** Action  :   populate phyto2d_balinc
96      !!
97      !! References :   Hemmings et al., 2008, J. Mar. Res.
98      !!                Ford et al., 2012, Ocean Sci.
99      !!                Skakala et al., 2018, JGR
100      !!---------------------------------------------------------------------------
101      !!
102      LOGICAL,  INTENT(in   )                               :: ld_chltot      ! Assim chltot y/n
103      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chltot    ! chltot increments
104      LOGICAL,  INTENT(in   )                               :: ld_chldia      ! Assim chldia y/n
105      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chldia    ! chldia 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
112      REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period
113      REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment
114      LOGICAL,  INTENT(in   )                               :: ld_phytobal    ! Balancing y/n
115      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld           ! Mixed layer depth
116      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth
117      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss
118      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto
119      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
121      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables
122      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments
123      !!
124      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
125      INTEGER                                               :: jkmax          ! Loop index
126      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)
130      REAL(wp)                                              :: n2be_p         ! N:biomass for total phy
131      REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo
132      REAL(wp)                                              :: n2be_d         ! N:biomass for detritus
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  !
189      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy
190      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
191      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters
192      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state
193      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments
194      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics
195      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics
196      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chltot_temp
197      !!---------------------------------------------------------------------------
198     
199      ! Set parameters, firstly molar mass of carbon and nitrogen
200      zmassc = 12.01
201      zmassn = 14.01
202     
203      ! Then mesozooplankton nitrogen(mmol):carbon(mg), ERSEM parameter Z4/qnc
204      ! Hardwire it for now due to difficulty of getting at parameters through FABM
205      ! I think the following should work:
206      ! z4qnc = model%state_variables(jp_fabm_z4c)%properties%get_property_by_name('qnc')%value
207      ! but, get_property_by_name is in the module fabm_properties, which we can't get at directly
208      ! We can do a "USE fabm" which itself does a "USE fabm_properties", but doing it indirectly
209      ! like that means get_property_by_name is treated as private so we can't use it
210      ! So hardwire to value used in SSB runs and v4/5 CMEMS reanalysis
211      z4qnc  = 0.0126
212
213      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
214      IF ( p_maxchlinc > 0.0 ) THEN
215         IF ( ld_chltot ) THEN
216            DO jj = 1, jpj
217               DO ji = 1, jpi
218                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
219               END DO
220            END DO
221         ELSE
222            DO jj = 1, jpj
223               DO ji = 1, jpi
224                  IF ( ld_chldia .AND. ld_chlnan .AND. ld_chlpic .AND. ld_chldin ) THEN
225                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
226                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
227                  ELSE IF ( ld_chldia .AND. ld_chlnan .AND. ld_chlpic ) THEN
228                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
229                        &                      pinc_chlpic(ji,jj)
230                  ELSE IF ( ld_chldia .AND. ld_chlnan .AND. ld_chldin ) THEN
231                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
232                        &                      pinc_chldin(ji,jj)
233                  ELSE IF ( ld_chldia .AND. ld_chlpic .AND. ld_chldin ) THEN
234                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + &
235                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
236                  ELSE IF ( ld_chlnan .AND. ld_chlpic .AND. ld_chldin ) THEN
237                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + &
238                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
239                  ELSE IF ( ld_chldia .AND. ld_chlnan ) THEN
240                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj)
241                  ELSE IF ( ld_chldia .AND. ld_chlpic ) THEN
242                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlpic(ji,jj)
243                  ELSE IF ( ld_chldia .AND. ld_chldin ) THEN
244                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chldin(ji,jj)
245                  ELSE IF ( ld_chlnan .AND. ld_chlpic ) THEN
246                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + pinc_chlpic(ji,jj)
247                  ELSE IF ( ld_chlnan .AND. ld_chldin ) THEN
248                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + pinc_chldin(ji,jj)
249                  ELSE IF ( ld_chlpic .AND. ld_chldin ) THEN
250                     pinc_chltot_temp(ji,jj) = pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
251                  ELSE IF ( ld_chldia ) THEN
252                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj)
253                  ELSE IF ( ld_chlnan ) THEN
254                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj)
255                  ELSE IF ( ld_chlpic ) THEN
256                     pinc_chltot_temp(ji,jj) = pinc_chlpic(ji,jj)
257                  ELSE IF ( ld_chldin ) THEN
258                     pinc_chltot_temp(ji,jj) = pinc_chldin(ji,jj)
259                  ENDIF
260                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot_temp(ji,jj), p_maxchlinc ) )
261                  IF ( pinc_chltot(ji,jj) .NE. pinc_chltot_temp(ji,jj) ) THEN
262                     zfrac = pinc_chltot(ji,jj) / pinc_chltot_temp(ji,jj)
263                     IF ( ld_chldia ) THEN
264                        pinc_chldia(ji,jj) = pinc_chldia(ji,jj) * zfrac
265                     ENDIF
266                     IF ( ld_chlnan ) THEN
267                        pinc_chlnan(ji,jj) = pinc_chlnan(ji,jj) * zfrac
268                     ENDIF
269                     IF ( ld_chlpic ) THEN
270                        pinc_chlpic(ji,jj) = pinc_chlpic(ji,jj) * zfrac
271                     ENDIF
272                     IF ( ld_chldin ) THEN
273                        pinc_chldin(ji,jj) = pinc_chldin(ji,jj) * zfrac
274                     ENDIF
275                  ENDIF
276               END DO
277            END DO
278         ENDIF
279      ENDIF
280     
281      ! Initialise balancing increments
282      phyto2d_balinc(:,:,:,:) = 0.0
283     
284      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
285
286         ! Set up model parameters to be passed into Hemmings balancing routine.
287         ! For now these are hardwired to the standard HadOCC parameter values
288         ! as this is what the scheme was developed for.
289         ! Obviously, HadOCC and ERSEM are rather different models, so this
290         ! isn't ideal, but there's not always direct analogues between the two
291         ! parameter sets, so it's the easiest way to get something running.
292         ! In the longer term, some serious MarMOT-based development is required.
293         modparm(1)  = 0.1                               ! grow_sat
294         modparm(2)  = 2.0                               ! psmax
295         modparm(3)  = 0.845                             ! par
296         modparm(4)  = 0.02                              ! alpha
297         modparm(5)  = 0.05                              ! resp_rate
298         modparm(6)  = 0.05                              ! pmort_rate
299         modparm(7)  = 0.01                              ! phyto_min
300         modparm(8)  = 0.05                              ! z_mort_1
301         modparm(9)  = 1.0                               ! z_mort_2
302         modparm(10) = 6.625                             ! c2n_p
303         modparm(11) = 5.625                             ! c2n_z
304         modparm(12) = 7.5                               ! c2n_d
305         modparm(13) = 0.01                              ! graze_threshold
306         modparm(14) = 2.0                               ! holling_coef
307         modparm(15) = 0.5                               ! graze_sat
308         modparm(16) = 2.0                               ! graze_max
309
310         ! Set up assimilation parameters to be passed into balancing routine
311         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
312         assimparm(2)  = balnutext
313         assimparm(3)  = balnutmin
314         assimparm(4)  = r
315         assimparm(5)  = beta_g
316         assimparm(6)  = beta_l
317         assimparm(7)  = beta_m
318         assimparm(8)  = a_g
319         assimparm(9)  = a_l
320         assimparm(10) = a_m
321         assimparm(11) = zfracb0
322         assimparm(12) = zfracb1
323         assimparm(13) = qrfmax
324         assimparm(14) = qafmax
325         assimparm(15) = zrfmax
326         assimparm(16) = zafmax
327         assimparm(17) = prfmax
328         assimparm(18) = incphymin
329         assimparm(19) = integnstep
330         assimparm(20) = pthreshold
331
332         ! Set up external tracer indices array bstate
333         i_tracer(1) = 1   ! nutrient
334         i_tracer(2) = 2   ! phytoplankton
335         i_tracer(3) = 3   ! zooplankton
336         i_tracer(4) = 4   ! detritus
337         i_tracer(5) = 5   ! DIC
338         i_tracer(6) = 6   ! Alkalinity
339
340         ! Set background state
341         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_n3n) + &
342            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_n4n)
343         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p1n) + &
344            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p2n) + &
345            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p3n) + &
346            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p4n)
347         bstate(:,:,:,i_tracer(3)) = (tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z4c) * z4qnc) + &
348            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z5n) + &
349            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z6n)
350         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_r4n) + &
351            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_r6n) + &
352            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_r8n)
353         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_o3c)
354         bstate(:,:,:,i_tracer(6)) = totalk_bkg(:,:,:)
355
356         ! Calculate carbon to chlorophyll ratio for combined phytoplankton
357         ! and nitrogen to biomass equivalent for PZD (hardwire as per HadOCC)
358         cchl_p(:,:) = 0.0
359         DO jj = 1, jpj
360            DO ji = 1, jpi
361               IF ( ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
362                  &   tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) ) .GT. 0.0 ) THEN
363                  cchl_p(ji,jj) = zmassc * ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p1c) +      &
364                     &                       tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) +      &
365                     &                       tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p3c) +      &
366                     &                       tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p4c)   ) /  &
367                     &            ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) +             &
368                     &              tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) +             &
369                     &              tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) +             &
370                     &              tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)   )
371               ENDIF
372            END DO
373         END DO
374         n2be_p = zmassn + ( zmassc * modparm(10) )
375         n2be_z = zmassn + ( zmassc * modparm(11) )
376         n2be_d = zmassn + ( zmassc * modparm(12) )
377
378         ! Call nitrogen balancing routine
379         CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   &
380            &               n2be_p, n2be_z, n2be_d, assimparm,                      &
381            &               INT(pincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       &
382            &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot(:,:), cchl_p(:,:), &
383            &               nbal_active, phyt_avg_bkg(:,:),                         &
384            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      &
385            &               subsurf_active, deepneg_active,                         &
386            &               deeppos_active, nutprof_active,                         &
387            &               bstate, outincs,                                        &
388            &               diag_active, diag,                                      &
389            &               diag_fulldepth_active, diag_fulldepth )
390         
391         ! Loop over each grid point partioning the increments
392         DO jk = 1, jpk
393            DO jj = 1, jpj
394               DO ji = 1, jpi
395
396                  ! Phytoplankton
397                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) > 0.0 ) .AND. &
398                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) > 0.0 ) .AND. &
399                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) > 0.0 ) .AND. &
400                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) > 0.0 ) .AND. &
401                     & ( pinc_chltot(ji,jj) /= 0.0 ) ) THEN
402                     IF ( ld_chltot ) THEN
403                        ! Phytoplankton nitrogen split up based on existing ratios
404                        zfrac_p1n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) / &
405                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
406                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
407                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
408                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) )
409                        zfrac_p2n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) / &
410                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
411                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
412                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
413                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) )
414                        zfrac_p3n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) / &
415                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
416                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
417                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
418                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) )
419                        zfrac_p4n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) / &
420                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
421                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
422                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
423                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) )
424                     ELSE
425                        ! Phytoplankton nitrogen split up based on assimilation increments
426                        zfrac_p1n = pinc_chldia(ji,jj) / pinc_chltot(ji,jj)
427                        zfrac_p2n = pinc_chlnan(ji,jj) / pinc_chltot(ji,jj)
428                        zfrac_p3n = pinc_chlpic(ji,jj) / pinc_chltot(ji,jj)
429                        zfrac_p4n = pinc_chldin(ji,jj) / pinc_chltot(ji,jj)
430                     ENDIF
431                     
432                     ! Other phytoplankton variables split up based on existing ratios with nitrogen
433                     zrat_chl1_p1n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
434                     zrat_p1c_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c)  / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
435                     zrat_p1p_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p)  / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
436                     zrat_p1s_p1n  = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s)  / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
437                     zrat_chl2_p2n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)
438                     zrat_p2c_p2n  = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c)  / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)
439                     zrat_p2p_p2n  = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p)  / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)
440                     zrat_chl3_p3n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)
441                     zrat_p3c_p3n  = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c)  / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)
442                     zrat_p3p_p3n  = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p)  / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)
443                     zrat_chl4_p4n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)
444                     zrat_p4c_p4n  = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c)  / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)
445                     zrat_p4p_p4n  = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p)  / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)
446                     
447                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p1n
448                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p2n
449                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p3n
450                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p4n
451                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_chl1_p1n
452                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1c_p1n
453                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1p_p1n
454                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s)  = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1s_p1n
455                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) * zrat_chl2_p2n
456                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) * zrat_p2c_p2n
457                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) * zrat_p2p_p2n
458                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) * zrat_chl3_p3n
459                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) * zrat_p3c_p3n
460                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) * zrat_p3p_p3n
461                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) * zrat_chl4_p4n
462                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c)  = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) * zrat_p4c_p4n
463                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p)  = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) * zrat_p4p_p4n
464                  ENDIF
465
466                  ! Zooplankton nitrogen split up based on existing ratios
467                  ! Update carbon and phosphorus according to existing ratios
468                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) > 0.0 ) .AND. &
469                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) > 0.0 ) .AND. &
470                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) > 0.0 ) ) THEN
471                     zfrac_z4n = ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * z4qnc ) / &
472                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * z4qnc ) + &
473                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
474                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) )
475                     zfrac_z5n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) / &
476                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * z4qnc ) + &
477                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
478                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) )
479                     zfrac_z6n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) / &
480                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * z4qnc ) + &
481                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
482                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) )
483                     zrat_z4c_z4n = 1.0 / z4qnc
484                     zrat_z5c_z5n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5c) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n)
485                     zrat_z5p_z5n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n)
486                     zrat_z6c_z6n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6c) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n)
487                     zrat_z6p_z6n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n)
488                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z5n
489                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z6n
490                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z4n * zrat_z4c_z4n
491                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5c) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) * zrat_z5c_z5n
492                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6c) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) * zrat_z6c_z6n
493                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) * zrat_z5p_z5n
494                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) * zrat_z6p_z6n
495                  ENDIF
496
497                  ! Nitrogen nutrient split between nitrate and ammonium based on existing ratios
498                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) > 0.0 ) .AND. &
499                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) > 0.0 ) ) THEN
500                     zfrac_n3n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) / &
501                        &        (tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) + tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n))
502                     zfrac_n4n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) / &
503                        &        (tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) + tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n))
504                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n3n
505                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n4n
506                  ENDIF
507
508                  ! Detritus nitrogen split up based on existing ratios
509                  ! Update carbon and phosphorus according to existing ratios
510                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) > 0.0 ) .AND. &
511                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) > 0.0 ) .AND. &
512                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) > 0.0 ) ) THEN
513                     zfrac_r4n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) / &
514                        &        (tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
515                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
516                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
517                     zfrac_r6n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) / &
518                        &        (tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
519                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
520                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
521                     zfrac_r8n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) / &
522                        &        (tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
523                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
524                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
525                     zrat_r4c_r4n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4c) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n)
526                     zrat_r4p_r4n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4p) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n)
527                     zrat_r6c_r6n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6c) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n)
528                     zrat_r6p_r6n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6p) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n)
529                     zrat_r6s_r6n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6s) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n)
530                     zrat_r8c_r8n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8c) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n)
531                     zrat_r8p_r8n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8p) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n)
532                     zrat_r8s_r8n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8s) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n)
533                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r4n
534                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r6n
535                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r8n
536                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4c) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) * zrat_r4c_r4n
537                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4p) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) * zrat_r4p_r4n
538                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6c) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) * zrat_r6c_r6n
539                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6p) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) * zrat_r6p_r6n
540                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6s) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) * zrat_r6s_r6n
541                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8c) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) * zrat_r8c_r8n
542                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8p) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) * zrat_r8p_r8n
543                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8s) = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) * zrat_r8s_r8n
544                  ENDIF
545
546                  ! DIC straight from balancing scheme
547                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3c) = outincs(ji,jj,jk,i_tracer(5))
548
549                  ! Alkalinity straight from balancing scheme
550                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3ba) = outincs(ji,jj,jk,i_tracer(6))
551
552                  ! Remove P/R silicon increments from silicate to conserve mass
553                  zfrac = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) + &
554                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6s) + &
555                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8s)
556                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n5s) - zfrac ) > 0.0 ) THEN
557                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n5s) = zfrac * (-1.0)
558                  ENDIF
559
560                  ! Remove P/Z/R phosphorus increments from phosphate to conserve mass
561                  zfrac = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) + &
562                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) + &
563                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) + &
564                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) + &
565                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) + &
566                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) + &
567                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4p) + &
568                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6p) + &
569                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8p)
570                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n1p) - zfrac ) > 0.0 ) THEN
571                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n1p) = zfrac * (-1.0)
572                  ENDIF
573                 
574               END DO
575            END DO
576         END DO
577     
578      ELSE   ! No nitrogen balancing - just update phytoplankton
579         
580         ! Split up total surface chlorophyll increments
581         DO jj = 1, jpj
582            DO ji = 1, jpi
583               IF ( ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
584                  & ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
585                  & ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
586                  & ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) ) THEN
587                  IF ( ld_chltot ) THEN
588                     ! Chlorophyll split up based on existing ratios
589                     zfrac_chl1 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) /  &
590                        &        ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
591                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
592                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
593                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) )
594                     zfrac_chl2 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) /  &
595                        &        ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
596                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
597                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
598                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) )
599                     zfrac_chl3 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) /  &
600                        &        ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
601                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
602                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
603                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) )
604                     zfrac_chl4 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) /  &
605                        &        ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
606                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
607                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
608                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) )
609                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = pinc_chltot(ji,jj) * zfrac_chl1
610                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = pinc_chltot(ji,jj) * zfrac_chl2
611                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = pinc_chltot(ji,jj) * zfrac_chl3
612                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = pinc_chltot(ji,jj) * zfrac_chl4
613                  ENDIF
614                  IF( ld_chldia ) THEN
615                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = pinc_chldia(ji,jj)
616                  ENDIF
617                  IF( ld_chlnan ) THEN
618                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = pinc_chlnan(ji,jj)
619                  ENDIF
620                  IF( ld_chlpic ) THEN
621                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = pinc_chlpic(ji,jj)
622                  ENDIF
623                  IF( ld_chldin ) THEN
624                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = pinc_chldin(ji,jj)
625                  ENDIF
626                 
627                  ! Maintain stoichiometric ratios of carbon, nitrogen, phosphorus and silicon
628                  IF ( ld_chltot .OR. ld_chldia ) THEN
629                     zrat_p1c_chl1 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p1c) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
630                     zrat_p1n_chl1 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p1n) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
631                     zrat_p1p_chl1 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p1p) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
632                     zrat_p1s_chl1 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p1s) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
633                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1c) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) * zrat_p1c_chl1
634                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1n) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) * zrat_p1n_chl1
635                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1p) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) * zrat_p1p_chl1
636                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1s) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) * zrat_p1s_chl1
637                  ENDIF
638                  IF ( ld_chltot .OR. ld_chlnan ) THEN
639                     zrat_p2c_chl2 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
640                     zrat_p2n_chl2 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p2n) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
641                     zrat_p2p_chl2 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p2p) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
642                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) * zrat_p2c_chl2
643                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p2n) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) * zrat_p2n_chl2
644                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p2p) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) * zrat_p2p_chl2
645                  ENDIF
646                  IF ( ld_chltot .OR. ld_chlpic ) THEN
647                     zrat_p3c_chl3 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p3c) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
648                     zrat_p3n_chl3 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p3n) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
649                     zrat_p3p_chl3 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p3p) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
650                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p3c) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) * zrat_p3c_chl3
651                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p3n) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) * zrat_p3n_chl3
652                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p3p) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) * zrat_p3p_chl3
653                  ENDIF
654                  IF ( ld_chltot .OR. ld_chldin ) THEN
655                     zrat_p4c_chl4 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p4c) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
656                     zrat_p4n_chl4 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p4n) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
657                     zrat_p4p_chl4 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p4p) / tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
658                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p4c) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) * zrat_p4c_chl4
659                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p4n) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) * zrat_p4n_chl4
660                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p4p) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) * zrat_p4p_chl4
661                  ENDIF
662               ENDIF
663            END DO
664         END DO
665         
666         ! Propagate through mixed layer
667         DO jj = 1, jpj
668            DO ji = 1, jpi
669               !
670               jkmax = jpk-1
671               DO jk = jpk-1, 1, -1
672                  IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
673                     & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
674                     pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
675                     jkmax = jk
676                  ENDIF
677               END DO
678               !
679               DO jk = 2, jkmax
680                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
681                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1c)
682                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1n)
683                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1p)
684                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1s)
685                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
686                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p2c)
687                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p2n)
688                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p2p)
689                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
690                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p3c)
691                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p3n)
692                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p3p)
693                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
694                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p4c)
695                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p4n)
696                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p4p)
697               END DO
698               !
699            END DO
700         END DO
701
702      ENDIF
703
704   END SUBROUTINE asm_phyto2d_bal_ersem
705
706#else
707   !!----------------------------------------------------------------------
708   !!   Default option : Empty routine
709   !!----------------------------------------------------------------------
710CONTAINS
711   SUBROUTINE asm_phyto2d_bal_ersem(  ld_chltot,                      &
712      &                               pinc_chltot,                    &
713      &                               ld_chldia,                      &
714      &                               pinc_chldia,                    &
715      &                               ld_chlnan,                      &
716      &                               pinc_chlnan,                    &
717      &                               ld_chlpic,                      &
718      &                               pinc_chlpic,                    &
719      &                               ld_chldin,                      &
720      &                               pinc_chldin,                    &
721      &                               pincper,                        &
722      &                               p_maxchlinc, ld_phytobal, pmld, &
723      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
724      &                               phyt_avg_bkg, mld_max_bkg,      &
725      &                               totalk_bkg,                     &
726      &                               tracer_bkg, phyto2d_balinc )
727      LOGICAL :: ld_chltot
728      REAL    :: pinc_chltot(:,:)
729      LOGICAL :: ld_chldia
730      REAL    :: pinc_chldia(:,:)
731      LOGICAL :: ld_chlnan
732      REAL    :: pinc_chlnan(:,:)
733      LOGICAL :: ld_chlpic
734      REAL    :: pinc_chlpic(:,:)
735      LOGICAL :: ld_chldin
736      REAL    :: pinc_chldin(:,:)
737      REAL    :: pincper
738      REAL    :: p_maxchlinc
739      LOGICAL :: ld_phytobal
740      REAL    :: pmld(:,:)
741      REAL    :: pgrow_avg_bkg(:,:)
742      REAL    :: ploss_avg_bkg(:,:)
743      REAL    :: phyt_avg_bkg(:,:)
744      REAL    :: mld_max_bkg(:,:)
745      REAL    :: totalk_bkg(:,:,:)
746      REAL    :: tracer_bkg(:,:,:,:)
747      REAL    :: phyto2d_balinc(:,:,:,:)
748      WRITE(*,*) 'asm_phyto2d_bal_ersem: You should not have seen this print! error?'
749   END SUBROUTINE asm_phyto2d_bal_ersem
750#endif
751
752   !!======================================================================
753END MODULE asmphyto2dbal_ersem
Note: See TracBrowser for help on using the repository browser.