source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/ASM/asmphytobal_medusa.F90 @ 13316

Last change on this file since 13316 was 13316, checked in by dford, 2 months ago

Allow nitrogen balancing when assimilating 3D chlorophyll data. See Met Office utils ticket 346.

File size: 33.2 KB
Line 
1MODULE asmphytobal_medusa
2   !!======================================================================
3   !!                       ***  MODULE asmphyto2dbal_medusa  ***
4   !! Calculate increments to MEDUSA 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  ! 2017-08 (D. Ford)  Adapted from asmphyto2dbal_hadocc
13   !!----------------------------------------------------------------------
14#if defined key_asminc && defined key_medusa
15   !!----------------------------------------------------------------------
16   !! 'key_asminc'          : assimilation increment interface
17   !! 'key_medusa'          : MEDUSA model
18   !!----------------------------------------------------------------------
19   !! asm_phyto2d_bal_medusa : routine to calculate increments to MEDUSA
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 zdftmx,        ONLY: ln_tmx_itf, &  ! Indonesian Throughflow
25      &                     mask_itf       ! tidal mixing mask
26   USE iom                                 ! i/o
27   USE sms_medusa                          ! MEDUSA parameters
28   USE par_medusa                          ! MEDUSA parameters
29   USE par_trc,       ONLY: jptra          ! Tracer parameters
30   USE bioanalysis                         ! Nitrogen balancing
31
32   IMPLICIT NONE
33   PRIVATE                   
34
35   PUBLIC asm_phyto_bal_medusa
36
37   ! Default values for biological assimilation parameters
38   ! Should match Hemmings et al. (2008)
39   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
40   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
41   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
42   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
43   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
44   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
45   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
46   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
47   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
48   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
49   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
50   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
51   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
52   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
53   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
54   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
55   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
56   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
57   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
58   !
59   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
60   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
61   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
62   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
63   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
64   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
65   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
66   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
67
68CONTAINS
69
70   SUBROUTINE asm_phyto_bal_medusa( kdeps,                              &
71      &                             ld_chltot,                          &
72      &                             pinc_chltot_3d,                     &
73      &                             ld_chldia,                          &
74      &                             pinc_chldia_3d,                     &
75      &                             ld_chlnon,                          &
76      &                             pinc_chlnon_3d,                     &
77      &                             ld_phytot,                          &
78      &                             pinc_phytot_3d,                     &
79      &                             ld_phydia,                          &
80      &                             pinc_phydia_3d,                     &
81      &                             ld_phynon,                          &
82      &                             pinc_phynon_3d,                     &
83      &                             pincper,                            &
84      &                             p_maxchlinc, ld_phytobal, pmld,     &
85      &                             pgrow_avg_bkg_3d, ploss_avg_bkg_3d, &
86      &                             phyt_avg_bkg_3d, mld_max_bkg,       &
87      &                             tracer_bkg, phyto_balinc )
88      !!---------------------------------------------------------------------------
89      !!                    ***  ROUTINE asm_phyto_bal_medusa  ***
90      !!
91      !! ** Purpose :   calculate increments to MEDUSA from phytoplankton increments
92      !!
93      !! ** Method  :   average up MEDUSA to look like HadOCC
94      !!                call nitrogen balancing scheme
95      !!                separate back out to MEDUSA
96      !!
97      !! ** Action  :   populate phyto_balinc
98      !!
99      !! References :   Hemmings et al., 2008, J. Mar. Res.
100      !!                Ford et al., 2012, Ocean Sci.
101      !!---------------------------------------------------------------------------
102      !!
103      INTEGER,  INTENT(in   )                               :: kdeps            ! No. inc deps 1 or jpk
104      LOGICAL,  INTENT(in   )                               :: ld_chltot        ! Assim chltot y/n
105      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_chltot_3d   ! chltot increments (3D)
106      LOGICAL,  INTENT(in   )                               :: ld_chldia        ! Assim chldia y/n
107      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_chldia_3d   ! chldia increments (3D)
108      LOGICAL,  INTENT(in   )                               :: ld_chlnon        ! Assim chlnon y/n
109      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_chlnon_3d   ! chlnon increments (3D)
110      LOGICAL,  INTENT(in   )                               :: ld_phytot        ! Assim phytot y/n
111      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_phytot_3d   ! phytot increments (3D)
112      LOGICAL,  INTENT(in   )                               :: ld_phydia        ! Assim phydia y/n
113      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_phydia_3d   ! phydia increments (3D)
114      LOGICAL,  INTENT(in   )                               :: ld_phynon        ! Assim phynon y/n
115      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_phynon_3d   ! phynon increments (3D)
116      REAL(wp), INTENT(in   )                               :: pincper          ! Assimilation period
117      REAL(wp), INTENT(in   )                               :: p_maxchlinc      ! Max chl increment
118      LOGICAL,  INTENT(in   )                               :: ld_phytobal      ! Balancing y/n
119      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld             ! Mixed layer depth
120      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: pgrow_avg_bkg_3d ! Avg phyto growth (3D)
121      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: ploss_avg_bkg_3d ! Avg phyto loss (3D)
122      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: phyt_avg_bkg_3d  ! Avg phyto (3D)
123      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg      ! Max MLD
124      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg       ! State variables
125      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto_balinc     ! Balancing increments
126      !!
127      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
128      INTEGER                                               :: jkmax          ! Loop index
129      INTEGER                                               :: jkinc          ! Loop index
130      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
131      REAL(wp)                                              :: n2be_p         ! N:biomass for total phy
132      REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo
133      REAL(wp)                                              :: n2be_d         ! N:biomass for detritus
134      REAL(wp)                                              :: zfrac          ! Fraction
135      REAL(wp)                                              :: zfrac_chn      ! Fraction of jpchn
136      REAL(wp)                                              :: zfrac_chd      ! Fraction of jpchd
137      REAL(wp)                                              :: zfrac_phn      ! Fraction of jpphn
138      REAL(wp)                                              :: zfrac_phd      ! Fraction of jpphd
139      REAL(wp)                                              :: zfrac_zmi      ! Fraction of jpzmi
140      REAL(wp)                                              :: zfrac_zme      ! Fraction of jpzme
141      REAL(wp)                                              :: zrat_pds_phd   ! Ratio of jppds:jpphd
142      REAL(wp)                                              :: zrat_chd_phd   ! Ratio of jpchd:jpphd
143      REAL(wp)                                              :: zrat_chn_phn   ! Ratio of jpchn:jpphn
144      REAL(wp)                                              :: zrat_phn_chn   ! Ratio of jpphn:jpchn
145      REAL(wp)                                              :: zrat_phd_chd   ! Ratio of jpphd:jpchd
146      REAL(wp)                                              :: zrat_pds_chd   ! Ratio of jppds:jpchd
147      REAL(wp)                                              :: zrat_dtc_det   ! Ratio of jpdtc:jpdet
148      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p_2d      ! C:Chl for total phy (2D)
149      REAL(wp),                DIMENSION(jpi,jpj,jpk)       :: cchl_p_3d      ! C:Chl for total phy (3D)
150      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
151      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters
152      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: bstate_2d      ! Background state (2D)
153      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate_3d      ! Background state (3D)
154      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: outincs_2d     ! Balancing increments (2D)
155      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs_3d     ! Balancing increments (3D)
156      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag              ! Depth-indep diagnostics
157      REAL(wp),                DIMENSION(jpi,jpj,1,22)      :: diag_fulldepth_2d ! Full-depth diagnostics (2D)
158      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth_3d ! Full-depth diagnostics (3D)
159      REAL(wp),                DIMENSION(jpi,jpj,1)         :: tmask_2d          ! Single-level tmask
160      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chltot_2d    ! chltot increments (2D)
161      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chldia_2d    ! chldia increments (2D)
162      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chlnon_2d    ! chlnon increments (2D)
163      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_phytot_2d    ! phytot increments (2D)
164      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_phydia_2d    ! phydia increments (2D)
165      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_phynon_2d    ! phynon increments (2D)
166      REAL(wp),                DIMENSION(jpi,jpj)           :: pgrow_avg_bkg_2d  ! Avg phyto growth (2D)
167      REAL(wp),                DIMENSION(jpi,jpj)           :: ploss_avg_bkg_2d  ! Avg phyto loss (2D)
168      REAL(wp),                DIMENSION(jpi,jpj)           :: phyt_avg_bkg_2d   ! Avg phyto (2D)
169      !!---------------------------------------------------------------------------
170
171      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
172      IF ( p_maxchlinc > 0.0 ) THEN
173         DO jk = 1, kdeps
174            DO jj = 1, jpj
175               DO ji = 1, jpi
176                  IF ( ld_chltot ) THEN
177                     pinc_chltot_3d(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot_3d(ji,jj,jk), p_maxchlinc ) )
178                  ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN
179                     pinc_chltot_3d(ji,jj,jk) = pinc_chldia_3d(ji,jj,jk) + pinc_chlnon_3d(ji,jj,jk)
180                     pinc_chltot_3d(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot_3d(ji,jj,jk), p_maxchlinc ) )
181                     IF ( pinc_chltot_3d(ji,jj,jk) .NE. ( pinc_chldia_3d(ji,jj,jk) + pinc_chlnon_3d(ji,jj,jk) ) ) THEN
182                        zfrac = pinc_chltot_3d(ji,jj,jk) / ( pinc_chldia_3d(ji,jj,jk) + pinc_chlnon_3d(ji,jj,jk) )
183                        pinc_chldia_3d(ji,jj,jk) = pinc_chldia_3d(ji,jj,jk) * zfrac
184                        pinc_chlnon_3d(ji,jj,jk) = pinc_chlnon_3d(ji,jj,jk) * zfrac
185                     ENDIF
186                  ELSE IF ( ld_chldia ) THEN
187                     pinc_chldia_3d(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chldia_3d(ji,jj,jk), p_maxchlinc ) )
188                     pinc_chltot_3d(ji,jj,jk) = pinc_chldia_3d(ji,jj,jk)
189                  ELSE IF ( ld_chlnon ) THEN
190                     pinc_chlnon_3d(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chlnon_3d(ji,jj,jk), p_maxchlinc ) )
191                     pinc_chltot_3d(ji,jj,jk) = pinc_chlnon_3d(ji,jj,jk)
192                  ENDIF
193               END DO
194            END DO
195         END DO
196      ENDIF
197
198      IF ( ld_phytot .OR. ld_phydia .OR. ld_phynon ) THEN
199         CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' )
200      ENDIF
201     
202      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
203
204         ! Set up model parameters to be passed into Hemmings balancing routine.
205         ! For now these are hardwired to the standard HadOCC parameter values
206         ! (except C:N ratios) as this is what the scheme was developed for.
207         ! Obviously, HadOCC and MEDUSA are rather different models, so this
208         ! isn't ideal, but there's not always direct analogues between the two
209         ! parameter sets, so it's the easiest way to get something running.
210         ! In the longer term, some serious MarMOT-based development is required.
211         modparm(1)  = 0.1                               ! grow_sat
212         modparm(2)  = 2.0                               ! psmax
213         modparm(3)  = 0.845                             ! par
214         modparm(4)  = 0.02                              ! alpha
215         modparm(5)  = 0.05                              ! resp_rate
216         modparm(6)  = 0.05                              ! pmort_rate
217         modparm(7)  = 0.01                              ! phyto_min
218         modparm(8)  = 0.05                              ! z_mort_1
219         modparm(9)  = 1.0                               ! z_mort_2
220         modparm(10) = ( xthetapn  + xthetapd  ) / 2.0   ! c2n_p
221         modparm(11) = ( xthetazmi + xthetazme ) / 2.0   ! c2n_z
222         modparm(12) = xthetad                           ! c2n_d
223         modparm(13) = 0.01                              ! graze_threshold
224         modparm(14) = 2.0                               ! holling_coef
225         modparm(15) = 0.5                               ! graze_sat
226         modparm(16) = 2.0                               ! graze_max
227
228         ! Set up assimilation parameters to be passed into balancing routine
229         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
230         assimparm(2)  = balnutext
231         assimparm(3)  = balnutmin
232         assimparm(4)  = r
233         assimparm(5)  = beta_g
234         assimparm(6)  = beta_l
235         assimparm(7)  = beta_m
236         assimparm(8)  = a_g
237         assimparm(9)  = a_l
238         assimparm(10) = a_m
239         assimparm(11) = zfracb0
240         assimparm(12) = zfracb1
241         assimparm(13) = qrfmax
242         assimparm(14) = qafmax
243         assimparm(15) = zrfmax
244         assimparm(16) = zafmax
245         assimparm(17) = prfmax
246         assimparm(18) = incphymin
247         assimparm(19) = integnstep
248         assimparm(20) = pthreshold
249
250         ! Set up external tracer indices array bstate
251         i_tracer(1) = 1   ! nutrient
252         i_tracer(2) = 2   ! phytoplankton
253         i_tracer(3) = 3   ! zooplankton
254         i_tracer(4) = 4   ! detritus
255         i_tracer(5) = 5   ! DIC
256         i_tracer(6) = 6   ! Alkalinity
257
258         ! Set background state
259         bstate_3d(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jpdin)
260         bstate_3d(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jpphn) + tracer_bkg(:,:,:,jpphd)
261         bstate_3d(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jpzmi) + tracer_bkg(:,:,:,jpzme)
262         bstate_3d(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jpdet)
263         bstate_3d(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jpdic)
264         bstate_3d(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jpalk)
265
266         ! Calculate carbon to chlorophyll ratio for combined phytoplankton
267         ! and nitrogen to biomass equivalent for PZD
268         ! Hardwire nitrogen mass to 14.01 for now as it doesn't seem to be set in MEDUSA
269         cchl_p_3d(:,:,:) = 0.0
270         DO jk = 1, jpk
271            DO jj = 1, jpj
272               DO ji = 1, jpi
273                  IF ( ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd ) ) .GT. 0.0 ) THEN
274                     cchl_p_3d(ji,jj,jk) = xmassc * ( ( tracer_bkg(ji,jj,jk,jpphn) * xthetapn ) +      &
275                        &                             ( tracer_bkg(ji,jj,jk,jpphd) * xthetapd )   ) /  &
276                        &                  ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd ) )
277                  ENDIF
278               END DO
279            END DO
280         END DO
281         n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) )
282         n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) )
283         n2be_d = 14.01 + ( xmassc * xthetad )
284
285         ! Call nitrogen balancing routine
286         IF (kdeps == 1) THEN
287            pinc_chltot_2d(:,:)   = pinc_chltot_3d(:,:,1)
288            cchl_p_2d(:,:)        = cchl_p_3d(:,:,1)
289            phyt_avg_bkg_2d(:,:)  = phyt_avg_bkg_3d(:,:,1)
290            pgrow_avg_bkg_2d(:,:) = pgrow_avg_bkg_3d(:,:,1)
291            ploss_avg_bkg_2d(:,:) = ploss_avg_bkg_3d(:,:,1)
292           
293            CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   &
294               &               n2be_p, n2be_z, n2be_d, assimparm,                      &
295               &               INT(pincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       &
296               &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot_2d(:,:), cchl_p_2d(:,:), &
297               &               nbal_active, phyt_avg_bkg_2d(:,:),                      &
298               &               gl_active, pgrow_avg_bkg_2d(:,:), ploss_avg_bkg_2d(:,:), &
299               &               subsurf_active, deepneg_active,                         &
300               &               deeppos_active, nutprof_active,                         &
301               &               bstate_3d, outincs_3d,                                  &
302               &               diag_active, diag,                                      &
303               &               diag_fulldepth_active, diag_fulldepth_3d )
304         ELSE
305            pmld(:,:) = 0.5
306           
307            DO jk = 1, kdeps
308               pinc_chltot_2d(:,:)   = pinc_chltot_3d(:,:,jk)
309               cchl_p_2d(:,:)        = cchl_p_3d(:,:,jk)
310               phyt_avg_bkg_2d(:,:)  = phyt_avg_bkg_3d(:,:,jk)
311               pgrow_avg_bkg_2d(:,:) = pgrow_avg_bkg_3d(:,:,jk)
312               ploss_avg_bkg_2d(:,:) = ploss_avg_bkg_3d(:,:,jk)
313               tmask_2d(:,:,1)       = tmask(:,:,jk)
314               bstate_2d(:,:,1,:)    = bstate_3d(:,:,jk,:)
315               outincs_2d(:,:,:,:)   = 0.0
316
317               CALL bio_analysis( jpi, jpj, 1, gdepw_n(:,:,2), i_tracer, modparm,            &
318                  &               n2be_p, n2be_z, n2be_d, assimparm,                         &
319                  &               INT(pincper), 1, INT(SUM(tmask_2d,3)), tmask_2d(:,:,:),    &
320                  &               pmld(:,:), pmld(:,:), pinc_chltot_2d(:,:), cchl_p_2d(:,:), &
321                  &               nbal_active, phyt_avg_bkg_2d(:,:),                         &
322                  &               gl_active, pgrow_avg_bkg_2d(:,:), ploss_avg_bkg_2d(:,:),   &
323                  &               subsurf_active, deepneg_active,                            &
324                  &               deeppos_active, nutprof_active,                            &
325                  &               bstate_2d, outincs_2d,                                     &
326                  &               diag_active, diag,                                         &
327                  &               diag_fulldepth_active, diag_fulldepth_2d )
328
329               outincs_3d(:,:,jk,:) = outincs_2d(:,:,1,:)
330            END DO
331         ENDIF
332         
333         ! Loop over each grid point partioning the increments
334         phyto_balinc(:,:,:,:) = 0.0
335         DO jk = 1, jpk
336            IF (kdeps == 1) THEN
337               jkinc = 1
338            ELSE
339               IF (jk > kdeps) THEN
340                  EXIT
341               ENDIF
342               jkinc = jk
343            ENDIF
344            DO jj = 1, jpj
345               DO ji = 1, jpi
346
347                  ! Phytoplankton
348                  IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. &
349                     & ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) .AND. &
350                     & ( pinc_chltot_3d(ji,jj,jkinc) /= 0.0 ) ) THEN
351                     IF ( ld_chltot ) THEN
352                        ! Phytoplankton nitrogen split up based on existing ratios
353                        zfrac_phn = tracer_bkg(ji,jj,jk,jpphn) / &
354                           &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd))
355                        zfrac_phd = tracer_bkg(ji,jj,jk,jpphd) / &
356                           &        (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd))
357                     ELSE IF ( ld_chldia .AND. ld_chlnon ) THEN
358                        ! Phytoplankton nitrogen split up based on assimilation increments
359                        zfrac_phn = pinc_chlnon_3d(ji,jj,jkinc) / pinc_chltot_3d(ji,jj,jkinc)
360                        zfrac_phd = pinc_chldia_3d(ji,jj,jkinc) / pinc_chltot_3d(ji,jj,jkinc)
361                     ENDIF
362
363                     ! Phytoplankton silicate split up based on existing ratios
364                     zrat_pds_phd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpphd)
365                     
366                     ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen
367                     ! Not using pinc_chltot directly as it's only 2D
368                     ! This method should give same results at surface as splitting pinc_chltot would
369                     zrat_chn_phn = tracer_bkg(ji,jj,jk,jpchn) / tracer_bkg(ji,jj,jk,jpphn)
370                     zrat_chd_phd = tracer_bkg(ji,jj,jk,jpchd) / tracer_bkg(ji,jj,jk,jpphd)
371                     
372                     phyto_balinc(ji,jj,jk,jpphn) = outincs_3d(ji,jj,jk,i_tracer(2)) * zfrac_phn
373                     phyto_balinc(ji,jj,jk,jpphd) = outincs_3d(ji,jj,jk,i_tracer(2)) * zfrac_phd
374                     phyto_balinc(ji,jj,jk,jppds) = phyto_balinc(ji,jj,jk,jpphd) * zrat_pds_phd
375                     phyto_balinc(ji,jj,jk,jpchn) = phyto_balinc(ji,jj,jk,jpphn) * zrat_chn_phn
376                     phyto_balinc(ji,jj,jk,jpchd) = phyto_balinc(ji,jj,jk,jpphd) * zrat_chd_phd
377                  ENDIF
378
379                  ! Zooplankton nitrogen split up based on existing ratios
380                  IF ( ( tracer_bkg(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpzme) > 0.0 ) ) THEN
381                     zfrac_zmi = tracer_bkg(ji,jj,jk,jpzmi) / &
382                        &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme))
383                     zfrac_zme = tracer_bkg(ji,jj,jk,jpzme) / &
384                        &        (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme))
385                     phyto_balinc(ji,jj,jk,jpzmi) = outincs_3d(ji,jj,jk,i_tracer(3)) * zfrac_zmi
386                     phyto_balinc(ji,jj,jk,jpzme) = outincs_3d(ji,jj,jk,i_tracer(3)) * zfrac_zme
387                  ENDIF
388
389                  ! Nitrogen nutrient straight from balancing scheme
390                  phyto_balinc(ji,jj,jk,jpdin) = outincs_3d(ji,jj,jk,i_tracer(1))
391
392                  ! Nitrogen detritus straight from balancing scheme
393                  phyto_balinc(ji,jj,jk,jpdet) = outincs_3d(ji,jj,jk,i_tracer(4))
394
395                  ! DIC straight from balancing scheme
396                  phyto_balinc(ji,jj,jk,jpdic) = outincs_3d(ji,jj,jk,i_tracer(5))
397
398                  ! Alkalinity straight from balancing scheme
399                  phyto_balinc(ji,jj,jk,jpalk) = outincs_3d(ji,jj,jk,i_tracer(6))
400
401                  ! Remove diatom silicate increment from nutrient silicate to conserve mass
402                  IF ( ( tracer_bkg(ji,jj,jk,jpsil) - phyto_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN
403                     phyto_balinc(ji,jj,jk,jpsil) = phyto_balinc(ji,jj,jk,jppds) * (-1.0)
404                  ENDIF
405
406                  ! Carbon detritus based on existing ratios
407                  IF ( ( tracer_bkg(ji,jj,jk,jpdet) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpdtc) > 0.0 ) ) THEN
408                     zrat_dtc_det = tracer_bkg(ji,jj,jk,jpdtc) / tracer_bkg(ji,jj,jk,jpdet)
409                     phyto_balinc(ji,jj,jk,jpdtc) = phyto_balinc(ji,jj,jk,jpdet) * zrat_dtc_det
410                  ENDIF
411
412                  ! Do nothing with iron or oxygen for the time being
413                  phyto_balinc(ji,jj,jk,jpfer) = 0.0
414                  phyto_balinc(ji,jj,jk,jpoxy) = 0.0
415                 
416               END DO
417            END DO
418         END DO
419     
420      ELSE   ! No nitrogen balancing
421     
422         ! Initialise individual chlorophyll increments to zero
423         phyto_balinc(:,:,:,jpchn) = 0.0
424         phyto_balinc(:,:,:,jpchd) = 0.0
425         
426         ! Split up total surface chlorophyll increments
427         DO jk = 1, kdeps
428            DO jj = 1, jpj
429               DO ji = 1, jpi
430                  IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. &
431                     & ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN
432                     IF ( ld_chltot ) THEN
433                        ! Chlorophyll split up based on existing ratios
434                        zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / &
435                           &        ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd) )
436                        zfrac_chd = tracer_bkg(ji,jj,jk,jpchd) / &
437                           &        ( tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd) )
438                        phyto_balinc(ji,jj,jk,jpchn) = pinc_chltot_3d(ji,jj,jk) * zfrac_chn
439                        phyto_balinc(ji,jj,jk,jpchd) = pinc_chltot_3d(ji,jj,jk) * zfrac_chd
440                     ENDIF
441                     IF( ld_chldia ) THEN
442                        phyto_balinc(ji,jj,jk,jpchd) = pinc_chldia_3d(ji,jj,jk)
443                     ENDIF
444                     IF( ld_chlnon ) THEN
445                        phyto_balinc(ji,jj,jk,jpchn) = pinc_chlnon_3d(ji,jj,jk)
446                     ENDIF
447
448                     ! Maintain stoichiometric ratios of nitrogen and silicate
449                     IF ( ld_chltot .OR. ld_chlnon ) THEN
450                        zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn)
451                        phyto_balinc(ji,jj,jk,jpphn) = phyto_balinc(ji,jj,jk,jpchn) * zrat_phn_chn
452                     ENDIF
453                     IF ( ld_chltot .OR. ld_chldia ) THEN
454                        zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd)
455                        phyto_balinc(ji,jj,jk,jpphd) = phyto_balinc(ji,jj,jk,jpchd) * zrat_phd_chd
456                        zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd)
457                        phyto_balinc(ji,jj,jk,jppds) = phyto_balinc(ji,jj,jk,jpchd) * zrat_pds_chd
458                     ENDIF
459                  ENDIF
460               END DO
461            END DO
462         END DO
463         
464         IF (kdeps == 1) THEN
465            ! Propagate through mixed layer
466            DO jj = 1, jpj
467               DO ji = 1, jpi
468                  !
469                  jkmax = jpk-1
470                  DO jk = jpk-1, 1, -1
471                     IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
472                        & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
473                        pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
474                        jkmax = jk
475                     ENDIF
476                  END DO
477                  !
478                  DO jk = 2, jkmax
479                     phyto_balinc(ji,jj,jk,jpchn) = phyto_balinc(ji,jj,1,jpchn)
480                     phyto_balinc(ji,jj,jk,jpchd) = phyto_balinc(ji,jj,1,jpchd)
481                     phyto_balinc(ji,jj,jk,jpphn) = phyto_balinc(ji,jj,1,jpphn)
482                     phyto_balinc(ji,jj,jk,jpphd) = phyto_balinc(ji,jj,1,jpphd)
483                     phyto_balinc(ji,jj,jk,jppds) = phyto_balinc(ji,jj,1,jppds)
484                  END DO
485                  !
486               END DO
487            END DO
488         ENDIF
489
490         ! Set other balancing increments to zero
491         phyto_balinc(:,:,:,jpzmi) = 0.0
492         phyto_balinc(:,:,:,jpzme) = 0.0
493         phyto_balinc(:,:,:,jpdin) = 0.0
494         phyto_balinc(:,:,:,jpsil) = 0.0
495         phyto_balinc(:,:,:,jpfer) = 0.0
496         phyto_balinc(:,:,:,jpdet) = 0.0
497         phyto_balinc(:,:,:,jpdtc) = 0.0
498         phyto_balinc(:,:,:,jpdic) = 0.0
499         phyto_balinc(:,:,:,jpalk) = 0.0
500         phyto_balinc(:,:,:,jpoxy) = 0.0
501
502      ENDIF
503     
504      ! If performing extra tidal mixing in the Indonesian Throughflow,
505      ! increments have been found to make the carbon cycle unstable
506      ! Therefore, mask these out
507      IF ( ln_tmx_itf ) THEN
508         DO jn = 1, jptra
509            DO jk = 1, jpk
510               phyto_balinc(:,:,jk,jn) = phyto_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
511            END DO
512         END DO
513      ENDIF
514
515   END SUBROUTINE asm_phyto_bal_medusa
516
517#else
518   !!----------------------------------------------------------------------
519   !!   Default option : Empty routine
520   !!----------------------------------------------------------------------
521CONTAINS
522   SUBROUTINE asm_phyto_bal_medusa( kdeps,                          &
523      &                             ld_chltot,                      &
524      &                             pinc_chltot_3d,                    &
525      &                             ld_chldia,                      &
526      &                             pinc_chldia_3d,                    &
527      &                             ld_chlnon,                      &
528      &                             pinc_chlnon_3d,                    &
529      &                             ld_phytot,                      &
530      &                             pinc_phytot_3d,                    &
531      &                             ld_phydia,                      &
532      &                             pinc_phydia_3d,                    &
533      &                             ld_phynon,                      &
534      &                             pinc_phynon_3d,                    &
535      &                             pincper,                        &
536      &                             p_maxchlinc, ld_phytobal, pmld, &
537      &                             pgrow_avg_bkg_3d, ploss_avg_bkg_3d,   &
538      &                             phyt_avg_bkg_3d, mld_max_bkg,      &
539      &                             tracer_bkg, phyto_balinc )
540      INTEGER :: kdeps
541      LOGICAL :: ld_chltot
542      REAL    :: pinc_chltot_3d(:,:,:)
543      LOGICAL :: ld_chldia
544      REAL    :: pinc_chldia_3d(:,:,:)
545      LOGICAL :: ld_chlnon
546      REAL    :: pinc_chlnon_3d(:,:,:)
547      LOGICAL :: ld_phytot
548      REAL    :: pinc_phytot_3d(:,:,:)
549      LOGICAL :: ld_phydia
550      REAL    :: pinc_phydia_3d(:,:,:)
551      LOGICAL :: ld_phynon
552      REAL    :: pinc_phynon_3d(:,:,:)
553      REAL    :: pincper
554      REAL    :: p_maxchlinc
555      LOGICAL :: ld_phytobal
556      REAL    :: pmld(:,:)
557      REAL    :: pgrow_avg_bkg_3d(:,:,:)
558      REAL    :: ploss_avg_bkg_3d(:,:,:)
559      REAL    :: phyt_avg_bkg_3d(:,:,:)
560      REAL    :: mld_max_bkg(:,:)
561      REAL    :: tracer_bkg(:,:,:,:)
562      REAL    :: phyto_balinc(:,:,:,:)
563      WRITE(*,*) 'asm_phyto_bal_medusa: You should not have seen this print! error?'
564   END SUBROUTINE asm_phyto_bal_medusa
565#endif
566
567   !!======================================================================
568END MODULE asmphytobal_medusa
Note: See TracBrowser for help on using the repository browser.