source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_medusa.F90 @ 8497

Last change on this file since 8497 was 8497, checked in by dford, 3 years ago

Include conversion from moles to grams in carbon to chlorophyll ratio calculation.

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