source: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmphytobal_hadocc.F90 @ 13354

Last change on this file since 13354 was 13354, checked in by jwhile, 8 weeks ago

Merged in trunk

File size: 18.8 KB
Line 
1MODULE asmphytobal_hadocc
2   !!======================================================================
3   !!                       ***  MODULE asmphytobal_hadocc  ***
4   !! Calculate increments to HadOCC based on surface phyto 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 bioanal.F90
13   !!----------------------------------------------------------------------
14#if defined key_asminc && defined key_hadocc
15   !!----------------------------------------------------------------------
16   !! 'key_asminc'          : assimilation increment interface
17   !! 'key_hadocc'          : HadOCC model
18   !!----------------------------------------------------------------------
19   !! asm_phyto_bal_hadocc : routine to calculate increments to HadOCC
20   !!----------------------------------------------------------------------
21   USE par_kind,      ONLY: wp             ! kind parameters
22   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
23   USE dom_oce,       ONLY: gdepw_n        ! domain information
24   USE iom                                 ! i/o
25   USE par_hadocc                          ! HadOCC parameters
26   USE had_bgc_stnd,  ONLY: kmt            ! HadOCC parameters
27   USE had_bgc_const                       ! HadOCC parameters
28   USE par_trc,       ONLY: jptra          ! Tracer parameters
29   USE bioanalysis                         ! Nitrogen balancing
30
31   IMPLICIT NONE
32   PRIVATE                   
33
34   PUBLIC asm_phyto_bal_hadocc
35
36   ! Default values for biological assimilation parameters
37   ! Should match Hemmings et al. (2008)
38   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
39   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
40   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
41   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
42   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
43   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
44   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
45   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
46   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
47   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
48   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
49   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
50   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
51   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
52   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
53   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
54   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
55   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
56   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
57   !
58   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
59   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
60   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
61   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
62   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
63   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
64   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
65   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
66
67CONTAINS
68
69   SUBROUTINE asm_phyto_bal_hadocc( kdeps,                              &
70      &                             ld_chltot,                          &
71      &                             pinc_chltot_3d,                     &
72      &                             ld_phytot,                          &
73      &                             pinc_phytot_3d,                     &
74      &                             pincper,                            &
75      &                             p_maxchlinc, ld_phytobal, pmld,     &
76      &                             pgrow_avg_bkg_3d, ploss_avg_bkg_3d, &
77      &                             phyt_avg_bkg_3d, mld_max_bkg,       &
78      &                             cchl_p_bkg_3d,                      &
79      &                             tracer_bkg, phyto_balinc )
80      !!---------------------------------------------------------------------------
81      !!                    ***  ROUTINE asm_phyto_bal_hadocc  ***
82      !!
83      !! ** Purpose :   calculate increments to HadOCC from 2d phytoplankton increments
84      !!
85      !! ** Method  :   call nitrogen balancing scheme
86      !!
87      !! ** Action  :   populate phyto_balinc
88      !!
89      !! References :   Hemmings et al., 2008, J. Mar. Res.
90      !!                Ford et al., 2012, Ocean Sci.
91      !!---------------------------------------------------------------------------
92      !!
93      INTEGER,  INTENT(in   )                               :: kdeps            ! No. inc deps 1 or jpk
94      LOGICAL,  INTENT(in   )                               :: ld_chltot        ! Assim chltot y/n
95      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_chltot_3d   ! chltot increments
96      LOGICAL,  INTENT(in   )                               :: ld_phytot        ! Assim phytot y/n
97      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,kdeps)     :: pinc_phytot_3d   ! phytot increments
98      REAL(wp), INTENT(in   )                               :: pincper          ! Assimilation period
99      REAL(wp), INTENT(in   )                               :: p_maxchlinc      ! Max chl increment
100      LOGICAL,  INTENT(in   )                               :: ld_phytobal      ! Balancing y/n
101      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld             ! Mixed layer depth
102      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: pgrow_avg_bkg_3d ! Avg phyto growth
103      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: ploss_avg_bkg_3d ! Avg phyto loss
104      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: phyt_avg_bkg_3d  ! Avg phyto
105      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg      ! Max MLD
106      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,kdeps)     :: cchl_p_bkg_3d    ! C:Chl
107      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg       ! State variables
108      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto_balinc     ! Balancing increments
109      !!
110      INTEGER                                               :: ji, jj, jk, jn    ! Loop counters
111      INTEGER                                               :: jkmax             ! Loop index
112      INTEGER,                 DIMENSION(6)                 :: i_tracer          ! Tracer indices
113      REAL(wp),                DIMENSION(16)                :: modparm           ! Model parameters
114      REAL(wp),                DIMENSION(20)                :: assimparm         ! Assimilation parameters
115      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: bstate_2d         ! Background state (2D)
116      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate_3d         ! Background state (3D)
117      REAL(wp),                DIMENSION(jpi,jpj,1,6)       :: outincs_2d        ! Balancing increments (2D)
118      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs_3d        ! Balancing increments (3D)
119      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag              ! Depth-indep diagnostics
120      REAL(wp),                DIMENSION(jpi,jpj,1,22)      :: diag_fulldepth_2d ! Full-depth diagnostics (2D)
121      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth_3d ! Full-depth diagnostics (3D)
122      REAL(wp),                DIMENSION(jpi,jpj)           :: cchl_p_bkg_2d     ! C:Chl for total phy (2D)
123      REAL(wp),                DIMENSION(jpi,jpj,1)         :: tmask_2d          ! Single-level tmask
124      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_chltot_2d    ! chltot increments (2D)
125      REAL(wp),                DIMENSION(jpi,jpj)           :: pinc_phytot_2d    ! phytot increments (2D)
126      REAL(wp),                DIMENSION(jpi,jpj)           :: pgrow_avg_bkg_2d  ! Avg phyto growth (2D)
127      REAL(wp),                DIMENSION(jpi,jpj)           :: ploss_avg_bkg_2d  ! Avg phyto loss (2D)
128      REAL(wp),                DIMENSION(jpi,jpj)           :: phyt_avg_bkg_2d   ! Avg phyto (2D)
129      !!---------------------------------------------------------------------------
130
131      IF ( ( .NOT. ld_chltot ) .AND. ( .NOT. ld_phytot ) ) THEN
132         CALL ctl_stop( ' Trying to do phyto balancing but nothing to assimilate' )
133      ENDIF
134     
135      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
136      IF ( p_maxchlinc > 0.0 ) THEN
137         IF ( ld_chltot ) THEN
138            DO jk = 1, kdeps
139               DO jj = 1, jpj
140                  DO ji = 1, jpi
141                     pinc_chltot_3d(ji,jj,jk) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot_3d(ji,jj,jk), p_maxchlinc ) )
142                  END DO
143               END DO
144            END DO
145         ENDIF
146      ENDIF
147
148      IF ( ld_phytot ) THEN
149         CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' )
150      ENDIF
151     
152      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
153
154         ! Set up model parameters to be passed into Hemmings balancing routine
155         modparm(1)  = grow_sat
156         modparm(2)  = psmax
157         modparm(3)  = par
158         modparm(4)  = alpha
159         modparm(5)  = resp_rate
160         modparm(6)  = pmort_rate
161         modparm(7)  = phyto_min
162         modparm(8)  = z_mort_1
163         modparm(9)  = z_mort_2
164         modparm(10) = c2n_p
165         modparm(11) = c2n_z
166         modparm(12) = c2n_d
167         modparm(13) = graze_threshold
168         modparm(14) = holling_coef
169         modparm(15) = graze_sat
170         modparm(16) = graze_max
171
172         ! Set up assimilation parameters to be passed into balancing routine
173         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
174         assimparm(2)  = balnutext
175         assimparm(3)  = balnutmin
176         assimparm(4)  = r
177         assimparm(5)  = beta_g
178         assimparm(6)  = beta_l
179         assimparm(7)  = beta_m
180         assimparm(8)  = a_g
181         assimparm(9)  = a_l
182         assimparm(10) = a_m
183         assimparm(11) = zfracb0
184         assimparm(12) = zfracb1
185         assimparm(13) = qrfmax
186         assimparm(14) = qafmax
187         assimparm(15) = zrfmax
188         assimparm(16) = zafmax
189         assimparm(17) = prfmax
190         assimparm(18) = incphymin
191         assimparm(19) = integnstep
192         assimparm(20) = pthreshold
193
194         ! Set up external tracer indices array bstate
195         i_tracer(1) = 1   ! nutrient
196         i_tracer(2) = 2   ! phytoplankton
197         i_tracer(3) = 3   ! zooplankton
198         i_tracer(4) = 4   ! detritus
199         i_tracer(5) = 5   ! DIC
200         i_tracer(6) = 6   ! Alkalinity
201
202         ! Set background state
203         bstate_3d(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jp_had_nut)
204         bstate_3d(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jp_had_phy)
205         bstate_3d(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jp_had_zoo)
206         bstate_3d(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jp_had_pdn)
207         bstate_3d(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jp_had_dic)
208         bstate_3d(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jp_had_alk)
209
210         ! Call nitrogen balancing routine
211         IF (kdeps == 1) THEN
212            pinc_chltot_2d(:,:)   = pinc_chltot_3d(:,:,1)
213            cchl_p_bkg_2d(:,:)    = cchl_p_bkg_3d(:,:,1)
214            phyt_avg_bkg_2d(:,:)  = phyt_avg_bkg_3d(:,:,1)
215            pgrow_avg_bkg_2d(:,:) = pgrow_avg_bkg_3d(:,:,1)
216            ploss_avg_bkg_2d(:,:) = ploss_avg_bkg_3d(:,:,1)
217           
218            CALL bio_analysis( jpi, jpj, jpk, ZDZ(:,:,:), i_tracer, modparm,               &
219               &               n2be_p, n2be_z, n2be_d, assimparm,                          &
220               &               INT(pincper), 1, kmt(:,:), tmask(:,:,:),                    &
221               &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot_2d(:,:), cchl_p_bkg_2d(:,:), &
222               &               nbal_active, phyt_avg_bkg_2d(:,:),                          &
223               &               gl_active, pgrow_avg_bkg_2d(:,:), ploss_avg_bkg_2d(:,:),    &
224               &               subsurf_active, deepneg_active,                             &
225               &               deeppos_active, nutprof_active,                             &
226               &               bstate_3d, outincs_3d,                                      &
227               &               diag_active, diag,                                          &
228               &               diag_fulldepth_active, diag_fulldepth_3d )
229         ELSE
230            pmld(:,:) = 0.5
231           
232            DO jk = 1, kdeps
233               pinc_chltot_2d(:,:)   = pinc_chltot_3d(:,:,jk)
234               cchl_p_bkg_2d(:,:)    = cchl_p_bkg_3d(:,:,jk)
235               phyt_avg_bkg_2d(:,:)  = phyt_avg_bkg_3d(:,:,jk)
236               pgrow_avg_bkg_2d(:,:) = pgrow_avg_bkg_3d(:,:,jk)
237               ploss_avg_bkg_2d(:,:) = ploss_avg_bkg_3d(:,:,jk)
238               tmask_2d(:,:,1)       = tmask(:,:,jk)
239               bstate_2d(:,:,1,:)    = bstate_3d(:,:,jk,:)
240               outincs_2d(:,:,:,:)   = 0.0
241
242               CALL bio_analysis( jpi, jpj, 1, gdepw_n(:,:,2), i_tracer, modparm,            &
243                  &               n2be_p, n2be_z, n2be_d, assimparm,                         &
244                  &               INT(pincper), 1, INT(SUM(tmask_2d,3)), tmask_2d(:,:,:),    &
245                  &               pmld(:,:), pmld(:,:), pinc_chltot_2d(:,:), cchl_p_bkg_2d(:,:), &
246                  &               nbal_active, phyt_avg_bkg_2d(:,:),                         &
247                  &               gl_active, pgrow_avg_bkg_2d(:,:), ploss_avg_bkg_2d(:,:),   &
248                  &               subsurf_active, deepneg_active,                            &
249                  &               deeppos_active, nutprof_active,                            &
250                  &               bstate_2d, outincs_2d,                                     &
251                  &               diag_active, diag,                                         &
252                  &               diag_fulldepth_active, diag_fulldepth_2d )
253
254               outincs_3d(:,:,jk,:) = outincs_2d(:,:,1,:)
255            END DO
256         ENDIF
257
258         ! Save balancing increments
259         phyto_balinc(:,:,:,jp_had_nut) = outincs_3d(:,:,:,i_tracer(1))
260         phyto_balinc(:,:,:,jp_had_phy) = outincs_3d(:,:,:,i_tracer(2))
261         phyto_balinc(:,:,:,jp_had_zoo) = outincs_3d(:,:,:,i_tracer(3))
262         phyto_balinc(:,:,:,jp_had_pdn) = outincs_3d(:,:,:,i_tracer(4))
263         phyto_balinc(:,:,:,jp_had_dic) = outincs_3d(:,:,:,i_tracer(5))
264         phyto_balinc(:,:,:,jp_had_alk) = outincs_3d(:,:,:,i_tracer(6))
265     
266      ELSE   ! No nitrogen balancing
267     
268         ! Initialise phytoplankton increment to zero
269         phyto_balinc(:,:,:,jp_had_phy) = 0.0
270         
271         ! Convert surface chlorophyll increment to phytoplankton nitrogen
272         DO jk = 1, kdeps
273            phyto_balinc(:,:,jk,jp_had_phy) = ( cchl_p_bkg_3d(:,:,jk) / (mw_carbon * c2n_p) ) * pinc_chltot_3d(:,:,jk)
274         END DO
275         
276         IF (kdeps == 1) THEN
277            ! Propagate through mixed layer
278            DO jj = 1, jpj
279               DO ji = 1, jpi
280                  !
281                  jkmax = jpk-1
282                  DO jk = jpk-1, 1, -1
283                     IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
284                        & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
285                        pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
286                        jkmax = jk
287                     ENDIF
288                  END DO
289                  !
290                  DO jk = 2, jkmax
291                     phyto_balinc(ji,jj,jk,jp_had_phy) = phyto_balinc(ji,jj,1,jp_had_phy)
292                  END DO
293                  !
294               END DO
295            END DO
296         ENDIF
297
298         ! Set other balancing increments to zero
299         phyto_balinc(:,:,:,jp_had_nut) = 0.0
300         phyto_balinc(:,:,:,jp_had_zoo) = 0.0
301         phyto_balinc(:,:,:,jp_had_pdn) = 0.0
302         phyto_balinc(:,:,:,jp_had_dic) = 0.0
303         phyto_balinc(:,:,:,jp_had_alk) = 0.0
304     
305      ENDIF
306
307   END SUBROUTINE asm_phyto_bal_hadocc
308
309#else
310   !!----------------------------------------------------------------------
311   !!   Default option : Empty routine
312   !!----------------------------------------------------------------------
313CONTAINS
314   SUBROUTINE asm_phyto_bal_hadocc( kdeps,                              &
315      &                             ld_chltot,                          &
316      &                             pinc_chltot_3d,                     &
317      &                             ld_phytot,                          &
318      &                             pinc_phytot,                        &
319      &                             pincper,                            &
320      &                             p_maxchlinc, ld_phytobal, pmld,     &
321      &                             pgrow_avg_bkg_3d, ploss_avg_bkg_3d, &
322      &                             phyt_avg_bkg_3d, mld_max_bkg,       &
323      &                             cchl_p_bkg_3d,                      &
324      &                             tracer_bkg, phyto_balinc )
325      INTEGER :: kdeps
326      LOGICAL :: ld_chltot
327      REAL    :: pinc_chltot_3d(:,:,:)
328      LOGICAL :: ld_phytot
329      REAL    :: pinc_phytot_3d(:,:,:)
330      REAL    :: pincper
331      REAL    :: p_maxchlinc
332      LOGICAL :: ld_phytobal
333      REAL    :: pmld(:,:)
334      REAL    :: pgrow_avg_bkg_3d(:,:,:)
335      REAL    :: ploss_avg_bkg_3d(:,:,:)
336      REAL    :: phyt_avg_bkg_3d(:,:,:)
337      REAL    :: mld_max_bkg(:,:)
338      REAL    :: cchl_p_bkg_3d(:,:,:)
339      REAL    :: tracer_bkg(:,:,:,:)
340      REAL    :: phyto_balinc(:,:,:,:)
341      WRITE(*,*) 'asm_phyto_bal_hadocc: You should not have seen this print! error?'
342   END SUBROUTINE asm_phyto_bal_hadocc
343#endif
344
345   !!======================================================================
346END MODULE asmphytobal_hadocc
Note: See TracBrowser for help on using the repository browser.