New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
asmlogchlbal_medusa.F90 in branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

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

Last change on this file since 8495 was 8495, checked in by dford, 7 years ago

Merge in changes from dev_r5518_GO6_package_asm_surf_bgc, and adapt to the updated MEDUSA structure.

File size: 23.7 KB
RevLine 
[8436]1MODULE asmlogchlbal_medusa
[8428]2   !!======================================================================
[8436]3   !!                       ***  MODULE asmlogchlbal_medusa  ***
4   !! Calculate increments to MEDUSA based on surface logchl increments
[8428]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   !!======================================================================
[8436]12   !! History :  3.6  ! 2017-08 (D. Ford)  Adapted from asmlogchlbal_hadocc
[8428]13   !!----------------------------------------------------------------------
[8436]14#if defined key_asminc && defined key_medusa && defined key_foam_medusa
[8428]15   !!----------------------------------------------------------------------
16   !! 'key_asminc'          : assimilation increment interface
[8436]17   !! 'key_medusa'          : MEDUSA model
18   !! 'key_foam_medusa'     : MEDUSA extras for FOAM OBS and ASM
[8428]19   !!----------------------------------------------------------------------
[8436]20   !! asm_logchl_bal_medusa : routine to calculate increments to MEDUSA
[8428]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
[8485]26   USE zdftmx,        ONLY: ln_tmx_itf, &  ! Indonesian Throughflow
27      &                     mask_itf       ! tidal mixing mask
[8428]28   USE iom                                 ! i/o
[8436]29   USE sms_medusa                          ! MEDUSA parameters
30   USE par_medusa                          ! MEDUSA parameters
[8428]31   USE par_trc,       ONLY: jptra          ! Tracer parameters
32   USE bioanalysis                         ! Nitrogen balancing
33
34   IMPLICIT NONE
35   PRIVATE                   
36
[8436]37   PUBLIC asm_logchl_bal_medusa
[8428]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
[8436]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,              &
[8440]76      &                              tracer_bkg, logchl_balinc )
[8428]77      !!---------------------------------------------------------------------------
[8436]78      !!                    ***  ROUTINE asm_logchl_bal_medusa  ***
[8428]79      !!
[8436]80      !! ** Purpose :   calculate increments to MEDUSA from logchl increments
[8428]81      !!
82      !! ** Method  :   convert logchl increments to chl increments
[8436]83      !!                average up MEDUSA to look like HadOCC
[8428]84      !!                call nitrogen balancing scheme
[8436]85      !!                separate back out to MEDUSA
[8428]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
[8436]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
[8440]102      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables
[8428]103      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc  ! Balancing increments
104      !!
[8436]105      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
[8428]106      INTEGER                                               :: jkmax          ! Loop index
107      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
[8436]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
[8428]121      REAL(wp),                DIMENSION(jpi,jpj)           :: chl_inc        ! Chlorophyll increments
[8436]122      REAL(wp),                DIMENSION(jpi,jpj)           :: medusa_chl     ! MEDUSA total chlorophyll
123      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy
[8428]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
[8440]140      medusa_chl(:,:) = tracer_bkg(:,:,1,jpchn) + tracer_bkg(:,:,1,jpchd)
[8428]141      DO jj = 1, jpj
142         DO ji = 1, jpi
[8436]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)
[8428]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
[8436]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
[8428]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
[8440]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)
[8428]243
[8436]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
[8440]250               IF ( ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) ) .GT. 0.0 ) THEN
251                  cchl_p(ji,jj) = ( ( tracer_bkg(ji,jj,1,jpphn) * xthetapn ) + ( tracer_bkg(ji,jj,1,jpphd) * xthetapd ) ) / &
252                     &            ( tracer_bkg(ji,jj,1,jpchn) + tracer_bkg(ji,jj,1,jpchd ) )
[8436]253               ENDIF
254            END DO
255         END DO
256         n2be_p = 14.01 + ( xmassc * ( ( xthetapn  + xthetapd  ) / 2.0 ) )
257         n2be_z = 14.01 + ( xmassc * ( ( xthetazmi + xthetazme ) / 2.0 ) )
258         n2be_d = 14.01 + ( xmassc * xthetad )
259
[8428]260         ! Call nitrogen balancing routine
[8436]261         CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   &
262            &               n2be_p, n2be_z, n2be_d, assimparm,                      &
263            &               INT(aincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       &
264            &               zmld(:,:), mld_max_bkg(:,:), chl_inc(:,:), cchl_p(:,:), &
265            &               nbal_active, phyt_avg_bkg(:,:),                         &
266            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      &
267            &               subsurf_active, deepneg_active,                         &
268            &               deeppos_active, nutprof_active,                         &
269            &               bstate, outincs,                                        &
270            &               diag_active, diag,                                      &
[8428]271            &               diag_fulldepth_active, diag_fulldepth )
[8436]272         
273         ! Loop over each grid point partioning the increments
274         logchl_balinc(:,:,:,:) = 0.0
275         DO jk = 1, jpk
276            DO jj = 1, jpj
277               DO ji = 1, jpi
278
[8440]279                  IF ( ( tracer_bkg(ji,jj,jk,jpphn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpphd) > 0.0 ) ) THEN
[8436]280                     ! Phytoplankton nitrogen and silicate split up based on existing ratios
[8440]281                     zfrac_phn = tracer_bkg(ji,jj,jk,jpphn) / (tracer_bkg(ji,jj,jk,jpphn) + tracer_bkg(ji,jj,jk,jpphd))
[8436]282                     zfrac_phd = 1.0 - zfrac_phn
[8440]283                     zrat_pds_phd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpphd)
[8436]284                     logchl_balinc(ji,jj,jk,jpphn) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phn
285                     logchl_balinc(ji,jj,jk,jpphd) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_phd
286                     logchl_balinc(ji,jj,jk,jppds) = logchl_balinc(ji,jj,jk,jpphd) * zrat_pds_phd
287
288                     ! Chlorophyll split up based on existing ratios to phytoplankton nitrogen
289                     ! Not using chl_inc directly as it's only 2D
290                     ! This method should give same results at surface as splitting chl_inc would
[8440]291                     zrat_chn_phn = tracer_bkg(ji,jj,jk,jpchn) / tracer_bkg(ji,jj,jk,jpphn)
292                     zrat_chd_phd = tracer_bkg(ji,jj,jk,jpchd) / tracer_bkg(ji,jj,jk,jpphd)
[8436]293                     logchl_balinc(ji,jj,jk,jpchn) = logchl_balinc(ji,jj,jk,jpphn) * zrat_chn_phn
294                     logchl_balinc(ji,jj,jk,jpchd) = logchl_balinc(ji,jj,jk,jpphd) * zrat_chd_phd
295                  ENDIF
296
[8440]297                  IF ( ( tracer_bkg(ji,jj,jk,jpzmi) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpzme) > 0.0 ) ) THEN
[8436]298                     ! Zooplankton nitrogen split up based on existing ratios
[8440]299                     zfrac_zmi = tracer_bkg(ji,jj,jk,jpzmi) / (tracer_bkg(ji,jj,jk,jpzmi) + tracer_bkg(ji,jj,jk,jpzme))
[8436]300                     zfrac_zme = 1.0 - zfrac_zmi
301                     logchl_balinc(ji,jj,jk,jpzmi) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zmi
302                     logchl_balinc(ji,jj,jk,jpzme) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_zme
303                  ENDIF
304
305                  ! Nitrogen nutrient straight from balancing scheme
306                  logchl_balinc(ji,jj,jk,jpdin) = outincs(ji,jj,jk,i_tracer(1))
307
308                  ! Nitrogen detritus straight from balancing scheme
309                  logchl_balinc(ji,jj,jk,jpdet) = outincs(ji,jj,jk,i_tracer(4))
310
311                  ! DIC straight from balancing scheme
312                  logchl_balinc(ji,jj,jk,jpdic) = outincs(ji,jj,jk,i_tracer(5))
313
314                  ! Alkalinity straight from balancing scheme
315                  logchl_balinc(ji,jj,jk,jpalk) = outincs(ji,jj,jk,i_tracer(6))
316
317                  ! Remove diatom silicate increment from nutrient silicate to conserve mass
[8440]318                  IF ( ( tracer_bkg(ji,jj,jk,jpsil) - logchl_balinc(ji,jj,jk,jppds) ) > 0.0 ) THEN
[8436]319                     logchl_balinc(ji,jj,jk,jpsil) = logchl_balinc(ji,jj,jk,jppds) * (-1.0)
320                  ENDIF
321
[8440]322                  IF ( ( tracer_bkg(ji,jj,jk,jpdet) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpdtc) > 0.0 ) ) THEN
[8436]323                     ! Carbon detritus based on existing ratios
[8440]324                     zrat_dtc_det = tracer_bkg(ji,jj,jk,jpdtc) / tracer_bkg(ji,jj,jk,jpdet)
[8436]325                     logchl_balinc(ji,jj,jk,jpdtc) = logchl_balinc(ji,jj,jk,jpdet) * zrat_dtc_det
326                  ENDIF
327
328                  ! Do nothing with iron or oxygen for the time being
329                  logchl_balinc(ji,jj,jk,jpfer) = 0.0
330                  logchl_balinc(ji,jj,jk,jpoxy) = 0.0
331                 
332               END DO
333            END DO
334         END DO
[8428]335     
336      ELSE   ! No nitrogen balancing
337     
[8436]338         ! Initialise individual chlorophyll increments to zero
339         logchl_balinc(:,:,:,jpchn) = 0.0
340         logchl_balinc(:,:,:,jpchd) = 0.0
[8428]341         
[8436]342         ! Split up total surface chlorophyll increments
343         DO jj = 1, jpj
344            DO ji = 1, jpi
345               IF ( medusa_chl(ji,jj) > 0.0 ) THEN
[8440]346                  zfrac_chn = tracer_bkg(ji,jj,1,jpchn) / medusa_chl(ji,jj)
[8436]347                  zfrac_chd = 1.0 - zfrac_chn
348                  logchl_balinc(ji,jj,1,jpchn) = chl_inc(ji,jj) * zfrac_chn
349                  logchl_balinc(ji,jj,1,jpchd) = chl_inc(ji,jj) * zfrac_chd
350               ENDIF
351            END DO
352         END DO
[8428]353         
354         ! Propagate through mixed layer
355         DO jj = 1, jpj
356            DO ji = 1, jpi
357               !
358               jkmax = jpk-1
359               DO jk = jpk-1, 1, -1
360                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
361                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
362                     zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
363                     jkmax = jk
364                  ENDIF
365               END DO
366               !
367               DO jk = 2, jkmax
[8436]368                  logchl_balinc(ji,jj,jk,jpchn) = logchl_balinc(ji,jj,1,jpchn)
369                  logchl_balinc(ji,jj,jk,jpchd) = logchl_balinc(ji,jj,1,jpchd)
[8428]370               END DO
371               !
372            END DO
373         END DO
374
375         ! Set other balancing increments to zero
[8436]376         logchl_balinc(:,:,:,jpphn) = 0.0
377         logchl_balinc(:,:,:,jpphd) = 0.0
378         logchl_balinc(:,:,:,jppds) = 0.0
379         logchl_balinc(:,:,:,jpzmi) = 0.0
380         logchl_balinc(:,:,:,jpzme) = 0.0
381         logchl_balinc(:,:,:,jpdin) = 0.0
382         logchl_balinc(:,:,:,jpsil) = 0.0
383         logchl_balinc(:,:,:,jpfer) = 0.0
384         logchl_balinc(:,:,:,jpdet) = 0.0
385         logchl_balinc(:,:,:,jpdtc) = 0.0
386         logchl_balinc(:,:,:,jpdic) = 0.0
387         logchl_balinc(:,:,:,jpalk) = 0.0
388         logchl_balinc(:,:,:,jpoxy) = 0.0
389
[8428]390      ENDIF
[8485]391     
392      ! If performing extra tidal mixing in the Indonesian Throughflow,
393      ! increments have been found to make the carbon cycle unstable
394      ! Therefore, mask these out
395      IF ( ln_tmx_itf ) THEN
396         DO jn = 1, jptra
397            DO jk = 1, jpk
398               logchl_balinc(:,:,jk,jn) = logchl_balinc(:,:,jk,jn) * ( 1.0 - mask_itf(:,:) )
399            END DO
400         END DO
401      ENDIF
[8428]402
[8436]403   END SUBROUTINE asm_logchl_bal_medusa
[8428]404
405#else
406   !!----------------------------------------------------------------------
407   !!   Default option : Empty routine
408   !!----------------------------------------------------------------------
409CONTAINS
[8436]410   SUBROUTINE asm_logchl_bal_medusa( logchl_bkginc, aincper, mld_choice_bgc, &
[8440]411      &                              k_maxchlinc, ld_logchlbal,              &
412      &                              pgrow_avg_bkg, ploss_avg_bkg,           &
413      &                              phyt_avg_bkg, mld_max_bkg,              &
414      &                              tracer_bkg, logchl_balinc )
[8428]415      REAL    :: logchl_bkginc(:,:)
416      REAL    :: aincper
417      INTEGER :: mld_choice_bgc
418      REAL    :: k_maxchlinc
[8440]419      LOGICAL :: ld_logchlbal
420      REAL    :: pgrow_avg_bkg(:,:)
421      REAL    :: ploss_avg_bkg(:,:)
422      REAL    :: phyt_avg_bkg(:,:)
423      REAL    :: mld_max_bkg(:,:)
424      REAL    :: tracer_bkg(:,:,:,:)
425      REAL    :: logchl_balinc(:,:,:,:)
[8436]426      WRITE(*,*) 'asm_logchl_bal_medusa: You should not have seen this print! error?'
427   END SUBROUTINE asm_logchl_bal_medusa
[8428]428#endif
429
430   !!======================================================================
[8436]431END MODULE asmlogchlbal_medusa
Note: See TracBrowser for help on using the repository browser.