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

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

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

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

File size: 15.2 KB
Line 
1MODULE asmphyto2dbal_hadocc
2   !!======================================================================
3   !!                       ***  MODULE asmphyto2dbal_hadocc  ***
4   !! Calculate increments to HadOCC based on surface phyto2d increments
5   !!
6   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
7   !! For licensing reasons this is kept in its own internal Met Office
8   !! branch (dev/frdf/vn3.6_nitrogen_balancing) rather than in the Paris
9   !! repository, and must be merged in when building.
10   !!
11   !!======================================================================
12   !! History :  3.6  ! 2017-08  (D. Ford)     Adapted from 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_phyto2d_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_phyto2d_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_phyto2d_bal_hadocc( ld_chltot,                      &
70      &                               pinc_chltot,                    &
71      &                               ld_phytot,                      &
72      &                               pinc_phytot,                    &
73      &                               pincper,                        &
74      &                               p_maxchlinc, ld_phytobal, pmld, &
75      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
76      &                               phyt_avg_bkg, mld_max_bkg,      &
77      &                               cchl_p_bkg,                     &
78      &                               tracer_bkg, phyto2d_balinc )
79      !!---------------------------------------------------------------------------
80      !!                    ***  ROUTINE asm_phyto2d_bal_hadocc  ***
81      !!
82      !! ** Purpose :   calculate increments to HadOCC from 2d phytoplankton increments
83      !!
84      !! ** Method  :   call nitrogen balancing scheme
85      !!
86      !! ** Action  :   populate phyto2d_balinc
87      !!
88      !! References :   Hemmings et al., 2008, J. Mar. Res.
89      !!                Ford et al., 2012, Ocean Sci.
90      !!---------------------------------------------------------------------------
91      !!
92      LOGICAL,  INTENT(in   )                               :: ld_chltot      ! Assim chltot y/n
93      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_chltot    ! chltot increments
94      LOGICAL,  INTENT(in   )                               :: ld_phytot      ! Assim phytot y/n
95      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pinc_phytot    ! phytot increments
96      REAL(wp), INTENT(in   )                               :: pincper        ! Assimilation period
97      REAL(wp), INTENT(in   )                               :: p_maxchlinc    ! Max chl increment
98      LOGICAL,  INTENT(in   )                               :: ld_phytobal    ! Balancing y/n
99      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj)           :: pmld           ! Mixed layer depth
100      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth
101      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss
102      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto
103      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD
104      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: cchl_p_bkg     ! Surface C:Chl
105      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables
106      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: phyto2d_balinc ! Balancing increments
107      !!
108      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
109      INTEGER                                               :: jkmax          ! Loop index
110      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
111      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
112      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters
113      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state
114      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments
115      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics
116      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics
117      !!---------------------------------------------------------------------------
118
119      IF ( ( .NOT. ld_chltot ) .AND. ( .NOT. ld_phytot ) ) THEN
120         CALL ctl_stop( ' Trying to do phyto2d balancing but nothing to assimilate' )
121      ENDIF
122     
123      ! If p_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
124      IF ( p_maxchlinc > 0.0 ) THEN
125         IF ( ld_chltot ) THEN
126            DO jj = 1, jpj
127               DO ji = 1, jpi
128                  pinc_chltot(ji,jj) = MAX( -1.0 * p_maxchlinc, MIN( pinc_chltot(ji,jj), p_maxchlinc ) )
129               END DO
130            END DO
131         ENDIF
132      ENDIF
133
134      IF ( ld_phytot ) THEN
135         CALL ctl_stop( ' No phytoplankton carbon assimilation quite yet' )
136      ENDIF
137     
138      IF ( ld_phytobal ) THEN   ! Nitrogen balancing
139
140         ! Set up model parameters to be passed into Hemmings balancing routine
141         modparm(1)  = grow_sat
142         modparm(2)  = psmax
143         modparm(3)  = par
144         modparm(4)  = alpha
145         modparm(5)  = resp_rate
146         modparm(6)  = pmort_rate
147         modparm(7)  = phyto_min
148         modparm(8)  = z_mort_1
149         modparm(9)  = z_mort_2
150         modparm(10) = c2n_p
151         modparm(11) = c2n_z
152         modparm(12) = c2n_d
153         modparm(13) = graze_threshold
154         modparm(14) = holling_coef
155         modparm(15) = graze_sat
156         modparm(16) = graze_max
157
158         ! Set up assimilation parameters to be passed into balancing routine
159         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
160         assimparm(2)  = balnutext
161         assimparm(3)  = balnutmin
162         assimparm(4)  = r
163         assimparm(5)  = beta_g
164         assimparm(6)  = beta_l
165         assimparm(7)  = beta_m
166         assimparm(8)  = a_g
167         assimparm(9)  = a_l
168         assimparm(10) = a_m
169         assimparm(11) = zfracb0
170         assimparm(12) = zfracb1
171         assimparm(13) = qrfmax
172         assimparm(14) = qafmax
173         assimparm(15) = zrfmax
174         assimparm(16) = zafmax
175         assimparm(17) = prfmax
176         assimparm(18) = incphymin
177         assimparm(19) = integnstep
178         assimparm(20) = pthreshold
179
180         ! Set up external tracer indices array bstate
181         i_tracer(1) = 1   ! nutrient
182         i_tracer(2) = 2   ! phytoplankton
183         i_tracer(3) = 3   ! zooplankton
184         i_tracer(4) = 4   ! detritus
185         i_tracer(5) = 5   ! DIC
186         i_tracer(6) = 6   ! Alkalinity
187
188         ! Set background state
189         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jp_had_nut)
190         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jp_had_phy)
191         bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jp_had_zoo)
192         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jp_had_pdn)
193         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jp_had_dic)
194         bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jp_had_alk)
195
196         ! Call nitrogen balancing routine
197         CALL bio_analysis( jpi, jpj, jpk, ZDZ(:,:,:), i_tracer, modparm,               &
198            &               n2be_p, n2be_z, n2be_d, assimparm,                          &
199            &               INT(pincper), 1, kmt(:,:), tmask(:,:,:),                    &
200            &               pmld(:,:), mld_max_bkg(:,:), pinc_chltot(:,:), cchl_p_bkg(:,:), &
201            &               nbal_active, phyt_avg_bkg(:,:),                             &
202            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),          &
203            &               subsurf_active, deepneg_active,                             &
204            &               deeppos_active, nutprof_active,                             &
205            &               bstate, outincs,                                            &
206            &               diag_active, diag,                                          &
207            &               diag_fulldepth_active, diag_fulldepth )
208
209         ! Save balancing increments
210         phyto2d_balinc(:,:,:,jp_had_nut) = outincs(:,:,:,i_tracer(1))
211         phyto2d_balinc(:,:,:,jp_had_phy) = outincs(:,:,:,i_tracer(2))
212         phyto2d_balinc(:,:,:,jp_had_zoo) = outincs(:,:,:,i_tracer(3))
213         phyto2d_balinc(:,:,:,jp_had_pdn) = outincs(:,:,:,i_tracer(4))
214         phyto2d_balinc(:,:,:,jp_had_dic) = outincs(:,:,:,i_tracer(5))
215         phyto2d_balinc(:,:,:,jp_had_alk) = outincs(:,:,:,i_tracer(6))
216     
217      ELSE   ! No nitrogen balancing
218     
219         ! Initialise phytoplankton increment to zero
220         phyto2d_balinc(:,:,:,jp_had_phy) = 0.0
221         
222         ! Convert surface chlorophyll increment to phytoplankton nitrogen
223         phyto2d_balinc(:,:,1,jp_had_phy) = ( cchl_p_bkg(:,:) / (mw_carbon * c2n_p) ) * pinc_chltot(:,:)
224         
225         ! Propagate through mixed layer
226         DO jj = 1, jpj
227            DO ji = 1, jpi
228               !
229               jkmax = jpk-1
230               DO jk = jpk-1, 1, -1
231                  IF ( ( pmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
232                     & ( pmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
233                     pmld(ji,jj) = gdepw_n(ji,jj,jk+1)
234                     jkmax = jk
235                  ENDIF
236               END DO
237               !
238               DO jk = 2, jkmax
239                  phyto2d_balinc(ji,jj,jk,jp_had_phy) = phyto2d_balinc(ji,jj,1,jp_had_phy)
240               END DO
241               !
242            END DO
243         END DO
244
245         ! Set other balancing increments to zero
246         phyto2d_balinc(:,:,:,jp_had_nut) = 0.0
247         phyto2d_balinc(:,:,:,jp_had_zoo) = 0.0
248         phyto2d_balinc(:,:,:,jp_had_pdn) = 0.0
249         phyto2d_balinc(:,:,:,jp_had_dic) = 0.0
250         phyto2d_balinc(:,:,:,jp_had_alk) = 0.0
251     
252      ENDIF
253
254   END SUBROUTINE asm_phyto2d_bal_hadocc
255
256#else
257   !!----------------------------------------------------------------------
258   !!   Default option : Empty routine
259   !!----------------------------------------------------------------------
260CONTAINS
261   SUBROUTINE asm_phyto2d_bal_hadocc( ld_chltot,                      &
262      &                               pinc_chltot,                    &
263      &                               ld_phytot,                      &
264      &                               pinc_phytot,                    &
265      &                               pincper,                        &
266      &                               p_maxchlinc, ld_phytobal, pmld, &
267      &                               pgrow_avg_bkg, ploss_avg_bkg,   &
268      &                               phyt_avg_bkg, mld_max_bkg,      &
269      &                               cchl_p_bkg,                     &
270      &                               tracer_bkg, phyto2d_balinc )
271      LOGICAL :: ld_chltot
272      REAL    :: pinc_chltot(:,:)
273      LOGICAL :: ld_phytot
274      REAL    :: pinc_phytot(:,:)
275      REAL    :: pincper
276      REAL    :: p_maxchlinc
277      LOGICAL :: ld_phytobal
278      REAL    :: pmld(:,:)
279      REAL    :: pgrow_avg_bkg(:,:)
280      REAL    :: ploss_avg_bkg(:,:)
281      REAL    :: phyt_avg_bkg(:,:)
282      REAL    :: mld_max_bkg(:,:)
283      REAL    :: cchl_p_bkg(:,:)
284      REAL    :: tracer_bkg(:,:,:,:)
285      REAL    :: phyto2d_balinc(:,:,:,:)
286      WRITE(*,*) 'asm_phyto2d_bal_hadocc: You should not have seen this print! error?'
287   END SUBROUTINE asm_phyto2d_bal_hadocc
288#endif
289
290   !!======================================================================
291END MODULE asmphyto2dbal_hadocc
Note: See TracBrowser for help on using the repository browser.