source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_medusa.F90 @ 9435

Last change on this file since 9435 was 9435, checked in by dford, 2 years ago

Allow assimilation of PFT chlorophyll.

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