source: branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_ersem_hem08/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90 @ 9331

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

Add balancing code.

File size: 36.9 KB
Line 
1MODULE asmlogchlbal_ersem
2   !!======================================================================
3   !!                       ***  MODULE asmlogchlbal_ersem  ***
4   !! Calculate increments to ERSEM based on surface logchl increments
5   !!======================================================================
6   !! History :  3.6  ! 2016-09  (D. Ford)     Original code
7   !!----------------------------------------------------------------------
8#if defined key_asminc && defined key_fabm
9   !!----------------------------------------------------------------------
10   !!   'key_asminc'         : assimilation increment interface
11   !!   'key_fabm'           : FABM-ERSEM coupling
12   !!----------------------------------------------------------------------
13   !!   asm_logchl_bal_ersem : routine to calculate increments to ERSEM
14   !!----------------------------------------------------------------------
15   USE par_kind,    ONLY: wp             ! kind parameters
16   USE par_oce,     ONLY: jpi, jpj, jpk  ! domain array sizes
17   USE dom_oce,     ONLY: gdepw_n        ! domain information
18   USE zdfmxl                            ! mixed layer depth
19   USE iom                               ! i/o
20   USE trc,         ONLY: trn            ! ERSEM state variables
21   USE par_fabm                          ! FABM parameters
22   USE par_trc,     ONLY: jptra          ! Tracer parameters
23   USE bioanalysis
24
25   IMPLICIT NONE
26   PRIVATE                   
27
28   PUBLIC asm_logchl_bal_ersem
29
30   ! Default values for biological assimilation parameters
31   ! Should match Hemmings et al. (2008)
32   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
33   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
34   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
35   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
36   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
37   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
38   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
39   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
40   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
41   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
42   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
43   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
44   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
45   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
46   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
47   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
48   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
49   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
50   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
51   !
52   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
53   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
54   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
55   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
56   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
57   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
58   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
59   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
60
61CONTAINS
62
63   SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, &
64      &                             k_maxchlinc, logchl_bkginc, logchl_balinc,   &
65      &                             pgrow_avg_bkg, ploss_avg_bkg,           &
66      &                             phyt_avg_bkg, mld_max_bkg,              &
67      &                             ld_logchlbal, aincper )
68      !!---------------------------------------------------------------------------
69      !!                    ***  ROUTINE asm_logchl_bal_ersem  ***
70      !!
71      !! ** Purpose :   calculate increments to ERSEM from logchl increments
72      !!
73      !! ** Method  :   convert logchl increments to chl increments
74      !!                split between the ERSEM PFTs
75      !!                spread through the mixed layer
76      !!                [forthcoming: calculate increments to nutrients and zooplankton]
77      !!
78      !! ** Action  :   populate logchl_balinc
79      !!
80      !! References :   forthcoming...
81      !!---------------------------------------------------------------------------
82      !!
83      LOGICAL,  INTENT(in   )                               :: ld_logchlpftinc
84      INTEGER,  INTENT(in   )                               :: npfts
85      INTEGER,  INTENT(in   )                               :: mld_choice_bgc
86      REAL(wp), INTENT(in   )                               :: k_maxchlinc
87      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,npfts)     :: logchl_bkginc
88      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc
89      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth
90      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss
91      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto
92      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD
93      LOGICAL,  INTENT(in   )                               :: ld_logchlbal   ! Balancing y/n
94      REAL(wp), INTENT(in   )                               :: aincper
95      !!
96      INTEGER                      :: ji, jj, jk
97      INTEGER                      :: jkmax
98      REAL(wp)                     :: chl_tot, chl_inc
99      REAL(wp), DIMENSION(jpi,jpj) :: zmld
100      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
101      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters
102      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state
103      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments
104      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics
105      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics
106      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
107      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p         ! C:Chl for total phy
108      REAL(wp),                DIMENSION(jpi,jpj)           :: chlinctot
109      REAL(wp)                                              :: n2be_p         ! N:biomass for total phy
110      REAL(wp)                                              :: n2be_z         ! N:biomass for total zoo
111      REAL(wp)                                              :: n2be_d         ! N:biomass for detritus
112      REAL(wp)                                              :: zfrac_n3n      ! Fraction of jp_fabm_n3n
113      REAL(wp)                                              :: zfrac_n4n      ! Fraction of jp_fabm_n4n
114      REAL(wp)                                              :: zfrac_r4n      ! Fraction of jp_fabm_r4n
115      REAL(wp)                                              :: zfrac_r6n      ! Fraction of jp_fabm_r6n
116      REAL(wp)                                              :: zfrac_r8n      ! Fraction of jp_fabm_r8n
117      REAL(wp)                                              :: zfrac_z4n      ! Fraction of jp_fabm_z4n
118      REAL(wp)                                              :: zfrac_z5n      ! Fraction of jp_fabm_z5n
119      REAL(wp)                                              :: zfrac_z6n      ! Fraction of jp_fabm_z6n
120      REAL(wp)                                              :: zfrac_p1n      ! Fraction of jp_fabm_p1n
121      REAL(wp)                                              :: zfrac_p2n      ! Fraction of jp_fabm_p2n
122      REAL(wp)                                              :: zfrac_p3n      ! Fraction of jp_fabm_p3n
123      REAL(wp)                                              :: zfrac_p4n      ! Fraction of jp_fabm_p4n
124      REAL(wp)                                              :: zrat_z4c_z4n
125      REAL(wp)                                              :: zrat_z5c_z5n
126      REAL(wp)                                              :: zrat_z5p_z5n
127      REAL(wp)                                              :: zrat_z6c_z6n
128      REAL(wp)                                              :: zrat_z6p_z6n
129      REAL(wp)                                              :: zrat_p1c_p1n
130      REAL(wp)                                              :: zrat_p1p_p1n
131      REAL(wp)                                              :: zrat_p1s_p1n
132      REAL(wp)                                              :: zrat_p2c_p2n
133      REAL(wp)                                              :: zrat_p2p_p2n
134      REAL(wp)                                              :: zrat_p3c_p3n
135      REAL(wp)                                              :: zrat_p3p_p3n
136      REAL(wp)                                              :: zrat_p4c_p4n
137      REAL(wp)                                              :: zrat_p4p_p4n
138      !!---------------------------------------------------------------------------
139
140      ! Split surface logchl incs into surface Chl1-4 incs
141      !
142      ! In order to transform logchl incs to chl incs, need to account for the background,
143      ! cannot simply do 10^logchl_bkginc. Need to:
144      ! 1) Add logchl inc to log10(background) to get log10(analysis)
145      ! 2) Take 10^log10(analysis) to get analysis
146      ! 3) Subtract background from analysis to get chl incs
147      ! If k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
148      !
149      ! Only apply increments if all of Chl1-4 background values are > 0
150      ! In theory, they always will be, and if any are not that's a sign
151      ! that something's going wrong which the assimilation might make worse
152      !
153      IF ( ld_logchlpftinc ) THEN
154         !
155         ! Assimilating separate PFTs, so separately transform each from LogChl to Chl
156         !
157         IF ( npfts /= 4 ) THEN
158            CALL ctl_stop( 'If assimilating PFTs into ERSEM, nn_asmpfts must be 4' )
159         ENDIF
160         DO jj = 1, jpj
161            DO ji = 1, jpi
162               IF ( ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
163                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
164                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
165                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) ) THEN
166                  IF ( logchl_bkginc(ji,jj,1) /= 0.0 ) THEN
167                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) ) + &
168                        &                                                   logchl_bkginc(ji,jj,1) )                      - &
169                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
170                  ENDIF
171                  IF ( logchl_bkginc(ji,jj,2) /= 0.0 ) THEN
172                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) ) + &
173                        &                                                   logchl_bkginc(ji,jj,2) )                      - &
174                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
175                  ENDIF
176                  IF ( logchl_bkginc(ji,jj,3) /= 0.0 ) THEN
177                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) ) + &
178                        &                                                   logchl_bkginc(ji,jj,3) )                      - &
179                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
180                  ENDIF
181                  IF ( logchl_bkginc(ji,jj,4) /= 0.0 ) THEN
182                     logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = 10**( LOG10( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) ) + &
183                        &                                                   logchl_bkginc(ji,jj,4) )                      - &
184                        &                                             trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
185                  ENDIF
186                  IF (k_maxchlinc > 0.0) THEN
187                     chl_inc = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
188                               logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
189                               logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
190                               logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
191                     IF ( ABS(chl_inc) > k_maxchlinc ) THEN
192                        chl_tot = ABS(chl_inc) / k_maxchlinc
193                        logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / chl_tot
194                        logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / chl_tot
195                        logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / chl_tot
196                        logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / chl_tot
197                     ENDIF
198                  ENDIF
199               ENDIF
200            END DO
201         END DO
202      ELSE
203         !
204         ! Assimilating total Chl, so transform total from LogChl to Chl
205         ! and split between PFTs according to the existing background ratios
206         !
207         IF ( npfts /= 1 ) THEN
208            CALL ctl_stop( 'If assimilating total chlorophyll, nn_asmpfts must be 1' )
209         ENDIF
210         DO jj = 1, jpj
211            DO ji = 1, jpi
212               IF ( ( logchl_bkginc(ji,jj,1)               /= 0.0 ) .AND. &
213                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) >  0.0 ) .AND. &
214                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) >  0.0 ) .AND. &
215                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) >  0.0 ) .AND. &
216                  & ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) >  0.0 ) ) THEN
217                  chl_tot = trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
218                     &      trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
219                  chl_inc = 10**( LOG10( chl_tot ) + logchl_bkginc(ji,jj,1) ) - chl_tot
220                  IF (k_maxchlinc > 0.0) THEN
221                     chl_inc = MAX( -1.0 * k_maxchlinc, MIN( chl_inc, k_maxchlinc ) )
222                  ENDIF
223                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / chl_tot
224                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / chl_tot
225                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / chl_tot
226                  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) = chl_inc * trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / chl_tot
227               ENDIF
228            END DO
229         END DO
230      ENDIF
231
232      ! Propagate surface Chl1-4 incs through mixed layer
233      ! First, choose mixed layer definition
234      !
235      SELECT CASE( mld_choice_bgc )
236      CASE ( 1 )                   ! Turbocline/mixing depth [W points]
237         zmld(:,:) = hmld(:,:)
238      CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
239         zmld(:,:) = hmlp(:,:)
240      CASE ( 3 )                   ! Kara MLD [Interpolated]
241#if defined key_karaml
242         IF ( ln_kara ) THEN
243            zmld(:,:) = hmld_kara(:,:)
244         ELSE
245            CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
246               &           ' but ln_kara=.false.' )
247         ENDIF
248#else
249         CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
250            &           ' but is not defined' )
251#endif
252      CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
253         zmld(:,:) = hmld_tref(:,:)
254      CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
255         zmld(:,:) = hmlpt(:,:)
256      END SELECT
257
258      IF ( ld_logchlbal ) THEN
259     
260         ! Set up model parameters to be passed into Hemmings balancing routine.
261         ! For now these are hardwired to the standard HadOCC parameter values
262         ! as this is what the scheme was developed for.
263         ! Obviously, HadOCC and ERSEM are rather different models, so this
264         ! isn't ideal, but there's not always direct analogues between the two
265         ! parameter sets, so it's the easiest way to get something running.
266         ! In the longer term, some serious MarMOT-based development is required.
267         modparm(1)  = 0.1                               ! grow_sat
268         modparm(2)  = 2.0                               ! psmax
269         modparm(3)  = 0.845                             ! par
270         modparm(4)  = 0.02                              ! alpha
271         modparm(5)  = 0.05                              ! resp_rate
272         modparm(6)  = 0.05                              ! pmort_rate
273         modparm(7)  = 0.01                              ! phyto_min
274         modparm(8)  = 0.05                              ! z_mort_1
275         modparm(9)  = 1.0                               ! z_mort_2
276         modparm(10) = 6.625                             ! c2n_p
277         modparm(11) = 5.625                             ! c2n_z
278         modparm(12) = 7.5                               ! c2n_d
279         modparm(13) = 0.01                              ! graze_threshold
280         modparm(14) = 2.0                               ! holling_coef
281         modparm(15) = 0.5                               ! graze_sat
282         modparm(16) = 2.0                               ! graze_max
283
284         ! Set up assimilation parameters to be passed into balancing routine
285         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
286         assimparm(2)  = balnutext
287         assimparm(3)  = balnutmin
288         assimparm(4)  = r
289         assimparm(5)  = beta_g
290         assimparm(6)  = beta_l
291         assimparm(7)  = beta_m
292         assimparm(8)  = a_g
293         assimparm(9)  = a_l
294         assimparm(10) = a_m
295         assimparm(11) = zfracb0
296         assimparm(12) = zfracb1
297         assimparm(13) = qrfmax
298         assimparm(14) = qafmax
299         assimparm(15) = zrfmax
300         assimparm(16) = zafmax
301         assimparm(17) = prfmax
302         assimparm(18) = incphymin
303         assimparm(19) = integnstep
304         assimparm(20) = pthreshold
305
306         ! Set up external tracer indices array bstate
307         i_tracer(1) = 1   ! nutrient
308         i_tracer(2) = 2   ! phytoplankton
309         i_tracer(3) = 3   ! zooplankton
310         i_tracer(4) = 4   ! detritus
311         i_tracer(5) = 5   ! DIC
312         i_tracer(6) = 6   ! Alkalinity
313
314         ! Set background state
315         bstate(:,:,:,i_tracer(1)) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) + trn(:,:,:,jp_fabm_m1+jp_fabm_n4n)
316         bstate(:,:,:,i_tracer(2)) = trn(:,:,:,jp_fabm_m1+jp_fabm_p1n) + trn(:,:,:,jp_fabm_m1+jp_fabm_p2n) + &
317            &                        trn(:,:,:,jp_fabm_m1+jp_fabm_p3n) + trn(:,:,:,jp_fabm_m1+jp_fabm_p4n)
318         ! Z4c needs converting by qnc, hardwire for now
319         bstate(:,:,:,i_tracer(3)) = (trn(:,:,:,jp_fabm_m1+jp_fabm_z4c) * 0.0126 ) + &
320            &                        trn(:,:,:,jp_fabm_m1+jp_fabm_z5n) + trn(:,:,:,jp_fabm_m1+jp_fabm_z6n)
321         bstate(:,:,:,i_tracer(4)) = trn(:,:,:,jp_fabm_m1+jp_fabm_r4n) + trn(:,:,:,jp_fabm_m1+jp_fabm_r6n) + &
322            &                        trn(:,:,:,jp_fabm_m1+jp_fabm_r8n)
323         bstate(:,:,:,i_tracer(5)) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c)
324         bstate(:,:,:,i_tracer(6)) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3a)
325
326         ! Calculate carbon to chlorophyll ratio for combined phytoplankton
327         ! and nitrogen to biomass equivalent for PZD
328         ! Need a single number, so base on HadOCC
329         cchl_p(:,:) = 0.0
330         DO jj = 1, jpj
331            DO ji = 1, jpi
332               IF ( ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
333                  &   trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) ) .GT. 0.0 ) THEN
334                  cchl_p(ji,jj) = ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_p1c)  + trn(ji,jj,1,jp_fabm_m1+jp_fabm_p2c) + &
335                     &              trn(ji,jj,1,jp_fabm_m1+jp_fabm_p3c)  + trn(ji,jj,1,jp_fabm_m1+jp_fabm_p4c) ) /  &
336                     &            ( trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
337                     &              trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + trn(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) )
338               ENDIF
339            END DO
340         END DO
341         n2be_p = ( 14.01 + ( 12.01 * 6.625 ) ) ! / ( 14.01 + ( 12.01 * 6.625 ) )
342         n2be_z = ( 14.01 + ( 12.01 * 5.625 ) ) ! / ( 14.01 + ( 12.01 * 6.625 ) )
343         n2be_d = ( 14.01 + ( 12.01 * 7.5 )   ) ! / ( 14.01 + ( 12.01 * 6.625 ) )
344
345         chlinctot(:,:) = logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl1) + logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
346            &             logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl3) + logchl_balinc(:,:,1,jp_fabm_m1+jp_fabm_chl4)
347
348         ! Call nitrogen balancing routine
349         CALL bio_analysis( jpi, jpj, jpk, gdepw_n(:,:,2:jpk), i_tracer, modparm,   &
350            &               n2be_p, n2be_z, n2be_d, assimparm,                      &
351            &               INT(aincper), 1, INT(SUM(tmask,3)), tmask(:,:,:),       &
352            &               zmld(:,:), mld_max_bkg(:,:), chlinctot(:,:), cchl_p(:,:), &
353            &               nbal_active, phyt_avg_bkg(:,:),                         &
354            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),      &
355            &               subsurf_active, deepneg_active,                         &
356            &               deeppos_active, nutprof_active,                         &
357            &               bstate, outincs,                                        &
358            &               diag_active, diag,                                      &
359            &               diag_fulldepth_active, diag_fulldepth )
360         
361         ! Loop over each grid point partioning the increments
362         DO jk = 1, jpk
363            DO jj = 1, jpj
364               DO ji = 1, jpi
365
366                  ! Nitrogen phytoplankton from balancing scheme
367                  ! Split according to current ratios [ChlTot] or assimilation [PFTs]
368                  ! Update carbon, phosphorus and silicon according to current ratios
369                  ! Already have chlorophyll
370                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) > 0.0 ) .AND. &
371                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) > 0.0 ) .AND. &
372                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) > 0.0 ) .AND. &
373                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) > 0.0 ) ) THEN
374
375                     IF ( ld_logchlpftinc ) THEN
376                        IF ( (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
377                           &  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
378                           &  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
379                           &  logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)) > 0.0 ) THEN
380                           zfrac_p1n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) / &
381                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
382                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
383                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
384                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
385                           zfrac_p2n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) / &
386                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
387                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
388                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
389                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
390                           zfrac_p3n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) / &
391                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
392                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
393                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
394                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
395                           zfrac_p4n = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) / &
396                              &        (logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) + &
397                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) + &
398                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) + &
399                              &         logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4))
400                        ELSE
401                           zfrac_p1n = 0.25
402                           zfrac_p2n = 0.25
403                           zfrac_p3n = 0.25
404                           zfrac_p4n = 0.25
405                        ENDIF
406                     ELSE
407                        zfrac_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) / &
408                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
409                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
410                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
411                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
412                        zfrac_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) / &
413                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
414                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
415                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
416                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
417                        zfrac_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) / &
418                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
419                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
420                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
421                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
422                        zfrac_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) / &
423                           &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) + &
424                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) + &
425                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) + &
426                           &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n))
427                     ENDIF
428                     
429                     zrat_p1c_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
430                     zrat_p1p_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
431                     zrat_p1s_p1n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n)
432                     zrat_p2c_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)
433                     zrat_p2p_p2n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n)
434                     zrat_p3c_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)
435                     zrat_p3p_p3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n)
436                     zrat_p4c_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)
437                     zrat_p4p_p4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n)
438                     
439                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p1n
440                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p2n
441                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p3n
442                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) = outincs(ji,jj,jk,i_tracer(2)) * zfrac_p4n
443                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1c_p1n
444                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1p_p1n
445                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) * zrat_p1s_p1n
446                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) * zrat_p2c_p2n
447                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) * zrat_p2p_p2n
448                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) * zrat_p3c_p3n
449                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) * zrat_p3p_p3n
450                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) * zrat_p4c_p4n
451                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) * zrat_p4p_p4n
452
453                  ENDIF
454
455                  ! Nitrogen nutrient from balancing scheme
456                  ! Split between nitrate and ammonium according to current ratios
457                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) > 0.0 ) .AND. &
458                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) > 0.0 ) ) THEN
459                     zfrac_n3n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) / &
460                        &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) + trn(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n))
461                     zfrac_n4n = 1.0 - zfrac_n3n
462                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n3n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n3n
463                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_n4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_n4n
464                  ENDIF
465
466                  ! Nitrogen zooplankton from balancing scheme
467                  ! Split according to current ratios
468                  ! Update carbon and phosphorus according to current ratios
469                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) > 0.0 ) .AND. &
470                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) > 0.0 ) .AND. &
471                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) > 0.0 ) ) THEN
472                     zfrac_z4n = (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) / &
473                        &        ((trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) + &
474                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
475                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n))
476                     zfrac_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) / &
477                        &        ((trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) * 0.0126) + &
478                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) + &
479                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n))
480                     zfrac_z6n = 1.0 - zfrac_z4n - zfrac_z5n
481                     zrat_z4c_z4n = 1.0 / 0.0126
482                     zrat_z5c_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n)
483                     zrat_z5p_z5n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n)
484                     zrat_z6c_z6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6c) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n)
485                     zrat_z6p_z6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) / trn(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n)
486                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z5n
487                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z6n
488                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z4c) = outincs(ji,jj,jk,i_tracer(3)) * zfrac_z4n * zrat_z4c_z4n
489                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) * zrat_z5c_z5n
490                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6c) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) * zrat_z6c_z6n
491                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z5n) * zrat_z5p_z5n
492                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6p) = logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_z6n) * zrat_z6p_z6n
493                  ENDIF
494
495                  ! Nitrogen detritus from balancing scheme
496                  ! Split according to current ratios
497                  IF ( ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) > 0.0 ) .AND. &
498                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) > 0.0 ) .AND. &
499                     & ( trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) > 0.0 ) ) THEN
500                     zfrac_r4n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) / &
501                        &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
502                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
503                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
504                     zfrac_r6n = trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) / &
505                        &        (trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) + &
506                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) + &
507                        &         trn(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n))
508                     zfrac_r8n = 1.0 - zfrac_r4n - zfrac_r6n
509                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r4n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r4n
510                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r6n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r6n
511                     logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_r8n) = outincs(ji,jj,jk,i_tracer(1)) * zfrac_r8n
512                  ENDIF
513
514                  ! DIC straight from balancing scheme
515                  logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3c) = outincs(ji,jj,jk,i_tracer(5))
516
517                  ! Alkalinity straight from balancing scheme
518                  logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_o3a) = outincs(ji,jj,jk,i_tracer(6))
519     
520               END DO
521            END DO
522         END DO
523     
524      ENDIF
525      !
526      ! Now set MLD to bottom of a level and propagate chlorophyll incs equally through mixed layer
527      ! If balancing, should really relate this back to phytoplankton, but stick with this for now
528      !
529      DO jj = 1, jpj
530         DO ji = 1, jpi
531            !
532            jkmax = jpk-1
533            DO jk = jpk-1, 1, -1
534               IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
535                  & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
536                  zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
537                  jkmax = jk
538               ENDIF
539            END DO
540            !
541            DO jk = 2, jkmax
542               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1)
543               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2)
544               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3)
545               logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4)
546            END DO
547            !
548         END DO
549      END DO
550
551   END SUBROUTINE asm_logchl_bal_ersem
552
553#else
554   !!----------------------------------------------------------------------
555   !!   Default option : Empty routine
556   !!----------------------------------------------------------------------
557CONTAINS
558   SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, &
559      &                             k_maxchlinc, logchl_bkginc, logchl_balinc )
560      LOGICAL :: ld_logchlpftinc
561      INTEGER :: npfts
562      INTEGER :: mld_choice_bgc
563      REAL    :: k_maxchlinc
564      REAL    :: logchl_bkginc(:,:,:)
565      REAL    :: logchl_balinc(:,:,:,:)
566      WRITE(*,*) 'asm_logchl_bal_ersem: You should not have seen this print! error?'
567   END SUBROUTINE asm_logchl_bal_ersem
568#endif
569
570   !!======================================================================
571END MODULE asmlogchlbal_ersem
Note: See TracBrowser for help on using the repository browser.