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 @ 10622

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

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