source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmphyto2dbal_medusa.F90 @ 10149

Last change on this file since 10149 was 10149, checked in by frrh, 23 months ago

Met Office GMED ticket 379: Merged David Ford's MEDUSA assimilation changes
using command:

svn merge -r 10054:10141 svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_v3

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