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/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_ersem.F90 @ 10728

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

Changes for assimilation of biogeochemical variables. See ​https://code.metoffice.gov.uk/trac/utils/ticket/174.
File size: 55.7 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 ! Temporary array
197      INTEGER,                 DIMENSION(1)                 :: zkmt             ! No. sea points in column
198      !!---------------------------------------------------------------------------
199     
200      ! Set parameters, firstly molar mass of carbon and nitrogen
201      zmassc = 12.01
202      zmassn = 14.01
203     
204      ! Then mesozooplankton nitrogen(mmol):carbon(mg), ERSEM parameter Z4/qnc
205      ! Hardwire it for now due to difficulty of getting at parameters through FABM
206      ! I think the following should work:
207      ! z4qnc = model%state_variables(jp_fabm_z4c)%properties%get_property_by_name('qnc')%value
208      ! but, get_property_by_name is in the module fabm_properties, which we can't get at directly
209      ! We can do a "USE fabm" which itself does a "USE fabm_properties", but doing it indirectly
210      ! like that means get_property_by_name is treated as private so we can't use it
211      ! So hardwire to value used in SSB runs and v4/5 CMEMS reanalysis
212      z4qnc  = 0.0126
213
214      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
215      IF ( p_maxchlinc > 0.0 ) THEN
216         IF ( ld_chltot ) THEN
217            DO jj = 1, jpj
218               DO ji = 1, jpi
219                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
220               END DO
221            END DO
222         ELSE
223            DO jj = 1, jpj
224               DO ji = 1, jpi
225                  IF ( ld_chldia .AND. ld_chlnan .AND. ld_chlpic .AND. ld_chldin ) THEN
226                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
227                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
228                  ELSE IF ( ld_chldia .AND. ld_chlnan .AND. ld_chlpic ) THEN
229                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
230                        &                      pinc_chlpic(ji,jj)
231                  ELSE IF ( ld_chldia .AND. ld_chlnan .AND. ld_chldin ) THEN
232                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj) + &
233                        &                      pinc_chldin(ji,jj)
234                  ELSE IF ( ld_chldia .AND. ld_chlpic .AND. ld_chldin ) THEN
235                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + &
236                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
237                  ELSE IF ( ld_chlnan .AND. ld_chlpic .AND. ld_chldin ) THEN
238                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + &
239                        &                      pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
240                  ELSE IF ( ld_chldia .AND. ld_chlnan ) THEN
241                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlnan(ji,jj)
242                  ELSE IF ( ld_chldia .AND. ld_chlpic ) THEN
243                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chlpic(ji,jj)
244                  ELSE IF ( ld_chldia .AND. ld_chldin ) THEN
245                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj) + pinc_chldin(ji,jj)
246                  ELSE IF ( ld_chlnan .AND. ld_chlpic ) THEN
247                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + pinc_chlpic(ji,jj)
248                  ELSE IF ( ld_chlnan .AND. ld_chldin ) THEN
249                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj) + pinc_chldin(ji,jj)
250                  ELSE IF ( ld_chlpic .AND. ld_chldin ) THEN
251                     pinc_chltot_temp(ji,jj) = pinc_chlpic(ji,jj) + pinc_chldin(ji,jj)
252                  ELSE IF ( ld_chldia ) THEN
253                     pinc_chltot_temp(ji,jj) = pinc_chldia(ji,jj)
254                  ELSE IF ( ld_chlnan ) THEN
255                     pinc_chltot_temp(ji,jj) = pinc_chlnan(ji,jj)
256                  ELSE IF ( ld_chlpic ) THEN
257                     pinc_chltot_temp(ji,jj) = pinc_chlpic(ji,jj)
258                  ELSE IF ( ld_chldin ) THEN
259                     pinc_chltot_temp(ji,jj) = pinc_chldin(ji,jj)
260                  ENDIF
261                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot_temp(ji,jj), p_maxchlinc ) )
262                  IF ( pinc_chltot(ji,jj) .NE. pinc_chltot_temp(ji,jj) ) THEN
263                     zfrac = pinc_chltot(ji,jj) / pinc_chltot_temp(ji,jj)
264                     IF ( ld_chldia ) THEN
265                        pinc_chldia(ji,jj) = pinc_chldia(ji,jj) * zfrac
266                     ENDIF
267                     IF ( ld_chlnan ) THEN
268                        pinc_chlnan(ji,jj) = pinc_chlnan(ji,jj) * zfrac
269                     ENDIF
270                     IF ( ld_chlpic ) THEN
271                        pinc_chlpic(ji,jj) = pinc_chlpic(ji,jj) * zfrac
272                     ENDIF
273                     IF ( ld_chldin ) THEN
274                        pinc_chldin(ji,jj) = pinc_chldin(ji,jj) * zfrac
275                     ENDIF
276                  ENDIF
277               END DO
278            END DO
279         ENDIF
280      ENDIF
281     
282      ! Initialise balancing increments
283      phyto2d_balinc(:,:,:,:) = 0.0
284     
285      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
286
287         ! Set up model parameters to be passed into Hemmings balancing routine.
288         ! For now these are hardwired to the standard HadOCC parameter values
289         ! as this is what the scheme was developed for.
290         ! Obviously, HadOCC and ERSEM are rather different models, so this
291         ! isn't ideal, but there's not always direct analogues between the two
292         ! parameter sets, so it's the easiest way to get something running.
293         ! In the longer term, some serious MarMOT-based development is required.
294         modparm(1)  = 0.1                               ! grow_sat
295         modparm(2)  = 2.0                               ! psmax
296         modparm(3)  = 0.845                             ! par
297         modparm(4)  = 0.02                              ! alpha
298         modparm(5)  = 0.05                              ! resp_rate
299         modparm(6)  = 0.05                              ! pmort_rate
300         modparm(7)  = 0.01                              ! phyto_min
301         modparm(8)  = 0.05                              ! z_mort_1
302         modparm(9)  = 1.0                               ! z_mort_2
303         !modparm(10) = 6.625                            ! c2n_p - set later per grid point
304         !modparm(11) = 5.625                            ! c2n_z - set later per grid point
305         !modparm(12) = 7.5                              ! c2n_d - set later per grid point
306         modparm(13) = 0.01                              ! graze_threshold
307         modparm(14) = 2.0                               ! holling_coef
308         modparm(15) = 0.5                               ! graze_sat
309         modparm(16) = 2.0                               ! graze_max
310
311         ! Set up assimilation parameters to be passed into balancing routine
312         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
313         assimparm(2)  = balnutext
314         assimparm(3)  = balnutmin
315         assimparm(4)  = r
316         assimparm(5)  = beta_g
317         assimparm(6)  = beta_l
318         assimparm(7)  = beta_m
319         assimparm(8)  = a_g
320         assimparm(9)  = a_l
321         assimparm(10) = a_m
322         assimparm(11) = zfracb0
323         assimparm(12) = zfracb1
324         assimparm(13) = qrfmax
325         assimparm(14) = qafmax
326         assimparm(15) = zrfmax
327         assimparm(16) = zafmax
328         assimparm(17) = prfmax
329         assimparm(18) = incphymin
330         assimparm(19) = integnstep
331         assimparm(20) = pthreshold
332
333         ! Set up external tracer indices array bstate
334         i_tracer(1) = 1   ! nutrient
335         i_tracer(2) = 2   ! phytoplankton
336         i_tracer(3) = 3   ! zooplankton
337         i_tracer(4) = 4   ! detritus
338         i_tracer(5) = 5   ! DIC
339         i_tracer(6) = 6   ! Alkalinity
340
341         ! Set background state
342         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_n3n) + &
343            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_n4n)
344         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p1n) + &
345            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p2n) + &
346            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p3n) + &
347            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p4n)
348         bstate(:,:,:,i_tracer(3)) = (tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z4c) * z4qnc) + &
349            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z5n) + &
350            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z6n)
351         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_r4n) + &
352            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_r6n) + &
353            &                        tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_r8n)
354         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_o3c)
355         bstate(:,:,:,i_tracer(6)) = totalk_bkg(:,:,:)
356
357         ! Calculate carbon to chlorophyll ratio for combined phytoplankton
358         ! and nitrogen to biomass equivalent for PZD (hardwire as per HadOCC)
359         cchl_p(:,:) = 0.0
360         DO jj = 1, jpj
361            DO ji = 1, jpi
362               IF ( ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
363                  &   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
364                  cchl_p(ji,jj) = zmassc * ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p1c) +      &
365                     &                       tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) +      &
366                     &                       tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p3c) +      &
367                     &                       tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p4c)   ) /  &
368                     &            ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) +             &
369                     &              tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) +             &
370                     &              tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) +             &
371                     &              tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)   )
372               ENDIF
373            END DO
374         END DO
375
376         ! Call nitrogen balancing routine - loop over grid points due to variable C:N ratios
377         DO jj = 1, jpj
378            DO ji = 1, jpi
379               ! Phytoplankton C:N
380               modparm(10) = ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p1c) + &
381                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) + &
382                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p3c) + &
383                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p4c)     ) / &
384                  &          ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p1n) + &
385                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p2n) + &
386                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p3n) + &
387                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_p4n)     )
388               ! Zooplankton C:N
389               modparm(11) = ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_z4c) + &
390                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_z5c) + &
391                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_z6c)     ) / &
392                  &          ( (tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_z4c) * z4qnc) + &
393                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_z5n) + &
394                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_z6n)     )
395               ! Detritus C:N
396               modparm(12) = ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_r4c) + &
397                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_r6c) + &
398                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_r8c)     ) / &
399                  &          ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_r4n) + &
400                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_r6n) + &
401                  &            tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_r8n)     )
402               ! Nitrogen to biomass equivalent for PZD
403               n2be_p = zmassn + ( zmassc * modparm(10) )
404               n2be_z = zmassn + ( zmassc * modparm(11) )
405               n2be_d = zmassn + ( zmassc * modparm(12) )
406               zkmt(:) = INT(SUM(tmask(ji,jj,:)))
407               CALL bio_analysis( 1, 1, jpk, gdepw_n(ji,jj,2:jpk), i_tracer, modparm,                 &
408                  &               n2be_p, n2be_z, n2be_d, assimparm,                                  &
409                  &               INT(pincper), 1, zkmt, tmask(ji,jj,:),                              &
410                  &               pmld(ji,jj), mld_max_bkg(ji,jj), pinc_chltot(ji,jj), cchl_p(ji,jj), &
411                  &               nbal_active, phyt_avg_bkg(ji,jj),                                   &
412                  &               gl_active, pgrow_avg_bkg(ji,jj), ploss_avg_bkg(ji,jj),              &
413                  &               subsurf_active, deepneg_active,                                     &
414                  &               deeppos_active, nutprof_active,                                     &
415                  &               bstate(ji,jj,:,:), outincs(ji,jj,:,:),                              &
416                  &               diag_active, diag(ji,jj,:),                                         &
417                  &               diag_fulldepth_active, diag_fulldepth(ji,jj,:,:) )
418            END DO
419         END DO
420         
421         ! Loop over each grid point partioning the increments
422         DO jk = 1, jpk
423            DO jj = 1, jpj
424               DO ji = 1, jpi
425
426                  ! Phytoplankton
427                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) > 0.0 ) .AND. &
428                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) > 0.0 ) .AND. &
429                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) > 0.0 ) .AND. &
430                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) > 0.0 ) .AND. &
431                     & ( pinc_chltot(ji,jj) /= 0.0 ) ) THEN
432                     IF ( ld_chltot ) THEN
433                        ! Phytoplankton nitrogen split up based on existing ratios
434                        zfrac_p1n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) / &
435                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
436                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
437                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
438                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) )
439                        zfrac_p2n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) / &
440                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
441                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
442                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
443                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) )
444                        zfrac_p3n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) / &
445                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
446                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
447                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
448                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) )
449                        zfrac_p4n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) / &
450                           &        ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
451                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
452                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
453                           &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) )
454                     ELSE
455                        ! Phytoplankton nitrogen split up based on assimilation increments
456                        zfrac_p1n = pinc_chldia(ji,jj) / pinc_chltot(ji,jj)
457                        zfrac_p2n = pinc_chlnan(ji,jj) / pinc_chltot(ji,jj)
458                        zfrac_p3n = pinc_chlpic(ji,jj) / pinc_chltot(ji,jj)
459                        zfrac_p4n = pinc_chldin(ji,jj) / pinc_chltot(ji,jj)
460                     ENDIF
461                     
462                     ! Other phytoplankton variables split up based on existing ratios with nitrogen
463                     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)
464                     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)
465                     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)
466                     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)
467                     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)
468                     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)
469                     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)
470                     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)
471                     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)
472                     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)
473                     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)
474                     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)
475                     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)
476                     
477                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p1n
478                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p2n
479                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p3n
480                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)  = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p4n
481                     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
482                     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
483                     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
484                     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
485                     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
486                     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
487                     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
488                     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
489                     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
490                     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
491                     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
492                     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
493                     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
494                  ENDIF
495
496                  ! Zooplankton nitrogen split up based on existing ratios
497                  ! Update carbon and phosphorus according to existing ratios
498                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) > 0.0 ) .AND. &
499                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) > 0.0 ) .AND. &
500                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) > 0.0 ) ) THEN
501                     zfrac_z4n = ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * z4qnc ) / &
502                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * z4qnc ) + &
503                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
504                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) )
505                     zfrac_z5n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) / &
506                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * z4qnc ) + &
507                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
508                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) )
509                     zfrac_z6n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) / &
510                        &        ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * z4qnc ) + &
511                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
512                        &          tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) )
513                     zrat_z4c_z4n = 1.0 / z4qnc
514                     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)
515                     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)
516                     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)
517                     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)
518                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z5n
519                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z6n
520                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z4n * zrat_z4c_z4n
521                     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
522                     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
523                     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
524                     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
525                  ENDIF
526
527                  ! Nitrogen nutrient split between nitrate and ammonium based on existing ratios
528                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) > 0.0 ) .AND. &
529                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) > 0.0 ) ) THEN
530                     zfrac_n3n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) / &
531                        &        (tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) + tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n))
532                     zfrac_n4n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) / &
533                        &        (tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) + tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n))
534                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n3n
535                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n4n
536                  ENDIF
537
538                  ! Detritus nitrogen split up based on existing ratios
539                  ! Update carbon and phosphorus according to existing ratios
540                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) > 0.0 ) .AND. &
541                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) > 0.0 ) .AND. &
542                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) > 0.0 ) ) THEN
543                     zfrac_r4n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) / &
544                        &        (tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
545                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
546                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
547                     zfrac_r6n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) / &
548                        &        (tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
549                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
550                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
551                     zfrac_r8n = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) / &
552                        &        (tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
553                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
554                        &         tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
555                     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)
556                     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)
557                     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)
558                     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)
559                     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)
560                     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)
561                     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)
562                     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)
563                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r4n
564                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r6n
565                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r8n
566                     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
567                     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
568                     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
569                     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
570                     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
571                     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
572                     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
573                     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
574                  ENDIF
575
576                  ! DIC straight from balancing scheme
577                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3c) = outincs(ji,jj,jk,i_tracer(5))
578
579                  ! Alkalinity straight from balancing scheme
580                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3ba) = outincs(ji,jj,jk,i_tracer(6))
581
582                  ! Remove P/R silicon increments from silicate to conserve mass
583                  zfrac = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) + &
584                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6s) + &
585                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8s)
586                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n5s) - zfrac ) > 0.0 ) THEN
587                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n5s) = zfrac * (-1.0)
588                  ENDIF
589
590                  ! Remove P/Z/R phosphorus increments from phosphate to conserve mass
591                  zfrac = phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) + &
592                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) + &
593                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) + &
594                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) + &
595                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) + &
596                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) + &
597                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4p) + &
598                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6p) + &
599                     &    phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8p)
600                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_n1p) - zfrac ) > 0.0 ) THEN
601                     phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n1p) = zfrac * (-1.0)
602                  ENDIF
603                 
604               END DO
605            END DO
606         END DO
607     
608      ELSE   ! No nitrogen balancing - just update phytoplankton
609         
610         ! Split up total surface chlorophyll increments
611         DO jj = 1, jpj
612            DO ji = 1, jpi
613               IF ( ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
614                  & ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
615                  & ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
616                  & ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) ) THEN
617                  IF ( ld_chltot ) THEN
618                     ! Chlorophyll split up based on existing ratios
619                     zfrac_chl1 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) /  &
620                        &        ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
621                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
622                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
623                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) )
624                     zfrac_chl2 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) /  &
625                        &        ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
626                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
627                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
628                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) )
629                     zfrac_chl3 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) /  &
630                        &        ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
631                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
632                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
633                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) )
634                     zfrac_chl4 = tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) /  &
635                        &        ( tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
636                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
637                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
638                        &          tracer_bkg(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) )
639                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = pinc_chltot(ji,jj) * zfrac_chl1
640                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = pinc_chltot(ji,jj) * zfrac_chl2
641                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = pinc_chltot(ji,jj) * zfrac_chl3
642                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = pinc_chltot(ji,jj) * zfrac_chl4
643                  ENDIF
644                  IF( ld_chldia ) THEN
645                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = pinc_chldia(ji,jj)
646                  ENDIF
647                  IF( ld_chlnan ) THEN
648                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = pinc_chlnan(ji,jj)
649                  ENDIF
650                  IF( ld_chlpic ) THEN
651                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = pinc_chlpic(ji,jj)
652                  ENDIF
653                  IF( ld_chldin ) THEN
654                     phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = pinc_chldin(ji,jj)
655                  ENDIF
656                 
657                  ! Maintain stoichiometric ratios of carbon, nitrogen, phosphorus and silicon
658                  IF ( ld_chltot .OR. ld_chldia ) THEN
659                     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)
660                     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)
661                     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)
662                     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)
663                     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
664                     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
665                     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
666                     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
667                  ENDIF
668                  IF ( ld_chltot .OR. ld_chlnan ) THEN
669                     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)
670                     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)
671                     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)
672                     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
673                     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
674                     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
675                  ENDIF
676                  IF ( ld_chltot .OR. ld_chlpic ) THEN
677                     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)
678                     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)
679                     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)
680                     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
681                     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
682                     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
683                  ENDIF
684                  IF ( ld_chltot .OR. ld_chldin ) THEN
685                     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)
686                     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)
687                     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)
688                     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
689                     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
690                     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
691                  ENDIF
692               ENDIF
693            END DO
694         END DO
695         
696         ! Propagate through mixed layer
697         DO jj = 1, jpj
698            DO ji = 1, jpi
699               !
700               jkmax = jpk-1
701               DO jk = jpk-1, 1, -1
702                  IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
703                     & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
704                     pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
705                     jkmax = jk
706                  ENDIF
707               END DO
708               !
709               DO jk = 2, jkmax
710                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
711                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1c)
712                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1n)
713                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1p)
714                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p1s)
715                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
716                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p2c)
717                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p2n)
718                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p2p)
719                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
720                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p3c)
721                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p3n)
722                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p3p)
723                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
724                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p4c)
725                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p4n)
726                  phyto2d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p)  = phyto2d_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_p4p)
727               END DO
728               !
729            END DO
730         END DO
731
732      ENDIF
733
734   END SUBROUTINE asm_phyto2d_bal_ersem
735
736#else
737   !!----------------------------------------------------------------------
738   !!   Default option : Empty routine
739   !!----------------------------------------------------------------------
740CONTAINS
741   SUBROUTINE asm_phyto2d_bal_ersem(  ld_chltot,                      &
742      &                               pinc_chltot,                    &
743      &                               ld_chldia,                      &
744      &                               pinc_chldia,                    &
745      &                               ld_chlnan,                      &
746      &                               pinc_chlnan,                    &
747      &                               ld_chlpic,                      &
748      &                               pinc_chlpic,                    &
749      &                               ld_chldin,                      &
750      &                               pinc_chldin,                    &
751      &                               pincper,                        &
752      &                               p_maxchlinc, ld_phytobal, pmld, &
753      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
754      &                               phyt_avg_bkg, mld_max_bkg,      &
755      &                               totalk_bkg,                     &
756      &                               tracer_bkg, phyto2d_balinc )
757      LOGICAL :: ld_chltot
758      REAL    :: pinc_chltot(:,:)
759      LOGICAL :: ld_chldia
760      REAL    :: pinc_chldia(:,:)
761      LOGICAL :: ld_chlnan
762      REAL    :: pinc_chlnan(:,:)
763      LOGICAL :: ld_chlpic
764      REAL    :: pinc_chlpic(:,:)
765      LOGICAL :: ld_chldin
766      REAL    :: pinc_chldin(:,:)
767      REAL    :: pincper
768      REAL    :: p_maxchlinc
769      LOGICAL :: ld_phytobal
770      REAL    :: pmld(:,:)
771      REAL    :: pgrow_avg_bkg(:,:)
772      REAL    :: ploss_avg_bkg(:,:)
773      REAL    :: phyt_avg_bkg(:,:)
774      REAL    :: mld_max_bkg(:,:)
775      REAL    :: totalk_bkg(:,:,:)
776      REAL    :: tracer_bkg(:,:,:,:)
777      REAL    :: phyto2d_balinc(:,:,:,:)
778      WRITE(*,*) 'asm_phyto2d_bal_ersem: You should not have seen this print! error?'
779   END SUBROUTINE asm_phyto2d_bal_ersem
780#endif
781
782   !!======================================================================
783END MODULE asmphyto2dbal_ersem
Note: See TracBrowser for help on using the repository browser.