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_hadocc.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_hadocc.F90 @ 8648

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

Change source of MLD for direct initialisation of logchl and pCO2 to work.
File size: 16.6 KB
Line 
1MODULE asmlogchlbal_hadocc
2   !!======================================================================
3   !!                       ***  MODULE asmlogchlbal_hadocc  ***
4   !! Calculate increments to HadOCC based on surface logchl increments
5   !!
6   !! IMPORTANT NOTE: This calls the bioanalysis routine of Hemmings et al.
7   !! For licensing reasons this is kept in its own internal Met Office
8   !! branch (dev/frdf/vn3.6_nitrogen_balancing) rather than in the Paris
9   !! repository, and must be merged in when building.
10   !!
11   !!======================================================================
12   !! History :  3.6  ! 2017-08  (D. Ford)     Adapted from 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_logchl_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 zdfmxl                              ! mixed layer depth
25   USE iom                                 ! i/o
26   USE par_hadocc                          ! HadOCC parameters
27   USE had_bgc_stnd,  ONLY: kmt            ! HadOCC parameters
28   USE had_bgc_const                       ! HadOCC parameters
29   USE par_trc,       ONLY: jptra          ! Tracer parameters
30   USE bioanalysis                         ! Nitrogen balancing
31
32   IMPLICIT NONE
33   PRIVATE                   
34
35   PUBLIC asm_logchl_bal_hadocc
36
37   ! Default values for biological assimilation parameters
38   ! Should match Hemmings et al. (2008)
39   REAL(wp), PARAMETER :: balnutext  =  0.6    !: Default nutrient balancing factor
40   REAL(wp), PARAMETER :: balnutmin  =  0.1    !: Fraction of phytoplankton loss to nutrient
41   REAL(wp), PARAMETER :: r          =  1      !: Reliability of model specific growth rate
42   REAL(wp), PARAMETER :: beta_g     =  0.05   !: Low rate bias correction for growth rate estimator
43   REAL(wp), PARAMETER :: beta_l     =  0.05   !: Low rate bias correction for primary loss rate estimator
44   REAL(wp), PARAMETER :: beta_m     =  0.05   !: Low rate bias correction for secondary loss rate estimator
45   REAL(wp), PARAMETER :: a_g        =  0.2    !: Error s.d. for log10 of growth rate estimator
46   REAL(wp), PARAMETER :: a_l        =  0.4    !: Error s.d. for log10 of primary loss rate estimator
47   REAL(wp), PARAMETER :: a_m        =  0.7    !: Error s.d. for log10 of secondary loss rate estimator
48   REAL(wp), PARAMETER :: zfracb0    =  0.7    !: Base zooplankton fraction of loss to Z & D
49   REAL(wp), PARAMETER :: zfracb1    =  0      !: Phytoplankton sensitivity of zooplankton fraction
50   REAL(wp), PARAMETER :: qrfmax     =  1.1    !: Maximum nutrient limitation reduction factor
51   REAL(wp), PARAMETER :: qafmax     =  1.1    !: Maximum nutrient limitation amplification factor
52   REAL(wp), PARAMETER :: zrfmax     =  2      !: Maximum zooplankton reduction factor
53   REAL(wp), PARAMETER :: zafmax     =  2      !: Maximum zooplankton amplification factor
54   REAL(wp), PARAMETER :: prfmax     =  10     !: Maximum phytoplankton reduction factor (secondary)
55   REAL(wp), PARAMETER :: incphymin  =  0.0001 !: Minimum size of non-zero phytoplankton increment
56   REAL(wp), PARAMETER :: integnstep =  20     !: Number of steps for p.d.f. integral evaluation
57   REAL(wp), PARAMETER :: pthreshold =  0.01   !: Fractional threshold level for setting p.d.f.
58   !
59   LOGICAL,  PARAMETER :: diag_active           = .TRUE.  !: Depth-independent diagnostics
60   LOGICAL,  PARAMETER :: diag_fulldepth_active = .TRUE.  !: Full-depth diagnostics
61   LOGICAL,  PARAMETER :: gl_active             = .TRUE.  !: Growth/loss-based balancing
62   LOGICAL,  PARAMETER :: nbal_active           = .TRUE.  !: Nitrogen balancing
63   LOGICAL,  PARAMETER :: subsurf_active        = .TRUE.  !: Increments below MLD
64   LOGICAL,  PARAMETER :: deepneg_active        = .FALSE. !: Negative primary increments below MLD
65   LOGICAL,  PARAMETER :: deeppos_active        = .FALSE. !: Positive primary increments below MLD
66   LOGICAL,  PARAMETER :: nutprof_active        = .TRUE.  !: Secondary increments
67
68CONTAINS
69
70   SUBROUTINE asm_logchl_bal_hadocc( logchl_bkginc, aincper, mld_choice_bgc, &
71      &                              k_maxchlinc, ld_logchlbal, ld_asmdin,   &
72      &                              pgrow_avg_bkg, ploss_avg_bkg,           &
73      &                              phyt_avg_bkg, mld_max_bkg,              &
74      &                              chl_bkg, cchl_p_bkg,                    &
75      &                              tracer_bkg, logchl_balinc )
76      !!---------------------------------------------------------------------------
77      !!                    ***  ROUTINE asm_logchl_bal_hadocc  ***
78      !!
79      !! ** Purpose :   calculate increments to HadOCC from logchl increments
80      !!
81      !! ** Method  :   convert logchl increments to chl increments
82      !!                call nitrogen balancing scheme
83      !!
84      !! ** Action  :   populate logchl_balinc
85      !!
86      !! References :   Hemmings et al., 2008, J. Mar. Res.
87      !!                Ford et al., 2012, Ocean Sci.
88      !!---------------------------------------------------------------------------
89      !!
90      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: logchl_bkginc  ! logchl increments
91      REAL(wp), INTENT(in   )                               :: aincper        ! Assimilation period
92      INTEGER,  INTENT(in   )                               :: mld_choice_bgc ! MLD criterion
93      REAL(wp), INTENT(in   )                               :: k_maxchlinc    ! Max chl increment
94      LOGICAL,  INTENT(in   )                               :: ld_logchlbal   ! Balancing y/n
95      LOGICAL,  INTENT(in   )                               :: ld_asmdin      ! Direct initialisation y/n
96      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: pgrow_avg_bkg  ! Avg phyto growth
97      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: ploss_avg_bkg  ! Avg phyto loss
98      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: phyt_avg_bkg   ! Avg phyto
99      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: mld_max_bkg    ! Max MLD
100      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: chl_bkg        ! Surface chlorophyll
101      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)           :: cchl_p_bkg     ! Surface C:Chl
102      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk,jptra) :: tracer_bkg     ! State variables
103      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc  ! Balancing increments
104      !!
105      INTEGER                                               :: ji, jj, jk, jn ! Loop counters
106      INTEGER                                               :: jkmax          ! Loop index
107      INTEGER,                 DIMENSION(6)                 :: i_tracer       ! Tracer indices
108      REAL(wp),                DIMENSION(jpi,jpj)           :: chl_inc        ! Chlorophyll increments
109      REAL(wp),                DIMENSION(jpi,jpj)           :: zmld           ! Mixed layer depth
110      REAL(wp),                DIMENSION(16)                :: modparm        ! Model parameters
111      REAL(wp),                DIMENSION(20)                :: assimparm      ! Assimilation parameters
112      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: bstate         ! Background state
113      REAL(wp),                DIMENSION(jpi,jpj,jpk,6)     :: outincs        ! Balancing increments
114      REAL(wp),                DIMENSION(jpi,jpj,22)        :: diag           ! Depth-indep diagnostics
115      REAL(wp),                DIMENSION(jpi,jpj,jpk,22)    :: diag_fulldepth ! Full-depth diagnostics
116      !!---------------------------------------------------------------------------
117
118      ! Convert log10(chlorophyll) increment back to a chlorophyll increment
119      ! In order to transform logchl incs to chl incs, need to account for model
120      ! background, cannot simply do 10^logchl_bkginc. Need to:
121      ! 1) Add logchl inc to log10(background) to get log10(analysis)
122      ! 2) Take 10^log10(analysis) to get analysis
123      ! 3) Subtract background from analysis to get chl incs
124      ! If k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
125      DO jj = 1, jpj
126         DO ji = 1, jpi
127            IF ( chl_bkg(ji,jj) > 0.0 ) THEN
128               chl_inc(ji,jj) = 10**( LOG10( chl_bkg(ji,jj) ) + logchl_bkginc(ji,jj) ) - chl_bkg(ji,jj)
129               IF ( k_maxchlinc > 0.0 ) THEN
130                  chl_inc(ji,jj) = MAX( -1.0 * k_maxchlinc, MIN( chl_inc(ji,jj), k_maxchlinc ) )
131               ENDIF
132            ELSE
133               chl_inc(ji,jj) = 0.0
134            ENDIF
135         END DO
136      END DO
137     
138      ! Select mixed layer
139      IF ( ld_asmdin ) THEN
140         CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', &
141            &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
142         zmld(:,:) = mld_max_bkg(:,:)
143      ELSE
144         SELECT CASE( mld_choice_bgc )
145         CASE ( 1 )                   ! Turbocline/mixing depth [W points]
146            zmld(:,:) = hmld(:,:)
147         CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
148            zmld(:,:) = hmlp(:,:)
149         CASE ( 3 )                   ! Kara MLD [Interpolated]
150#if defined key_karaml
151            IF ( ln_kara ) THEN
152               zmld(:,:) = hmld_kara(:,:)
153            ELSE
154               CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
155                  &           ' but ln_kara=.false.' )
156            ENDIF
157#else
158            CALL ctl_stop( ' Kara mixed layer requested for LogChl assimilation,', &
159               &           ' but is not defined' )
160#endif
161         CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
162            !zmld(:,:) = hmld_tref(:,:)
163            CALL ctl_stop( ' hmld_tref mixed layer requested for LogChl assimilation,', &
164               &           ' but is not available in this version' )
165         CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
166            zmld(:,:) = hmlpt(:,:)
167         END SELECT
168      ENDIF
169     
170      IF ( ld_logchlbal ) THEN   ! Nitrogen balancing
171
172         ! Set up model parameters to be passed into Hemmings balancing routine
173         modparm(1)  = grow_sat
174         modparm(2)  = psmax
175         modparm(3)  = par
176         modparm(4)  = alpha
177         modparm(5)  = resp_rate
178         modparm(6)  = pmort_rate
179         modparm(7)  = phyto_min
180         modparm(8)  = z_mort_1
181         modparm(9)  = z_mort_2
182         modparm(10) = c2n_p
183         modparm(11) = c2n_z
184         modparm(12) = c2n_d
185         modparm(13) = graze_threshold
186         modparm(14) = holling_coef
187         modparm(15) = graze_sat
188         modparm(16) = graze_max
189
190         ! Set up assimilation parameters to be passed into balancing routine
191         ! Not sure what assimparm(1) is meant to be, but it doesn't get used
192         assimparm(2)  = balnutext
193         assimparm(3)  = balnutmin
194         assimparm(4)  = r
195         assimparm(5)  = beta_g
196         assimparm(6)  = beta_l
197         assimparm(7)  = beta_m
198         assimparm(8)  = a_g
199         assimparm(9)  = a_l
200         assimparm(10) = a_m
201         assimparm(11) = zfracb0
202         assimparm(12) = zfracb1
203         assimparm(13) = qrfmax
204         assimparm(14) = qafmax
205         assimparm(15) = zrfmax
206         assimparm(16) = zafmax
207         assimparm(17) = prfmax
208         assimparm(18) = incphymin
209         assimparm(19) = integnstep
210         assimparm(20) = pthreshold
211
212         ! Set up external tracer indices array bstate
213         i_tracer(1) = 1   ! nutrient
214         i_tracer(2) = 2   ! phytoplankton
215         i_tracer(3) = 3   ! zooplankton
216         i_tracer(4) = 4   ! detritus
217         i_tracer(5) = 5   ! DIC
218         i_tracer(6) = 6   ! Alkalinity
219
220         ! Set background state
221         bstate(:,:,:,i_tracer(1)) = tracer_bkg(:,:,:,jp_had_nut)
222         bstate(:,:,:,i_tracer(2)) = tracer_bkg(:,:,:,jp_had_phy)
223         bstate(:,:,:,i_tracer(3)) = tracer_bkg(:,:,:,jp_had_zoo)
224         bstate(:,:,:,i_tracer(4)) = tracer_bkg(:,:,:,jp_had_pdn)
225         bstate(:,:,:,i_tracer(5)) = tracer_bkg(:,:,:,jp_had_dic)
226         bstate(:,:,:,i_tracer(6)) = tracer_bkg(:,:,:,jp_had_alk)
227
228         ! Call nitrogen balancing routine
229         CALL bio_analysis( jpi, jpj, jpk, ZDZ(:,:,:), i_tracer, modparm,               &
230            &               n2be_p, n2be_z, n2be_d, assimparm,                          &
231            &               INT(aincper), 1, kmt(:,:), tmask(:,:,:),                    &
232            &               zmld(:,:), mld_max_bkg(:,:), chl_inc(:,:), cchl_p_bkg(:,:), &
233            &               nbal_active, phyt_avg_bkg(:,:),                             &
234            &               gl_active, pgrow_avg_bkg(:,:), ploss_avg_bkg(:,:),          &
235            &               subsurf_active, deepneg_active,                             &
236            &               deeppos_active, nutprof_active,                             &
237            &               bstate, outincs,                                            &
238            &               diag_active, diag,                                          &
239            &               diag_fulldepth_active, diag_fulldepth )
240
241         ! Save balancing increments
242         logchl_balinc(:,:,:,jp_had_nut) = outincs(:,:,:,i_tracer(1))
243         logchl_balinc(:,:,:,jp_had_phy) = outincs(:,:,:,i_tracer(2))
244         logchl_balinc(:,:,:,jp_had_zoo) = outincs(:,:,:,i_tracer(3))
245         logchl_balinc(:,:,:,jp_had_pdn) = outincs(:,:,:,i_tracer(4))
246         logchl_balinc(:,:,:,jp_had_dic) = outincs(:,:,:,i_tracer(5))
247         logchl_balinc(:,:,:,jp_had_alk) = outincs(:,:,:,i_tracer(6))
248     
249      ELSE   ! No nitrogen balancing
250     
251         ! Initialise phytoplankton increment to zero
252         logchl_balinc(:,:,:,jp_had_phy) = 0.0
253         
254         ! Convert surface chlorophyll increment to phytoplankton nitrogen
255         logchl_balinc(:,:,1,jp_had_phy) = ( cchl_p_bkg(:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:)
256         
257         ! Propagate through mixed layer
258         DO jj = 1, jpj
259            DO ji = 1, jpi
260               !
261               jkmax = jpk-1
262               DO jk = jpk-1, 1, -1
263                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
264                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
265                     zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
266                     jkmax = jk
267                  ENDIF
268               END DO
269               !
270               DO jk = 2, jkmax
271                  logchl_balinc(ji,jj,jk,jp_had_phy) = logchl_balinc(ji,jj,1,jp_had_phy)
272               END DO
273               !
274            END DO
275         END DO
276
277         ! Set other balancing increments to zero
278         logchl_balinc(:,:,:,jp_had_nut) = 0.0
279         logchl_balinc(:,:,:,jp_had_zoo) = 0.0
280         logchl_balinc(:,:,:,jp_had_pdn) = 0.0
281         logchl_balinc(:,:,:,jp_had_dic) = 0.0
282         logchl_balinc(:,:,:,jp_had_alk) = 0.0
283     
284      ENDIF
285
286   END SUBROUTINE asm_logchl_bal_hadocc
287
288#else
289   !!----------------------------------------------------------------------
290   !!   Default option : Empty routine
291   !!----------------------------------------------------------------------
292CONTAINS
293   SUBROUTINE asm_logchl_bal_hadocc( logchl_bkginc, aincper, mld_choice_bgc, &
294      &                              k_maxchlinc, ld_logchlbal,              &
295      &                              pgrow_avg_bkg, ploss_avg_bkg,           &
296      &                              phyt_avg_bkg, mld_max_bkg,              &
297      &                              chl_bkg, cchl_p_bkg,                    &
298      &                              tracer_bkg, logchl_balinc )
299      REAL    :: logchl_bkginc(:,:)
300      REAL    :: aincper
301      INTEGER :: mld_choice_bgc
302      REAL    :: k_maxchlinc
303      LOGICAL :: ld_logchlbal
304      REAL    :: pgrow_avg_bkg(:,:)
305      REAL    :: ploss_avg_bkg(:,:)
306      REAL    :: phyt_avg_bkg(:,:)
307      REAL    :: mld_max_bkg(:,:)
308      REAL    :: chl_bkg(:,:)
309      REAL    :: cchl_p_bkg(:,:)
310      REAL    :: tracer_bkg(:,:,:,:)
311      REAL    :: logchl_balinc(:,:,:,:)
312      WRITE(*,*) 'asm_logchl_bal_hadocc: You should not have seen this print! error?'
313   END SUBROUTINE asm_logchl_bal_hadocc
314#endif
315
316   !!======================================================================
317END MODULE asmlogchlbal_hadocc
Note: See TracBrowser for help on using the repository browser.